!-----------------------------------------------------------------------------!
!           Group on Data Assimilation Development - GDAD/CPTEC/INPE          !
!-----------------------------------------------------------------------------!
!BOP
!
! ! TITLE: GSI Conventional Diagnostic Documentation \\ Version 0.01
!
! !AUTHORS: Jo\~a Gerd Z. de Mattos
!
! !AFFILIATION: Group on Data Assimilation Development, CPTEC/INPE, Cachoeira Paulista - SP
! 
! !DATE: May 15, 2012
!
! !INTRODUCTION: Package Overview
!
!    diag_GSI_conv is a Fortran 90 collection of routines/functions for 
!    .
!    .
!    .
!
!EOP
!---------------------------------------------------------------------
!BOP
!---------------------------------------------------------------------
! !ROUTINE: diag_gsi_cov.f90
!---------------------------------------------------------------------
!               INPE/CPTEC Data Assimilation Group                   
!---------------------------------------------------------------------
!
! !REVISION HISTORY:
!  15 May 2012 - J. G. de Mattos - Initial Version
!
!
!
!---------------------------------------------------------------------
!

PROGRAM diag_GSI_conv

  !
  ! !USES:
  !
  use kinds, only: r_kind,r_single,i_kind
  use m_utils
  use m_string
  IMPLICIT NONE


  !
  ! Auxiliary variables for diagnostic file reading
  !

  character(len=3)                            :: var
  character(8),allocatable,dimension(:)       :: cdiagbuf
  real(r_single),allocatable,dimension(:,:)   :: rdiagbuf
  real(r_single),allocatable,dimension(:,:)   :: diag 
  integer(i_kind)                             :: idate
  integer(i_kind)                             :: nchar,ninfo,nobs,mype

  !
  !  misc.
  !

  real(r_single)      :: zero_single 
  real(r_single)      :: tiny_single
  real                :: small_num
  character(8)        :: stationID
  character           :: ch
  character(len=1024) :: fname, ctlfname, mapfname
  integer             :: lu
  integer             :: idx
  integer             :: iflg
  integer             :: i,j,k, P, V, ic,ttt
  integer             :: HSN(2)
  integer             :: ios
  character(len=50)   :: rowfmt
  character(len=3)    :: DataType
  character(len=10)   :: TimeFile
  character(len=1024) :: tmp

  INTEGER   :: nymd
  INTEGER   :: nhms

  real :: Tim
  integer :: Lev
  integer :: Flag

  real, allocatable :: rtmp(:)
  logical :: OK

  !
  ! Variaveis
  !

  type(data), dimension(nvars+1) :: conv

  !
  !EOP
  !-------------------------------------------------------------------
  !BOC
  !

  if ( iargc() .lt. 2) then
     write(*,'(a)')'error: diag_gsi.x [ YYYYMMDDHH ] [ 01/02 ]'
     write(*,'(a)')'       01 - Process Diag FGS files'
     write(*,'(a)')'       02 - Process Diag ANL files'

     call exit()
  endif

  call getarg(1,TimeFile) ! Data YYYYMMDDHH
  read(TimeFile,'(I8.8,I8.8)') nymd,nhms
  IF(nhms.LT.100) nhms = nhms * 10000
  IF(nhms.GT.100.AND.nhms.LT.10000) nhms = nhms * 100

  call getarg(2,DataType) ! 01- Guess, 02- Analysis



!  year   = INT ( ymd / 10000 )
!  month  = MOD ( ymd,  10000 ) / 100
!  day    = MOD ( ymd,    100 )
!
!  select case(hms)
!  case( : 99 )
!     hour   = hms
!     minute = 0
!     second = 0
!  case(100:9999)
!     hour   = INT ( hms / 100 )
!     minute = MOD ( hms,  10 )
!     second = 0 
!  case default
!     hour   = INT ( hms / 10000 )
!     minute = MOD ( hms,  10000 ) / 100
!     second = MOD ( hms,    100 )
!  end select

  call readcard()

  !
  ! lim_qm: Valores baseados na descricao contida no manual do GSI. 
  !        Estes valores são utilizados para controlar os dados aceitos ou
  !        rejeitados.
  !        The parameter lim_qm is a threshold set in GSI read in procedure. The current values of
  !        lim_qm in GSI are listed in the following table:
  !
  !        +----------------------+---------------+---------------+
  !        |The value of namelist | lim_qm for Ps | lim_qm others |
  !        |option noiqc          |               |               |
  !        +----------------------+---------------+---------------+
  !        |True (without OI QC)  |       7       |       8       |
  !        +----------------------+---------------+---------------+
  !        |False (with OI QC)    |       4       |       4       |
  !        +----------------------+---------------+---------------+
  !

  do i=1,nvars
     conv(i)%name = vname(i)
     if (noiqc) then
        conv(i)%lim_qm = 8
        if(vname(i).eq." ps")conv(i)%lim_qm = 7
     else
        conv(i)%lim_qm = 4
     endif
  enddo

  !
  !
  !

  tiny_single = tiny(zero_single)
  small_num   = 1.e3 * tiny_single


  !
  !  Abrindo e lendo arquivos de diagnostico gerados pelo GSI
  ! Se o GSI foi executado em paralelo são abertos mais do que um arquivo. Isto
  ! é controlado pela variaveis nproc.

  lu = 10
  
  print*, '>>>>>>>> NPROC: ',nproc
  
  DO I= 0,nproc-1

     FName = TRIM(DirIN)//'/'//TRIM(FInp)
     CALL str_template(FName, nymd,nhms)
     call replace_(FName, "%e", int2str(I,'(I4.4)'))
     call replace_(FName, "%type", TRIM(DataType))

