#!/bin/ksh
#----------------------------------------------------------
#-- NCL Doc Example script:  NUG_use_Fortran_subroutines.ksh
#--
#--     Write a short Fortran subroutine --> ex01.f
#--     Run the wrapper                 --> ex01.so
#--     Write NCL script                --> NUG_use_Fortran_subroutines.ncl
#--     Run NCL script                  --> write results to stdout
#----------------------------------------------------------
example=${0%.*}
ncl_script=${example}.ncl
#----------------------------------------------------------
#-- write the NCL script
#----------------------------------------------------------
cat << EOF > ${ncl_script}
external EX01 "./ex01.so"
 
begin
   print("")
;-- Calculate three values of a quadratic equation
   nump = 3
   x    = (/ -1., 0.0, 1.0 /)
   qval = new(nump,float)
;-- Call the NCL version of your Fortran subroutine.     
   EX01::cquad(-1., 2., 3., nump, x, qval)
   print("Polynomial value = " + qval)     ;-- should be (/0,3,4/)
;-- Calculate an arc length.
   xc = (/ 0., 1., 2. /)
   yc = (/ 0., 1., 0. /)
;-- Call the NCL version of your Fortran function.
   arclen = EX01::arcln(nump,xc,yc)
   print("Arc length = " + arclen)         ;-- should be 2.82843
   print("")
end
EOF
#----------------------------------------------------------
#-- fortran77 code
#----------------------------------------------------------
cat << EOF > ex01.f
C NCLFORTSTART
      subroutine cquad (a, b, c, nq, x, quad)
      real x(nq), quad(nq)
C NCLEND
C
C  Calculate quadratic polynomial values.
C
      do 10 i=1,nq
        quad(i) = a*x(i)**2 + b*x(i) + c
   10 continue

      return
      end

C NCLFORTSTART
      function arcln (numpnt, pointx, pointy)
      dimension pointx(numpnt),pointy(numpnt)
C NCLEND

C
C  Calculate arc lengths.
C
      if (numpnt .lt. 2) then
        print *, 'arcln: number of points must be at least 2'
        stop
      endif
      arcln = 0.
      do 10 i=2,numpnt
        pdist = sqrt((pointx(i)-pointx(i-1))**2 + 
     +                        (pointy(i)-pointy(i-1))**2)
        arcln = arcln + pdist
   10 continue

      return
      end
EOF

#----------------------------------------------------------
#-- run the NCL wrapper. Generates the ex01.so shared object file.
#----------------------------------------------------------
WRAPIT ex01.f

#----------------------------------------------------------
#-- run ncl
#----------------------------------------------------------
ncl -n ${ncl_script}

exit
