! This program reads wind on grid points, transform them back to
! spectrum space, and then computes psi, chi, vorticity, and divergence
! Author: Saha
! History: Chuang: modified to be used by ncep post and to add dynamical allocation
! April 2010   Moorthi : Updated to use less memory, speedup and implicit none
!
      implicit none
!
      integer, parameter :: komax=100
      character*200 grbwnd,indwnd
      character*200 file,psichifile
      character*10 kyr,kmth
!
      real, allocatable ::  buff(:)
      real, allocatable ::  ui(:,:,:), vi(:,:,:), uo(:,:,:), vo(:,:,:)
      real, allocatable ::  div(:,:,:),  zo(:,:,:)
      real, allocatable ::  psio(:,:,:), so(:,:,:)
!
      integer KPDS(200), KGDS(22), JPDS(200), JGDS(22)
      integer, allocatable :: npress(:)
!
      logical*1, allocatable :: lbms(:,:), lb(:)
      real    :: po(komax), th(komax), pv(komax)
!
      integer :: idim, jdim, kdim, ldim, nt, nl, kl, klb, num_threads
     &,          n, lev, klen, jcap, idrt, ntimes, k, ierr, iret, kskp
     &,          ndata, kpo, kth, kpv, maxdim
      integer :: imax, jmax
      integer :: num_parthds
      data       imax/1760/, jmax/880/
      namelist/nampgb/kpo,po,kth,th,kpv,pv,imax,jmax
!
       read(5,nampgb,iostat=ierr)
       if(ierr == 0)then
         kdim = kpo
	 allocate (npress(kdim))
	 do k=1,kdim
	   npress(k) = nint(po(k))
	 end do  
         if (npress(kdim) == 0 ) then 
           print*,'kdim, npress = ',kdim,npress
           print *,' May be you should use Pa instead of mb for'
     &,            ' pressure?'
         endif
       else
         call abort
       end if	 
       maxdim = imax*jmax
       allocate (lb(maxdim))
       allocate (buff(maxdim))
       call getenv("PGBOUT",grbwnd)
       write(*,*) "grbwnd= ",grbwnd
!
       call getenv("PGIOUT",indwnd)
       write(*,*) "indwnd= ",indwnd
!
       call getenv("psichifile",psichifile)
       write(*,*) "psichifile= ",psichifile
!
!       call getenv("YY",kyr)
!       read(kyr,'(i4)') iyr
!       write(*,*) "year= ",iyr
!
!       call getenv("MM",kmth)
!       read(kmth,'(i2)') imth
!       write(*,*) "month= ",imth
!
       call baopenr(11,grbwnd,ierr)
         if(ierr.ne.0) then
         print *,'error opening file ',grbwnd
         call abort
       endif
       call baopenr(12,indwnd,ierr)
         if(ierr.ne.0) then
         print *,'error opening file ',indwnd
         call abort
       endif
!
       call baopenwt(51,psichifile,ierr)
         if(ierr.ne.0) then
         print *,'error opening file ',psichifile
         call abort
       endif
! get grid dimensions
       jpds = -1
       jgds = -1
       CALL GETGB(11,12,maxdim,0,JPDS,JGDS,
     *            NDATA,KSKP,KPDS,KGDS,LB,buff,IRET) 
       if (iret == 0)then
         idim = kgds(2)
	 jdim = kgds(3)
	 idrt = kgds(1)
	 if(idrt == 0)then
	   jcap = (jdim-3)/2
	 else
	   jcap = jdim-1
	 end if    
	 print*,'idim,jdim idrt,jcap= ',idim,jdim,idrt,jcap
       else
         call abort	 
       end if 	     
       num_threads = num_parthds()
       ldim        = min(kdim,num_threads)