#ifdef SUMMARY
     print*,trim(fname)
#endif

     OPEN ( UNIT   = lu+I,          &
          FILE   = trim(fname),   &
          STATUS = 'OLD',         &
          IOSTAT = ios,           &
          ACCESS = 'SEQUENTIAL',  &
          FORM   = 'UNFORMATTED', CONVERT='SWAP' )

     if(ios > 0 ) then
        write(*,*) ' file is unavailabe: ', trim(fname)
        stop 123
     endif

  ENDDO

  !
  ! Reading observation diagnostics
  !

  !
  read(lu, ERR=997) idate ! time of observation analisys 
  !
  !
  ! zero fill to counters
  !
  do i=0,nlev
     conv(:)%use  ( i ) = 0
     conv(:)%nuse ( i ) = 0
     conv(:)%rej  ( i ) = 0
     conv(:)%mon  ( i ) = 0

     Do j=1,2
        conv(:)%rmse ( i,j ) = 0
        conv(:)%vies ( i,j ) = 0
     enddo

  enddo
  conv(:)%count(1)=0
  conv(:)%count(2)=0

  conv(:)%nobs = 0
  ttt=0
  idx=1
  var='tt'
  mype=0
  FILES : DO P=0,nproc-1
     !     if(P.eq.93.or.&
     !        P.eq.99.or.&
     !        P.eq.100.or.&
     !        P.ge.101)cycle
#ifdef DEBUG
     write(*,)'=======================',p,'======================='
#endif
     VARIABLES : DO

        !        print*,conv(idx)%nobs,var,mype
        read(lu+P, iostat=ios,end=110) var, nchar,ninfo,nobs,mype
        !        if(ios/=0)read(lu+P, iostat=ios,end=110)
        !        if(ios/=0)then
                 print*,var, nchar,ninfo,nobs,mype,ios
        !         cycle
        !        endif
        idx = maxloc(index(vname(:),var),1)
#ifdef DEBUG
        write(*,'(2I,A,5I)')P,conv(idx)%nobs,var, nchar,ninfo,nobs,mype, ios
