Skip to content

copy relevant code from fm/feature/UNST-8560_nest_02 #580

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mihai-pricope wants to merge 1 commit into
base: main
Choose a base branch
from
Open

copy relevant code from fm/feature/UNST-8560_nest_02 #580

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 12 additions & 8 deletions src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_io/unc_write_his.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
246 changes: 199 additions & 47 deletions src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/fill_valobs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Loading