!-----------------------------------------------------------------------------!
!             Modeling and Development Division - DMD/CPTEC/INPE              !
!-----------------------------------------------------------------------------!
!BOI
! !TITLE: Conversor de arquivo fct para unf
!
! !AUTHORS: João Gerd Zell de Mattos
!
! !AFFILIATION: Modeling and Development Division, CPTEC/INPE
!
! !DATE:  20 de Fevereiro de 2018
!
! !DESCRIPTION:
!
! !USES:
!     * sigio_BAMMod - Modulo contendo as rotinas necessarias para a leitura e 
!                      escrita dos arquivos do BAM
!EOI
!-----------------------------------------------------------------------------!
!

program PostAlt

#ifdef useMPI
    use MPI
#endif
    use Output, only : writefield, sendfield
    use utilmod, only: writectl, fldInfo

    use sigio_BAMMod
    use transformTools
    use LegendreTransform, only: getxMaxyMax
    use ModConstants
    use typeKinds, only: r4 => Single, r8=>Double
    use scan_MetForm, only: hmxr1, svap, RdRv

    implicit none
 
    real,  parameter   :: ps = 1000.00
    integer, parameter :: fileOutUnit = 99

!    integer, parameter :: stderr = 0
!    integer, parameter :: stdinp = 5
!    integer, parameter :: stdout = 6


    type(bamFile) :: bam, tmp

    integer :: i, j, k
    integer :: iField, iLev
    integer :: iMax, jMax, kMax, ofs
    integer :: nFields, nFieldsOut, NLevs
    integer :: nargs

    integer :: myPe, myLev
    integer :: peRoot
    integer :: nPe

    integer :: iret, ierror
    integer :: icount, iPe
    integer :: total,wrt
    integer :: jcap, imax_b
    integer :: mEnd
    type(specTrans) :: sp 

    real(r8), allocatable :: ak(:), bk(:)
    real(r4), allocatable :: clat(:)

    real(r4), allocatable :: field(:), prsl(:)
    real(r4), allocatable :: gq(:), gt(:), gtv(:), gw(:), gws(:), ges(:),gur(:) 
    real(r4), allocatable :: gurst(:),gursv(:),gurnt(:),gurnv(:)
    real(r4), allocatable :: lon(:), lat(:), lon1, londel
    real(r8), allocatable :: gu(:,:,:), gv(:,:,:)

    character(len=   2) :: first
    character(len=1024) :: arg
    character(len=1024) :: ftype
    character(len=1024) :: fileDir
    character(len=1024) :: fileUnf
    character(len=1024) :: fileOut
    character(len=1024) :: fileDirTmp
    character(len=1024) :: varName
    character(len=1024) :: pathOut, bufr
    character(len=  60) :: fieldName
    character(len=  60), allocatable :: fldNames(:), fldNamesOut(:)
    character(len=   4), allocatable :: shortNames(:)

    type(fldInfo), pointer :: rootVars => null()
    type(fldInfo), pointer :: vars => null()

    logical :: extra, calcWind, writePsLevels, calcUmrl, umrlMask
    logical :: onlyCtl, onlyOneVar
    logical :: initSpec



    !
    ! Initialize MPI
    !
#ifdef useMPI
    call MPI_INIT(ierror)
    call MPI_COMM_RANK(MPI_COMM_WORLD,MyPe,ierror)
    call MPI_COMM_SIZE(MPI_COMM_WORLD,NPe,ierror)
#else
    myPe = 0
    nPe  = 1
