      module gfs_transforms
!  GFS_TRANSFORMS - this module gathers together all the spectral transform
!  code from the 2001/2004 GFS.
!  Adapted from 2001 GFS code by Eric Kostelich.
!  Summary of changes:
!  (1) converted from fixed to free-form source
!  (2) explicitly declares double-precision arithmetic where needed.
!  (3) fixes illegal aliases and eliminates deprecated calls to
!      Schwartztrauber FFT routines.
!  (4) USEs the appropriate modules instead of include files.
!  (5) Avoids implicit array rank remapping except in calls to FFT routines.
!  (6) All dummy arguments have explicit INTENT attributes.
!  (7) IMPLICIT NONE added throughout with explicit declarations added as
!      necessary.
!
      use resolution
      use precision,only:D
      use constant,only:REARTH
      use gfs_init
      implicit none
      contains
!
!***********************************************************************
!
      SUBROUTINE FLN2FOUR(num,latd,nvar,FLN,GRID,PLN, &
                          nvect,ivect,FTOP,PLNTOP,nlats,nlons)
!
!***********************************************************************
!
      use resolution
      use precision,only:D
      implicit none
      integer,intent(in)::num,latd,nvar,nvect,ivect(nvect),nlats,nlons(latd)
      real(D),intent(in)::fln(LNTD,2*num),pln(LNTD,latd)
      real(D),intent(out)::GRID(LONFD,2,nvar,latd)
      real(D),intent(in):: FTOP(0:JCAP,2,nvect),PLNTOP(0:JCAP,latd)
!
!  Local variables
!
      integer l,ireal,iimag,j,lmax,i,iv,latmin
      real(D) a(latd,2*num),b(latd,2*num),ui,ur
      integer ii,ir,lenp,leni,indp,indi,i1,i2,n
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
      DO l=0,JCAP-1
        DO j=nlats,1,-1
          IF ((nlons(j)/2) .ge. l) latmin=j
        ENDDO
        n=nlats-latmin+1
        lenp=(JCAP-l+2)/2
        leni=(JCAP-l+1)/2
        indp=JCAP*l+l*(1-l)/2+1+l
        indi=indp+lenp

        CALL DGEMM('t','n',n,2*num,lenp,1._D,PLN(indp,latmin),LNTD, &
                   FLN(indp,1),LNTD,0._D,a(latmin,1),latd)
        CALL DGEMM('t','n',n,2*num,leni,1._D,PLN(indi,latmin),LNTD, &
                   FLN(indi,1),LNTD,0._D,b(latmin,1),latd)

        ir=2*l+1
        ii=2*l+2
        DO i=1,num
          i1=2*i-1
          i2=2*i
          DO j=latmin,nlats
            grid(ir,1,i,j)=a(j,i1)+b(j,i1)
            grid(ir,2,i,j)=a(j,i1)-b(j,i1)
            grid(ii,1,i,j)=a(j,i2)+b(j,i2)
            grid(ii,2,i,j)=a(j,i2)-b(j,i2)
          ENDDO
        ENDDO
      ENDDO