#endif


        !        if (idx.eq.5)conv(idx)%nobs=ttt

        !        if(idx.eq.5)print*,conv(idx)%nobs, nobs, conv(idx)%nobs+nobs

        if (nobs > 0) then
           print*, "1",nobs
           allocate(cdiagbuf(nobs),rdiagbuf(ninfo,nobs))
           print*, "2",nobs

           read(lu+P,ERR=999,end=109) cdiagbuf, rdiagbuf

           STATION : do i=1,nobs

              j = conv(idx)%nobs + i
              k = minloc(rdiagbuf(6,i)-levs,mask=(rdiagbuf(6,i)-levs).GE.0,DIM=1)
              !if (mype.eq.94)print*,i,j, ii, conv(idx)%nobs, idx, conv(:)%nobs, conv(5)%nobs

              if (var .eq. " ps") k = 1
              conv(idx)%kx   (j) = nint(rdiagbuf(1,i)) ! observation type
              conv(idx)%lat  (j) = rdiagbuf(3,i)       ! observation latitude (degrees)
              conv(idx)%lon  (j) = rdiagbuf(4,i)       ! observation longitude (degrees)
              conv(idx)%lev  (j) = levs(k)             ! observation level reference
              conv(idx)%prs  (j) = rdiagbuf(6,i)       ! observation pressure (hPa)
              conv(idx)%time (j) = rdiagbuf(8,i) * 60  ! obs time (minutes relative to analysis time)
              conv(idx)%pbqc (j) = rdiagbuf(9,i)       ! input prepbufr qc or event mark
              conv(idx)%iuse (j) = int(rdiagbuf(12,i)) ! analysis usage flag (1=use, -1=monitoring )
              conv(idx)%iusev(j) = int(rdiagbuf(11,i)) ! analysis usage flag ( value )
              conv(idx)%robs (j) = rdiagbuf(17,i)      ! observation
              conv(idx)%diff (j) = rdiagbuf(18,i)      ! obs-ges used in analysis (K)
              conv(idx)%rmod (j) = conv(idx)%robs(j)-conv(idx)%diff(j)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
              !!
              !! Este passo causa falha de segmentacao quando rdiagbuf muito proximo a zero
              !! Necessita Investigacao
              !
              !              conv(idx)%sigo (j) = udef
              ! 
              !             if (rdiagbuf(16,i) > small_num) then! final inverse observation error (K**-1)
              !                conv(idx)%sigo(j) = 1.0/rdiagbuf(16,i)
              !             end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

              !
              ! get station ID
              !

              stationID = cdiagbuf(i)
              iflg = 0
              do ic=8,1,-1
                 ch = stationID(ic:ic)
                 if (ch > ' ' .and. ch <= 'z') then
                    iflg = 1
                 else
                    stationID(ic:ic) = ' '
                 end if
                 if (ch == ' '  .and. iflg == 1) then
                    stationID(ic:ic) = '_'
                 endif
              enddo

              conv(idx)%ID(j) = stationID


              !
              ! Counting number of observations accepted, reject and monitored
              !

              if (conv(idx)%iuse(j).eq. 1) conv(idx)%use(k)  = conv(idx)%use(k)  + 1
              if (conv(idx)%iuse(j).eq.-1) conv(idx)%nuse(k) = conv(idx)%nuse(k) + 1

              !
              !   The QC process creates a number indicating the data quality for each observation.
              ! These numbers are called QC markers in PrepBUFR files and are important as parts of the
              ! observation information. GSI uses QC markers to decide how to use the data. A brief
              ! summary of the meaning of the QC markers is as follows:
              ! 
              !    +-----------------+-----------------------------------------------------------+
              !    | QC markes range | Data Process in GSI                                       |
              !    +-----------------+-----------------------------------------------------------+
              !    |  > 15 or        |GSI skips these observations during reading procedure. That|
              !    |  <= 0           |means these observations are tossed                        | 
              !    +-----------------+-----------------------------------------------------------+
              !    |  >= lim_qm      |These observations will be in monitoring status. That means|
              !    |  and            |these observations will be read in and be processed through|
              !    |  < = 15         |GSI QC process (gross check) and innovation calculation    | 
              !    |                 |stage but will not be used in inner iteration.             |
              !    +-----------------+-----------------------------------------------------------+
              !    |  > 0            |Observations will be used in further gross check (failure  |
              !    |  and            |observation will be list in rejection), innovation         |
              !    |  < lim_qm       |caalculation, and the analysis (inner iteration).          |
              !    +-----------------+-----------------------------------------------------------+
              !

              if (conv(idx)%iuse(j).eq.-1 .and.&
                   (conv(idx)%pbqc(j) > 15.0 .or.&
                   conv(idx)%pbqc(j) <=  0.0)) conv(idx)%rej(k) = conv(idx)%rej(k) + 1

              if (conv(idx)%iuse(j).eq.-1  .and.&
                   (conv(idx)%pbqc(j) >  0.0 .and.&
                   conv(idx)%pbqc(j) <  conv(idx)%lim_qm)) conv(idx)%rej(k) = conv(idx)%rej(k) + 1

              if (conv(idx)%iuse(j).eq.-1.and.&
                   (conv(idx)%pbqc(j) >= conv(idx)%lim_qm .and.&
                   conv(idx)%pbqc(j) <= 15.0)) conv(idx)%mon(k) = conv(idx)%mon(k) + 1
              !
              ! Data adjustments
              !

              if (var == "gps")then
                  conv(idx)%robs (j) = rdiagbuf(17,i)      ! observation (data(igps,i))
                  conv(idx)%incref (j) = rdiagbuf(5,i)     ! incremental refractivity (x100 %)---- (data(igps,i)-nrefges)/data(igps,i)
                  conv(idx)%diff (j) = rdiagbuf(17,i)*rdiagbuf(5,i)      ! obs-ges used in analysis
                  conv(idx)%diffpt (j) = conv(idx)%diff(j)*100      ! obs-ges used in analysis *100%
                  conv(idx)%rmod (j) = conv(idx)%robs(j)-conv(idx)%diff(j)   ! ges refrativity (nrefges)
                  conv(idx)%elevgps (j) = rdiagbuf(7,i)    ! height in meters/ impact height in meters
                  conv(idx)%kmelev (j) = conv(idx)%elevgps(j)/1000 ! height in km/ impact height in km
                  conv(idx)%hobl (j) = rdiagbuf(19,i)      ! model vertical grid  (midpoint)
                  conv(idx)%errf (j) = rdiagbuf(16,i)      ! final inverse observation error due to superob factor (N**-1)
                  conv(idx)%normdiff (j) = conv(idx)%diff(j)*conv(idx)%errf(j)  ! normalized bending angles differences 

              endif


              !
              ! When the data is q, unit convert kg/kg -> g/kg **/
              !

              if (var == "  q") then
                 conv(idx)%robs(j) = conv(idx)%robs(j) * 1000.0
                 conv(idx)%rmod(j) = conv(idx)%rmod(j) * 1000.0
                 conv(idx)%diff(j) = conv(idx)%diff(j) * 1000.0
                 !                 conv(idx)%sigo(j)  = conv(idx)%sigo(j) * 1000.0
              end if

              !
              ! When the data is pw, replase the rprs to udef **/
              !

              if (var .eq. " pw") conv(idx)%time(j) = udef

              !
              !
              !

              !
              ! When the data is uv U return on uv idx and V return nvars + 1
              !

              if ( var .eq. " uv")then

                 !
                 ! UWIND
                 !

                 conv(idx)%robs(j) = rdiagbuf(17,i)
                 conv(idx)%diff(j) = rdiagbuf(18,i)
                 conv(idx)%rmod(j) = conv(idx)%robs(j)-conv(idx)%diff(j)

                 !
                 ! VWIND
                 !

                 conv(nvars+1)%robs(j) = rdiagbuf(20,i)
                 conv(nvars+1)%diff(j) = rdiagbuf(21,i)
                 conv(nvars+1)%rmod(j) = conv(nvars+1)%robs(j)-conv(nvars+1)%diff(j)


                 !
                 ! SPD
                 !

                 conv(idx)%robs(j) = sqrt(conv(idx)%robs(j)**2 + conv(nvars+1)%robs(j)**2)
                 conv(idx)%rmod(j) = sqrt(conv(idx)%rmod(j)**2 + conv(nvars+1)%rmod(j)**2)
                 conv(idx)%diff(j) = conv(idx)%robs(j)-conv(idx)%rmod(j)


              endif


              if(conv(idx)%robs(j) .gt. 1.0e8) then
                 conv(idx)%robs(j) = udef
                 conv(idx)%diff(j) = udef 
              endif

              !
              ! Statistical calculations
              !

              if (conv(idx)%diff(j).ne.udef)then
                 if (conv(idx)%lat(j) .gt. 0.0)then
                    conv(idx)%vies(k,1) = conv(idx)%vies(k,1) + conv(idx)%diff(j)
                    conv(idx)%rmse(k,1) = conv(idx)%rmse(k,1) + conv(idx)%diff(j)*conv(idx)%diff(j)
                    conv(idx)%count(1)=conv(idx)%count(1)+1
                 elseif(conv(idx)%lat(j) .lt. 0.0)then
                    conv(idx)%vies(k,2) = conv(idx)%vies(k,2) + conv(idx)%diff(j)
                    conv(idx)%rmse(k,2) = conv(idx)%rmse(k,2) + conv(idx)%diff(j)*conv(idx)%diff(j)
                    conv(idx)%count(2)=conv(idx)%count(2)+1
                 endif
              endif

              !
              !
              !


           enddo STATION
