program rms_SWM
!**********************************************************************
! *              .      .    .
! * PROGRAM:    rms_SWM
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2003-09-09
! *
! * ABSTRACT:   Prepare formated 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
! * 10/24/2003    ..... M. ZUPANSKI: Root-Mean-Squared error calculation
! * 05/21/2004    ..... M. ZUPANSKI: Add shallow-water model (CSU)
! *
! **********************************************************************
!-----
     integer,parameter::in_file1=20     ! input file1 #
     integer,parameter::in_file2=21     ! input file2 #
     integer,parameter::out_rms=51      ! formatted RMS file #
!-----
     integer, dimension(1) :: NNXP,NNYP,NNZP,NNSOIL,NNSNOW
     character(len=30)     :: filename

     real :: rms_h,rms_psi,rms_chi

     real,dimension(:,:),allocatable::h_1,h_2 
     real,dimension(:,:),allocatable::psi_1,psi_2 
     real,dimension(:,:),allocatable::chi_1,chi_2 
     integer :: im,jm
     integer :: imx,jmx
     character(len=10) :: namex
!-----
      NAMELIST /MODEL_DIMENSION/ NNXP,NNYP,NNZP,NNSOIL,NNSNOW
!==============start calculation===================

!--- get model dimensions
      write(filename,1000)
 1000 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)
!-----------------------------------------
         im=NNXP(1)
         jm=NNYP(1)

     write(*,*) "im,jm=",im,jm
!-----------------------------------------

     allocate(h_1(1:im,1:jm))
     allocate(h_2(1:im,1:jm))
     allocate(psi_1(1:im,1:jm))
     allocate(psi_2(1:im,1:jm))
     allocate(chi_1(1:im,1:jm))
     allocate(chi_2(1:im,1:jm))

 100  format(E20.10)
 110  format(2i7)
 200  format(3E20.10)
 300 format(a10)

!-- read in first file
     rewind in_file1
     read(in_file1,110) imx,jmx
      if(imx.ne.im.or.jmx.ne.jm) then
      write(*,*) "RMS: Bad indexes im,jm=",im,jm," imx,jmx=",imx,jmx
      end if

!-- height
     read(in_file1,300) namex
     write(*,*) "READ variable ",namex
        do i=1,im
         do j=1,jm
          read(in_file1,100) h_1(i,j)
         end do
        end do

!-- Stream function
     read(in_file1,300) namex
     write(*,*) "READ variable ",namex
        do i=1,im
         do j=1,jm
          read(in_file1,100) psi_1(i,j)
         end do
        end do

!-- velocity potential
     read(in_file1,300) namex
     write(*,*) "READ variable ",namex
        do i=1,im
         do j=1,jm
          read(in_file1,100) chi_1(i,j)
         end do
        end do
!-------------------------------------------------------

!-- read in second file
     rewind in_file2
     read(in_file2,110) imx,jmx
      if(imx.ne.im.or.jmx.ne.jm) then
      write(*,*) "RMS: Bad indexes im,jm=",im,jm," imx,jmx=",imx,jmx
      end if

!-- height
     read(in_file2,300) namex
     write(*,*) "READ variable ",namex
        do i=1,im
         do j=1,jm
          read(in_file2,100) h_2(i,j)
         end do
        end do

!-- Stream function
     read(in_file2,300) namex
     write(*,*) "READ variable ",namex
        do i=1,im
         do j=1,jm
          read(in_file2,100) psi_2(i,j)
         end do
        end do

!-- velocity potential
     read(in_file2,300) namex
     write(*,*) "READ variable ",namex
        do i=1,im
         do j=1,jm
          read(in_file2,100) chi_2(i,j)
         end do
        end do
!-------------------------------------------------------


!-- calculate rms error

      rms_h=0.0
      rms_psi=0.0
      rms_chi=0.0
      do i=1,im
       do j=1,jm
         rms_h  =rms_h  +(h_1(i,j)  -h_2(i,j)  )**2
         rms_psi=rms_psi+(psi_1(i,j)-psi_2(i,j))**2
         rms_chi=rms_chi+(chi_1(i,j)-chi_2(i,j))**2
       end do
      end do
      rms_h=sqrt(rms_h/float(im*jm))
      rms_psi=sqrt(rms_psi/float(im*jm))
      rms_chi=sqrt(rms_chi/float(im*jm))

!-- write formatted file
     rewind out_rms
     write(out_rms,200) rms_h,rms_psi,rms_chi

            stop
            end