program subroutine_und_function
  implicit none
  integer :: l, q, quadratsumme
  double precision :: x, ll
  !quadratsumme ist als Funktion ganzzahlig definiert.
  !Trotzdem wird der Name "quadratsumme" ohne Deklaration
  !als "real" interpretiert, und das Ergebnis ist immer Null.
  !Es ist zu beachten, dass ganzzahlige Variablen nur 
  !bis 2^31 (etwa 2.15e9) gehen. Die Quadratsumme wird also falsch
  !fuer groessere l (ab l=1861)
  call biswieweit(l)  !<=== Aufruf der Subroutine
  q=quadratsumme(l)   !<=== Aufruf des Funktions-Unterprogramms
  write(*,fmt=*) 'Quadratsumme ist: ',q
  ll=l; if (ll > 0.0d0) then
    x=ll*(ll+1)*(2*ll+1)/6
    else; x=0.0d0
  end if
  write(*,fmt=*) 'Quadratsumme doppelt genau zum Vergleich: ',x
end program subroutine_und_function

subroutine biswieweit(n)
  implicit none
  integer, intent(out) :: n
  write(*,fmt="('Bitte ganzzahliges n angeben: ')",advance="no")
  read(*,fmt=*) n
end subroutine biswieweit

function quadratsumme(x) result(y)
  implicit none
  integer :: j,x,y,s 
  !Die Summe der Quadrate y=1+2^2+...+x^2 ist 
  !bis zu einem vorgegebenen x auszurechnen
  !Das x wird ueber die Subroutine "biswieweit" eingeschleust.
  s=0
  do j=1,x
    s=s+j*j
  end do
  y=s;
end function quadratsumme