109        continue
           deallocate(cdiagbuf,rdiagbuf)
        else
           read(lu+P)
        endif

        conv(idx)%nobs = conv(idx)%nobs + nobs

     END DO  VARIABLES

110  continue
     CLOSE(lu+P)

  END DO  FILES


  !
  ! Adding observations at all levels
  !

  do i=1,nvars
     conv(i)%use (0) = sum(conv(i)%use (1:nlev))
     conv(i)%nuse(0) = sum(conv(i)%nuse(1:nlev))
     conv(i)%rej (0) = sum(conv(i)%rej (1:nlev))
     conv(i)%mon (0) = sum(conv(i)%mon (1:nlev))

     do j=1,2

        conv(i)%vies(0,j) = sum(conv(i)%vies(1:nlev,j))
        conv(i)%rmse(0,j) = sum(conv(i)%rmse(1:nlev,j))

        conv(i)%vies(     0,j) = conv(i)%vies(     0,j) / (conv(i)%count(j)*nlev)
        conv(i)%vies(1:nlev,j) = conv(i)%vies(1:nlev,j) / conv(i)%count(j)

        conv(i)%rmse(     0,j) = sqrt( conv(i)%rmse(     0,j) / (conv(i)%count(j)*nlev) )
        conv(i)%rmse(1:nlev,j) = sqrt( conv(i)%rmse(1:nlev,j) / conv(i)%count(j) )
     enddo

  enddo

  !
  ! writing output to plot on GMT and gnuplot
  !

  !
  ! Line and bar Graphics
  !

  !
  ! write to file
  !

  IF(DataType.EQ."01") DataType="fgs"
  ! It'll work till 6 outerloops
  IF(DataType.EQ."02" .OR. DataType.EQ."03" .OR. DataType.EQ."04" .OR. DataType.EQ."05" .OR. DataType.EQ."06") DataType="anl"

  !
  ! Spatial Output
  !
  DO v=1,nvars

     FName = TRIm(DirOu)//'/'//TRIM(FDatSSO)
     call replace_(FName, "%obstype", trim(adjustl(vname(v))))
     call replace_(FName, "%type", TRIM(DataType))
     call str_template(FName,nymd,nhms)

     Open(Unit   = 10, &
          File   = trim(FName),&
          Status = 'unknown',&
          Form   = 'Unformatted',&
          access = 'sequential'&
!          recordtype = 'stream'&
          )
