program obs_format
!**********************************************************************
! *              .      .    .
! * PROGRAM:    obs_format
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2003-09-09
! *
! * ABSTRACT:   Prepare formated x1d file, suitable for plotting
! *
! * PROGRAM LOG:
! *
! * 09/09/2003    ..... M. ZUPANSKI: Original 'transform.F'
! * 10/17/2003    ..... M. ZUPANSKI: Innovation normalization
! * 10/22/2003    ..... M. ZUPANSKI: Innovation histogram
! * 10/22/2003    ..... M. ZUPANSKI: plotting
! * 12/05/2003    ..... M. ZUPANSKI: observation plotting
! *
! **********************************************************************
!-----
     integer,parameter::idim=20        ! model dimensions
     integer,parameter::in_file=21     ! unformatted x1d file #
     integer,parameter::iaddress=22    ! address file
     integer,parameter::out_file=51    ! formatted x1d file #
!-----
     integer :: N_obs
     integer :: N1_model,N2_model
     integer :: iobs,nnobs

     real,dimension(:),allocatable :: x_obs
     real,dimension(:,:),allocatable :: x1_plot,x2_plot,x3_plot
     integer,dimension(:),allocatable :: i1_loc,j1_loc
     integer,dimension(:),allocatable :: i2_loc,j2_loc
     integer,dimension(:),allocatable :: i3_loc,j3_loc
!-----
     integer, dimension(1) :: NNXP,NNYP,NNZP,NNSOIL,NNSNOW
     character(len=30)     :: filename
     character(len=9)      :: fgname1,fgname2,fgname3
!-----
      NAMELIST /MODEL_DIMENSION/ NNXP,NNYP,NNZP,NNSOIL,NNSNOW
!==============start calculation===================

!==============start calculation===================
!--- 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(*,*)' 500 OPEN UNIT ERROR IER=',IER
      READ(500,MODEL_DIMENSION)
      CLOSE(500,STATUS='KEEP')
      write(*,*) "NNXP=",NNXP(1)," NNYP=",NNYP(1)," NNZP=",NNZP(1)
!-----------------------------------------
         N1_model=NNXP(1)
         N2_model=NNYP(1)

 100  format(E20.10)
 110  format(i7)
 300  format(a9)

!-- read in unformatted obs file
     rewind in_file
     read(in_file) N_obs
     allocate(x_obs(1:N_obs))
     read(in_file) x_obs

      do i=1,min(10,N_obs)
      write(*,*) i," input obs=",x_obs(i)
      end do
    
!-- read in address file
     rewind iaddress
!- 1st
     read(iaddress) fgname1
     read(iaddress) icount1
     allocate(i1_loc(1:icount1)) ; allocate(j1_loc(1:icount1))
     read(iaddress) i1_loc(1:icount1),j1_loc(1:icount1)
     write(*,*) "READ: name=",fgname1," icount1=",icount1

!-- 2nd
     read(iaddress) fgname2
     read(iaddress) icount2
     allocate(i2_loc(1:icount2)) ; allocate(j2_loc(1:icount2))
     read(iaddress) i2_loc(1:icount2),j2_loc(1:icount2)
     write(*,*) "READ: name=",fgname2," icount2=",icount2

!-- 3rd
     read(iaddress) fgname3
     read(iaddress) icount3
     allocate(i3_loc(1:icount3)) ; allocate(j3_loc(1:icount3))
     read(iaddress) i3_loc(1:icount3),j3_loc(1:icount3)
     write(*,*) "READ: name=",fgname3," icount3=",icount3
    
!---- get obs at model points (h, psi, chi) -----

     allocate(x1_plot(1:N1_model,1:N2_model))
     allocate(x2_plot(1:N1_model,1:N2_model))
     allocate(x3_plot(1:N1_model,1:N2_model))

      x1_plot(:,:)=0.0
      x2_plot(:,:)=0.0
      x3_plot(:,:)=0.0

      nnobs=0
       do iobs=1,N_obs
       nnobs=nnobs+1
        if(nnobs.le.icount1) then
        icount_obs=nnobs
        x1_plot(i1_loc(icount_obs),j1_loc(icount_obs))=x_obs(iobs)
          else if (nnobs.le.icount1+icount2) then
          icount_obs=nnobs-icount1
          x2_plot(i2_loc(icount_obs),j2_loc(icount_obs))=x_obs(iobs)
            else 
            icount_obs=nnobs-(icount1+icount2)
            x3_plot(i3_loc(icount_obs),j3_loc(icount_obs))=x_obs(iobs)
         end if

        end do

      deallocate(x_obs) 
      deallocate(i1_loc) ; deallocate(j1_loc)
      deallocate(i2_loc) ; deallocate(j2_loc)
      deallocate(i3_loc) ; deallocate(j3_loc)

!----------------------------------
!-- write formatted file ('write_history_file' from ../ensda/src/model/ensda_programs_...F)

      REWIND out_file

      write(out_file,110) N1_model,N2_model
      write(out_file,300) fgname1
!-- height
        do i=1,N1_model
         do j=1,N2_model
          write(out_file,100) x1_plot(i,j)
         end do
        end do

      write(out_file,300) fgname2
!-- stream function
        do i=1,N1_model
         do j=1,N2_model
          write(out_file,100) x2_plot(i,j)
         end do
        end do

      write(out_file,300) fgname3
!-- velocity potential
        do i=1,N1_model
         do j=1,N2_model
          write(out_file,100) x3_plot(i,j)
         end do
        end do

       deallocate(x1_plot) ; deallocate(x2_plot) ; deallocate(x3_plot)

            stop
            end