program createobs
! **********************************************************************
! *              .      .    .
! * PROGRAM: createobs
! *
! * PRGMMR:    D. ZUPANSKI,M. Zupanski          ORG: CIRA/CSU   DATE: 2003-10-06
! *
! * ABSTRACT:  THIS PROGRAM READS FIRST GUESS VARIABLES AND 
! *            CREATES OBSERVATIONS BY ADDING RANDOM PERTURBATIONS
! *
! * INPUT FILES:
! * fguess_${OBSTYPE}.${rdate}
! * first_guess_list.h (include file)
! * first_guess
!
! * OUTPUT FILES:
! * obs_${OBSTYPE}.${rdate}
! *
! * PROGRAM LOG:
! *
! * 10/06/2003    ..... D. ZUPANSKI
! *
! **********************************************************************

      integer,parameter :: i_fguess=21         ! first guess file unit
      integer,parameter :: idev_input=25       ! input standard deviation
      integer,parameter :: i_obs=51            ! obs file unit
      integer,parameter :: i_loc=52            ! obs interpolation addresses file unit

      INCLUDE 'define_var_list_type.h' 
      INCLUDE 'define_fgvar_type.h'

      TYPE (var_list_type), dimension(:), allocatable :: fgvar_list
      TYPE (fgvar_type), dimension(:), allocatable :: fgvar
  
      integer, dimension(1) :: NNXP,NNYP,NNZP,NNSOIL,NNSNOW
      integer               :: fgvar_max
      character*20          :: filename
      character*9           :: fg_name,obs_name
      integer               :: N_obs
      integer               :: icycle
      real, dimension(:), allocatable :: obs,obs_err,stddev_input
      real, dimension(:), allocatable :: obs_x,obs_y
      real, dimension(:), allocatable :: fg1d,obs1d
      real, dimension(:,:), allocatable :: fg2d,obs2d
      real, dimension(:,:,:), allocatable :: fg3d,obs3d
      real, dimension(:,:,:,:), allocatable :: fg4d,obs4d
      character*3 H_obs

       integer,dimension(:),allocatable :: seed
       integer,parameter :: iseed=900

!------------------------------
      NAMELIST /MODEL_DIMENSION/ NNXP,NNYP,NNZP,NNSOIL,NNSNOW
      NAMELIST /HOBS/ H_obs
      NAMELIST /CURRENT_CYCLE/ icycle
!==============================================================
      write(*,*) " start createobs"