!     l=JCAP
!     ------
      ir=2*JCAP+1
      ii=2*JCAP+2
      DO j=nlats,1,-1
        IF ((nlons(j)/2) .ge. JCAP) latmin=j
      ENDDO
      DO j=latmin,nlats
        DO i=1,num
          i1=2*i-1
          i2=2*i
          grid(ir,1,i,j)=PLN(lnt,j)*FLN(lnt,i1)
          grid(ir,2,i,j)=PLN(lnt,j)*FLN(lnt,i1)
          grid(ii,1,i,j)=PLN(lnt,j)*FLN(lnt,i2)
          grid(ii,2,i,j)=PLN(lnt,j)*FLN(lnt,i2)
        ENDDO
      ENDDO

      DO j=1,nlats
        lmax=min(JCAP,nlons(j)/2)
        DO l=0,lmax
          ireal=2*l+1
          iimag=2*l+2
          DO i=1,nvect
            iv=ivect(i)
            ur=FTOP(l,1,i)*PLNTOP(l,j)
            ui=FTOP(l,2,i)*PLNTOP(l,j)
            grid(ireal,1,iv,j)=grid(ireal,1,iv,j)+ur
            grid(iimag,1,iv,j)=grid(iimag,1,iv,j)+ui
            IF (mod(JCAP1-l,2).eq.0) THEN
              grid(ireal,2,iv,j)=grid(ireal,2,iv,j)+ur
              grid(iimag,2,iv,j)=grid(iimag,2,iv,j)+ui
            ELSE
              grid(ireal,2,iv,j)=grid(ireal,2,iv,j)-ur
              grid(iimag,2,iv,j)=grid(iimag,2,iv,j)-ui
            ENDIF
          ENDDO
        ENDDO
      ENDDO

      RETURN
      END SUBROUTINE FLN2FOUR
!
!***********************************************************************
!
      subroutine four2grid(num,latd,nvar,grid,nlats,nlons)
!
!***********************************************************************
!
      implicit none
      integer,intent(in)::num,latd,nvar,nlats,nlons(latd)
      real(D),intent(inout)::grid(LONFD,2,nvar,latd)
!
!  Local variables
!
      real(D) scale
      real(D),allocatable::table(:),temp(:,:,:,:)
      integer::ifac(15)
      integer j,ldx,ldy,maxlons
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
      scale=1.0
      ldx=lonfd/2
      ldy=lonfd
      maxlons=maxval(nlons)
      allocate(table(maxlons))  ! per FFTPACK docs
      allocate(temp(LONFD,2,nvar,latd))   ! so we don't exhaust stack
      temp=grid  ! original aliasing is illegal
      do j=1,nlats
        call rffti1(nlons(j),table,ifac)
        call dcrft(temp(:,:,:,j),ldx,grid(:,:,:,j),ldy,nlons(j),2*num,scale, &
           table,ifac)
      enddo
      return
      end subroutine four2grid

!***********************************************************************
!  Adapted from 2001 GFS code by Eric Kostelich.
!  Summary of changes:
!  (1) from fixed to free-form source
!  (2) explicitly declares double-precision arithmetic where needed.
!  (3) fixes illegal aliases and eliminates deprecated calls to
!      Schwartztrauber FFT routines.
!
!***********************************************************************
!
      SUBROUTINE GRID2FOUR(num,latd,nvar,grid,nlats,wgt,nlons)
!
!***********************************************************************
!
      implicit none
      integer,intent(in)::num,latd,nvar,nlats,nlons(latd)
      real(D),intent(inout)::grid(LONFD,2,nvar,latd)
      real(D),intent(in)::wgt(latd)
!
!  Local variables
!
      real(D)::scale
      real(D),allocatable::table(:),temp(:,:,:,:)
      integer::ifac(15)
      integer j,ldx,ldy,maxlons
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
      ldx=lonfd
      ldy=lonfd/2
      maxlons=maxval(nlons)
      allocate(temp(LONFD,2,nvar,latd))  ! so we don't exhaust stack
      allocate(table(maxlons))       ! per FFTPACK docs
      temp=grid                      ! aliasing in original code is illegal
      do j=1,nlats
        call rffti1(nlons(j),table,ifac)
        scale=wgt(j)/nlons(j)
        call drcft(temp(:,:,:,j),ldx,grid(:,:,:,j),ldy,nlons(j),2*num, &
                   scale,table,ifac)
      enddo
      return
      end subroutine grid2four
!
!***********************************************************************
!
      SUBROUTINE UVSUMS(latd,DUDT,DVDT,UTOP,VTOP,PLNTOP, &
                        nlats,nlons)
!
!***********************************************************************
!
      implicit none
      integer,intent(in)::latd,nlats,nlons(latd)
      real(D),intent(in)::DUDT(LONFD,2,latd),DVDT(LONFD,2,latd), &
           PLNTOP(0:JCAP,latd)
      real(D),intent(inout):: UTOP(0:JCAP,2),VTOP(0:JCAP,2)
