!
      program sim_obs_error
!
!  Main program for creation of simulated conventioanl data for an OSSE 
!
!  Initial program written by Ronald Errico January 2008
!  Addition to perturb surface pressure     September 2008
!
      use m_kinds       ! precision specification module for real variables
      use m_obs_pert    ! perturbation routines
      use m_relhum      ! transforms between specific and relative humidity
      use m_bufr_rw     ! Read/Write BUFR data
!
!  link to libNCEP_bufr_r4i4.a library
!  Also requires a routine to compute eigenvectors/values for symm. matrix
!
      implicit none
!
! If ltest_all is true, then the program will use fake, low resolution NR 
! data and fake observations and will print out many intermediate results
! for testing purposes.  Alternatively, ltest=.true. or .false. can be 
! specified at desired points in the code to print out only those portions. 
      logical, parameter :: ltest_all=.false. 
!
! Specify maximum dimensions for obs data arrays
      integer, parameter :: obs_ndim1=305  ! >= # of p-levels or channels 
      integer, parameter :: obs_ndim2=2    ! >= # of fields of data (excl ps)
      integer, parameter :: obs_ndim3=2    ! >= # of p-lev info 
      integer, parameter :: obs_ndim4=12*3+25 ! >= # of hdr info for each obs 
! obs_ndim1: maximum number of p-levels or channels for observations 
! obs_ndim2: maximum number of observation fields in buffer (This only refers 
!            to obs derived from 3d-fields; e.g; 2 for u,v, but don't count ps)
! obs_ndim3: number of for each level or channel, besides the obs values
!            themselves; e.g., p-levels, channels, and QC marks, etc.
! obs_ndim4: maximum number of single value report info (lat, lon etc.)
!            12*3+25 is required value for AIRS file reports
!
      integer, parameter :: n_errors=11    ! number of types of counters
! n_errors: number of types of counters for errors and observations processed
!

      logical :: l_comp_z   ! true if height values present
      logical :: ltest
      integer :: ierrors(n_errors)
      integer :: obs_nlevs
      integer :: obs_nfields
      integer :: nobs
      real(rkind3) :: obs_info(obs_ndim4)  ! id,lon,lat,time,type,height,ps,...
      real(rkind3) :: obs_levels(obs_ndim1,obs_ndim3) ! plevs, QCflags, ...
      real(rkind3) :: obs_values(obs_ndim1,obs_ndim2) ! u,v or t,q, or channels
      real(rkind3) :: corr_dist(obs_ndim2) ! correlation distance parameter
!  correlation distance here is the ratio of two pressure levels over which
!  relative "distance" the error correlation drops to 0.1 
!
! Variables realted to the data bufr file
      logical lsmg1
      integer luin, luout, luin_tab  !  unit numbers
      integer argc
      integer(4) iargc
      integer ireadmg
      integer ireadsb
      integer n_mesg                 ! counter of messages in file
      character(len=120) file_rc     ! resource file name
      character(len=120) inputfile   ! bufr input file of simulated obs
      character(len=120) outputfile  ! bufr output file of pvert. sim. obs.
      character(len=120) cdatetime   ! synoptic time yyyymmddhh
      character(len=5)   d_type      ! Type of observation
      character(len=8)  subset       ! name of current BUFR subset
      character(len=8)  psubset      ! name of previous BUFR subset
      integer      idate             ! synoptic date/time YYYYMMDD
      logical check_type
      logical      modtype  ! flag indicating whether type of report to modify
      real(rkind3)   bmiss            ! value of flag to indicate missing value
      data luin /8/, luout /9/, luin_tab /7/  ! do not use units 10-12
!
! Variables concerning errors and error table  
      integer :: icolumns(obs_ndim2)
      integer :: luin_err, luin_rc
      integer :: i, j, k
      integer, parameter :: random_type=1111 ! part of seed for random 
      integer :: err_index
      integer :: itype
      integer*8 :: i_random(1)
      integer*8 :: idatetime
      integer, allocatable :: err_itype(:)
      integer :: n_err1, n_err2, n_err3
      logical :: lrelhum,lmsg1
      logical, parameter :: lsim=.false.  ! false since only perts added
      real(rkind3)  :: pert_fac
      real(rkind3), allocatable :: err_tab(:,:,:)
      character(len=2) :: d_type_m  ! type of data
      data luin_err /20/, luin_rc /21/
      real(rkind3) :: pert_fac_wind,  corr_dist_wind
      real(rkind3) :: pert_fac_mass,  corr_dist_mass, corr_dist_q
      real(rkind3) :: pert_fac_airs,  corr_dist_airs
      integer      :: random_type_wind, random_type_mass, random_type_airs
      integer      :: random_case  
      character(len=120) :: file_err_table
!
!
      print *,' '
      print *,'BEGIN PROGRAM sim_obs_error'    
!
! Read and check arguments
      argc = iargc()
      if (argc .lt. 5) then
         print *,'usage must be: sim_error.x d_type datetime file_rc ', &
                 'inputbufr outputbufr'
         stop
      endif
      call GetArg( 1_4, d_type)
      call GetArg( 2_4, cdatetime)
      call GetArg( 3_4, file_rc)
      call GetArg( 4_4, inputfile)
      call GetArg( 5_4, outputfile)
      if ( (d_type == 'WIND_') .or. (d_type == 'MASS_') .or. &
           (d_type == 'HIRS2') .or. (d_type == 'HIRS3') .or. &
           (d_type == 'AMSUA') .or. (d_type == 'AMSUB') .or. &
           (d_type == 'AIRS_') .or. (d_type == '_MSU_') ) then
        print *,'Processing for d_type=',d_type
      else
        print *,' Specification of d_type=',d_type,' neither ', &
              'WIND_, MASS_, HIRS2, ','HIRS3, AMSUA, AMSUB, ',  & 
              'AIRS_ or _MSU_ as required'
        stop
      endif
!
! Set some variables
      ltest=ltest_all
      bmiss = 10.e10  
      ierrors(:)=0
      n_mesg=0
!
! Allocate array for table of errors
      n_err1=33      ! max number of p-levels in table
      n_err2=6       ! number of fields in table
      n_err3=35      ! max number of sub-types of observations in table
      if ( .not. ((d_type == 'WIND_').or.(d_type == 'MASS_'))) then
        n_err1=40    ! max number of channels (AIRS needs more)
        n_err2=2     
        n_err3=6     ! max number of sub-types of observations in table
      endif
      if (d_type == 'AIRS_') n_err1=301
      allocate (err_itype(n_err3))
      allocate (err_tab(n_err1,n_err2,n_err3))
!
! Read resource file
      call read_error_rc (obs_ndim2,luin_rc,file_rc,d_type,i_random, &
                          pert_fac,corr_dist,file_err_table)
!
! Set some variables
      if (.not. (d_type == 'AIRS_')) then 
        luin_tab=luin             ! bufr table is in data file 
      endif
      if (d_type == 'MASS_') then ! q error in terms of rh
        lrelhum=.true.
        call relhum_setup         ! setup table for saturation vapor press.
      else
        lrelhum=.false.
      endif
!
! Set seed for random number generator to datetime
      read (cdatetime(1:10),'(i10)') idatetime
      i_random(1)=idatetime+i_random(1)
      call random_seed (put=i_random(1:1))
      print ('(a,i12,a,i12)'),'Seed for random number generator = ' &
             ,i_random(1),'  idatetime=',idatetime 
!
! Read error table
     call read_error_table (n_err1,n_err2,n_err3,luin_err,d_type, &
                            file_err_table,err_itype,err_tab)
!
! Open AIRS bufr table if required
      if (d_type=='AIRS_') then 
        open(unit=luin_tab,file ='airs_bufr_table', form='formatted')
        print *,' '
        print ('(3a,i4)'),' input file=','airs_bufr_table', &
                          ' opened on unit=',luin_tab       
      endif
! 
! Setup files for reading and writing observations
      open(unit=luin,file =trim(inputfile), form='unformatted')
      print *,' '
      print ('(3a,i4)'),' input file=',trim(inputfile), &
                        ' opened on unit=',luin
      open(unit=luout,file=trim(outputfile),form='unformatted')
      print *,' '
      print ('(3a,i4)'),' output file=',trim(outputfile), &
                        ' opened on unit=',luout
      call openbf(luin, 'IN ',luin_tab)
      call openbf(luout,'OUT',luin_tab)
!
! Compress messages in BUFR files written (so several records fit in 1 msg)
      if ((d_type /= 'WIND_') .and. (d_type /= 'MASS_')) call cmpmsg('Y')
!
! Begin loop to read input observations
      nobs=0  
      psubset =''
      if (d_type == 'AIRS_' .or. d_type == 'WIND_' .or. d_type == 'MASS_') then
        lmsg1=.false.   ! do not need to create first message with 0 reports
      else
        lmsg1=.true.    
      endif
!
      do while (ireadmg(luin,subset,idate).eq. 0)
! If new data source (new  'subset')
        if (subset .ne. psubset ) then
          print *,'Processing subset ',subset,' for idate ',idate
! Create first message as just date/time stamp with no reports
! GSI is designed to skip the first message on file, so this must be message
! filled if the first set of reports are to be read by GSI
          if (lmsg1) then
            call openmg(luout,subset,idate)  ! idate here is yyyymmddhh
            call minimg(luout,0)             ! add minutes=00 to idate write
            call closmg(luout)
            lmsg1=.false. 
          endif
          psubset = subset
        endif
        modtype = check_type (subset, d_type)
        if ( modtype ) then
          call openmb(luout,subset,idate)
          do while ( ireadsb(luin) .eq. 0 )
!
            if (d_type == 'WIND_' ) then 
              call read_write_obs_uv (obs_ndim1,obs_ndim2,obs_ndim3, &
                            obs_ndim4,obs_nlevs,obs_nfields,nobs,    &
                            ierrors(9),luin,luout,ltest,lsim,bmiss,  & 
                            obs_info,obs_levels,obs_values,          &
                            l_comp_z,.true.)  
              if (obs_nlevs > 0) then
                icolumns(1:2)=4                      
                itype=nint(obs_info(5)) 
                call pert_find_itype (n_err3,itype,err_itype,err_index)
                if (err_index == 0) then
                  ierrors(1)=ierrors(1)+1
                else   
                  call pert_obs (n_err1,n_err2,obs_ndim1,obs_ndim2,nobs, &
                        obs_nlevs,obs_nfields,icolumns,bmiss,pert_fac,   &
                        err_tab(:,:,err_index),obs_levels,obs_values,    &
                        corr_dist,ltest,lrelhum,.false.)  
                  call read_write_obs_uv (obs_ndim1,obs_ndim2,obs_ndim3, &
                                obs_ndim4,obs_nlevs,obs_nfields,nobs,    &
                                ierrors(9),luin,luout,ltest,lsim,bmiss,  & 
                                obs_info,obs_levels,obs_values,          &
                                l_comp_z,.false.)  
                endif  ! test on err_index
              endif    ! test on obs_nlevs
!
            elseif (d_type == 'MASS_' ) then 
              call read_write_obs_tq (obs_ndim1,obs_ndim2,obs_ndim3,  &
                            obs_ndim4,obs_nlevs,obs_nfields,nobs,     &
                            ierrors(9),luin,luout,ltest,.true.,bmiss, & 
                            obs_info,obs_levels,obs_values,           &
                            l_comp_z,.true.)  
              if (obs_nlevs > 0) then
                icolumns(1)=2                                    
                icolumns(2)=3                                    
                itype=nint(obs_info(5)) 
                call pert_find_itype (n_err3,itype,err_itype,err_index)
                if (err_index == 0) then
                  ierrors(1)=ierrors(1)+1
                else   
                  call pert_ps (n_err1,n_err2,obs_ndim1,nobs,obs_nlevs,  &
                        pert_fac,err_tab(:,:,err_index),obs_levels(1,1), &
                        ltest)
                  call pert_obs (n_err1,n_err2,obs_ndim1,obs_ndim2,nobs, &
                        obs_nlevs,obs_nfields,icolumns,bmiss,pert_fac,   &
                        err_tab(:,:,err_index),obs_levels,obs_values,    &
                        corr_dist,ltest,lrelhum,.false.)  
                  call read_write_obs_tq (obs_ndim1,obs_ndim2,obs_ndim3, &
                        obs_ndim4,obs_nlevs,obs_nfields,nobs,            &
                        ierrors(9),luin,luout,ltest,lsim,bmiss,          & 
                        obs_info,obs_levels,obs_values,                  &
                        l_comp_z,.false.)  

                endif  ! test on err_index
              endif    ! test on obs_nlevs
!
            elseif (d_type == 'AIRS_' ) then ! (both AIRS and AMSU-A on AQUA) 
              call read_write_obs_airs (obs_ndim1,obs_ndim2,         &
                     obs_ndim3,obs_ndim4,obs_nlevs,obs_nfields,nobs, &
                     ierrors(9),luin,luout,ltest,lsim,bmiss,         & 
                     obs_info,obs_levels,obs_values,.true.) 
              itype=nint(obs_info(1))  ! instrument ID number 
              call pert_find_itype (n_err3,itype,err_itype,err_index)
              if (err_index == 0) then
                ierrors(1)=ierrors(1)+1
              else   
                icolumns(1)=2
                call pert_obs (n_err1,n_err2,obs_ndim1,obs_ndim2,nobs, &
                      obs_nlevs,obs_nfields,icolumns,bmiss,pert_fac,   &
                      err_tab(:,:,err_index),obs_levels,obs_values,    &
                      corr_dist,ltest,lrelhum,.true.)
! Set HSB data to missing
                obs_values(297:301,1)=bmiss                
                call read_write_obs_airs (obs_ndim1,obs_ndim2,obs_ndim3, &
                          obs_ndim4,obs_nlevs,obs_nfields,nobs,          &
                          ierrors(9),luin,luout,ltest,lsim,bmiss,        & 
                          obs_info,obs_levels,obs_values,.false.)  
              endif
!
            else   ! HIRS, AMSU, or MSU                  
               call read_write_obs_amsu_hirs (obs_ndim1,obs_ndim2,   &
                     obs_ndim3,obs_ndim4,obs_nlevs,obs_nfields,nobs, &
                     ierrors(9),luin,luout,ltest,lsim,bmiss,         & 
                     obs_info,obs_levels,obs_values,.true.)  
              itype=nint(obs_info(14))  ! satellite ID number
              call pert_find_itype (n_err3,itype,err_itype,err_index)
              if (err_index == 0) then
                ierrors(1)=ierrors(1)+1
              else   
                icolumns(1)=2 
                call pert_obs (n_err1,n_err2,obs_ndim1,obs_ndim2,nobs, &
                      obs_nlevs,obs_nfields,icolumns,bmiss,pert_fac,   &
                      err_tab(:,:,err_index),obs_levels,obs_values,    &
                      corr_dist,ltest,lrelhum,.true.)
                call read_write_obs_amsu_hirs (obs_ndim1,obs_ndim2,    &
                       obs_ndim3,obs_ndim4,obs_nlevs,obs_nfields,nobs, &
                       ierrors(9),luin,luout,ltest,lsim,bmiss,         & 
                       obs_info,obs_levels,obs_values,.false.)  
              endif  ! test on err_index    
!
            endif    ! test on d_type
!
! End loop over buffer read
          enddo               ! loop over do while ( ireadsb(luin) .eq. 0 )
          call closmg(luout)  
        endif                 ! check on modtype
        n_mesg=n_mesg+1       ! count messages read in original file
      enddo                   ! loop over (ireadmg(luin,subset,idate).eq. 0)
      call closbf(luout)        
!
! Print out number of obs records processed
      print *,' '
      print ('(i7,2a)'),n_mesg,' observation message-groups read all types'
      print ('(i7,2a)'),ierrors(1),' observation reports found with no',  &
                         ' corresponding entry in error table'
      print ('(i7,a)' ),ierrors(9),' read/write errors detected'
      print ('(i7,2a)'),nobs,' observation reports processed for type = ' &
                       ,d_type
      print *,' '
      print *,'program completed' 
!
      end program sim_obs_error
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
     subroutine fix_error_table (n_err1,n_err2,n_err3,err_tab)
!
     use m_kinds
     implicit none
     integer :: n_err1   
     integer :: n_err2
     integer :: n_err3
     real(rkind3) :: err_tab(n_err1,n_err2,n_err3)
!
! The following are set to signify "large" errors
     real(rkind3), parameter :: T_err=4.
     real(rkind3), parameter :: q_err=0.05   ! % rel. humidity
     real(rkind3), parameter :: w_err=5.     ! each wind component
     real(rkind3), parameter :: p_err=200.   ! hPa
     integer :: i,j
!
! Change units of p from hPa to Pa
     do i=1,n_err3
       do j=1,n_err1
         err_tab(j,1,i)=100.*err_tab(j,1,i)    ! p levels     
         if (err_tab(j,5,i) > 1.e8) then       ! replace large p errors
           err_tab(j,5,i)= p_err
         else
           err_tab(j,5,i)=100.*err_tab(j,5,i)  ! stdv of p error  
         endif 
       enddo
     enddo
!
! replace values of 0.1e10 read in table by reasonably large values 
     do i=1,n_err3
       do j=1,n_err1
         if (err_tab(j,2,i) > 1.e8) then   ! replace large T errors
           err_tab(j,2,i)= T_err
         endif
         if (err_tab(j,4,i) > 1.e8) then   ! replace large wind errors
           err_tab(j,4,i)= w_err
         endif
       enddo
     enddo
!
! Change q error to fraction of rh
! Also, don't allow large errors at very low p where rh computed poorly
     do i=1,n_err3
       do j=1,n_err1
         if (err_tab(j,3,i) > 1.e8) then       ! replace large q errors
            err_tab(j,3,i)= q_err
         else
            err_tab(j,3,i)=0.1*err_tab(j,3,i)  ! rescale error to fraction rh  
         endif
         if (err_tab(j,1,i) < 3.e4) then   ! reduce values at low pressure
           err_tab(j,3,i)=q_err*(err_tab(j,1,i)/3.e4)**3
         endif            
       enddo
     enddo
! 
     end subroutine fix_error_table 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine read_error_rc (obs_ndim2,luin_rc,file_rc,d_type,i_random, &
                                pert_fac,corr_dist,file_err_table)
!
      use m_kinds 
      implicit none      
      integer      :: obs_ndim2
      integer      :: luin_rc
      integer(8)   :: i_random(1)
      real(rkind3) :: pert_fac
      real(rkind3) :: corr_dist(obs_ndim2)
      character(len=*) :: file_rc
      character(len=*) :: d_type
      character(len=*) :: file_err_table
!
      integer :: i
      integer(8) :: n_lines
      integer(8) :: random_case, random_type
      character(len=16)  :: name
!
! Read resource file for constructing errors
      open(unit=luin_rc,file=trim(file_rc), form='formatted')
      print *,' '
      print ('(3a,i4)'),' input file=',trim(file_rc), &
              ' opened on unit=',luin_rc
!
      read (luin_rc,'(a16)') name      ! this read acts to skips this line
      read (luin_rc,'(a16,2x,i5)') name,n_lines  ! number of lines in file
      read (luin_rc,'(a16,2x,i5)') name,random_case
      do i=4,n_lines
        read (luin_rc,'(a5)') name(1:5)
        if (name(1:5) == d_type) exit
        if (i==n_lines) then
          print *,' d_type=',d_type,' not found on rc file'
          stop
        endif 
      enddo
!      
      read (luin_rc,'(a16,2x,f5.2)') name,pert_fac
      read (luin_rc,'(a16,2x,i5)'  ) name,random_type
      read (luin_rc,'(a16,2x,f5.2)') name,corr_dist(1)
      read (luin_rc,'(a16,2x,f5.2)') name,corr_dist(2)
      read (luin_rc,'(a16,2x,a)'   ) name,file_err_table
!
      close (luin_rc)
      i_random(1)=random_type+random_case
! 
      end subroutine read_error_rc 
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine read_error_table (n_err1,n_err2,n_err3,luin_err,d_type, &
                                  file_err_table,err_itype,err_tab)
!
! Note that if this routine requires changes, likely also does the
! subroutine "find_rad_subset" used with the module "m_rdata_boxes"
!
      use m_kinds
      implicit none
      integer :: n_err1,n_err2,n_err3
      integer :: luin_err
      integer :: err_itype(n_err3)
      real(rkind3) :: err_tab(n_err1,n_err2,n_err3)     
      character(len=*) :: file_err_table
      character(len=*) :: d_type
!
      integer :: i, j
      integer :: ic, inum, i1, i2
      integer :: ihead   ! number of header records to skip in sat_error file
      integer :: igroups ! number of sat platforms or instruments in table
      integer :: n_channels
      character(len=18) :: adum
      character(len=5)  :: i_type
      character(len=4)  :: sat_name
      character(len=1)  :: adum1
!
! Read table of obs error standard deviations for conventional obs.
      open(unit=luin_err,file=trim(file_err_table), form='formatted')
      print *,' '
      print ('(3a,i4)'),' input file=',trim(file_err_table), &
              ' opened on unit=',luin_err
!
      if ((d_type=='WIND_') .or. (d_type=='MASS_') ) then 
        do i=1,n_err3  ! loop over observation types
          read (luin_err,'(i4)') err_itype(i)
          do j=1,n_err1
            read (luin_err,'(1x,6e12.5)') err_tab(j,1:6,i)
          enddo
        enddo
!
      else  ! prepare table for satellites
!        
        ic=0
        err_itype(:)=0
!
! Read header: 
! igroups=number of distinct
        read (luin_err,'(7x,i2,30x,i2)') ihead,igroups  
        do i=1,ihead                     ! skip these header records
          read (luin_err,'(a1)') adum1
        enddo
!
! Check the following records and skip or read as required
        do i=1,igroups   
          read (luin_err,'(i2)') inum
          read (luin_err,'(1x,a5,1x,a4,1x,i3)') i_type,sat_name,n_channels
          if (i_type == d_type) then  ! interpret these required records
            ic=ic+1
            if ( (ic > n_err3 ) .or. (n_channels > n_err1) ) then
              print *,' '
              print *,'Problem reading sat_error file'
              print ('(i3,a,i3)'),ic,'=ic must not be > n_err3=',n_err3
              print ('(i3,a,i3)'),n_channels,'=n_channels must not be', &
                     ' > n_err1=', n_err1
              stop
            endif   
            do j=1,n_channels
              read (luin_err,'(1x,a18,2i5,f7.3)') adum,i1,i2,err_tab(j,2,ic)
              err_tab(j,1,ic)=j
            enddo
!
! Set err_itype to either sat ID or to instrument ID numbers
            if (sat_name == ' N14') err_itype(ic)=205         
            if (sat_name == ' N15') err_itype(ic)=206         
            if (sat_name == ' N16') err_itype(ic)=207         
            if (sat_name == ' N17') err_itype(ic)=208         
            if (sat_name == ' N18') err_itype(ic)=209         
            if ((d_type == 'AIRS_') .and. (n_channels == 301)) &
                 err_itype(ic)=420  ! instrument type ID as AIRS+AMSUA+HSB
            if ((d_type == 'AIRS_') .and. (n_channels == 15)) &
                 err_itype(ic)=570  ! instrument type ID as AMSUA
          else  ! skip these records
            do j=1,n_channels
              read (luin_err,'(a1)') adum1  
            enddo
          endif  ! test on d_type  
        enddo    ! loop over all groups of entries in table      
!
      endif      ! test on whther d_type is WIND_/MASS_ or Satellite 
      close (luin_err)
!
! Change some values in the error table
      if ( (d_type=='WIND_') .or. (d_type=='MASS_') ) then
        call fix_error_table (n_err1,n_err2,n_err3,err_tab)
      endif
!
     end subroutine read_error_table