#endif

    peRoot = 0
    
    !
    ! reading command line options
    !
    if (myPe .eq. peRoot)then
       nargs = iargc()
       call getarg(1,first) 
       if(iargc() .eq. 0 .or. trim(first) .eq. '-h')then
          write(*,'(A)')''
          write(*,'(A)')'readBAM.x: le e converte arquivo espectral do modelo BAM'
          write(*,'(A)')'           em campos em ponto de grade.'
          write(*,'(A)')''
          write(*,'(A)')'     postAlt.x [opções] -d <dir file> -f <unf file>'
          write(*,'(A)')''
          write(*,'(A)')'     opções:'
          write(*,'(A)')'             -d : indica nome do arquivo dir'
          write(*,'(A)')'             -f : indica nome do arquivo espectral'
          write(*,'(A)')'             -t : indica o tipo de arquivo: '
          write(*,'(A)')'                    * unf: arquivo de análise (chopping ou GSI)'
          write(*,'(A)')'                    * icn: arquivo de análise gerado pelo BAM (padrão)'
          write(*,'(A)')'                    * inz: arquivo de análise gerado pelo BAM (inicializado)'
          write(*,'(A)')'                    * fct: arquivo de previsão gerado pelo BAM (padrão)'
          
          write(*,'(A)')'             -w : calcula componentes do vento'
          write(*,'(A)')'             -l : calcula níveis de pressão'
          write(*,'(A)')'             -q : calcula umidade relativa)'
          write(*,'(A)')'             -Q : calcula umidade relativa (com valores não físicos)'
          write(*,'(A)')'             -c : produz somente o arquivo ctl'
          write(*,'(A)')'             -v : indica a variável a ser extraída'
          write(*,'(A)')'             -p : indica o diretório de saída'
          write(*,'(A)')'             -T : indica o truncamento de saída'

          write(*,'(A)')'             -r : indica um arquivo Dir extra do qual serão'
          write(*,'(A)')'                  obtidos os nomes das variáveis que devem se'
          write(*,'(A)')'                  lidas do arquivo espectral. Por padrão os nomes'
          write(*,'(A)')'                  são lidos do arquivo dir especificado pela opção'
          write(*,'(A)')'                  -d <dir file>'
          write(*,'(A)')''
          stop
       endif
   
       fileDir       = ''
       fileUnf       = ''
       fileDirTmp    = ''
       ftype         = ''
       pathOut       = './'
       extra         = .false.
       onlyCtl       = .false.
       onlyOneVar    = .false.
       calcWind      = .false.
       calcUmrl      = .false.
       umrlMask      = .false.
       writePsLevels = .false.
       mEnd          = 0
       initSpec      =.true.
   
       i = 1
       do while( i<= iargc())
          call getarg(i,arg)
          select case (trim(arg))
             case ('-d')
                call getarg(i+1,arg)
                fileDir=trim(arg)
                i = i + 1
             case ('-f')
                call getarg(i+1,arg)
                fileUnf=trim(arg)
                i = i + 1
             case ('-w')
                call getarg(i+1,arg)
                fileDirTmp=trim(arg)
                extra=.true.
                i = i + 1
             case ('-t')
                call getarg(i+1,arg)
                ftype=trim(arg)
                i = i + 1
             case ('-v')
                call getarg(i+1,arg)
                varName=trim(arg)
                onlyOneVar=.true.
                i = i + 1
             case ('-p')
                call getarg(i+1,arg)
                pathOut=trim(arg)
                i = i + 1
             case ('-T')
                call getarg(i+1,arg)
                bufr=trim(arg)
                read(bufr,*)mEnd
                initSpec = .false.
                i = i + 1
             case ('-r')
                call getarg(i+1,arg)
                fileDirTmp=trim(arg)
                i = i + 1
                extra=.true.
             case ('-u')
                calcWind=.true.
             case ('-l')
                writePsLevels=.true.
             case ('-c')
                onlyCtl=.true.
             case ('-q', '-Q')
                calcUmrl = .true.
                if (trim(arg).eq.'-q') umrlMask = .true.
             case default
                print*,'WARNING: Unknow option:',trim(arg)
                print*,''
                print*,'        use postAlt.x -h para ajuda'
                print*,''
                stop
          end select
          i = i + 1
       enddo
   
       !
       ! sanity check
       !
       if (fileDir .eq. '')then
          write(*,'(A)')'ERROR: missing set dir file!'
          write(*,'(A)')'       use -d <dir file>'
          stop
       endif
   
       if (fileUnf .eq. '')then
          write(*,'(A)')'ERROR: missing set Unf (spectral) file!'
          write(*,'(A)')'       use -f <dir file>'
          stop
       endif
   
   
       if (ftype .eq. '')then
          write(*,'(A)')'WARNING: missing set file type (anl or fct) !'
          write(*,'(A)')'         setting: ftype = fct'
          ftype = 'fct'
       endif
   
       !
       ! Abrindo arquivo de saída
       !
       if (.not.onlyCtl)then
          i=index(fileUnf,'/',.true.) + 1
          fileOut = trim(pathOut)//'/'//trim(fileUnf(i:len_trim(fileUnf)))//'.gs4r'
          open(unit=fileOutUnit,file=trim(fileOut),form='unformatted')
       endif

    endif