!     print*,trim(adjustl(vname(v))), conv(v)%nobs

     Tim  = 0.0
!
! O Lev = 2 identifica que havera niveis na vertical
!
     Lev  = 2 
     Flag = 1



     DO I=1,conv(v)%nobs
        write(StationID,'(I8)')I
        WRITE(10)&
             StationID ,&
             conv(v)%lat(I),&
             conv(v)%lon(I),&
             Tim,           &
             Lev,           &
             Flag

!             print*, "LAT: ", conv(v)%lat(I)
!             print*, "LON: ", conv(v)%lon(I)
!
! A primeira variavel escrita corresponde ao nivel da observacao e não conta no
! ctl.
!
        WRITE(10)&
             conv(v)%lev  (I),&
             REAL(I) ,&
             conv(v)%kx   (I),&
             conv(v)%lat  (I),&
             conv(v)%lon  (I),&
             conv(v)%lev  (I),&
             conv(v)%prs  (I),&
             conv(v)%time (I),&
             conv(v)%pbqc (I),&
             conv(v)%iuse (I),&
             conv(v)%iusev(I),&
             conv(v)%robs (I),&
             conv(v)%rmod (I),&
             conv(v)%diff (I)

     ENDDO



     WRITE(10)' ENDFILE',0.0,0.0,0.0,0,1
!     Lev = 0
!     WRITE(10)&
!          conv(v)%ID(I) ,&
!          conv(v)%lat(I),&
!          conv(v)%lon(I),&
!          Tim,           &
!          Lev,           &
!          Flag
     Close(10)


     FName =  TRIM(FDatSSO)
     call replace_(FName, "%obstype", trim(adjustl(vname(v))))
     call replace_(FName, "%type", TRIM(DataType))
     call str_template(FName,nymd,nhms)

     MapFName = TRIM(FMAP)
     call replace_(MapFName, "%obstype", trim(adjustl(vname(v))))
     call replace_(MapFName, "%type", TRIM(DataType))
     call str_template(MapFName,nymd,nhms)

     ctlFName = TRIm(DirOu)//'/'//TRIM(FCTL)
     call replace_(CtlFName, "%obstype", trim(adjustl(vname(v))))
     call replace_(CtlFName, "%type", TRIM(DataType))
     call str_template(CtlFName,nymd,nhms)

     tmp='%h2Z%d2%MC%y4'
     call str_template(tmp,nymd,nhms)
     call writectl(FName,ctlFname,MapFName,tmp)

     ! A execucao do stnmap foi incluida no shell script run_stat.sh
     !write(tmp,'(2A)')'stnmap -i ',TRIM(ctlFName)
     !call system(tmp)


  ENDDO


  !dset ps_spatial_summary_obs_01.bin
  !DTYPE  station 
  !STNMAP psspatial_summary_obs_01.map
  !UNDEF  -999.0
  !TITLE  Station Data Sample
  !TDEF   1 linear 12z18jan1992 12hr
  !VARS 12
  !ID    1 99 Station ID
  !kx    1 99 observation type
  !rlat  1 99 observation latitude (degrees)
  !rlon  1 99 observation longitude (degrees)
  !rlev  1 99 observation level reference
  !prs   1 99 observation pressure (hPa)
  !rtime 1 99 obs time (minutes relative to analysis time)
  !pbqc  1 99 input prepbufr qc or event mark
  !iuse  1 99 analysis usage flag (1=use, -1=monitoring )
  !iusev 1 99 analysis usage flag ( value ) 
  !robs  1 99 observation               
  !rmod  1 99 model 
  !diff  1 99 OMF
  !ENDVARS
                                                           

  !
  ! Spatial Output
  !

  allocate(rtmp(nlev))
  DO v=1,nvars

     !     WRITE(FName,'(A,A,A)')'spatial_summary_obs_count_',trim(adjustl(vname(v))),'.dat'
     !     Open(Unit=10, File=Trim(FName))

     FName = TRIm(DirOu)//'/'//TRIM(FTxtSSO)
     call replace_(FName, "%obstype", trim(adjustl(vname(v))))
     call replace_(FName, "%type", TRIM(DataType))
     call str_template(FName,nymd,nhms)

     Open(Unit=10,File=trim(FName))

     DO I=1,conv(v)%nobs

 !       if(conv(v)%iuse(I).ne.1)cycle

        WRITE(10,'(13(F15.6,2x))')&
             conv(v)%kx   (I),&
             conv(v)%iuse   (I),&
             conv(v)%robs (I),&
             conv(v)%lev  (I),&
             conv(v)%lat  (I),&
             conv(v)%lon  (I),&
             conv(v)%diff (I),&
             conv(v)%diffpt (I),&
             conv(v)%rmod (I),&
             conv(v)%kmelev (I),&
             conv(v)%incref (I),&
             conv(v)%hobl (I),&
             conv(v)%normdiff (I)
