program projectRMM
      implicit none
 
c Originally written by Matthew Wheeler in ~2003
c Modified in February 2008 to work on ASCII input data (for the EOF sturctures)

c This program computes the RMM1 and RMM2 indices for any day (or series
c of days) that we have *all* sources of data for.
c
c The required input data is daily-averaged anomalies of u850, u200, and olr.
c These fields must have already had the seasonal cycle removed as well as 
c the mean of the previous 120 days (as an estimate of interannual variability).
c They should be averages from 15S-15N on a 2.5 degree grid (i.e. 144 points).
c The input data should be centred on the longitudes 0deg, 2.5deg, ....., 357.5deg.
c 
c We must also read-in the file that contains the EOFs so that we can make
c the projection onto those structures.
c
c No missing data is allowed.

c  Variables for input EOF file
      integer NX,TNX
      parameter (NX=144,TNX=3*NX)
      real eigval(2),eigvec(TNX,2),norm(TNX)
      integer imd,isp,it
      
c  Input data fields
      real olr(NX),u850(NX),u200(NX)
c  Other variables
      integer fstyr,fstmo,fstda      !day to start producing RMM1 and RMM2
      integer nyr,nmo,nda1,nda2,nda3,i
      real datmat(TNX),pc(2)
      character*(*) inpute,input1,input2,input3

c INPUT EOF FILE    (for the spatial structure to project onto)
      parameter(inpute='WH04_EOFstruc.txt')

c INPUT DATA FILES  (daily anomalies with interannual removed)
      parameter(input1='olr15av.anom-120dm.98toRealtime.1x.b')
      parameter(input2=
     # 'u85015av.anom-120dm.02toRealtime.1x.b')
      parameter(input3=
     # 'u20015av.anom-120dm.02toRealtime.1x.b')

      open(unit=11,file=inpute,status='old')

      open(unit=21,file='/bm/gdata/mwheeler/maproom/OLR/'//input1,
     #     form='unformatted',status='old')
      open(unit=22,file='/bm/gdata/mwheeler/maproom/NCEP+GASP/'//input2,
     #     form='unformatted',status='old')
      open(unit=23,file='/bm/gdata/mwheeler/maproom/NCEP+GASP/'//input3,
     #     form='unformatted',status='old')

c OUTPUT FILE
      open(15,file='RMMprojections.02toRealtime.txt')

c   Read in the EOF spatial structures
c   ----------------------------------

      read(11,*)
      read(11,*)
      read(11,*)
      read(11,*)
      read(11,*) (eigval(imd),imd=1,2)  ! eigenvalues
      read(11,*)
      read(11,*)
      read(11,*)
      read(11,*)
      do isp=1,TNX
       read(11,*) (eigvec(isp,imd),imd=1,2)  ! eigenvectors
      enddo
      read(11,*)
      read(11,*)
      do isp=1,TNX
       read(11,*) norm(isp)     ! normalization factor 
      enddo

c --------------------------------------------------------------------
c    Advance field data to first day
c    -------------------------------
c As wind data starts later, get first day from that file.....
 61   read(22) nyr,nmo,nda2    ! year, month, day of u850 data
      read(22) (u850(i),i=1,NX)
      fstyr=nyr
      fstmo=nmo
      fstda=nda2

c    get to the starting year in each other file.
 62   read(21) nyr,nmo,nda1    ! year month day of OLR data.
      read(21) (olr(i),i=1,NX)
      if(nyr.ne.fstyr) goto 62
      if(nmo.ne.fstmo) goto 62
      if(nda1.ne.fstda) goto 62

 63   read(23) nyr,nmo,nda3    ! u200 data
      read(23) (u200(i),i=1,NX)
      if(nyr.ne.fstyr) goto 63
      if(nmo.ne.fstmo) goto 63
      if(nda3.ne.fstda) goto 63

c --------------------------------------------------------------------
      it=0

 66   CONTINUE     !!!!!!!!!!!!!!!!!!!!!!!!!!! BIG LOOP TO HERE !!!!!!!
        it=it+1       

c       read-in the OLR, u850 and u200.
        if (it.ne.1) then
         read(21,end=888) nyr,nmo,nda1
         read(21) (olr(i),i=1,NX)
         read(22,end=888) nyr,nmo,nda2
         read(22) (u850(i),i=1,NX)
         read(23,end=888) nyr,nmo,nda3
         read(23) (u200(i),i=1,NX)
        endif

c       check that the dates match
        if(nda1.ne.nda2.or.nda1.ne.nda3) then
         print*,'Something wrong with dates between fields'
         print*,'nda1=',nda1,' nda2=',nda2,' nda3=',nda3
         STOP
        endif

c      -------
c      Create extended vector of OLR, u850, and u200 data.
        DO i=1,TNX
          if(i.le.NX) then
            datmat(i)=olr(i)
          elseif(i.le.(NX*2)) then
            datmat(i)=u850(i-NX)
          elseif(i.le.(NX*3)) then
            datmat(i)=u200(i-(NX*2))
          endif
c        Divide datmat by the "norm" as calculated in the EOF program.
         datmat(i)=datmat(i)/norm(i) ! I normalize by the "global" s.d. of each field.
        ENDDO
c      -------

c  Compute RMM1 and RMM2  (the first two normalized PCs)
c  -----------------------------------------------------
       do 300 imd=1,2
        pc(imd)=0.
         do isp=1,TNX
          pc(imd)=pc(imd)+(eigvec(isp,imd)*datmat(isp))
         enddo
  300  continue
c   Now normalize (by EOF-calcualted s.d.) the newly calculated PCs
       do imd=1,2
        pc(imd)=pc(imd)/sqrt(eigval(imd))
       enddo


c   WRITE OUT THE DESIRED RMM VALUES
c   --------------------------------
       print*,'Writing to file',nyr,nmo,nda1
       write(15,*) nyr,nmo,nda1,pc(1),pc(2)

      GOTO 66 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 888  print*,'last input/output on ',nyr,nmo,' and day',nda1,
     # ' or ',nda2,' or ',nda3

      END   !main program