program rms_GEOS
!**********************************************************************
! *              .      .    .
! * PROGRAM:    rms_GEOS
! * PRGMMR:     D. ZUPANSKI     ORG: CIRA/CSU   DATE: 2004-01-16
! *
! * 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
! * 10/24/2003    ..... M. ZUPANSKI: Root-Mean-Squared error calculation
! * 01/16/2004    ..... D. ZUPANSKI: Root-Mean-Squared error calculation
!                       for GEOS column model
! *
! **********************************************************************
!-----
     integer,parameter::in_file1=20     ! input file1 #
     integer,parameter::in_file2=21     ! input file2 #
     integer,parameter::out_rms=51      ! formatted RMS file #
!-----
     integer :: im,jm,nlayr
     integer, dimension(1) :: NNXP,NNYP,NNZP,NNSOIL,NNSNOW
     character(len=30)     :: filename

     real :: rms_pt,rms_q
     real,dimension(:,:,:),allocatable::PT_hist_1,Q_hist_1
     real,dimension(:,:,:),allocatable::PT_hist_2,Q_hist_2
     logical :: global
!-----
      NAMELIST /MODEL_DIMENSION/ NNXP,NNYP,NNZP,NNSOIL,NNSNOW
!==============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)
!-----------------------------------------
         im=NNXP(1)
         jm=NNYP(1)
         nlayr=NNZP(1)
         GLOBAL=.false.

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

 100  format(E20.10)

              if(GLOBAL)then
           allocate(PT_hist_1(im,jm,nlayr))
           allocate(Q_hist_1(im,jm,nlayr))
           allocate(PT_hist_2(im,jm,nlayr))
           allocate(Q_hist_2(im,jm,nlayr))
              else
           allocate(PT_hist_1(1,1,nlayr))
           allocate(Q_hist_1(1,1,nlayr))
           allocate(PT_hist_2(1,1,nlayr))
           allocate(Q_hist_2(1,1,nlayr))
              endif

!-- read in first file
     rewind in_file1
     read(in_file1,100) PT_hist_1
     read(in_file1,100) Q_hist_1

!-- read in second file
     rewind in_file2
     read(in_file2,100) PT_hist_2
     read(in_file2,100) Q_hist_2

!-- calculate rms error

      rms_pt=0.0
      rms_q=0.0
      do l=1,nlayr
       do j=1,jm
        do i=1,im
      rms_pt=rms_pt+(PT_hist_1(i,j,l)-PT_hist_2(i,j,l))**2
      rms_q=rms_q+(Q_hist_1(i,j,l)-Q_hist_2(i,j,l))**2
        end do
       end do
      end do
      rms_pt=sqrt(rms_pt/float(im*jm*nlayr))
      rms_q=sqrt(rms_q/float(im*jm*nlayr))

!-- write formatted file
     rewind out_rms
     write(out_rms,100) rms_pt
     write(out_rms,100) rms_q

            stop
            end