! conv(v)%errf (I),&
!conv(v)%time (I),&


     ENDDO
     Close(10)


!    FName = TRIm(DirOu)//'/KX_%obstype_%type.txt'
!    call replace_(FName, "%obstype", trim(adjustl(vname(v))))
!    call replace_(FName, "%type", TRIM(DataType))
!    call str_template(FName,nymd,nhms)

!    Open(Unit=10,File=trim(FName),form='formatted')

!    DO I=1,conv(v)%nobs

!       if(conv(v)%iuse(I).ne.1)cycle

!       WRITE(10,'(F9.3)')conv(v)%kx (I)
!    ENDDO
!    Close(10)

     !
     ! salvando valor médio sobre o globo de OMF/OMF para os dados assimilados
     !


     DO I=1,nkx

        Flag = count(conv(v)%kx.eq.kx(I).and.conv(v)%iuse.eq.1)
        if (Flag.eq.0)cycle


        Open(Unit=11,File=trim(FName),position='append')

        

        DO J=1,nlev

           rtmp(J) = udef
           Flag    = count(conv(v)%lev.eq.levs(J))
           if (Flag.eq.0)cycle

           rtmp(J) = sum(conv(v)%diff,&
                     mask=conv(v)%iuse.eq.1.and.conv(v)%lev.eq.levs(J).and.nint(conv(v)%kx).eq.kx(I))/&
                     count(conv(v)%iuse.eq.1.and.conv(v)%lev.eq.levs(J).and.nint(conv(v)%kx).eq.kx(I))

           if (isnan(rtmp(J))) rtmp(J) = udef ! funcao do gfortran e ifort (fortran padrao)
        ENDDO

        FName = TRIM(FTxtOMF)
        call replace_(FName, "%obstype", trim(adjustl(vname(v))))
        call replace_(FName, "%type", TRIM(DataType))
        call replace_(FName, "%RepType", int2str(kx(I),'(I4.4)'))
        call str_template(FName,nymd,nhms)

        Open(Unit=11,File=trim(FName),position='append')
        write(tmp,'(A,I3.2,A)')'(A,1x',nlev,'F9.3)'
        WRITE(11,tmp)Trim(TimeFile),(rtmp(J),j=1,nlev)
        close(11)


     ENDDO
