!-----------------------------------------------------------------------------!
!           Group on Data Assimilation Development - GDAD/CPTEC/INPE          !
!-----------------------------------------------------------------------------!
!BOP
module m_utils

  ! Modules from INPAK90 lib

  USE m_inpak90
  USE m_die, only : perr

  implicit none
  public


  LOGICAL, PARAMETER :: noiqc  = .true.   ! Logical Flag to OI QC (See GSI Manual)
  INTEGER, PARAMETER :: maxobs = 2000000  ! Numero maximo de observacoes a serem lidas

  !
  ! Type of Observation to be read
  !
  INTEGER, PARAMETER  :: nvars  = 7 ! number of conventional variables
  CHARACTER(len=3), DIMENSION(nvars) :: vname = (/     &
       "sst",&
       " ps",& ! Surface Pressure
       " pw",& ! 
       " uv",& ! Wind Speed
       "  t",& ! Air Temperature
       "  q",& ! Specific Humidity
       "gps"/) ! GPS RO- Nao fique perdido

  !
  ! Levels to be Analyzed.
  ! These levels will be analyzed by layers
  !
  ! ex. 1000.00 - 1200.0
  !      900.00 -  999.9
  !        .         . 
  !        .         .
  !        .         .
  !        0.00 - 2000.0 this level will contain information about all levels
  !


  integer, parameter :: nlev = 14
  real, parameter, dimension(nlev) :: levs = (/       &
       1000.0,&
       900.0,&
       800.0,&
       700.0,&
       600.0,&
       500.0,&
       400.0,&
       300.0,&
       250.0,&
       200.0,&
       150.0,&
       100.0,&
       50.0,&
       0.0  &
       /)

