# calculate S/N eofs with fortran code.
# very SILLY to use when number of time points * number of ensemble member 
# is less than number of grid points.
export MP_SET_NUMTHREADS=8
ingrid<<EOF> data.r4
\beginingrid
defHTMLwords
HTMLwords 
(deviations.cdf) readfile .TAUX
0 replaceNaN 
data.r4
EOF
nt=204 #number of time points
nens=4 #number of ensemble members
cat<<EOF>prog.f90
program driver
  !program to calulate inter-ensemble eofs
  !number of time-steps times number of ensemble members
  parameter(m=$nens*$nt)
  !number grid points
  parameter(n=128*64)
  !number of unmasked grid points
  parameter(ngg=n)
  !lg2n is log base 2 of n 
  parameter(lg2n=13)
  parameter(lwork=1+5*n+2*n*lg2n+3*n**2)
  parameter(liwork=3 + 5*n)
  real, allocatable, dimension(:,:) :: xx
  double precision, allocatable, dimension(:,:) :: dxx
  double precision, allocatable, dimension(:,:) :: ts
  double precision, allocatable, dimension(:,:) :: c
  double precision w(n),tmp(n),total
  double precision, allocatable, dimension(:) :: work
  integer iwork(liwork)
  integer info,error
  
  !allocate space for and read data
  allocate(xx(n,m),stat=error)
  if  (error /= 0) then
     print *,"failed"
     stop
  end if
  
  tic = second()
  
  open(1,file='data.r4',form='unformatted',status='unknown',&
       &access='direct',RECL=m*n*4)
  read(1,REC=1) xx
  toc = second()
  print *,"Read time  ", toc-tic
  
  !allocate, convert to double and deallocate real data
  allocate(dxx(n,m),stat=error)
  if  (error /= 0) then
     print *,"failed"
     stop
  end if
  
  dxx = xx
  deallocate(xx)

  !Form covariance
  
  allocate(c(n,n),stat=error)
  if  (error /= 0) then
     print *,"failed"
     stop
  end if
  
  tic = second()
  CALL DSYRK ('U','N',n,m,1d0/(m-1d0),dxx,n,0d0,c,n)
  toc = second()
  print *,"Form covariance",toc-tic
  
  !Compute eigenvalues and eigenvectors
  
  allocate(work(lwork),stat=error)
  if  (error /= 0) then
     print *,"failed"
     stop
  end if
  tic = second()
  call DSYEVD('V','U',n,c,n,w,work,lwork,iwork,liwork,info)
  toc = second()
  deallocate(work)
  print *,"Compute eigendecompostion", toc-tic

  !write eigenvalues normalized by trace
  open(2,file='eigenvalues',form='unformatted',status='unknown',&
       &access='direct',RECL=4)
  total = sum(w)
  do ii=n,1,-1
     write(2,rec=n-ii+1) real(w(ii)/total)
  end do
  close(2)

  !write sqrt of eigenvalues normalized by sqrt ngg
  open(2,file='sv',form='unformatted',status='unknown',&
       &access='direct',RECL=4)
  do ii=n,1,-1
     write(2,rec=n-ii+1) real(sqrt(w(ii)/dble(ngg)))
  end do
  close(2)

  !write eigenvectors multiplied by sqrt ngg
  open(2,file='eigenvectors',form='unformatted',status='unknown',&
       &access='direct',RECL=n*4)

  do ii=1,n
     write(2,rec=ii) real(c(:,n-ii+1)*sqrt(real(ngg)))
  end do
  close(2)

  !compute time-series  
  allocate(ts(n,m),stat=error)
  if  (error /= 0) then
     print *,"failed"
     stop
  end if

  tic = second()
  CALL DGEMM ('T','N', n, m, n, 1d0, c, n, dxx, n, 0d0, ts, n)
  toc = second()

  open(2,file='time-series',form='unformatted',status='unknown',&
       &access='direct',RECL=n*4)

  print *,"Compute time-series ", toc-tic

  !write time-series normalized by sqrt(ngg)
  do ii=1,m
     do j=n,1,-1
        tmp(j) = ts(n-j+1,ii)
     end do
     write(2,rec=ii) real(tmp)/sqrt(real(ngg))
  end do
  close(2)

end program driver
EOF
f90 -64 -freeform -O2 -o eofs.x prog.f90 -lcomplib.sgimath_mp 
./eofs.x