!
!  Local variables
!
      integer lev,ireal,iimag,lmax,j
      real(D) url,uil,vrl,vil
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

      DO j=1,nlats
        lmax=min(JCAP,nlons(j)/2)
        DO lev=0,lmax
          ireal=2*lev+1
          iimag=2*lev+2
          IF (mod(JCAP1-lev,2).eq.0) THEN
            url=DUDT(ireal,1,j)+DUDT(ireal,2,j)
            uil=DUDT(iimag,1,j)+DUDT(iimag,2,j)
            vrl=DVDT(ireal,1,j)+DVDT(ireal,2,j)
            vil=DVDT(iimag,1,j)+DVDT(iimag,2,j)
          ELSE
            url=DUDT(ireal,1,j)-DUDT(ireal,2,j)
            uil=DUDT(iimag,1,j)-DUDT(iimag,2,j)
            vrl=DVDT(ireal,1,j)-DVDT(ireal,2,j)
            vil=DVDT(iimag,1,j)-DVDT(iimag,2,j)
          ENDIF
          UTOP(lev,1)=UTOP(lev,1)+url*PLNTOP(lev,j)
          UTOP(lev,2)=UTOP(lev,2)+uil*PLNTOP(lev,j)
          VTOP(lev,1)=VTOP(lev,1)+vrl*PLNTOP(lev,j)
          VTOP(lev,2)=VTOP(lev,2)+vil*PLNTOP(lev,j)
        ENDDO
      ENDDO
      RETURN
      END SUBROUTINE UVSUMS
!
!***********************************************************************
!
      SUBROUTINE FOUR2FLN(num,latd,nvar,GRID,FLN,PLN,nlats,nlons)
!
!***********************************************************************
!
      implicit none
      integer,intent(in)::latd,num,nvar,nlats,nlons(latd)
      real(D),intent(out)::FLN(LNTD,2*num)
      real(D),intent(in)::GRID(LONFD,2,nvar,latd),PLN(LNTD,latd)
      integer k,lenp,leni,indp,indi,ii,ir,n,latmin,l,i,j
      real(D) a(latd,2*num),b(latd,2*num)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
      DO l=0,JCAP-1
        DO j=nlats,1,-1
          IF ((nlons(j)/2) .ge. l) latmin=j
        ENDDO
        n=nlats-latmin+1
        ir=2*l+1
        ii=2*l+2
        lenp=(JCAP-l+2)/2
        leni=(JCAP-l+1)/2
        indp=JCAP*l+l*(1-l)/2+1+l
        indi=indp+lenp

        DO j=latmin,nlats
          DO k=1,num
            i=2*k-1
            a(j,i  )=grid(ir,1,k,j)+grid(ir,2,k,j)
            b(j,i  )=grid(ir,1,k,j)-grid(ir,2,k,j)
            a(j,i+1)=grid(ii,1,k,j)+grid(ii,2,k,j)
            b(j,i+1)=grid(ii,1,k,j)-grid(ii,2,k,j)
          ENDDO
        ENDDO

        CALL DGEMM('n','n',lenp,2*num,n,1._D,PLN(indp,latmin),LNTD, &
                   a(latmin,1),latd,0._D,FLN(indp,1),LNTD)
        CALL DGEMM('n','n',leni,2*num,n,1._D,PLN(indi,latmin),LNTD, &
                   b(latmin,1),latd,0._D,FLN(indi,1),LNTD)

      ENDDO