!--- get current cycle
      CLOSE(152)
      OPEN(UNIT=152,FILE='cycle.name',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(3,*)' 152 OPEN UNIT ERROR IER=',IER
      READ(152,CURRENT_CYCLE)
      CLOSE(152,STATUS='KEEP')
                                                                                
   if(icycle.eq.1) then
                                                                                
      write(*,*) "START-UP CYCLE: icycle=",icycle
      write(*,*) "NO NEED TO READ put_seed file"
                                                                                
   else

      write(*,*) "OTHER CYCLES: icycle=",icycle
      write(*,*) "READ put_seed file"
                                                                                
      CLOSE(iseed)
      OPEN(iseed,FILE='put_seed',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*) iseed,' OPEN UNIT ERROR IER=',IER
      READ(iseed,*) k
      allocate(seed(1:k))
      READ(iseed,*) seed(1:k)
      CLOSE(iseed,status='keep')
                                                                                
      call random_seed(put=seed(1:k))
      deallocate(seed)
                                                                                
   end if

!--- get observation operator type
      CLOSE(201)
      OPEN(UNIT=201,FILE='H_obs.parm',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(3,*)' 201 OPEN UNIT ERROR IER=',IER
      READ(201,HOBS)
      CLOSE(201,STATUS='KEEP')
      write(*,*) "H_obs=",H_obs

!------------------------------
!--- get model dimensions 
      write(filename,200)
 200  format('model_dimension.name')
      CLOSE(500)
      OPEN(UNIT=500,FILE=filename,FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(3,*)' 500 OPEN UNIT ERROR IER=',IER
      READ(500,MODEL_DIMENSION)  
      CLOSE(500,STATUS='KEEP')
      write(*,*) "NNXP=",NNXP(1)," NNYP=",NNYP(1)," NNZP=",NNZP(1)

      INCLUDE 'first_guess_list.h' 

      write(*,*) "max_num_of_first_guess=",max_num_of_first_guess
       do nfg=1,max_num_of_first_guess
      write(*,*) "ndim ",fgvar_list(nfg)%ndim
      write(*,*) "start1 ",fgvar_list(nfg)%start_index(1)
      write(*,*) "start2 ",fgvar_list(nfg)%start_index(2)
      write(*,*) "start3 ",fgvar_list(nfg)%start_index(3)
      write(*,*) "start4 ",fgvar_list(nfg)%start_index(4)
      write(*,*) "end1 ",fgvar_list(nfg)%end_index(1)
      write(*,*) "end2 ",fgvar_list(nfg)%end_index(2)
      write(*,*) "end3 ",fgvar_list(nfg)%end_index(3)
      write(*,*) "end4 ",fgvar_list(nfg)%end_index(4)
      write(*,*) "name ",fgvar_list(nfg)%name
      write(*,*) "description ",fgvar_list(nfg)%description
       end do

      CLOSE(105)
      OPEN(105,FILE='first_guess_anlsis',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' 105 OPEN UNIT ERROR IER=',IER
      READ(105,*) fgvar_max
      allocate(fgvar(1:fgvar_max))
       do i=1,fgvar_max
      READ(105,*) fg_name
        fgvar(i)%name=fg_name
      write(*,*) "fg_name=",fgvar(i)%name
       end do
      CLOSE(105,STATUS='KEEP')

        write(filename,1020)
 1020 format('obs_std_dev_input')
      write(*,*) "input: obs_std_dev_input filename=",filename
      CLOSE(idev_input)
      OPEN(UNIT=idev_input,FILE=filename,FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) idev_input,' OPEN UNIT ERROR IER=',IER
!----------------------------------------------
        N_obs=0
       do i=1,fgvar_max
        do j=1,max_num_of_first_guess
        if(fgvar(i)%name.eq.fgvar_list(j)%name) then

         i1s=fgvar_list(j)%start_index(1)
         i2s=fgvar_list(j)%start_index(2)
         i3s=fgvar_list(j)%start_index(3)
         i4s=fgvar_list(j)%start_index(4)
         i1e=fgvar_list(j)%end_index(1)  
         i2e=fgvar_list(j)%end_index(2)  
         i3e=fgvar_list(j)%end_index(3)  
         i4e=fgvar_list(j)%end_index(4)  

         i4dim=(i1e-i1s+1)*(i2e-i2s+1)*(i3e-i3s+1)*(i4e-i4s+1)
         i3dim=(i1e-i1s+1)*(i2e-i2s+1)*(i3e-i3s+1)
         i2dim=(i1e-i1s+1)*(i2e-i2s+1)
         i1dim= i1e-i1s+1
        write(*,*) " i4dim,i3dim,i2dim,i1dim=",i4dim,i3dim,i2dim,i1dim

        if(fgvar_list(j)%ndim.eq.4) then
        N_obs=N_obs + i4dim
          else if(fgvar_list(j)%ndim.eq.3) then
        N_obs=N_obs + i3dim
        write(*,*) " i=",i," j=",j," N_obs=",N_obs
            else if(fgvar_list(j)%ndim.eq.2) then
        N_obs=N_obs + i2dim
             else if(fgvar_list(j)%ndim.eq.1) then
        N_obs=N_obs + i1dim
          else
     write(*,*) "FGUESS DOES NOT HAVE APPROPRIATE ndim=",fgvar_list(j)%ndim
        end if
!----------------------------------------------
        end if      !!!if(fgvar(i)%name.eq.fgvar_list(j)%name) then
        end do
       end do
      write(*,*) "N_obs=",N_obs
!===================================================
!======================================
     allocate(obs(1:N_obs))
     allocate(obs_err(1:N_obs))
     allocate(obs_x(1:N_obs))
     allocate(obs_y(1:N_obs))
!======================================
!======================================
!======================================

     rewind i_fguess
     rewind i_obs
     rewind i_loc

! 100 format(e20.10)

        iobs=0
        nvar=0
        icount=0
       do i=1,fgvar_max
       write(*,*) "do i=",i
        do j=1,max_num_of_first_guess
        write(*,*) "do j=",j
        if(fgvar(i)%name.eq.fgvar_list(j)%name) then
        nvar=nvar+1

         i1s=fgvar_list(j)%start_index(1)
         i2s=fgvar_list(j)%start_index(2)
         i3s=fgvar_list(j)%start_index(3)
         i4s=fgvar_list(j)%start_index(4)
         i1e=fgvar_list(j)%end_index(1)  
         i2e=fgvar_list(j)%end_index(2)  
         i3e=fgvar_list(j)%end_index(3)  
         i4e=fgvar_list(j)%end_index(4)  

!-- write first guess on unit i_fguess

        fg_name=fgvar(i)%name
        obs_err0=(fgvar_list(j)%stddev)
        if(fgvar_list(j)%ndim.eq.4) then

       call H_X_4d &
      (N_obs,iobs,i_loc,i_fguess,i1s,i1e,i2s,i2e,i3s,i3e,i4s,i4e, &
       H_obs,obs_err0,obs,obs_err,obs_x,obs_y)

          else if(fgvar_list(j)%ndim.eq.3) then

          allocate(stddev_input(i3s:i3e))
           rewind idev_input
                      do 
           read(idev_input,*,end=220) fg_name
           read(idev_input,*) stddev_input
               if(fg_name.eq.fgvar(i)%name) then
           exit
               else
           cycle
               end if
                      end do
           if(fg_name.eq.'u        ') then
            stddev_input(:)=obs_err0
            stddev_input(:)=stddev_input(:)*1.8
             else if(fg_name.eq.'v        ') then
            stddev_input(:)=obs_err0
            stddev_input(:)=stddev_input(:)*1.8
               else if(fg_name.eq.'w        ') then
            stddev_input(:)=obs_err0
            stddev_input(:)=stddev_input(:)*1.8
                 else if(fg_name.eq.'exnr     ') then
            stddev_input(:)=obs_err0
            stddev_input(:)=stddev_input(:)*1.8
                   else if(fg_name.eq.'thetail  ') then
            stddev_input(:)=obs_err0
            stddev_input(:)=stddev_input(:)*1.8
                   else if(fg_name.eq.'r_total  ') then
            stddev_input(:)=stddev_input(:)*1.0
            stddev_input(:)=stddev_input(:)*1.8
                   endif
           write(*,*) " fg_name=",fg_name
           write(*,*) " stddev_input=",stddev_input
       call H_X_3d &
      (N_obs,iobs,i_loc,i_fguess,i1s,i1e,i2s,i2e,i3s,i3e, &
       H_obs,obs_err0,obs,obs_err,stddev_input,obs_x,obs_y)

          deallocate(stddev_input)
            else if(fgvar_list(j)%ndim.eq.2) then

       call H_X_2d &
      (N_obs,iobs,i_loc,i_fguess,i1s,i1e,i2s,i2e, &
       H_obs,obs_err0,obs,obs_err,obs_x,obs_y)

             else if(fgvar_list(j)%ndim.eq.1) then

       call H_X_1d &
      (N_obs,iobs,i_loc,i_fguess,i1s,i1e, &
       H_obs,obs_err0,obs,obs_err,obs_x,obs_y)

          else
     write(*,*) "FGUESS DOES NOT HAVE APPROPRIATE ndim=",fgvar_list(j)%ndim
          end if
!----------------------------------------------
        end if      !!!if(fgvar(i)%name.eq.fgvar_list(j)%name) then
        end do
       end do
!----------------------------------------------
!-- write out residuals

     write(*,*) "BEFORE WRITE OBS: iobs=",iobs," N_obs=",N_obs

     rewind i_obs
     WRITE(i_obs) iobs
!    WRITE(i_obs) obs(1:iobs),obs_err(1:iobs)
     WRITE(i_obs) obs(1:iobs),obs_err(1:iobs),obs_x(1:iobs),obs_y(1:iobs)

     deallocate(obs)   ; deallocate(obs_err)
     deallocate(obs_x) ; deallocate(obs_y)

!--- save seed---------------------
                                                                                
      call random_seed(size=k)
      allocate(seed(1:k))
      call random_seed(get=seed(1:k))
                                                                                
      CLOSE(iseed)
      OPEN(iseed,FILE='put_seed',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*) iseed,' OPEN UNIT ERROR IER=',IER
      WRITE(iseed,*) k
      WRITE(iseed,*) seed(1:k)
      CLOSE(iseed,status='keep')
                                                                                
      deallocate(seed)

     write(*,*) " end createobs"
     stop

   220  write(*,*) " PROBLEM: fg_name on idev_input is not appropriate"

     write(*,*) " end createobs"
     stop

      end program createobs