pure function expScaleFunction(hdiag) result(Si) real(dp), intent(in) :: hdiag real(dp) :: Si Si = 1.0 / (sfBeta * exp(sFAlpha * hdiag)) end function expScaleFunction