#ifdef useMPI
!    if (myPe .eq. peRoot) print*,len(fileDir),trim(fileDir)
    call MPI_Bcast(      fileDir,   len(fileDir), MPI_CHARACTER, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(      fileUnf,   len(fileUnf), MPI_CHARACTER, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(   fileDirTmp,len(fileDirTmp), MPI_CHARACTER, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(        fType,     len(fType), MPI_CHARACTER, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(        extra,              1,   MPI_LOGICAL, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(      onlyCtl,              1,   MPI_LOGICAL, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(   onlyOneVar,              1,   MPI_LOGICAL, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(     calcWind,              1,   MPI_LOGICAL, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(     calcUmrl,              1,   MPI_LOGICAL, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(     umrlMask,              1,   MPI_LOGICAL, peRoot, MPI_COMM_WORLD, ierror)
    call MPI_Bcast(writePsLevels,              1,   MPI_LOGICAL, peRoot, MPI_COMM_WORLD, ierror)
#endif

    !
    ! Abrindo arquivo dir template, de onde são obtidas os
    ! nomes das variáveis de saída
    !

    if (extra) call tmp%open(trim(fileDirTmp),ftype='dir',initSpec=.false.)

    !
    ! Abrindo o arquivo do BAM de onde serão lidos os campos
    !

    
    call bam%open(hFile    = trim(fileDir), &
                  bFile    = trim(fileUnf), &
                  initSpec = initSpec,      &
                  ftype    = ftype,         &
                  mode     = 'r')

    jcap = bam%getOneDim('mend')
    kMax = bam%getOneDim('kmax')

    if( .not.initSpec)then
       call getxMaxyMax(mEnd, iMax, jMax)
       iMax_b = ((2*jcap+1)/iMax+1)*iMax
       call bam%initTransform(jcap, iMax_b, jMax)
    else
       mEnd = jcap
       iMax = bam%getOneDim('xmax')
       jMax = bam%getOneDim('ymax')
    endif
    
    if (extra)then
       nFields = tmp%getNFields()
       allocate(fldNames(nFields))
       call tmp%getVarNames(fldNames)
       call tmp%close()
    else
       nFields = bam%getNFields()
       allocate(fldNames(nFields))
       call bam%getVarNames(fldNames)
    endif
    !
    ! get total fields to be write
    !
    total = 0
    do iField = 1, nFields
       total = total + bam%getNLevels(trim(fldNames(iField))) 
    enddo

    if (.not.onlyCtl)then
       !
       ! Iniciando a leitura e conversão
       !
   
       if (onlyOneVar)then

          if (myPe.eq.peRoot)then

             allocate(field(imax*jmax))
             nlevs=bam%getNLevels(trim(varName))
             do j=1,nLevs
                call bam%getField(trim(varName),j,field,mEnd)
                write(fileOutUnit)field
                print*,j, minval(field),maxval(field)
             enddo
             deallocate(field)

             !
             ! include name to output ctl
             !
             allocate(vars)
             vars%fName  = varName
             vars%nLevs  = nLevs

             call writeCtl(bam, mEnd, pathOut, vars)
             call bam%close()
#ifdef useMPI
             call MPI_FINALIZE(ierror)
             write(stdout,'(A,1x,I4)')'Program ends', ierror
#endif
          endif

          stop 0
       endif
   

       allocate(field(imax*jmax))

       icount = 1 
       do iField = 1, nFields

          nLevs = bam%getNLevels(trim(fldNames(iField)))

          !
          ! include field name to output ctl
          !

          if(.not.associated(rootVars))then
             allocate(rootVars)
             vars => rootVars
          else
             allocate(vars%next)
             vars => vars%next
          endif
          
          vars%fName  = trim(fldNames(iField))
          vars%nLevs  = nLevs

          
          do iLev = 1, nLevs

             iPe = icount - 1

             if (myPe .eq. iPe)then
                call bam%getField(trim(fldNames(iField)),iLev,field,mEnd)
                myLev     = iLev
                fieldName = fldNames(iField)
                write(*,'(A40,2I6.1,2F16.9)')fldNames(iField), myPe, iLev, minval(field),maxval(field)

                if (myPe .eq.peRoot)then
                   wrt = min(nPe,total)
                   call writefield(fileOutUnit, field, wrt, peRoot)
                   total = total - nPe
                else
                   call sendfield(field,peRoot)
                endif

             endif

             if (icount .eq. nPe ) icount = 0

             icount  = icount + 1
          enddo
       enddo
       deallocate(field)
       if (writePsLevels)then

          !
          ! include name to output ctl
          !

          if(.not.associated(rootVars))then
             allocate(rootVars)
             vars => rootVars
          else
             allocate(vars%next)
             vars => vars%next
          endif

          vars%fName  = 'PRESSURE LEVELS'
          vars%nLevs  = kMax

          if (myPe .eq. peRoot)then
             ! get pressure levels
             allocate(field(imax*jmax))
             call bam%getField('LN SURFACE PRESSURE',1,field,mEnd,iret)
             if (iret.eq.0)then
                allocate(ak(kmax+1))
                allocate(bk(kmax+1))
                allocate(prsl(imax*jmax))
      
                call BAM%GetVerticalCoord('ak', ak)
                call BAM%GetVerticalCoord('bk', bk)

                if (bam%isHybrid.and.ftype.ne.'unf')then
                   field = exp(field)*1e-3
                else
                   field = exp(field)
                endif

                ak    = ak*1e-3
                
                do k=1,kmax
                   do i=1,imax*jmax
                      prsl(i) = ak(k) + bk(k)*field(i)
                   enddo
                   write(*,'(A40,2I6.1,2F16.9)')'PRESSURE LEVLS', myPe, iLev, minval(prsl),maxval(prsl)
                   call writefield(fileOutUnit, prsl, 1, peRoot)
                enddo
      
                deallocate(ak)
                deallocate(bk)
                deallocate(prsl)
                deallocate(field)
                
             else
                print*,'Surface pressure field not found, skip ....'
             endif
          endif
       end if

 
       if (calcWind)then

          !-----------------------------------------!
          ! include name to output ctl
          !

          if(.not.associated(rootVars))then
             allocate(rootVars)
             vars => rootVars
          else
             allocate(vars%next)
             vars => vars%next
          endif

          vars%fName  = 'ZONAL WIND (U)'
          vars%nLevs  = kMax
          allocate(vars%next)
          vars => vars%next
          vars%fName  = 'MERIDIONAL WIND (V)'
          vars%nLevs  = kMax
          !
          !-----------------------------------------!
          !

          allocate(clat(JMax))
          call BAM%GetWCoord('clat', clat, iret)
          if(iret.ne.0) stop
      
          allocate(gu(imax,jmax,kmax))
          allocate(gv(imax,jmax,kmax))
          total  = kMax * 2
          icount = 1
          do k = 1, kmax
             iPe = icount - 1
             if (myPe .eq. iPe)then
                print*,'getting wind for lev:',k
                call bam%getUV(k, gu(:,:,k), gv(:,:,k))
                do j=1,jmax
                   gu(1:imax,j,k) = gu(1:imax,j,k)/clat(j)
                   gv(1:imax,j,k) = gv(1:imax,j,k)/clat(j)
                enddo
                write(*,'(A40,4F16.9)')'GU<min,max>, GV<min,max>',minval(gu(:,:,k)),maxval(gu(:,:,k)),minval(gv(:,:,k)),maxval(gv(:,:,k))
             endif
             if(icount .eq. nPe)icount = 0
             icount = icount + 1
          enddo

          total  = kMax
          icount = 1
          do k=1,kmax
            iPe = icount - 1
            if (myPe .eq. iPe)then
                allocate(field(iMax*jMax))
                if (myPe .eq. peRoot)then
                   wrt = min(nPe,total)
                   field = real(reshape(gu(:,:,k),[iMax*jMax]),r4)
                   call writefield(fileOutUnit, field, wrt, peRoot, mpiTag=101)
                   total = total - nPe
                else
                   field = real(reshape(gu(:,:,k),[iMax*jMax]),r4)
                   call sendfield(field,peRoot, mpiTag=101)
                endif
                deallocate(field)
            endif
            if (icount .eq. nPe) icount = 0
            icount = icount + 1
          enddo

          total  = kMax
          icount = 1
          do k=1,kmax
            iPe = icount - 1
            if (myPe .eq. iPe)then
                allocate(field(iMax*jMax))
                if (myPe .eq. peRoot)then
                   wrt = min(nPe,total)
                   field = real(reshape(gv(:,:,k),[iMax*jMax]),r4)
                   call writefield(fileOutUnit, field, wrt, peRoot, mpiTag=102)
                   total = total - nPe
                else
                   field = real(reshape(gv(:,:,k),[iMax*jMax]),r4)
                   call sendfield(field,peRoot, mpiTag=102)
                endif
                deallocate(field)
            endif
            if (icount .eq. nPe) icount = 0
            icount = icount + 1
          enddo

            
          deallocate(clat)
          deallocate(gu)
          deallocate(gv)
       endif
    endif

    if (calcUmrl)then

          if(.not.associated(rootVars))then
             allocate(rootVars)
             vars => rootVars
          else
             allocate(vars%next)
             vars => vars%next
          endif

          vars%fName  = 'RELATIVE HUMIDITY'
          vars%nLevs  = kMax
          if(.not.umrlMask)then
             allocate(vars%next)
             vars => vars%next            
             vars%fName  = 'TOTAL COUNT OF SATURATED SPECIFIC HUMIDITY'
             vars%nLevs  = 1
             allocate(vars%next)
             vars => vars%next            
             vars%fName  = 'VERTICAL SUM OF SATURATED SPEFIFIC HUMIDITY'
             vars%nLevs  = 1
             allocate(vars%next)
             vars => vars%next            
             vars%fName  = 'TOTAL COUNT OF NEGATIVE SPECIFIC HUMIDITY'
             vars%nLevs  = 1
             allocate(vars%next)
             vars => vars%next            
             vars%fName  = 'VERTICAL SUM OF NEGATIVE SPECIFIC HUMIDITY'
             vars%nLevs  = 1
          endif

          allocate(gq(imax*jmax))
          allocate(gt(imax*jmax))
          allocate(gtv(imax*jmax))
          allocate(gw(imax*jmax))
          allocate(gws(imax*jmax))
          allocate(ges(imax*jmax))
          allocate(field(imax*jmax))
          allocate(prsl(imax*jmax))
          allocate(gur(imax*jmax))
          if(.not.umrlMask)then
             allocate(gurst(imax*jmax))
             allocate(gursv(imax*jmax))
             allocate(gurnt(imax*jmax))
             allocate(gurnv(imax*jmax))
          endif
          allocate(ak(kmax+1))
          allocate(bk(kmax+1))

          if(myPe .eq. peRoot)then
             call bam%getField('LN SURFACE PRESSURE', 1, field, mEnd,iret)
             call BAM%GetVerticalCoord('ak', ak)
             call BAM%GetVerticalCoord('bk', bk)
      
                if (bam%isHybrid.and.ftype.ne.'unf')then
                   field = exp(field)*1e-3
                else
                   field = exp(field)
                endif
                ak    = ak*1e-3
          endif

#ifdef useMPI
          call MPI_Bcast (field, imax*jmax, MPI_FLOAT, peRoot, MPI_COMM_WORLD, ierror)
          call MPI_Bcast (   ak,    kmax+1, MPI_FLOAT, peRoot, MPI_COMM_WORLD, ierror)
          call MPI_Bcast (   bk,    kmax+1, MPI_FLOAT, peRoot, MPI_COMM_WORLD, ierror)
#endif

          gurst = 0.0
          gursv = 0.0
          gurnt = 0.0
          gurnv = 0.0

          total = kMax
          icount = 1
          do k = 1, kmax
             iPe = icount - 1
             if (myPe .eq. iPe)then
                print*,'getting umrl for lev:',k
                call bam%getField('SPECIFIC HUMIDITY',k, gq, mEnd,iret)
                call bam%getField('VIRTUAL TEMPERATURE',k, gtv, mEnd,iret)
                gw   =  hmxr1(gq)
                gt   = gtv / (( 1 + 0.61 * ( gq / ( 1 - gq ) ) ) )
                ges  = svap(gt-273.16)
                prsl = (ak(k) + bk(k)*field)*1e3 ! convert to Pah
                gws  = (RdRv*ges) / (prsl - ges);
                gur  = (gw/gws)*100.0

                if (umrlMask)then
                   where(gur > 100.0) gur = 100.0
                   where(gur < 0.0) gur = 0.0
                else
#ifdef useMPI
                   where(gur > 100.0)
                       gurst = 1.0
                       gursv = gur
                   endwhere

                   where(gur < 0.0)
                       gurnt = 1.0
                       gurnv = gur
                   endwhere

#else
                   where(gur > 100.0)
                      gurst = gurst + 1.0
                      gursv = gursv + gur
                   end where

                   where(gur < 0.0)
                      gurnt = gurnt + 1.0
                      gurnv = gurnv + gur
                   end where
#endif          
                endif

                !write gur
                if (myPe .eq. peRoot)then
                   wrt = min(nPe,total)
                   call writefield(fileOutUnit, gur, wrt, peRoot, mpiTag=101)
                   total = total - nPe
                else
                   call sendfield(gur,peRoot, mpiTag=101)
                endif
             endif
             if(icount .eq. nPe)icount = 0
             icount = icount + 1
          enddo
          
          if (myPe .eq. peRoot.and..not.umrlMask)then
             call writefield(fileOutUnit, gurst, wrt, peRoot, mpiTag=101)
             call writefield(fileOutUnit, gursv, wrt, peRoot, mpiTag=101)
             call writefield(fileOutUnit, gurnt, wrt, peRoot, mpiTag=101)
             call writefield(fileOutUnit, gurnv, wrt, peRoot, mpiTag=101)
          endif


          deallocate(gq)
          deallocate(gt)
          deallocate(gtv)
          deallocate(gw)
          deallocate(gws)
          deallocate(ges)
          deallocate(field)
          deallocate(prsl)
          deallocate(gur)
          if(.not.umrlMask)then
             deallocate(gurst)
             deallocate(gursv)
             deallocate(gurnt)
             deallocate(gurnv)
          endif
          deallocate(ak)
          deallocate(bk)
   
       
    endif
    

    if(myPe .eq. peRoot)then
       call writeCtl(bam, mEnd, pathOut, rootVars)
    endif
    
    call bam%close()
#ifdef useMPI
    call MPI_FINALIZE(ierror)
    if(MyPe.eq.peRoot)    write(stdout,'(A,1x,I4)')'Program ends', ierror
#endif
end program