27 public :: swfdfwtype, dfw_cr
32 integer(I4B),
pointer :: is2d => null()
33 integer(I4B),
pointer :: icentral => null()
34 integer(I4B),
pointer :: iswrcond => null()
35 real(DP),
pointer :: unitconv
36 real(DP),
pointer :: timeconv
37 real(DP),
pointer :: lengthconv
38 real(DP),
dimension(:),
pointer,
contiguous :: hnew => null()
39 real(DP),
dimension(:),
pointer,
contiguous :: manningsn => null()
40 integer(I4B),
dimension(:),
pointer,
contiguous :: idcxs => null()
41 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
42 integer(I4B),
dimension(:),
pointer,
contiguous :: icelltype => null()
45 integer(I4B),
pointer :: icalcvelocity => null()
46 integer(I4B),
pointer :: isavvelocity => null()
47 real(DP),
dimension(:, :),
pointer,
contiguous :: vcomp => null()
48 real(DP),
dimension(:),
pointer,
contiguous :: vmag => null()
49 integer(I4B),
pointer :: nedges => null()
50 integer(I4B),
pointer :: lastedge => null()
51 integer(I4B),
dimension(:),
pointer,
contiguous :: nodedge => null()
52 integer(I4B),
dimension(:),
pointer,
contiguous :: ihcedge => null()
53 real(DP),
dimension(:, :),
pointer,
contiguous :: propsedge => null()
54 real(DP),
dimension(:),
pointer,
contiguous :: grad_dhds_mag => null()
55 real(DP),
dimension(:),
pointer,
contiguous :: dhdsja => null()
58 integer(I4B),
pointer :: inobspkg => null()
59 type(ObsType),
pointer :: obs => null()
62 type(SwfCxsType),
pointer :: cxs
67 procedure :: allocate_scalars
68 procedure :: allocate_arrays
70 procedure :: source_options
71 procedure :: log_options
72 procedure :: source_griddata
73 procedure :: log_griddata
78 procedure :: dfw_qnm_fc_nr
84 procedure :: dfw_save_model_flows
85 procedure :: dfw_print_model_flows
87 procedure :: dfw_df_obs
88 procedure :: dfw_rp_obs
89 procedure :: dfw_bd_obs
92 procedure :: get_cond_swr
93 procedure :: get_cond_n
94 procedure :: get_flow_area_nm
95 procedure :: calc_velocity
96 procedure :: sav_velocity
97 procedure,
public :: increase_edge_count
98 procedure,
public :: set_edge_properties
99 procedure :: calc_dhds
100 procedure :: write_cxs_tables
108 subroutine dfw_cr(dfwobj, name_model, input_mempath, inunit, iout, &
113 type(SwfDfwType),
pointer :: dfwobj
114 character(len=*),
intent(in) :: name_model
115 character(len=*),
intent(in) :: input_mempath
116 integer(I4B),
intent(in) :: inunit
117 integer(I4B),
intent(in) :: iout
118 type(SwfCxsType),
pointer,
intent(in) :: cxs
120 logical(LGP) :: found_fname
122 character(len=*),
parameter :: fmtheader = &
123 "(1x, /1x, 'DFW -- DIFFUSIVE WAVE (DFW) PACKAGE, VERSION 1, 9/25/2023', &
124 &' INPUT READ FROM MEMPATH: ', A, /)"
130 call dfwobj%set_names(1, name_model,
'DFW',
'DFW')
133 call dfwobj%allocate_scalars()
136 dfwobj%input_mempath = input_mempath
137 dfwobj%inunit = inunit
141 call mem_set_value(dfwobj%input_fname,
'INPUT_FNAME', dfwobj%input_mempath, &
148 call obs_cr(dfwobj%obs, dfwobj%inobspkg)
154 write (iout, fmtheader) input_mempath
158 end subroutine dfw_cr
162 subroutine dfw_df(this, dis)
164 class(SwfDfwType) :: this
165 class(DisBaseType),
pointer,
intent(inout) :: dis
171 if (this%dis%is_2d())
then
180 call this%allocate_arrays()
187 end subroutine dfw_df
195 subroutine allocate_scalars(this)
198 class(SwfDfwtype) :: this
201 call this%NumericalPackageType%allocate_scalars()
205 call mem_allocate(this%icentral,
'ICENTRAL', this%memoryPath)
206 call mem_allocate(this%iswrcond,
'ISWRCOND', this%memoryPath)
207 call mem_allocate(this%unitconv,
'UNITCONV', this%memoryPath)
208 call mem_allocate(this%lengthconv,
'LENGTHCONV', this%memoryPath)
209 call mem_allocate(this%timeconv,
'TIMECONV', this%memoryPath)
210 call mem_allocate(this%inobspkg,
'INOBSPKG', this%memoryPath)
211 call mem_allocate(this%icalcvelocity,
'ICALCVELOCITY', this%memoryPath)
212 call mem_allocate(this%isavvelocity,
'ISAVVELOCITY', this%memoryPath)
213 call mem_allocate(this%nedges,
'NEDGES', this%memoryPath)
214 call mem_allocate(this%lastedge,
'LASTEDGE', this%memoryPath)
220 this%lengthconv =
done
223 this%icalcvelocity = 0
224 this%isavvelocity = 0
228 end subroutine allocate_scalars
232 subroutine allocate_arrays(this)
234 class(SwfDfwType) :: this
240 'MANNINGSN', this%memoryPath)
242 'IDCXS', this%memoryPath)
244 'ICELLTYPE', this%memoryPath)
247 call mem_allocate(this%nodedge, 0,
'NODEDGE', this%memoryPath)
248 call mem_allocate(this%ihcedge, 0,
'IHCEDGE', this%memoryPath)
249 call mem_allocate(this%propsedge, 0, 0,
'PROPSEDGE', this%memoryPath)
252 call mem_allocate(this%vcomp, 3, 0,
'VCOMP', this%memoryPath)
253 call mem_allocate(this%vmag, 0,
'VMAG', this%memoryPath)
255 do n = 1, this%dis%nodes
256 this%manningsn(n) =
dzero
258 this%icelltype(n) = 1
262 if (this%is2d == 1)
then
264 'GRAD_DHDS_MAG', this%memoryPath)
266 'DHDSJA', this%memoryPath)
267 do n = 1, this%dis%nodes
268 this%grad_dhds_mag(n) =
dzero
270 do n = 1, this%dis%njas
271 this%dhdsja(n) =
dzero
275 end subroutine allocate_arrays
279 subroutine dfw_load(this)
281 class(SwfDfwType) :: this
284 call this%source_options()
285 call this%source_griddata()
287 end subroutine dfw_load
291 subroutine source_options(this)
299 class(SwfDfwType) :: this
301 integer(I4B) :: isize
302 type(SwfDfwParamFoundType) :: found
303 type(CharacterStringType),
dimension(:),
pointer, &
304 contiguous :: obs6_fnames
308 this%input_mempath, found%icentral)
310 this%input_mempath, found%iswrcond)
312 this%input_mempath, found%lengthconv)
314 this%input_mempath, found%timeconv)
316 this%input_mempath, found%iprflow)
318 this%input_mempath, found%ipakcb)
320 this%input_mempath, found%isavvelocity)
323 if (found%icentral) this%icentral = 1
324 if (found%ipakcb) this%ipakcb = -1
327 this%unitconv = this%lengthconv**
donethird
328 this%unitconv = this%unitconv * this%timeconv
331 if (found%isavvelocity) this%icalcvelocity = this%isavvelocity
334 call get_isize(
'OBS6_FILENAME', this%input_mempath, isize)
338 errmsg =
'Multiple OBS6 keywords detected in OPTIONS block.'// &
339 ' Only one OBS6 entry allowed.'
344 call mem_setptr(obs6_fnames,
'OBS6_FILENAME', this%input_mempath)
346 found%obs6_filename = .true.
347 this%obs%inputFilename = obs6_fnames(1)
348 this%obs%active = .true.
350 this%obs%inUnitObs = this%inobspkg
351 call openfile(this%inobspkg, this%iout, this%obs%inputFilename,
'OBS')
352 call this%obs%obs_df(this%iout, this%packName, this%filtyp, this%dis)
353 call this%dfw_df_obs()
357 if (this%iout > 0)
then
358 call this%log_options(found)
361 end subroutine source_options
365 subroutine log_options(this, found)
367 class(SwfDfwType) :: this
368 type(SwfDfwParamFoundType),
intent(in) :: found
370 write (this%iout,
'(1x,a)')
'Setting DFW Options'
372 if (found%lengthconv)
then
373 write (this%iout,
'(4x,a, G0)')
'Mannings length conversion value &
374 &specified as ', this%lengthconv
377 if (found%timeconv)
then
378 write (this%iout,
'(4x,a, G0)')
'Mannings time conversion value &
379 &specified as ', this%timeconv
382 if (found%lengthconv .or. found%timeconv)
then
383 write (this%iout,
'(4x,a, G0)')
'Mannings conversion value calculated &
384 &from user-provided length_conversion and &
385 &time_conversion is ', this%unitconv
388 if (found%iprflow)
then
389 write (this%iout,
'(4x,a)')
'Cell-by-cell flow information will be printed &
390 &to listing file whenever ICBCFL is not zero.'
393 if (found%ipakcb)
then
394 write (this%iout,
'(4x,a)')
'Cell-by-cell flow information will be printed &
395 &to listing file whenever ICBCFL is not zero.'
398 if (found%obs6_filename)
then
399 write (this%iout,
'(4x,a)')
'Observation package is active.'
402 if (found%isavvelocity) &
403 write (this%iout,
'(4x,a)')
'Velocity will be calculated at cell &
404 ¢ers and written to DATA-VCOMP in budget &
405 &file when requested.'
407 if (found%iswrcond)
then
408 write (this%iout,
'(4x,a, G0)')
'Conductance will be calculated using &
409 &the SWR development option.'
412 write (this%iout,
'(1x,a,/)')
'End Setting DFW Options'
414 end subroutine log_options
418 subroutine source_griddata(this)
426 class(SwfDfwType) :: this
428 character(len=LENMEMPATH) :: idmMemoryPath
429 type(SwfDfwParamFoundType) :: found
430 integer(I4B),
dimension(:),
pointer,
contiguous :: map
437 if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
441 idmmemorypath, map, found%manningsn)
442 call mem_set_value(this%idcxs,
'IDCXS', idmmemorypath, map, found%idcxs)
445 if (.not. found%manningsn)
then
446 write (errmsg,
'(a)')
'Error in GRIDDATA block: MANNINGSN not found.'
451 call store_error_filename(this%input_fname)
455 if (this%iout > 0)
then
456 call this%log_griddata(found)
459 end subroutine source_griddata
463 subroutine log_griddata(this, found)
465 class(SwfDfwType) :: this
466 type(SwfDfwParamFoundType),
intent(in) :: found
468 write (this%iout,
'(1x,a)')
'Setting DFW Griddata'
470 if (found%manningsn)
then
471 write (this%iout,
'(4x,a)')
'MANNINGSN set from input file'
474 if (found%idcxs)
then
475 write (this%iout,
'(4x,a)')
'IDCXS set from input file'
478 call this%write_cxs_tables()
480 write (this%iout,
'(1x,a,/)')
'End Setting DFW Griddata'
482 end subroutine log_griddata
484 subroutine write_cxs_tables(this)
487 class(SwfDfwType) :: this
500 end subroutine write_cxs_tables
504 subroutine dfw_ar(this, ibound, hnew)
507 class(SwfDfwType) :: this
508 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
509 real(DP),
dimension(:),
pointer,
contiguous,
intent(inout) :: hnew
514 this%ibound => ibound
517 if (this%icalcvelocity == 1)
then
518 call mem_reallocate(this%vcomp, 3, this%dis%nodes,
'VCOMP', this%memoryPath)
519 call mem_reallocate(this%vmag, this%dis%nodes,
'VMAG', this%memoryPath)
520 call mem_reallocate(this%nodedge, this%nedges,
'NODEDGE', this%memoryPath)
521 call mem_reallocate(this%ihcedge, this%nedges,
'IHCEDGE', this%memoryPath)
522 call mem_reallocate(this%propsedge, 5, this%nedges,
'PROPSEDGE', &
524 do n = 1, this%dis%nodes
525 this%vcomp(:, n) =
dzero
531 call this%obs%obs_ar()
533 end subroutine dfw_ar
537 subroutine dfw_rp(this)
540 class(SwfDfwType) :: this
543 call this%dfw_rp_obs()
545 end subroutine dfw_rp
549 subroutine dfw_ad(this, irestore)
550 class(SwfDfwType) :: this
551 integer(I4B),
intent(in) :: irestore
554 call this%obs%obs_ad()
556 end subroutine dfw_ad
564 subroutine dfw_fc(this, kiter, matrix_sln, idxglo, rhs, stage, stage_old)
567 class(SwfDfwType) :: this
568 integer(I4B) :: kiter
569 class(MatrixBaseType),
pointer :: matrix_sln
570 integer(I4B),
intent(in),
dimension(:) :: idxglo
571 real(DP),
intent(inout),
dimension(:) :: rhs
572 real(DP),
intent(inout),
dimension(:) :: stage
573 real(DP),
intent(inout),
dimension(:) :: stage_old
577 if (this%is2d == 1)
then
578 call this%calc_dhds()
582 call this%dfw_qnm_fc_nr(kiter, matrix_sln, idxglo, rhs, stage, stage_old)
584 end subroutine dfw_fc
589 subroutine dfw_qnm_fc_nr(this, kiter, matrix_sln, idxglo, rhs, stage, stage_old)
593 class(SwfDfwType) :: this
594 integer(I4B) :: kiter
595 class(MatrixBaseType),
pointer :: matrix_sln
596 integer(I4B),
intent(in),
dimension(:) :: idxglo
597 real(DP),
intent(inout),
dimension(:) :: rhs
598 real(DP),
intent(inout),
dimension(:) :: stage
599 real(DP),
intent(inout),
dimension(:) :: stage_old
601 integer(I4B) :: n, m, ii, idiag
608 do n = 1, this%dis%nodes
611 idiag = this%dis%con%ia(n)
614 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
617 if (this%dis%con%mask(ii) == 0) cycle
620 m = this%dis%con%ja(ii)
623 qnm = this%qcalc(n, m, stage(n), stage(m), ii)
624 rhs(n) = rhs(n) - qnm
628 qeps = this%qcalc(n, m, stage(n) + eps, stage(m), ii)
629 derv = (qeps - qnm) / eps
630 call matrix_sln%add_value_pos(idxglo(idiag), derv)
631 rhs(n) = rhs(n) + derv * stage(n)
635 qeps = this%qcalc(n, m, stage(n), stage(m) + eps, ii)
636 derv = (qeps - qnm) / eps
637 call matrix_sln%add_value_pos(idxglo(ii), derv)
638 rhs(n) = rhs(n) + derv * stage(m)
643 end subroutine dfw_qnm_fc_nr
647 subroutine dfw_fn(this, kiter, matrix_sln, idxglo, rhs, stage)
649 class(SwfDfwType) :: this
650 integer(I4B) :: kiter
651 class(MatrixBaseType),
pointer :: matrix_sln
652 integer(I4B),
intent(in),
dimension(:) :: idxglo
653 real(DP),
intent(inout),
dimension(:) :: rhs
654 real(DP),
intent(inout),
dimension(:) :: stage
661 end subroutine dfw_fn
665 function qcalc(this, n, m, stage_n, stage_m, ipos)
result(qnm)
667 class(SwfDfwType) :: this
668 integer(I4B),
intent(in) :: n
669 integer(I4B),
intent(in) :: m
670 real(DP),
intent(in) :: stage_n
671 real(DP),
intent(in) :: stage_m
672 integer(I4B),
intent(in) :: ipos
674 integer(I4B) :: isympos
681 isympos = this%dis%con%jas(ipos)
683 cl1 = this%dis%con%cl1(isympos)
684 cl2 = this%dis%con%cl2(isympos)
686 cl1 = this%dis%con%cl2(isympos)
687 cl2 = this%dis%con%cl1(isympos)
691 if (this%iswrcond == 0)
then
692 cond = this%get_cond(n, m, ipos, stage_n, stage_m, cl1, cl2)
693 else if (this%iswrcond == 1)
then
694 cond = this%get_cond_swr(n, m, ipos, stage_n, stage_m, cl1, cl2)
698 qnm = cond * (stage_m - stage_n)
707 function get_cond(this, n, m, ipos, stage_n, stage_m, cln, clm)
result(cond)
711 class(SwfDfwType) :: this
712 integer(I4B),
intent(in) :: n
713 integer(I4B),
intent(in) :: m
714 integer(I4B),
intent(in) :: ipos
715 real(DP),
intent(in) :: stage_n
716 real(DP),
intent(in) :: stage_m
717 real(DP),
intent(in) :: cln
718 real(DP),
intent(in) :: clm
726 real(DP) :: range = 1.d-6
728 real(DP) :: smooth_factor
729 real(DP) :: length_nm
737 length_nm = cln + clm
739 if (length_nm >
dprec)
then
742 depth_n = stage_n - this%dis%bot(n)
743 depth_m = stage_m - this%dis%bot(m)
746 if (this%is2d == 0)
then
747 dhds_n = abs(stage_m - stage_n) / (cln + clm)
750 dhds_n = this%grad_dhds_mag(n)
751 dhds_m = this%grad_dhds_mag(m)
755 if (this%icentral == 0)
then
757 if (stage_n > stage_m)
then
766 call squadratic(depth_n, range, dydx, smooth_factor)
767 depth_n = depth_n * smooth_factor
768 call squadratic(depth_m, range, dydx, smooth_factor)
769 depth_m = depth_m * smooth_factor
772 call this%dis%get_flow_width(n, m, ipos, width_n, width_m)
776 cn = this%get_cond_n(n, depth_n, cln, width_n, dhds_n)
777 cm = this%get_cond_n(m, depth_m, clm, width_m, dhds_m)
782 if ((cn + cm) >
dprec)
then
783 cond = cn * cm / (cn + cm)
790 end function get_cond
796 function get_cond_n(this, n, depth, dx, width, dhds)
result(c)
799 class(SwfDfwType) :: this
800 integer(I4B),
intent(in) :: n
801 real(DP),
intent(in) :: depth
802 real(DP),
intent(in) :: dx
803 real(DP),
intent(in) :: width
804 real(DP),
intent(in) :: dhds
810 real(DP) :: conveyance
813 rough = this%manningsn(n)
814 conveyance = this%cxs%get_conveyance(this%idcxs(n), width, depth, rough)
815 dhds_sqr = dhds**
dhalf
816 if (dhds_sqr <
dem10)
then
821 c = this%unitconv * conveyance / dx / dhds_sqr
823 end function get_cond_n
831 function get_cond_swr(this, n, m, ipos, stage_n, stage_m, cln, clm)
result(cond)
835 class(SwfDfwType) :: this
836 integer(I4B),
intent(in) :: n
837 integer(I4B),
intent(in) :: m
838 integer(I4B),
intent(in) :: ipos
839 real(DP),
intent(in) :: stage_n
840 real(DP),
intent(in) :: stage_m
841 real(DP),
intent(in) :: cln
842 real(DP),
intent(in) :: clm
852 real(DP) :: range = 1.d-6
854 real(DP) :: smooth_factor
855 real(DP) :: length_nm
859 real(DP) :: area_n, area_m, area_avg
860 real(DP) :: rhn, rhm, rhavg
868 length_nm = cln + clm
870 if (length_nm >
dprec)
then
873 depth_n = stage_n - this%dis%bot(n)
874 depth_m = stage_m - this%dis%bot(m)
877 if (this%icentral == 0)
then
879 if (stage_n > stage_m)
then
888 call squadratic(depth_n, range, dydx, smooth_factor)
889 depth_n = depth_n * smooth_factor
890 call squadratic(depth_m, range, dydx, smooth_factor)
891 depth_m = depth_m * smooth_factor
894 call this%dis%get_flow_width(n, m, ipos, width_n, width_m)
897 weight_n = clm / length_nm
898 weight_m =
done - weight_n
901 area_n = this%cxs%get_area(this%idcxs(n), width_n, depth_n)
902 area_m = this%cxs%get_area(this%idcxs(m), width_m, depth_m)
903 area_avg = weight_n * area_n + weight_m * area_m
906 if (this%is2d == 0)
then
907 rhn = this%cxs%get_hydraulic_radius(this%idcxs(n), width_n, &
909 rhm = this%cxs%get_hydraulic_radius(this%idcxs(m), width_m, &
911 rhavg = weight_n * rhn + weight_m * rhm
913 rhavg = area_avg / width_n
918 if (this%is2d == 0)
then
919 dhds_nm = abs(stage_m - stage_n) / (length_nm)
921 dhds_n = this%grad_dhds_mag(n)
922 dhds_m = this%grad_dhds_mag(m)
923 dhds_nm = weight_n * dhds_n + weight_m * dhds_m
925 dhds_sqr = dhds_nm**
dhalf
926 if (dhds_sqr <
dem10)
then
931 weight_n = cln / length_nm
932 weight_m =
done - weight_n
933 rough_n = this%cxs%get_roughness(this%idcxs(n), width_n, depth_n, &
934 this%manningsn(n), dhds_nm)
935 rough_m = this%cxs%get_roughness(this%idcxs(m), width_m, depth_m, &
936 this%manningsn(m), dhds_nm)
937 ravg = (weight_n + weight_m) / &
938 (weight_n / rough_n + weight_m / rough_m)
939 rinv_avg =
done / ravg
942 cond = this%unitconv * rinv_avg * area_avg * rhavg / dhds_sqr / length_nm
946 end function get_cond_swr
954 function get_flow_area_nm(this, n, m, stage_n, stage_m, cln, clm, &
955 ipos)
result(area_avg)
959 class(SwfDfwType) :: this
960 integer(I4B),
intent(in) :: n
961 integer(I4B),
intent(in) :: m
962 real(DP),
intent(in) :: stage_n
963 real(DP),
intent(in) :: stage_m
964 real(DP),
intent(in) :: cln
965 real(DP),
intent(in) :: clm
966 integer(I4B),
intent(in) :: ipos
976 real(DP) :: length_nm
977 real(DP) :: range = 1.d-6
979 real(DP) :: smooth_factor
984 depth_n = stage_n - this%dis%bot(n)
985 depth_m = stage_m - this%dis%bot(m)
988 if (this%icentral == 0)
then
990 if (stage_n > stage_m)
then
999 call squadratic(depth_n, range, dydx, smooth_factor)
1000 depth_n = depth_n * smooth_factor
1001 call squadratic(depth_m, range, dydx, smooth_factor)
1002 depth_m = depth_m * smooth_factor
1005 call this%dis%get_flow_width(n, m, ipos, width_n, width_m)
1008 length_nm = cln + clm
1009 weight_n = clm / length_nm
1010 weight_m =
done - weight_n
1013 area_n = this%cxs%get_area(this%idcxs(n), width_n, depth_n)
1014 area_m = this%cxs%get_area(this%idcxs(m), width_m, depth_m)
1015 area_avg = weight_n * area_n + weight_m * area_m
1017 end function get_flow_area_nm
1025 subroutine calc_dhds(this)
1029 class(SwfDfwType) :: this
1033 integer(I4B) :: ipos
1034 integer(I4B) :: isympos
1038 do n = 1, this%dis%nodes
1039 this%grad_dhds_mag(n) =
dzero
1040 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1041 m = this%dis%con%ja(ipos)
1042 isympos = this%dis%con%jas(ipos)
1046 cl1 = this%dis%con%cl1(isympos)
1047 cl2 = this%dis%con%cl2(isympos)
1049 cl1 = this%dis%con%cl2(isympos)
1050 cl2 = this%dis%con%cl1(isympos)
1055 if (cl1 + cl2 >
dprec)
then
1056 this%dhdsja(isympos) = (this%hnew(m) - this%hnew(n)) / (cl1 + cl2)
1058 this%dhdsja(isympos) =
dzero
1068 end subroutine calc_dhds
1073 subroutine dfw_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax)
1075 class(SwfDfwType) :: this
1076 integer(I4B),
intent(in) :: neqmod
1077 real(DP),
dimension(neqmod),
intent(inout) :: x
1078 real(DP),
dimension(neqmod),
intent(in) :: xtemp
1079 real(DP),
dimension(neqmod),
intent(inout) :: dx
1080 integer(I4B),
intent(inout) :: inewtonur
1081 real(DP),
intent(inout) :: dxmax
1082 integer(I4B),
intent(inout) :: locmax
1090 do n = 1, this%dis%nodes
1091 if (this%ibound(n) < 1) cycle
1092 if (this%icelltype(n) > 0)
then
1093 botm = this%dis%bot(n)
1096 if (x(n) < botm)
then
1100 if (abs(dxx) > abs(dxmax))
then
1110 end subroutine dfw_nur
1114 subroutine dfw_cq(this, stage, stage_old, flowja)
1116 class(SwfDfwType) :: this
1117 real(DP),
intent(inout),
dimension(:) :: stage
1118 real(DP),
intent(inout),
dimension(:) :: stage_old
1119 real(DP),
intent(inout),
dimension(:) :: flowja
1121 integer(I4B) :: n, ipos, m
1124 do n = 1, this%dis%nodes
1125 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1126 m = this%dis%con%ja(ipos)
1128 qnm = this%qcalc(n, m, stage(n), stage(m), ipos)
1130 flowja(this%dis%con%isym(ipos)) = -qnm
1134 end subroutine dfw_cq
1138 subroutine dfw_bd(this, isuppress_output, model_budget)
1142 class(SwfDfwType) :: this
1143 integer(I4B),
intent(in) :: isuppress_output
1144 type(BudgetType),
intent(inout) :: model_budget
1149 end subroutine dfw_bd
1153 subroutine dfw_save_model_flows(this, flowja, icbcfl, icbcun)
1155 class(SwfDfwType) :: this
1156 real(DP),
dimension(:),
intent(in) :: flowja
1157 integer(I4B),
intent(in) :: icbcfl
1158 integer(I4B),
intent(in) :: icbcun
1160 integer(I4B) :: ibinun
1163 if (this%ipakcb < 0)
then
1165 elseif (this%ipakcb == 0)
then
1168 ibinun = this%ipakcb
1170 if (icbcfl == 0) ibinun = 0
1173 if (ibinun /= 0)
then
1175 call this%dis%record_connection_array(flowja, ibinun, this%iout)
1179 if (this%isavvelocity /= 0)
then
1180 if (ibinun /= 0)
call this%sav_velocity(ibinun)
1183 end subroutine dfw_save_model_flows
1187 subroutine dfw_print_model_flows(this, ibudfl, flowja)
1192 class(SwfDfwType) :: this
1193 integer(I4B),
intent(in) :: ibudfl
1194 real(DP),
intent(inout),
dimension(:) :: flowja
1196 character(len=LENBIGLINE) :: line
1197 character(len=30) :: tempstr
1198 integer(I4B) :: n, ipos, m
1201 character(len=*),
parameter :: fmtiprflow = &
1202 &
"(/,4x,'CALCULATED INTERCELL FLOW FOR PERIOD ', i0, ' STEP ', i0)"
1205 if (ibudfl /= 0 .and. this%iprflow > 0)
then
1206 write (this%iout, fmtiprflow)
kper,
kstp
1207 do n = 1, this%dis%nodes
1209 call this%dis%noder_to_string(n, tempstr)
1210 line = trim(tempstr)//
':'
1211 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1212 m = this%dis%con%ja(ipos)
1213 call this%dis%noder_to_string(m, tempstr)
1214 line = trim(line)//
' '//trim(tempstr)
1216 write (tempstr,
'(1pg15.6)') qnm
1217 line = trim(line)//
' '//trim(adjustl(tempstr))
1219 write (this%iout,
'(a)') trim(line)
1223 end subroutine dfw_print_model_flows
1227 subroutine dfw_da(this)
1233 class(SwfDfwType) :: this
1247 if (this%is2d == 1)
then
1266 call this%obs%obs_da()
1267 deallocate (this%obs)
1272 call this%NumericalPackageType%da()
1277 end subroutine dfw_da
1283 subroutine calc_velocity(this, flowja)
1286 class(SwfDfwType) :: this
1287 real(DP),
intent(in),
dimension(:) :: flowja
1291 integer(I4B) :: ipos
1292 integer(I4B) :: isympos
1318 real(DP),
allocatable,
dimension(:) :: vi
1319 real(DP),
allocatable,
dimension(:) :: di
1320 real(DP),
allocatable,
dimension(:) :: viz
1321 real(DP),
allocatable,
dimension(:) :: diz
1322 real(DP),
allocatable,
dimension(:) :: nix
1323 real(DP),
allocatable,
dimension(:) :: niy
1324 real(DP),
allocatable,
dimension(:) :: wix
1325 real(DP),
allocatable,
dimension(:) :: wiy
1326 real(DP),
allocatable,
dimension(:) :: wiz
1327 real(DP),
allocatable,
dimension(:) :: bix
1328 real(DP),
allocatable,
dimension(:) :: biy
1329 logical :: nozee = .true.
1333 if (this%icalcvelocity /= 0 .and. this%dis%con%ianglex == 0)
then
1334 call store_error(
'Error. ANGLDEGX not provided in '// &
1335 'discretization file. ANGLDEGX required for '// &
1336 'calculation of velocity.', terminate=.true.)
1341 do n = 1, this%dis%nodes
1344 ic = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
1347 do m = 1, this%nedges
1348 if (this%nodedge(m) == n)
then
1354 if (ic > nc) nc = ic
1371 do n = 1, this%dis%nodes
1383 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1384 m = this%dis%con%ja(ipos)
1385 isympos = this%dis%con%jas(ipos)
1386 ihc = this%dis%con%ihc(isympos)
1388 call this%dis%connection_normal(n, m, ihc, xn, yn, zn, ipos)
1389 call this%dis%connection_vector(n, m, nozee,
done,
done, &
1390 ihc, xc, yc, zc, dltot)
1391 cl1 = this%dis%con%cl1(isympos)
1392 cl2 = this%dis%con%cl2(isympos)
1394 cl1 = this%dis%con%cl2(isympos)
1395 cl2 = this%dis%con%cl1(isympos)
1397 ooclsum =
done / (cl1 + cl2)
1400 di(ic) = dltot * cl1 * ooclsum
1401 area = this%get_flow_area_nm(n, m, this%hnew(n), this%hnew(m), &
1403 if (area >
dzero)
then
1404 vi(ic) = flowja(ipos) / area
1413 do m = 1, this%nedges
1414 if (this%nodedge(m) == n)
then
1417 ihc = this%ihcedge(m)
1418 area = this%propsedge(2, m)
1421 nix(ic) = -this%propsedge(3, m)
1422 niy(ic) = -this%propsedge(4, m)
1423 di(ic) = this%propsedge(5, m)
1424 if (area >
dzero)
then
1425 vi(ic) = this%propsedge(1, m) / area
1443 dsumz = dsumz + diz(iz)
1445 denom = (ncz -
done)
1447 dsumz = dsumz +
dem10 * dsumz
1449 if (dsumz >
dzero) wiz(iz) =
done - diz(iz) / dsumz
1451 wiz(iz) = wiz(iz) / denom
1459 vz = vz + wiz(iz) * viz(iz)
1468 wix(ic) = di(ic) * abs(nix(ic))
1469 wiy(ic) = di(ic) * abs(niy(ic))
1470 dsumx = dsumx + wix(ic)
1471 dsumy = dsumy + wiy(ic)
1478 dsumx = dsumx +
dem10 * dsumx
1479 dsumy = dsumy +
dem10 * dsumy
1481 wix(ic) = (dsumx - wix(ic)) * abs(nix(ic))
1482 wiy(ic) = (dsumy - wiy(ic)) * abs(niy(ic))
1489 bix(ic) = wix(ic) * sign(
done, nix(ic))
1490 biy(ic) = wiy(ic) * sign(
done, niy(ic))
1491 dsumx = dsumx + wix(ic) * abs(nix(ic))
1492 dsumy = dsumy + wiy(ic) * abs(niy(ic))
1499 bix(ic) = bix(ic) * dsumx
1500 biy(ic) = biy(ic) * dsumy
1501 axy = axy + bix(ic) * niy(ic)
1502 ayx = ayx + biy(ic) * nix(ic)
1515 vx = vx + (bix(ic) - axy * biy(ic)) * vi(ic)
1516 vy = vy + (biy(ic) - ayx * bix(ic)) * vi(ic)
1518 denom =
done - axy * ayx
1519 if (denom /=
dzero)
then
1524 this%vcomp(1, n) = vx
1525 this%vcomp(2, n) = vy
1526 this%vcomp(3, n) = vz
1527 this%vmag(n) = sqrt(vx**2 + vy**2 + vz**2)
1542 end subroutine calc_velocity
1549 subroutine increase_edge_count(this, nedges)
1551 class(SwfDfwType) :: this
1552 integer(I4B),
intent(in) :: nedges
1554 this%nedges = this%nedges + nedges
1556 end subroutine increase_edge_count
1562 subroutine set_edge_properties(this, nodedge, ihcedge, q, area, nx, ny, &
1565 class(SwfDfwType) :: this
1566 integer(I4B),
intent(in) :: nodedge
1567 integer(I4B),
intent(in) :: ihcedge
1568 real(DP),
intent(in) :: q
1569 real(DP),
intent(in) :: area
1570 real(DP),
intent(in) :: nx
1571 real(DP),
intent(in) :: ny
1572 real(DP),
intent(in) :: distance
1574 integer(I4B) :: lastedge
1576 this%lastedge = this%lastedge + 1
1577 lastedge = this%lastedge
1578 this%nodedge(lastedge) = nodedge
1579 this%ihcedge(lastedge) = ihcedge
1580 this%propsedge(1, lastedge) = q
1581 this%propsedge(2, lastedge) = area
1582 this%propsedge(3, lastedge) = nx
1583 this%propsedge(4, lastedge) = ny
1584 this%propsedge(5, lastedge) = distance
1588 if (this%lastedge == this%nedges) this%lastedge = 0
1590 end subroutine set_edge_properties
1596 subroutine sav_velocity(this, ibinun)
1598 class(SwfDfwType) :: this
1599 integer(I4B),
intent(in) :: ibinun
1601 character(len=16) :: text
1602 character(len=16),
dimension(3) :: auxtxt
1604 integer(I4B) :: naux
1607 text =
' DATA-VCOMP'
1609 auxtxt(:) = [
' vx',
' vy',
' vz']
1610 call this%dis%record_srcdst_list_header(text, this%name_model, &
1611 this%packName, this%name_model, &
1612 this%packName, naux, auxtxt, ibinun, &
1613 this%dis%nodes, this%iout)
1616 do n = 1, this%dis%nodes
1617 call this%dis%record_mf6_list_entry(ibinun, n, n,
dzero, naux, &
1621 end subroutine sav_velocity
1626 subroutine dfw_df_obs(this)
1628 class(SwfDfwType) :: this
1630 integer(I4B) :: indx
1634 call this%obs%StoreObsType(
'ext-outflow', .true., indx)
1635 this%obs%obsData(indx)%ProcessIdPtr => dfwobsidprocessor
1637 end subroutine dfw_df_obs
1639 subroutine dfwobsidprocessor(obsrv, dis, inunitobs, iout)
1641 type(ObserveType),
intent(inout) :: obsrv
1642 class(DisBaseType),
intent(in) :: dis
1643 integer(I4B),
intent(in) :: inunitobs
1644 integer(I4B),
intent(in) :: iout
1647 character(len=LINELENGTH) :: string
1650 string = obsrv%IDstring
1654 obsrv%NodeNumber = n
1656 errmsg =
'Error reading data from ID string'
1661 end subroutine dfwobsidprocessor
1666 subroutine dfw_bd_obs(this)
1668 class(SwfDfwType) :: this
1674 character(len=100) :: msg
1675 type(ObserveType),
pointer :: obsrv => null()
1678 if (this%obs%npakobs > 0)
then
1679 call this%obs%obs_bd_clear()
1680 do i = 1, this%obs%npakobs
1681 obsrv => this%obs%pakobs(i)%obsrv
1682 do j = 1, obsrv%indxbnds_count
1683 n = obsrv%indxbnds(j)
1685 select case (obsrv%ObsTypeId)
1687 msg =
'Unrecognized observation type: '//trim(obsrv%ObsTypeId)
1690 call this%obs%SaveOneSimval(obsrv, v)
1700 end subroutine dfw_bd_obs
1705 subroutine dfw_rp_obs(this)
1709 class(SwfDfwType),
intent(inout) :: this
1714 class(ObserveType),
pointer :: obsrv => null()
1720 do i = 1, this%obs%npakobs
1721 obsrv => this%obs%pakobs(i)%obsrv
1724 nn1 = obsrv%NodeNumber
1725 if (nn1 < 1 .or. nn1 > this%dis%nodes)
then
1726 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
1727 trim(adjustl(obsrv%ObsTypeId)), &
1728 'reach must be greater than 0 and less than or equal to', &
1729 this%dis%nodes,
'(specified value is ', nn1,
')'
1732 if (obsrv%indxbnds_count == 0)
then
1733 call obsrv%AddObsIndex(nn1)
1735 errmsg =
'Programming error in dfw_rp_obs'
1741 do j = 1, obsrv%indxbnds_count
1742 nn1 = obsrv%indxbnds(j)
1743 if (nn1 < 1 .or. nn1 > this%dis%nodes)
then
1744 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
1745 trim(adjustl(obsrv%ObsTypeId)), &
1746 'reach must be greater than 0 and less than or equal to', &
1747 this%dis%nodes,
'(specified value is ', nn1,
')'
1760 end subroutine dfw_rp_obs
1762 end module swfdfwmodule
This module contains the BudgetModule.
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
real(dp), parameter dtwothirds
real constant 2/3
real(dp), parameter dp9
real constant 9/10
real(dp), parameter dem10
real constant 1e-10
real(dp), parameter donethird
real constant 1/3
integer(i4b), parameter lenbigline
maximum length of a big line
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dzero
real constant zero
real(dp), parameter dprec
real constant machine precision
real(dp), parameter dtwo
real constant 2
integer(i4b), parameter lenmempath
maximum length of the memory path
real(dp), parameter done
real constant 1
This module defines variable data types.
real(dp) function, public get_perturbation(x)
Calculate a numerical perturbation given the value of x.
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
subroutine, public memorystore_remove(component, subcomponent, context)
subroutine, public get_isize(name, mem_path, isize)
@ brief Get the number of elements for this variable
This module contains the base numerical package type.
This module contains the derived types ObserveType and ObsDataType.
This module contains the derived type ObsType.
subroutine, public obs_cr(obs, inobs)
@ brief Create a new ObsType object
This module contains simulation methods.
subroutine, public store_error(msg, terminate)
Store an error message.
integer(i4b) function, public count_errors()
Return number of errors.
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
subroutine, public store_error_unit(iunit, terminate)
Store the file unit number.
This module contains simulation variables.
character(len=maxcharlen) errmsg
error message string
character(len=linelength) idm_context
subroutine squadratic(x, range, dydx, y)
@ brief sQuadratic
integer(i4b), pointer, public kstp
current time step number
integer(i4b), pointer, public kper
current stress period number
subroutine, public vector_interpolation_2d(dis, flowja, nedges, nodedge, propsedge, vcomp, vmag, flowareaja)
Interpolate 2D vector components at cell center.
Derived type for the Budget object.
This class is used to store a single deferred-length character string. It was designed to work in an ...