!     close(11)
 
     !
     ! salvando valor médio sobre o globo de OMF/OMF dos dados monitorados
     ! Apenas para os dados 153 GPSIPW e 126 RASSDA vento sensor acustico adicionar outros
     !


     DO I=1,nkx

        Flag = count((((kx(I).eq.153).and.vname(v).eq.' pw').or.((kx(I).eq.126).and.vname(v).eq.'  t')).and.conv(v)%iuse.eq.-1)
        if (Flag.eq.0)cycle


        Open(Unit=11,File=trim(FName),position='append')

        

        DO J=1,nlev

           rtmp(J) = udef
           Flag    = count(conv(v)%lev.eq.levs(J))
           if (Flag.eq.0)cycle

           rtmp(J) = sum(conv(v)%diff,&
                     mask=conv(v)%iuse.eq.-1.and.conv(v)%lev.eq.levs(J).and.nint(conv(v)%kx).eq.kx(I))/&
                     count(conv(v)%iuse.eq.-1.and.conv(v)%lev.eq.levs(J).and.nint(conv(v)%kx).eq.kx(I))

           if (isnan(rtmp(J))) rtmp(J) = udef ! funcao do gfortran e ifort (fortran padrao)
        ENDDO

        FName = TRIM(FTxtOMF)
        call replace_(FName, "%obstype", "%obstype_mon")
        call replace_(FName, "%obstype", trim(adjustl(vname(v))))
        call replace_(FName, "%type", TRIM(DataType))
        call replace_(FName, "%RepType", int2str(kx(I),'(I4.4)'))
        call str_template(FName,nymd,nhms)

        Open(Unit=11,File=trim(FName),position='append')
        write(tmp,'(A,I3.2,A)')'(A,1x',nlev,'F9.3)'
        WRITE(11,tmp)Trim(TimeFile),(rtmp(J),j=1,nlev)
        close(11)


     ENDDO 
 

  ENDDO


  !
  ! Count Summary
  !

  FName = TRIm(DirOu)//'/'//TRIM(FSOC)
  call replace_(FName, "%type", TRIM(DataType))
  call str_template(FName,nymd,nhms)

  Open(Unit=10, File=trim(FName))

  WRITE(10,'(  6A13)')'Obs Name','# of OBS','Used','Not Used','Rejeited','Monitored'  
  DO I=1,nvars
     WRITE(10,'(A12,5I12)') vname(i), conv(i)%nobs, conv(i)%use(0), conv(i)%nuse(0), conv(i)%rej(0), conv(i)%mon(0)
  ENDDO
  Close(10)

  !
  ! used
  !

  FName = TRIm(DirOu)//'/'//TRIM(FSOCL)
  call replace_(FName, "%stat", "tused")
  call replace_(FName, "%type", TRIM(DataType))
  call str_template(FName,nymd,nhms)

  Open(Unit=11, File=trim(FName))

  write(rowfmt,'(A,I4,A)') '(A9,',nvars,'A9)'
  write(11,rowfmt)'Levels',vname(1:nvars)
  write(rowfmt,'(A,I4,A)') '(F9.1,',nvars,'I9)'
  Do i=1,nlev
     write(11,rowfmt)levs(i),(conv(j)%use(i),j=1,nvars)
  Enddo

  !
  ! Not Used
  !

  FName = TRIm(DirOu)//'/'//TRIM(FSOCL)
  call replace_(FName, "%stat", "nused")
  call replace_(FName, "%type", TRIM(DataType))
  call str_template(FName,nymd,nhms)

  Open(Unit=11, File=trim(FName))

  write(rowfmt,'(A,I4,A)') '(A9,',nvars,'A9)'
  write(11,rowfmt)'Levels',vname(1:nvars)
  write(rowfmt,'(A,I4,A)') '(F9.1,',nvars,'I9)'
  Do i=1,nlev
     write(11,rowfmt)levs(i),(conv(j)%nuse(i),j=1,nvars)
  Enddo

  !
  ! Rejeited
  !

  FName = TRIm(DirOu)//'/'//TRIM(FSOCL)
  call replace_(FName, "%stat", "rejec")
  call replace_(FName, "%type", TRIM(DataType))
  call str_template(FName,nymd,nhms)

  Open(Unit=11, File=trim(FName))

  write(rowfmt,'(A,I4,A)') '(A9,',nvars,'A9)'
  write(11,rowfmt)'Levels',vname(1:nvars)
  write(rowfmt,'(A,I4,A)') '(F9.1,',nvars,'I9)'
  Do i=1,nlev
     write(11,rowfmt)levs(i),(conv(j)%rej(i),j=1,nvars)
  Enddo

  !
  ! Monitored
  !

  FName = TRIm(DirOu)//'/'//TRIM(FSOCL)
  call replace_(FName, "%stat", "monit")
  call replace_(FName, "%type", TRIM(DataType))
  call str_template(FName,nymd,nhms)

  Open(Unit=11, File=trim(FName))


  write(rowfmt,'(A,I4,A)') '(A9,',nvars,'A9)'
  write(11,rowfmt)'Levels',vname(1:nvars)
  write(rowfmt,'(A,I4,A)') '(F9.1,',nvars,'I9)'
  Do i=1,nlev
     write(11,rowfmt)levs(i),(conv(j)%mon(i),j=1,nvars)
  Enddo

  !
  ! Vies and RMS Summary
  !


  !
  ! HN
  !

