calculate_wrf_dvgch Module Function

    module function calculate_wrf_dvgch(self, h) result(theta_w)
        implicit none
        class(type_wrf_dvgch), intent(in) :: self
        real(real64), intent(in) :: h
        real(real64) :: theta_w

        if (h < 0.0d0) then
            theta_w = self%theta_r + (self%theta_s - self%theta_r) * &
                      (self%w1 * (1.0d0 + abs(self%alpha1 * h)**self%n1)**(-self%m1) &
                       + self%w2 * (1.0d0 + abs(self%alpha1 * h)**self%n2)**(-self%m2))
        else
            theta_w = self%theta_s
        end if

    end function calculate_wrf_dvgch