!     l=JCAP
!     ------
      ir=2*JCAP+1
      ii=2*JCAP+2
      DO j=nlats,1,-1
        IF ((nlons(j)/2) .ge. JCAP) latmin=j
      ENDDO
      DO j=latmin,nlats
        DO k=1,num
          i=2*k-1
          a(j,i  )=grid(ir,1,k,j)+grid(ir,2,k,j)
          a(j,i+1)=grid(ii,1,k,j)+grid(ii,2,k,j)
        ENDDO
      ENDDO
      n=nlats-latmin+1
      CALL DGEMM('n','n',1,2*num,n,1._D,PLN(lnt,latmin),LNTD, &
                 a(latmin,1),latd,0._D,FLN(lnt,1),LNTD)

      CALL FromEvenOdd(2*num,fln)

      RETURN
      END SUBROUTINE FOUR2FLN
!
!***********************************************************************
!
      SUBROUTINE DZ2UV(Div,Z,U,V,UTOP,VTOP)
!
!***********************************************************************
!
      implicit none

      real(D),intent(in)::Div(LNTD,2),Z(LNTD,2)
      real(D),intent(out)::U(LNTD,2),V(LNTD,2),UTOP(0:JCAP,2),VTOP(0:JCAP,2)
!
!  Local variables
!
      real(D) epsan(LNTD),alnnp(LNTD),epsantop(0:JCAP)
      integer first,n,l,ind,i

      SAVE first,epsan,epsantop,alnnp
      DATA first/1/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
      IF  (first.eq.1) THEN
        ind=1
        DO l=0,JCAP
          epsan(ind)=0.
          ind=ind+JCAP1-l
        ENDDO
        ind=0
        DO l=0,JCAP
          ind=ind+1
          DO n=l+1,JCAP
            ind=ind+1
            epsan(ind)=sqrt( (n*n-l*l)/(4._D*n*n-1._D) ) * REARTH / n
          ENDDO
          epsantop(l)=sqrt((JCAP1*JCAP1-l*l)/(4.*JCAP1*JCAP1-1._D)) &
              * REARTH/JCAP1
        ENDDO

        DO i=1,JCAP1
          alnnp(i)=0.
        ENDDO
        ind=JCAP1
        DO l=1,JCAP
        DO n=l,JCAP
          ind=ind+1
          alnnp(ind)=l/(n*(n+1.)) * REARTH 
        ENDDO
        ENDDO

        first=0
      ENDIF

      DO i=1,LNT
        U(i,1)= alnnp(i)*Div(i,2)
        U(i,2)=-alnnp(i)*Div(i,1)
        V(i,1)= alnnp(i)*Z(i,2)
        V(i,2)=-alnnp(i)*Z(i,1)
      ENDDO

      ind=0
      DO l=0,JCAP
        ind=ind+1
        DO n=l+1,JCAP
          ind=ind+1
          U(ind,1)=U(ind,1)-epsan(ind)*Z(ind-1,1)
          U(ind,2)=U(ind,2)-epsan(ind)*Z(ind-1,2)
          V(ind,1)=V(ind,1)+epsan(ind)*Div(ind-1,1)
          V(ind,2)=V(ind,2)+epsan(ind)*Div(ind-1,2)
        ENDDO
      ENDDO

      ind=0
      DO l=0,JCAP
        DO n=l,JCAP-1
          ind=ind+1
          U(ind,1)=U(ind,1)+epsan(ind+1)*Z(ind+1,1)
          U(ind,2)=U(ind,2)+epsan(ind+1)*Z(ind+1,2)
          V(ind,1)=V(ind,1)-epsan(ind+1)*Div(ind+1,1)
          V(ind,2)=V(ind,2)-epsan(ind+1)*Div(ind+1,2)
        ENDDO
        ind=ind+1
      ENDDO

!     The case n=JCAP
      ind=0
      DO l=0,JCAP
        ind=ind+(JCAP1-l)
        UTOP(l,1)=-epsantop(l)*Z(ind,1)
        UTOP(l,2)=-epsantop(l)*Z(ind,2)
        VTOP(l,1)=+epsantop(l)*Div(ind,1)
        VTOP(l,2)=+epsantop(l)*Div(ind,2)
      ENDDO

      RETURN
      END SUBROUTINE DZ2UV