!  Open(Unit=10, File=trim(dirou)//'/RMSE_HN_'//bufr//'.txt')
!  write(rowfmt,'(A,I4,A)') '(A9,',nlev+1,'F15.9)'
!  DO I=1,nvars
!     WRITE(10,rowfmt) vname(i), (conv(i)%rmse(J,1),J=0,nlev)
!  ENDDO
!  Close(10)


!  Open(Unit=10, File=trim(dirou)//'/VIES_HN_'//bufr//'.txt')
!  write(rowfmt,'(A,I4,A)') '(A9,',nlev+1,'F15.9)'
!  DO I=1,nvars
!     WRITE(10,rowfmt) vname(i), (conv(i)%vies(J,1),J=0,nlev)
!  ENDDO
!  Close(10)

  !
  ! HS
  !

!  Open(Unit=10, File=trim(dirou)//'/RMSE_HS_'//bufr//'.txt')
!  write(rowfmt,'(A,I4,A)') '(A9,',nlev+1,'F15.9)'
!  DO I=1,nvars
!     WRITE(10,rowfmt) vname(i), (conv(i)%rmse(J,2),J=0,nlev)
!  ENDDO
!  Close(10)


!  Open(Unit=10, File=trim(dirou)//'/VIES_HS_'//bufr//'.txt')
!  write(rowfmt,'(A,I4,A)') '(A9,',nlev+1,'F15.9)'
!  DO I=1,nvars
!     WRITE(10,rowfmt) vname(i), (conv(i)%vies(J,2),J=0,nlev)
!  ENDDO
!  Close(10)

#ifdef SUMMARY
  !
  ! Writing summary of count
  !

  WRITE(6,'(1x,A75)')'---------------------------------------------------------------------------'
  WRITE(6,'(1x,A75)')'                          S  U  M  M  A  R  Y                              '
  WRITE(6,'(1x,A75)')'---------------------------------------------------------------------------'
  WRITE(6,'(  6A12)')'Obs Name','# of OBS','Used','Not Used','Rejected','Monitored'
  WRITE(6,'(1X,A75)')'---------------------------------------------------------------------------'
  DO I=1,nvars
     WRITE(6,'(A12,5I12)') vname(i), conv(i)%nobs, conv(i)%use(0), conv(i)%nuse(0), conv(i)%rej(0), conv(i)%mon(0)
  ENDDO
  WRITE(6,'(1X,A75)')'---------------------------------------------------------------------------'
  WRITE(6,'(A12,5I12)')'Total',sum(conv(:)%nobs),sum(conv(:)%use(0)),sum(conv(:)%nuse(0)), sum(conv(:)%rej(0)), sum(conv(:)%mon(0))
  WRITE(6,'(1X,A75)')'---------------------------------------------------------------------------'

  write(rowfmt,'(A,I4,A)') '(A6,',nlev,'F9.1)'  
  WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------'
  write(6,'(1x,A120)')'                                            U S E D   B Y   L E V E L                                                   '
  write(6,rowfmt)'Top',(levs(i),i=1,nlev)
  write(6,rowfmt)'Bot',1200.0,(levs(i)-0.1,i=1,nlev-1)
  WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------'
  write(rowfmt,'(A,I4,A)') '(A6,',nlev,'I9)'
  do i=1,nvars
     write(6,rowfmt)vname(i), conv(i)%use(1:nlev)
  enddo
  write(rowfmt,'(A,I4,A)') '(A6,',nlev,'F9.1)'
  WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------'
  write(6,'(1x,A120)')'                                        N O T    U S E D   B Y   L E V E L                                              '
  write(6,rowfmt)'Top',(levs(i),i=1,nlev)
  write(6,rowfmt)'Bot',1200.0,(levs(i)-0.1,i=1,nlev-1)
  WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------'
  write(rowfmt,'(A,I4,A)') '(A6,',nlev,'I9)'
  do i=1,nvars
     write(6,rowfmt)vname(i), conv(i)%nuse(1:nlev)
  enddo
  write(rowfmt,'(A,I4,A)') '(A6,',nlev,'F9.1)'
  WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------'
  write(6,'(1x,A120)')'                                        R E J E I T E D    B Y   L E V E L                                              '
  write(6,rowfmt)'Top',(levs(i),i=1,nlev)
  write(6,rowfmt)'Bot',1200.0,(levs(i)-0.1,i=1,nlev-1)
  WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------'
  write(rowfmt,'(A,I4,A)') '(A6,',nlev,'I9)'
  do i=1,nvars
     write(6,rowfmt)vname(i), conv(i)%rej(1:nlev)
  enddo
  write(rowfmt,'(A,I4,A)') '(A6,',nlev,'F9.1)'
  WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------'
  write(6,'(1x,A120)')'                                        M O N I T O R E D   B Y   L E V E L                                             '
  write(6,rowfmt)'Top',(levs(i),i=1,nlev)
  write(6,rowfmt)'Bot',1200.0,(levs(i)-0.1,i=1,nlev-1)
  WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------'
  write(rowfmt,'(A,I4,A)') '(A6,',nlev,'I9)'
  do i=1,nvars
     write(6,rowfmt)vname(i), conv(i)%mon(1:nlev)
  enddo
#endif

  STOP 'ENDS NORMALLY'
997 PRINT *,'error read in diag file 997'
  stop 1
998 PRINT *,'error read in diag file 998'
  stop 2
999 PRINT *,'error read in diag file 999'
  stop 3


END PROGRAM diag_GSI_conv
!
!EOCdiag_gsi_conv.f90
!---------------------------------------------------------------------