!                                       allocate arrays now
      allocate (ui(idim,jdim,ldim))
      allocate (vi(idim,jdim,ldim))
      allocate (uo(idim,jdim,ldim))
      allocate (vo(idim,jdim,ldim))
      allocate (div(idim,jdim,ldim))
      allocate (zo(idim,jdim,ldim))
      allocate (psio(idim,jdim,ldim))
      allocate (so(idim,jdim,ldim))
      allocate (lbms(idim,jdim))
!     allocate deg(jdim)
                  
!      do j=1,jdim
!      deg(j)=90.0-float(j-1)*2.5
!      enddo
!
      ntimes = (kdim-1)/ldim + 1
      do nt=1,ntimes
        klb = (nt-1)*ldim
        klen = ldim
        if (klb+klen > kdim) klen = kdim - klb
        do kl=1,klen
          nl  = kl + klb
          lev = npress(nl)
!
          N        = -1
          JPDS     = -1
          JPDS(5)  = 33
          JPDS(6)  = 100
          JPDS(7)  = lev
!         JPDS(8)  = mod(IYR-1,100) + 1
!         JPDS(9)  = IMTH
!         JPDS(21) = ((IYR-1)/100) + 1
          JGDS     = -1
          CALL GETGB(11,12,idim*jdim,N,JPDS,JGDS,NDATA,KSKP,KPDS,KGDS,
     &               LBMS,UI(1,1,kl),IRET)
          if(iret.ne.0) then
            print *,' error in GETGB for rc = ',iret,jpds
            call abort
          endif
!
          N        = -1
          JPDS     = -1
          JPDS(5)  = 34
          JPDS(6)  = 100
          JPDS(7)  = lev
!         JPDS(8)  = mod(IYR-1,100) + 1
!         JPDS(9)  = IMTH
!         JPDS(21) = ((IYR-1)/100) + 1
          JGDS     = -1
          CALL GETGB(11,12,idim*jdim,N,JPDS,JGDS,NDATA,KSKP,KPDS,KGDS,
     &               LBMS,VI(1,1,kl),IRET)
          if(iret.ne.0) then
            print *,' error in GETGB for rc = ',iret,jpds
            call abort
          endif
        enddo             ! end of kl loop
!
!       print*,'sample ui,vi=',ui(idim/2,jdim/2),vi(idim/2,jdim/2)
        call sptrunv(0,jcap,idrt,idim,jdim,idrt,idim,jdim,klen,
     &               0,0,0,0,0,0,0,0,ui(1,1,1),vi(1,1,1),
     &              .false.,uo(1,1,1),vo(1,1,1),.false.,div,zo,.true.
     &              ,psio(1,1,1),so(1,1,1))
!
!        call gridav(ui(1:idim,1:jdim,nl),idim,jdim,deg,globu)
!        call gridav(vi(1:idim,1:jdim,nl),idim,jdim,deg,globv)
!        call gridav(po(1:idim,1:jdim,nl),idim,jdim,deg,globp)
!        call gridav(so(1:idim,1:jdim,nl),idim,jdim,deg,globs)
!        print *,iyr,imth,lev,globu,globv,globp,globs
!
        do kl=1,klen
          nl  = kl + klb
          KPDS(5)  = 35
          KPDS(22) = -4
	  KPDS(7)  = npress(nl)
          CALL PUTGB(51,NDATA,KPDS,KGDS,LBMS,SO(1:idim,1:jdim,kl)
     &	  ,IRET)
          if(iret.ne.0) then
            print *,' error in PUTGB for psi for iret ',iret,
     *        (kpds(k),k=5,7)
            call abort
          endif
!
          KPDS(5)  = 36
          KPDS(22) = -4
	  KPDS(7)  = npress(nl)
          CALL PUTGB(51,NDATA,KPDS,KGDS,LBMS,PsiO(1:idim,1:jdim,kl)
     &	  ,IRET)
          if(iret.ne.0) then
            print *,' error in PUTGB for chi for iret ',iret,
     *        (kpds(k),k=5,7)
            call abort
          endif
!
        enddo               ! end of kl loop
      enddo                 ! end of nt loop
!
      stop
      end