module beta_dist_mod use fms_mod,only: stdlog, write_version_number, & error_mesg, FATAL use mpp_mod,only: get_unit implicit none private ! Provide values of the beta distribtion as a function of the CDF (the incomplete beta ! function). Returns a value as a function of two beta distribution parameters p, q ! (here they must be integers) and the value x of the CDF between 0 and 1. ! In this version we build tables using the NAG library function nag_beta_deviate, then ! look up values from a table. The table can be built at run time or read from ! a file (this version uses netcdf format). ! betaDeviateTable is a 3D table with dimensions ! x, p, q. The range of P and Q are specified when the tables are built. ! The arrays bounds are from 0 to nSteps + 1, just in case we draw exactly 0 or 1. ! character(len=128) :: version = '$Id: betaDistribution.F90,v 16.0 2008/07/30 22:06:18 fms Exp $' character(len=128) :: tagname = '$Name: tikal $' logical :: module_is_initialized = .false. real, parameter :: failureValue = -1. real, parameter :: Xmin = 0., Xmax = 1. integer :: Pmax, Qmax, numXSteps real, dimension(:, :, :), allocatable & :: betaDeviateTable, incompleteBetaTable character(len = 32) & :: fileName = "INPUT/BetaDistributionTable.txt" interface interpolateFromTable module procedure interpolateFromTable_s, interpolateFromTable_1D, interpolateFromTable_2D, & interpolateFromTable_3D, interpolateFromTable_4D end interface ! interpolateFromTable interface beta_deviate module procedure betaDeviate_s, betaDeviate_1D, betaDeviate_2D, & betaDeviate_3D, betaDeviate_4D end interface ! betaDeviate interface incomplete_beta module procedure incompleteBeta_s, incompleteBeta_1D, incompleteBeta_2D, & incompleteBeta_3D, incompleteBeta_4D end interface ! incompleteBeta public :: beta_dist_init, beta_deviate, incomplete_beta, beta_dist_end contains ! --------------------------------------------------------- subroutine test_beta integer :: i real :: x, inc_x, inv_inc_x, inv_x, inc_inv_x print *, "TESTING BETA" print *, "x inc(x) inv(inc(x)) inv(x) inc(inv(x))" do i = 0, numXSteps x = real(i)/real(numXSteps) inc_x = incomplete_beta( x, 5, 5) inv_inc_x = beta_deviate( inc_x, 5, 5) inv_x = beta_deviate( x, 5, 5) inc_inv_x = incomplete_beta(inv_x, 5, 5) write(*, "(5(f10.7, 1x))") x, inc_x, inv_inc_x, inv_x, inc_inv_x end do end subroutine test_beta ! --------------------------------------------------------- subroutine beta_dist_init ! Initialize the tables containing the incomplete beta function ! and its inverse (beta deviate). ! If the table parameters are supplied we use the NAG libraries to ! compute a new table and write it to a file; if just ! the file name is supplied we read the table from the file. ! !--------------------------------------------------------------------- ! if routine has already been executed, exit. !--------------------------------------------------------------------- if (module_is_initialized) return !--------------------------------------------------------------------- ! mark the module as initialized if we're able to read the tables !--------------------------------------------------------------------- call write_version_number(version, tagname) module_is_initialized = readFromFile(fileName) end subroutine beta_dist_init !--------------------------------------------------------------------- subroutine beta_dist_end !--------------------------------------------------------------------- ! be sure module has been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized ) return if(allocated(betaDeviateTable)) deallocate(betaDeviateTable) if(allocated(incompleteBetaTable)) deallocate(incompleteBetaTable) module_is_initialized = .false. end subroutine beta_dist_end !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! SEMI-PRIVATE PROCEDURES ! Not accessed directly but through generic interface ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! --------------------------------------------------------- ! Functions to look up the beta deviate (inverse incomplete beta distribution) ! from a table ! Overloaded, to allow for input arguments from 0 to 4 dimensions ! It might be more efficient to loop over dimensions higher than 1 to ! avoid using reshape. ! --------------------------------------------------------- function betaDeviate_s(x, p, q) result (betaDeviate) real, intent( in) :: x integer, intent( in) :: p, q real :: betaDeviate if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax, x < Xmin, x > Xmax /) )) then betaDeviate = FailureValue else betaDeviate = interpolateFromTable(x, p, q, betaDeviateTable) end if end function betaDeviate_s ! --------------------------------------------------------- function betaDeviate_1D(x, p, q) result (betaDeviate) real, dimension(:), intent( in) :: x integer, intent( in) :: p, q real, dimension(size(x)) :: betaDeviate if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax /) )) then betaDeviate(:) = FailureValue else betaDeviate(:) = interpolateFromTable(x, p, q, betaDeviateTable) end if end function betaDeviate_1D ! --------------------------------------------------------- function betaDeviate_2D(x, p, q) result (betaDeviate) real, dimension(:, :), intent( in) :: x integer, intent( in) :: p, q real, dimension(size(x, 1), & size(x, 2)) :: betaDeviate if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax /) )) then betaDeviate(:, :) = FailureValue else betaDeviate(:, :) = interpolateFromTable(x, p, q, betaDeviateTable) end if end function betaDeviate_2D ! --------------------------------------------------------- function betaDeviate_3D(x, p, q) result (betaDeviate) real, dimension(:, :, :), intent( in) :: x integer, intent( in) :: p, q real, dimension(size(x, 1), & size(x, 2), & size(x, 3)) :: betaDeviate if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax /) )) then betaDeviate(:, :, :) = FailureValue else betaDeviate(:, :, :) = interpolateFromTable(x, p, q, betaDeviateTable) end if end function betaDeviate_3D ! --------------------------------------------------------- function betaDeviate_4D(x, p, q) result (betaDeviate) real, dimension(:, :, :, :), intent( in) :: x integer, intent( in) :: p, q real, dimension(size(x, 1), size(x, 2), & size(x, 3), size(x, 4)) :: betaDeviate if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax /) )) then betaDeviate(:, :, :, :) = FailureValue else betaDeviate(:, :, :, :) = interpolateFromTable(x, p, q, betaDeviateTable) end if end function betaDeviate_4D ! --------------------------------------------------------- ! --------------------------------------------------------- ! Functions to look up the incomplete beta function from a table. ! Overloaded, to allow for input arguments from 0 to 4 dimensions ! It might be more efficient to loop over dimensions higher than 1 to ! avoid using reshape. ! --------------------------------------------------------- function incompleteBeta_s(x, p, q) result (incompleteBeta) real, intent( in) :: x integer, intent( in) :: p, q real :: incompleteBeta if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax, x < Xmin, x > Xmax /) )) then incompleteBeta = FailureValue else incompleteBeta = interpolateFromTable(x, p, q, incompleteBetaTable) end if end function incompleteBeta_s ! --------------------------------------------------------- function incompleteBeta_1D(x, p, q) result (incompleteBeta) real, dimension(:), intent( in) :: x integer, intent( in) :: p, q real, dimension(size(x)) :: incompleteBeta if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax /) )) then incompleteBeta(:) = FailureValue else incompleteBeta(:) = interpolateFromTable(x, p, q, incompleteBetaTable) end if end function incompleteBeta_1D ! --------------------------------------------------------- function incompleteBeta_2D(x, p, q) result (incompleteBeta) real, dimension(:, :), intent( in) :: x integer, intent( in) :: p, q real, dimension(size(x, 1), & size(x, 2)) :: incompleteBeta if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax /) )) then incompleteBeta(:, :) = FailureValue else incompleteBeta(:, :) = interpolateFromTable(x, p, q, incompleteBetaTable) end if end function incompleteBeta_2D ! --------------------------------------------------------- function incompleteBeta_3D(x, p, q) result (incompleteBeta) real, dimension(:, :, :), intent( in) :: x integer, intent( in) :: p, q real, dimension(size(x, 1), & size(x, 2), & size(x, 3)) :: incompleteBeta if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax /) )) then incompleteBeta(:, :, :) = FailureValue else incompleteBeta(:, :, :) = interpolateFromTable(x, p, q, incompleteBetaTable) end if end function incompleteBeta_3D ! --------------------------------------------------------- function incompleteBeta_4D(x, p, q) result (incompleteBeta) real, dimension(:, :, :, :), intent( in) :: x integer, intent( in) :: p, q real, dimension(size(x, 1), size(x, 2), & size(x, 3), size(x, 4)) :: incompleteBeta if (.not. module_is_initialized ) then call error_mesg ('beta_dist_mod', 'module has not been initialized', FATAL ) endif if(any( (/ p < 1, p > pMax, q < 1, q > qMax /) )) then incompleteBeta(:, :, :, :) = FailureValue else incompleteBeta(:, :, :, :) = interpolateFromTable(x, p, q, incompleteBetaTable) end if end function incompleteBeta_4D ! --------------------------------------------------------- !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! PRIVATE PROCEDURES ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! --------------------------------------------------------- function interpolateFromTable_s(x, p, q, table) result (values) real, intent( in) :: x integer, intent( in) :: p, q real, dimension(0:, :, :), intent( in) :: table real :: values ! Local variables integer :: xIndex real :: xWeight ! Check to see if we're at an endpoint of the distribution if (abs(x - xMax) < spacing(Xmax)) then values = 1. else if (abs(x - XMin) < spacing(Xmin)) then values = 0. else ! ! Linear interpolation in x ! xIndex = int( x * numXSteps) xWeight = numXSteps * x - xIndex values = (1. - xWeight) * table(xIndex , p, q) + & ( xWeight) * table(xIndex + 1, p, q) end if end function interpolateFromTable_s ! --------------------------------------------------------- function interpolateFromTable_1D(x, p, q, table) result(values) real, dimension(:), intent( in) :: x integer, intent( in) :: p, q real, dimension(0:, :, :), intent( in) :: table real, dimension(size(x)) :: values ! Local variables integer, dimension(size(x)) :: xIndex real, dimension(size(x)) :: xWeight ! Check for parameters out of bounds, and be sure the table has been initialized. ! where(x(:) < Xmin .or. x(:) > Xmax) values(:) = FailureValue else where ! ! Linear interpolation in x ! xIndex(:) = int( x(:) * numXSteps) xWeight(:) = numXSteps * x(:) - xIndex(:) values(:) = (1. - xWeight(:)) * table(xIndex(:) , p, q) + & ( xWeight(:)) * table(xIndex(:) + 1, p, q) end where ! Check to see if we're at an endpoint of the distribution where (abs(x(:) - xMax) < spacing(Xmax)) values(:) = 1. where (abs(x(:) - XMin) < spacing(Xmin)) values(:) = 0. end function interpolateFromTable_1D ! --------------------------------------------------------- function interpolateFromTable_2D(x, p, q, table) result(values) real, dimension(:, :), intent( in) :: x integer, intent( in) :: p, q real, dimension(0:, :, :), intent( in) :: table real, dimension(size(x, 1), & size(x, 2)) :: values ! Local variables integer :: i do i = 1, size(x, 2) values(:, i) = interpolateFromTable(x(:, i), p, q, table) end do end function interpolateFromTable_2D ! --------------------------------------------------------- function interpolateFromTable_3D(x, p, q, table) result(values) real, dimension(:, :, :), intent( in) :: x integer, intent( in) :: p, q real, dimension(0:, :, :), intent( in) :: table real, dimension(size(x, 1), & size(x, 2), & size(x, 3)) :: values ! Local variables integer :: i do i = 1, size(x, 3) values(:, :, i) = interpolateFromTable(x(:, :, i), p, q, table) end do end function interpolateFromTable_3D ! --------------------------------------------------------- function interpolateFromTable_4D(x, p, q, table) result(values) real, dimension(:, :, :, :), intent( in) :: x integer, intent( in) :: p, q real, dimension(0:, :, :), intent( in) :: table real, dimension(size(x, 1), & size(x, 2), & size(x, 3), & size(x, 4)) :: values ! Local variables integer :: i do i = 1, size(x, 4) values(:, :, :, i) = interpolateFromTable(x(:, :, :, i), p, q, table) end do end function interpolateFromTable_4D ! --------------------------------------------------------- ! Reading and writing procedures ! --------------------------------------------------------- function readFromFile(fileName) character(len = *), intent( in) :: fileName logical :: readFromFile ! Local variables ! integer, parameter :: unit = 909 integer :: unit unit = get_unit() open(unit = unit, file = trim(fileName), status = 'old') read(unit, '(3(i5, 1x))') Pmax, Qmax, numXSteps allocate( betaDeviateTable(0:numXSteps + 1, Pmax, Qmax), & incompleteBetaTable(0:numXSteps + 1, Pmax, Qmax)) read(unit, '(8(f10.8, 1x))') betaDeviateTable read(unit, '(8(f10.8, 1x))') incompleteBetaTable close(unit) ! print *, "Reading beta distribution tables..." ! write(*, '(8(f10.8, 1x))') betaDeviateTable(:8, 5, 5) ! write(*, '(8(f10.8, 1x))') incompleteBetaTable(:8, 5, 5) readFromFile = .true. end function readFromFile ! function readFromFile(fileName) ! use netcdf ! character(len = *), intent( in) :: fileName ! logical :: readFromFile ! ! ! Local variables - all related to netcdf ! integer, dimension(16) :: status ! integer :: fileId, ncDimId, ncVarId ! ! status( :) = nf90_NoErr ! status( 1) = nf90_open(trim(fileName), nf90_NoWrite, fileID) ! status( 2) = nf90_inq_dimid(fileId, "x", ncDimId) ! status( 3) = nf90_Inquire_Dimension(fileId, ncDimId, len = numXSteps) ! ! if(allocated(betaDeviateTable)) deallocate(betaDeviateTable) ! allocate(betaDeviateTable(0:numXSteps-1, Pmax, Qmax)) ! status( 4) = nf90_inq_varId(fileId, "betaDeviateTable", ncVarId) ! status( 5) = nf90_get_var(fileId, ncVarId, betaDeviateTable) ! ! if(allocated(incompleteBetaTable)) deallocate(incompleteBetaTable) ! allocate(incompleteBetaTable(0:numXSteps-1, Pmax, Qmax)) ! status( 7) = nf90_inq_varId(fileId, "incompleteBetaTable", ncVarId) ! status( 8) = nf90_get_var(fileId, ncVarId, betaDeviateTable) ! ! status( 9) = nf90_close(fileId) ! ! ! ! ! Make the definitions of numXSteps, consistent with what's in routine ! ! initializeIncBetaTables, which has a zero at one end an one extra element in each ! ! direction, just in case we hit the last element. ! ! ! numXSteps = numXSteps - 2 ! ! readFromFile = all(status(:) == nf90_NoErr) ! end function readFromFile end module beta_dist_mod