Self compiling fortran happy number program in a bash script

#Author Matthew Watson. This file is in the public domain.
#!/bin/bash
#Self compiling fortran program in a bash script. Compiles and runs happyab
#Call it sh happy.sh 3 1 for max of 10^3 and output to file on.
#sh happy.sh 8 0 for max of 10^8 and no output
if [ -e "./.happy$1$2" ]
  then ./.happy$1$2>gzip -f --fast>happy.txt.gz && exit
  else //usr/bin/tail -n +10 $0 | gfortran -O3 -cpp -DPTEN=$1 -DOUT=$2 -o .happy$1$2 -x f95 - && ./.happy$1$2>gzip -f --fast>happy.txt.gz && exit
fi
#if PTEN < 3 || !(defined PTEN)
#define POWTEN 3
#else
#define POWTEN PTEN
#endif
#if OUT==1 || !(defined OUT)
#define DO_OUTPUT .TRUE.
#else
#define DO_OUTPUT .FALSE.
#endif
!Calculate all of the happy numbers up to 10^POWTEN and write the
!result to happy.txt if DO_OUTPUT is set
program happy
  implicit none
  logical,parameter :: writeoutput = DO_OUTPUT
  integer,parameter :: maxnum = 10**POWTEN

  !Happy numbers are not very dense, for any value we can possibly
  !calculate, there are going to be less than 0.15*n happy numbers.
  integer,parameter :: numhap = 15*10**(POWTEN-2)

  !We only need to cache a number if we're going to look it up.
  !So if our largest number is 99..n..9 or smaller, we only need n*9*9
  !cache values.
  integer,parameter :: numcache = (POWTEN) * 90

  integer :: idx = 1, w = 2, num = 2
  integer :: happyI(numhap)

  !Cache unhappy numbers rather than happy ones. They're ~10 times as
  !common and short cirtuiting the if block speeds things up marginally
  !and implies that our number is happy iff the first check is false and
  !w<num.

  logical :: unHappyArr(numcache) = .FALSE.

  unHappyArr(4) = .TRUE.        !We need one unhappy number to get started
  happyI(1) = 1

  do                            !In this block we calculate the sum of squares
    if (num > numcache) then    !of digits, check if it's in our unhappy cache.
      exit                      !if so, cache it. If not, check if we've
    end if                      !checked that number already.
                                !If we have, it must be happy, so num must be
    w = sum_sq(w)               !happy. Add it to our list.

    if(unHappyArr(w)) then

      if(num < numcache) then
        unHappyArr(num) = .TRUE.
      end if
      num = num + 1
      w = num

    else if(w < num) then       !Take advantage of the fact that if a number
                                !is less than num, we must have already
      idx = idx + 1             !checked it.
      happyI(idx) = num
      num = num + 1
      w = num
    end if

  end do
  do                            !Same as the above loop, but we're done caching
    if (num > maxnum) then      !numbers.
      exit                      !Also, once we are past 162, the possible
    end if                      !values for w are all less than num
                                !so we no longer need to check w<num
    w = sum_sq(w)

    if(unHappyArr(w)) then
      num = num + 1
      w = num
    else
      idx = idx + 1
      happyI(idx) = num
      num = num + 1
      w = num
    end if

  end do
  if(writeoutput) then
    print *,happyI(1:idx)
  elseif(idx < 1000) then         !Don't print if we have lots, too slow.
    print *, happyI(1:idx)
  end if

contains

  !Calculate the sum of the squares of the digits
  !of a number, then return that sum.


  function sum_sq (number) result (result)

    implicit none
    integer, intent (in) :: number
    integer :: result,digit,rest,work
    result = 0
    work = number
    do
      if (work == 0) then             !exit if the remaining number is 0
        exit                          !then we have nothing to do. return
      end if
      rest = work / 10                !Integer divide by 10 to shift right
      digit = work - 10 * rest        !Take mod 10 (one digit)
      result = result + digit * digit !add square
      work = rest                     !Set number to shifted and repeat
    end do

  end function sum_sq

end program happy