Skip to content
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

revise Gamma distribution in develop branch #205

Merged
merged 1 commit into from
Aug 25, 2021
Merged

revise Gamma distribution in develop branch #205

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
14 changes: 6 additions & 8 deletions route/build/src/process_param.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,13 +35,12 @@ SUBROUTINE basinUH(dt, fshape, tscale, IERR, MESSAGE)
IMPLICIT NONE
! input
REAL(DP), INTENT(IN) :: dt ! model time step
REAL(DP), INTENT(IN) :: fshape ! shapef parameter in gamma distribution
REAL(DP), INTENT(IN) :: tscale ! time scale parameter
REAL(DP), INTENT(IN) :: fshape ! shape parameter in gamma distribution
REAL(DP), INTENT(IN) :: tscale ! time scale parameter in gamma distribution
! output
INTEGER(I4B), INTENT(OUT) :: IERR ! error code
CHARACTER(*), INTENT(OUT) :: MESSAGE ! error message
! locals
REAL(DP) :: alamb ! scale parameter
REAL(DP) :: ntdh_min ! minimum number of time delay points
REAL(DP) :: ntdh_max ! maximum number of time delay points
REAL(DP) :: ntdh_try ! trial number of time delay points
Expand All @@ -56,11 +55,10 @@ SUBROUTINE basinUH(dt, fshape, tscale, IERR, MESSAGE)
! ---------------------------------------------------------------------------------------

ierr=0; message='basinUH/'
! use a Gamma distribution with shape parameter, fshape = 2.5, and time parameter, tscale, input
alamb = fshape/tscale ! scale parameter
! use a Gamma distribution with shape parameter, fshape, and time parameter, tscale, input
! find the desired number of future time steps
! check if the cummulative Gamma distribution is close to 1.00 for given model time step, tscale and fsahpe.
X_VALUE = alamb*dt
X_VALUE = dt/tscale
cumprob = gammp(fshape, X_VALUE)
if(cumprob > 0.999_dp) then ! in case if the cumprob is close to 1 in one model time step
ntdh_try = 1.999_dp
Expand All @@ -69,7 +67,7 @@ SUBROUTINE basinUH(dt, fshape, tscale, IERR, MESSAGE)
ntdh_max = 1000._dp
ntdh_try = 0.5_dp*(ntdh_min + ntdh_max)
do itry=1,maxtry
x_value = alamb*dt*ntdh_try
x_value = dt*ntdh_try/tscale
cumprob = gammp(fshape, x_value)
!print*, tscale, ntdh_try, cumprob, x_value, itry
if(cumprob < 0.99_dp) ntdh_min = ntdh_try
Expand All @@ -89,7 +87,7 @@ SUBROUTINE basinUH(dt, fshape, tscale, IERR, MESSAGE)
PSAVE = 0. ! cumulative probability at JTIM-1
DO JTIM=1,NTDH
TFUTURE = REAL(JTIM, kind(dp))*DT ! future time
CUMPROB = gammp(fshape,alamb*TFUTURE) ! cumulative probability at JTIM
CUMPROB = gammp(fshape,TFUTURE/tscale) ! cumulative probability at JTIM
FRAC_FUTURE(JTIM) = MAX(0._DP, CUMPROB-PSAVE) ! probability between JTIM-1 and JTIM
PSAVE = CUMPROB ! cumulative probability at JTIM-1
!WRITE(*,'(I5,1X,F20.5,1X,2(F11.5))') JTIM, TFUTURE, FRAC_FUTURE(JTIM), CUMPROB
Expand Down