! TYPE kx
!    integer          :: code
!    character(len=5) :: name
! END TYPE kx
! ! ver tabela 2 prepbr (http://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/table_2.htm)
  integer, parameter      :: nkx=87
  integer, dimension(nkx) :: kx = (/&
       111,112,120,122,126,130,131,132,&
       133,134,135,150,151,152,153,156,&
       157,158,159,164,165,174,175,180,&
       181,182,183,187,188,191,210,220,&
       221,222,223,224,227,228,229,230,&
       231,232,233,234,235,241,242,243,&
       244,245,246,247,248,249,250,251,&
       252,253,254,255,256,257,258,259,&
       280,281,282,283,284,285,286,287,&
       288,289,290,003,004,722,740,741,&
       742,743,744,745,820,042,786/)

            
  !
  ! Prefixo utilizado no arquivo de saida para identificar os arquivos de saida
  !

  character(len=3) :: NameFGS='fgs'
  character(len=3) :: NameANL='anl'

  !
  ! Informacoes que vem do namelist. Aqui alguns ficam com valores padrao
  !

  INTEGER :: nproc = 71      ! Numero de processadores utilizados pelo
                              ! sistema de assimilacao de dados
  real    :: udef   = -999.9  ! Valor Indefinido

  character(len=1024) :: DirIn ! Diretorio de entrada
  character(len=1024) :: DirOu ! Diretorio de saida

  !
  ! Arquivos de Entrada/Saida
  !

  character(len=1024) :: FInp    ! Arquivo pe GSI de entrada

  character(len=1024) :: FDatSSO ! Arquivo Binario com Distribuicao Espacial
  character(len=1024) :: FTxtSSO ! Arquivo ASCII com Distribuicao Espacial
                                 ! [somente valores de observacao e modelo]

  character(len=1024) :: FTxTOMF ! Arquivo ASCII com o valor do OMF ou OMF Medio no Globo

  character(len=1024) :: FSOC    ! Arquivo ASCII com tabela resumindo a contagem
                                 ! dos dados
  character(len=1024) :: FSOCL   ! Arquivo ASCII com tabela resumindo a contagem
                                 ! dos dados por nivel 

  character(len=1024) :: FCTL    ! Arquivo Ctl para o Arquio Binario do GrADS
                                 ! com a distribuicao Espacial

  character(len=1024) :: FMAP    ! Arquivo map para o Arquio Binario do GrADS
                                 ! com a distribuicao Espacial


  TYPE data
     character(len=3) :: name                  ! observation name
     integer          :: lim_qm                ! GSI/PreBufr quality control marker
     integer          :: nobs                  ! total # of observation
     integer          :: use  ( 0:nlev )       ! total # of used observation
     integer          :: nuse ( 0:nlev )       ! total # of unused observation
     integer          :: rej  ( 0:nlev )       ! total # of rejeited observation by GSI Quality Control
     integer          :: mon  ( 0:nlev )       ! total # of monitored observation
     integer          :: count(   2    )
     real             :: vies ( 0:nlev, 2 )
     real             :: rmse ( 0:nlev, 2 )
     character(len=8) :: ID   ( maxobs )       ! stattion ID (Name or Number)
     real             :: kx   ( maxobs )       ! observation type
     real             :: lat  ( maxobs )       ! observation latitude (degrees)
     real             :: lon  ( maxobs )       ! observation longitude (degrees)
     real             :: lev  ( maxobs )       ! observation level reference
     real             :: prs  ( maxobs )       ! observation pressure (hPa)
     real             :: time ( maxobs )       ! obs time (minutes relative to analysis time)
     real             :: pbqc ( maxobs )       ! input prepbufr qc or event mark
     real             :: iuse ( maxobs )       ! analysis usage flag (1=use, -1=monitoring ) 
     real             :: iusev( maxobs )       ! analysis usage flag ( value ) 
     real             :: robs ( maxobs )       ! observation
     real             :: rmod ( maxobs )       ! model
     real             :: diff ( maxobs )       ! obs-ges used in analysis
     real             :: sigo ( maxobs )       ! observation error      
     real             :: elevgps ( maxobs )      ! rdiagbuf(7,i) = elev !height in meters (setupref.f90)
     					       ! rdiagbuf(7,i) = tpdpres(i)-rocprof ! impact height in meters (setupbend.f90) (responsavel Ivette Banos)
     real             :: incref ( maxobs )       ! incremental refractivity (x100 %)
     real             :: hobl ( maxobs )         ! model vertical grid  (midpoint)
     real             :: errf ( maxobs )         ! final inverse observation error
     real             :: normdiff ( maxobs )     ! normalized bending angles differences
     real             :: kmelev ( maxobs )       ! elevgps !height in km
     real             :: diffpt ( maxobs )       ! obs-ges used in analysis

  END TYPE data


  INTERFACE get_value
     MODULE PROCEDURE get_real, get_m_real,&
          get_int, get_m_int,&
          get_char, get_m_char
  END INTERFACE get_value


contains


  !
  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
  !
  ! !ROUTINE: readcard
  !  \label{readcard}
  !
  ! !REVISION HISTORY:
  !  ## ### ####   J. G. de Mattos - Initial Code
  !  22 Oct 2012 - R. Mello        - Correcao na aquisição dos indices 
  !                                  dos modelos.
  !
  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
  !
  ! !INTERFACE:
  !
  SUBROUTINE readcard(istat)
    IMPLICIT NONE
    CHARACTER(len=*),PARAMETER    :: myname_="readcard"
    INTEGER                       :: I, J, iret, narg, ndom
    INTEGER, OPTIONAL,INTENT(OUT) :: istat
    CHARACTER(len=300)            :: config, formato
    LOGICAL                       :: exists
    INTEGER, ALLOCATABLE          :: tmp(:)
    INTEGER, EXTERNAL             :: iargc
    INTEGER                       :: len_strg
    INTEGER                       :: idx

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
    ! !DESCRIPTION:
    !
    !  The code in this file initializes the SCAMTec configuration management 
    !  utility.
    !  The generic model-independent runtime specifications of a SCAMTec
    !  are read by this routine. The routine also initializes the SCAMTeC
    !  log buffers, and calls some of the registries that set up the 
    !  component plugin definitions.
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
    !
    !EOP
#ifdef DEBUG
    WRITE(*,'(a72)')'!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!'
    WRITE(*,'(1x,a10,a8,a8,44x,a1)')'! INSIDE :', TRIM(myname_),' Routine','!'
    WRITE(*,'(a72)')'!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!'
#endif
    if(present(istat)) istat=0

    !         call i90_Release ( )
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    ! Load descritor file

!    narg=iargc()
!    IF( iargc() .EQ. 1)then
!       call getarg(1,config)
!    ELSE
       config='../util/conf/diag_gsi.conf'
!    END IF

    INQUIRE(FILE=trim(config),EXIST=exists)
    IF(.NOT. exists)THEN
       iret = -1
       call perr(myname_,'config file not found',iret)
       if(present(istat))istat=iret
       return
    END IF

    call i90_LoadF ( TRIM(config), iret )

    if(iret /= 0) then
       call perr(myname_,'i90_LoadF("'//trim(config)//'")',iret)
       if(present(istat))istat=iret
       return
    endif
    !
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    ! Reading time parameters
    !
    call get_value( 'Numero de Processadores:' , NProc   , iret )
    call get_value( 'Arquivo pe GSI:'          , FInp    , iret ) 
    call get_value( 'Diretorio de Entrada:'    , DirIn   , iret )
    call get_value( 'Diretorio de Saida:'      , DirOu   , iret )
    call get_value( 'Arquivo Espacial Binario:', FDatSSO , iret )
    call get_value( 'Arquivo Espacial ASCII:'  , FTxtSSO , iret )
    call get_value( 'Arquivo OMF/OMA ASCII:'   , FTxtOMF , iret )
    call get_value( 'Tabela Contagem Total:'   , FSOC    , iret )
    call get_value( 'Tabela Contagem Nivel:'   , FSOCL   , iret )

    !
    ! Gerando Nome dos Arquivos CTL e MAP para serem utilizados no GrADS.
    ! Os nomes serão baseados nos nomes passados pelo FDatSSO
    !

    len_strg = LEN_TRIM(FDatSSO)
    idx      = index(FDatSSO(len_strg-5:len_strg),'.')
    FCTL     = TRIM(FDatSSO(1:len_strg-6+idx))//'ctl'
    FMAP     = TRIM(FDatSSO(1:len_strg-6+idx))//'map'




  end subroutine readcard

  subroutine writectl(BinFName,CtlFName,MapFName,TimeFile)
    implicit none

    character(len=*), INTENT(IN) :: BinFName
    character(len=*), INTENT(IN) :: CtlFName
    character(len=*), INTENT(IN) :: MapFName
    character(len=*), INTENT(IN) :: TimeFile

    Open(Unit = 10,&
         File = TRIM(CtlFName),&
         Form = 'Formatted'&
         )

    write(10,'(2A)')'dset ^',TRIM(BinFName)
    write(10,'(2A)')'options sequential'
    write(10,'( A)')'DTYPE  station' 
    write(10,'(2A)')'STNMAP ^',TRIM(MapFName)
    write(10,'( A)')'UNDEF  -999.0'
    write(10,'( A)')'TITLE  Station Data Sample'
    write(10,'(3A)')'TDEF   1 linear ',TRIM(TimeFile),' 6hr'
    write(10,'( A)')'VARS 13'
    write(10,'( A)')'ID    1 99 Station ID'
    write(10,'( A)')'kx    1 99 observation type'
    write(10,'( A)')'rlat  1 99 observation latitude (degrees)'
    write(10,'( A)')'rlon  1 99 observation longitude (degrees)'
    write(10,'( A)')'rlev  1 99 observation level reference'
    write(10,'( A)')'prs   1 99 observation pressure (hPa)'
    write(10,'( A)')'rtime 1 99 obs time (minutes relative to analysis time)'
    write(10,'( A)')'pbqc  1 99 input prepbufr qc or event mark'
    write(10,'( A)')'iuse  1 99 analysis usage flag (1=use, -1=monitoring )'
    write(10,'( A)')'iusev 1 99 analysis usage flag ( value )'
    write(10,'( A)')'robs  1 99 observation'
    write(10,'( A)')'rmod  1 99 model'
    write(10,'( A)')'diff  1 99 OMF'
    write(10,'( A)')'ENDVARS'


  end subroutine writectl

  ! 
  ! rotinas para ler dados do namelist
  ! usam o inpack90: libmpeu

  subroutine get_real(label,value,iret)
    implicit none
    character(len=*), parameter    :: myname_="get_real"

    character(len=*),  intent(in)  :: label
    real,              intent(out) :: value
    integer,           intent(out) :: iret

    character(len=36) :: msg

    iret = 0

    call i90_label(label, iret)
    if (iret .eq. 0 )then
       value = i90_gfloat(iret)
    else
       write(msg,'(3A)')'i90_label("',trim(label),'")' 
       call perr(myname_,msg,iret)
       return
    endif

  endsubroutine get_real

  subroutine get_m_real(label,value,iret)
    implicit none
    character(len=*), parameter    :: myname_="get_real"

    character(len=*),  intent(in)  :: label
    real, dimension(:),intent(out) :: value
    integer,           intent(out) :: iret

    integer           :: i, n
    character(len=36) :: msg

    n=size(value,dim=1)

    iret = 0

    call i90_label(label, iret)
    if (iret .eq. 0 )then
       do i=1,n
          value(i) = i90_gfloat(iret)
       enddo
    else
       write(msg,'(3A)')'i90_label("',trim(label),'")' 
       call perr(myname_,msg,iret)
       return
    endif


  endsubroutine get_m_real


  subroutine get_int(label,value,iret)
    implicit none
    character(len=*), parameter    :: myname_="get_int"

    character(len=*),  intent(in)  :: label
    integer,           intent(out) :: value
    integer,           intent(out) :: iret

    character(len=36) :: msg

    iret = 0

    call i90_label(label, iret)
    if (iret .eq. 0 )then
       value = i90_gint(iret)
    else
       write(msg,'(3A)')'i90_label("',trim(label),'")' 
       call perr(myname_,msg,iret)
       return
    endif

  endsubroutine get_int

  subroutine get_m_int(label,value,iret)
    implicit none
    character(len=*), parameter    :: myname_="get_int"

    character(len=*),      intent(in)  :: label
    integer, dimension(:), intent(out) :: value
    integer,               intent(out) :: iret

    integer           :: i, n
    character(len=36) :: msg

    n=size(value)
    iret = 0

    call i90_label(label, iret)
    if (iret .eq. 0 )then
       do i=1,n
          value(i) = i90_gint(iret)
       enddo
    else
       write(msg,'(3A)')'i90_label("',trim(label),'")' 
       call perr(myname_,msg,iret)
       return
    endif

  endsubroutine get_m_int


  subroutine get_char(label,value,iret)
    implicit none
    character(len=*), parameter    :: myname_="get_char"

    character(len=*),  intent(in)  :: label
    character(len=*),  intent(out) :: value
    integer,           intent(out) :: iret


    character(len=36) :: msg

    iret = 0

    call i90_label(label, iret)
    if (iret .eq. 0 )then
       call i90_Gtoken(value,iret)
    else
       write(msg,'(3A)')'i90_label("',trim(label),'")' 
       call perr(myname_,msg,iret)
       return
    endif

  endsubroutine get_char

  subroutine get_m_char(label,value,iret)
    implicit none
    character(len=*), parameter    :: myname_="get_char"

    character(len=*),                intent(in)  :: label
    character(len=*), dimension(:),  intent(out) :: value
    integer,                         intent(out) :: iret

    integer           :: i, n
    character(len=36) :: msg

    n=size(value)

    iret = 0

    call i90_label(label, iret)
    if (iret .eq. 0 )then
       do i=1,n
          call i90_Gtoken(value(i),iret)
       enddo
    else
       write(msg,'(3A)')'i90_label("',trim(label),'")' 
       call perr(myname_,msg,iret)
       return
    endif

  endsubroutine get_m_char

end module m_utils
