module pumamod
  implicit none
#ifndef TRUNC
#define TRUNC 21
#endif
#ifndef PARTITION
#define PARTITION 1
#endif
#if TRUNC == 21
  integer, parameter :: NLAT=32
#elif TRUNC == 42
  integer, parameter :: NLAT=64
#elif TRUNC == 85
  integer, parameter :: NLAT=128
#endif
  integer, parameter :: NTRU=TRUNC, NTP1=NTRU+1, NCSP=(NTRU+1)*(NTRU+2)/2, NESP=2*NCSP
  integer, parameter :: NLON=2*NLAT, NLPP=NLAT/PARTITION, NHOR=NLON*NLPP, NLEV=10, NROOT=0
  integer :: nfilter=0, mypid=1, nud=6, ngptfilter=1, nspvfilter=1, nfilterexp=4, nqspec=1
  real :: gwd(NLAT), sid(NLAT), plavor=0.17, landhoskn0=15.0, filterkappa=12.0
  real :: sd(NESP,NLEV), st(NESP,NLEV), sz(NESP,NLEV), sq(NESP,NLEV), sp(NESP)
  real :: gd(NHOR,NLEV), gt(NHOR,NLEV), gz(NHOR,NLEV), gu(NHOR,NLEV), gv(NHOR,NLEV)
  real :: gq(NHOR,NLEV), gp(NHOR), gpj(NHOR)
end module pumamod

program transform_benchmark
  use pumamod
  use legmod
  implicit none
  real :: fc(2,NLON/2,NLPP), result(2,NCSP)
  real :: pd(2,NCSP,NLEV), pz(2,NCSP,NLEV), savedpz(2,NCSP,NLEV)
  real :: pu(2,NLON/2,NLPP,NLEV), pv(2,NLON/2,NLPP,NLEV)
  real :: start, finish, fcseconds, dvseconds
  integer :: i,j,k,v,r,repeats,scenario,iterations
  character(32) :: arg
  call get_command_argument(1,arg)
  read(arg,*) repeats
  do k=1,NLPP
    do j=1,NCSP
      qi(j,k)=sin(real(j+3*k))*0.03
      qc(j,k)=cos(real(2*j+k))*0.01
      qu(j,k)=sin(real(3*j+2*k))*0.02
      qv(j,k)=cos(real(j+5*k))*0.02
    enddo
  enddo
  open(20,file='output.bin',access='stream',form='unformatted',status='replace')
  do scenario=1,3
    do k=1,NLPP
      do j=1,NLON/2
        do i=1,2
          fc(i,j,k)=sin(real(17*i+7*j+3*k))*0.01
        enddo
      enddo
    enddo
    do v=1,NLEV
      do k=1,NCSP
        do i=1,2
          pd(i,k,v)=cos(real(i+5*k+v))*0.02
          pz(i,k,v)=sin(real(3*i+11*k+2*v))*0.03
        enddo
      enddo
    enddo
    skgpsp=1.0
    skspgp=1.0
    plavor=0.17
    if (scenario==2) then
      fc=0.0
      pd=0.0
      pz=0.0
      plavor=0.0
    else if (scenario==3) then
      do i=1,NTP1
        skgpsp(i)=exp(-real(i-1)/real(NTP1))
        skspgp(i)=exp(-2.0*real(i-1)/real(NTP1))
      enddo
      plavor=-0.31
    endif
    savedpz=pz
    call fc2sp(fc,result)
    call dv2uv(pd,pz,pu,pv)
    if (any(pz/=savedpz)) error stop 'dv2uv changed its vorticity input'
    write(20) result,pu,pv,pz
    if (scenario==1) then
      iterations=0
      call cpu_time(start)
      do
        do r=1,repeats
          call fc2sp(fc,result)
        enddo
        iterations=iterations+repeats
        call cpu_time(finish)
        if (finish-start>=0.3) exit
      enddo
      fcseconds=(finish-start)/iterations
      iterations=0
      call cpu_time(start)
      do
        do r=1,repeats
          call dv2uv(pd,pz,pu,pv)
        enddo
        iterations=iterations+repeats
        call cpu_time(finish)
        if (finish-start>=0.3) exit
      enddo
      dvseconds=(finish-start)/iterations
      print '(2(es24.16,1x))',fcseconds,dvseconds
    endif
  enddo
  close(20)
end program transform_benchmark