!
!***********************************************************************
!
      SUBROUTINE UV2DZ(U,V,UTOP,VTOP,Div,Z)
!
!***********************************************************************
!
      implicit none

      real(D),intent(out)::Div(LNTD,2),Z(LNTD,2)
      real(D),intent(in)::U(LNTD,2),V(LNTD,2),UTOP(0:JCAP,2),VTOP(0:JCAP,2)
!
!  Local variables
!
      real(D) epsl(LNTD),epstop(0:JCAP),ra,lra
      integer first,n,l,ind

      SAVE first,epsl,epstop
      DATA first/1/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
      IF  (first.eq.1) THEN
        ind=1
        DO l=0,JCAP
          epsl(ind)=0.
          ind=ind+JCAP1-l
        ENDDO
        ind=0
        DO l=0,JCAP
          ind=ind+1
          DO n=l+1,JCAP
            ind=ind+1
            epsl(ind)=sqrt( (n*n-l*l)/(4.*n*n-1.) )
          ENDDO
          epstop(l)=sqrt((JCAP1*JCAP1-l*l)/(4.*JCAP1*JCAP1-1.))
        ENDDO
        first=0
      ENDIF

      ra = 1./REARTH

      ind=0
!     For l=JCAP-1 and l=JCAP the loop is empty
      DO l=0,JCAP-2
        ind=ind+1
        DO n=l+1,JCAP-1
          ind=ind+1
          Z(ind,1)=( -l*V(ind,2) - n*epsl(ind+1)*U(ind+1,1) &
                     + (n+1.)*epsl(ind)*U(ind-1,1) )*ra
          Z(ind,2)=(  l*V(ind,1) - n*epsl(ind+1)*U(ind+1,2) &
                     + (n+1.)*epsl(ind)*U(ind-1,2) )*ra
          Div(ind,1)=( -l*U(ind,2) + n*epsl(ind+1)*V(ind+1,1) &
                     - (n+1.)*epsl(ind)*V(ind-1,1) )*ra
          Div(ind,2)=(  l*U(ind,1) + n*epsl(ind+1)*V(ind+1,2) &
                     - (n+1.)*epsl(ind)*V(ind-1,2) )*ra
        ENDDO
        ind=ind+1
      ENDDO

!     The case n=l
      Z(1,1)=0.
      Z(1,2)=0.
      Div(1,1)=0.
      Div(1,2)=0.
      ind=1
      DO l=1,JCAP-1
        ind=ind+(JCAP2-l)
        lra=l*ra
        Z(ind,1)=( -V(ind,2) - epsl(ind+1)*U(ind+1,1) )*lra
        Z(ind,2)=(  V(ind,1) - epsl(ind+1)*U(ind+1,2) )*lra
        Div(ind,1)=( -U(ind,2) + epsl(ind+1)*V(ind+1,1) )*lra
        Div(ind,2)=(  U(ind,1) + epsl(ind+1)*V(ind+1,2) )*lra
      ENDDO

!     The case n=JCAP
      ind=0
      DO l=0,JCAP
        ind=ind+(JCAP1-l)
        Z(ind,1)=( -l*V(ind,2) - JCAP*epstop(l)*UTOP(l,1) &
                    + JCAP1*epsl(ind)*U(ind-1,1) )*ra
        Z(ind,2)=(  l*V(ind,1) - JCAP*epstop(l)*UTOP(l,2) &
                    + JCAP1*epsl(ind)*U(ind-1,2) )*ra
        Div(ind,1)=( -l*U(ind,2) + JCAP*epstop(l)*VTOP(l,1) &
                    - JCAP1*epsl(ind)*V(ind-1,1) )*ra
        Div(ind,2)=(  l*U(ind,1) + JCAP*epstop(l)*VTOP(l,2) &
                    - JCAP1*epsl(ind)*V(ind-1,2) )*ra
      ENDDO

      RETURN
      END SUBROUTINE UV2DZ
      end module gfs_transforms
