kupdate.f90

! ----------------------------------------------------------------------
!
  SUBROUTINE KUPDATE(Y,Z_IN,ZUSE,CV_IN,ZMEAS_IN,V_OUT,N_RL,N_Y,N_Z,N_R,&
    MY_COMM,RANK,ERR,M,MEAS_SWITCH,RANGE_U,RANGE_A,N_ZR,N_ZP,N_YP,N_YR,&
    INNOV_2MOM,TR_ERR,ABS_TR_ERR,REL_L1_NORM,NU_BAR)
!
! ----------------------------------------------------------------------
!
!  COMPUTES THE POSTERIOR STATE VARIABLES BASED ON THE PRIOR STATE, 
!  THE PREDICTED OBSERVATIONS AND THE ACTUAL OBSERVATION.  THIS VERSION
!  OF THE CODE OVERWRITES THE Y MATRIX DURING THE UPDATE; THE PRIOR STATE
!  MUST BE STORED BEFORE THE CALL TO KUPDATE.  THIS CODE CALLS A MATRIX 
!  INVERSION ROUTINE THAT IS BASED ON AN L-U FACTORIZATION.  IN THIS MPI
!  VERSION OF THIS CODE, THE MEANS ARE FIRST COMPUTED BY REDUCING ACROSS
!  ALL PROCESSES USING THE MPI_SUM FUNCTION.  THE MEANS ARE THEN USED ON 
!  EACH INDIVIDUAL PROCESS IN ORDER TO COMPUTE THE LOCAL DEVIANCE FROM THE
!  MEAN (TILDA VARIABLES).  THESE TILDA VARIABLES ARE USED TO COMPUTE THE
!  LOCAL SUM OF THE SQUARE DEVIANCE (SSE) FROM THE MEAN.  THESE ARE THEN 
!  REDUCED ACROSS ALL PROCESSES IN ORDER TO OBTAIN THE GLOBAL SSE FOR THE 
!  CALCULATTION OF THE COVARIANCE OF Z WITH Z AND OF Y WITH Z.  THIS 
!  CALCULATION IS CARRIED OUT ON THE 0 NODE, AND THE KALMAN GAIN IS COMPUTED
!  THERE, THEN BROADCAST TO ALL PROCESSES.  THE UPDATE TO THE LOCAL STATES
!  IS THEN CARRIED OUT ON EACH INDIVIDUAL PROCESSOR.
!
!  Y -IN/OUTPUT OF PRIOR/POSTERIOR STATE VARIABLES OVERWRITTEN DURING THE UPDATE
!  Z_IN - VECTOR OF PREDICTED OBSERVATIONS FOR EACH REPLICATE.  THIS IS REDUCED 
!    TO Z BASED ON THE ZUSE VECTOR 
!  ZUSE - SPECIFIES WHICH OF THE MEASUREMENTS ARE BEING USED DURING THIS UPDATE
!  CV_IN - MEASUREMENT ERROR COVARIANCE.  THIS IS REDUCED TO CV BASED ON THE 
!    ZUSE VECTOR
!  ZMEAS_IN - THE VECTOR OF OBSERVATIONS.  THIS IS REDUCED TO Z BASED ON 
!  V_OUT - OUTPUT RANDOM ERROR WHICH IS GENERATED TO KEEP REPLICATES FROM
!    DEVELOPING CORRELATION AMONG THEMSELVES
!  N_RL - NUMBER OF REPLICATES ON EACH PROCESS
!  N_Y - TOTAL NUMBER OF STATES
!  N_Z - TOTAL NUMBER OF MEASUREMENTS
!  N_R - TOTAL NUMBER OF REPLICATES
!  MY_COMM - LOCAL COPY OF GLOBAL MPI_COMM_WORLD COMMUNICATOR
!  RANK - RANK OF PROCESS RUNNING THIS SUBROUTINE
!  ERR - GLOBAL ERROR VARIABLE
!  M - MEASUREMENT INTERVAL NUMBER, NEEDED TO SEED RANDOM NUMBER GENERATOR
!
!  CALLS: INVERT_MATRIX,MVNRND
!  CALLED FROM: ENKF
!
!  VERSION HISTORY:
!  WRITTEN - MD - MAY 2005
!  UPDATED FOR MPI - MD - AUG 2005

IMPLICIT NONE
INCLUDE 'mpif.h'

INTEGER,INTENT(IN)::N_Y,N_RL,N_Z,N_R,M,MEAS_SWITCH,N_ZR,N_YP,N_YR
INTEGER,INTENT(IN) ::  ZUSE(N_Z),N_ZP(N_ZR)
INTEGER,INTENT(IN)::MY_COMM,RANK,ERR
REAL,INTENT(IN)::Z_IN(N_Z,N_RL),CV_IN(N_Z,N_Z),ZMEAS_IN(N_Z),RANGE_U,RANGE_A
REAL,INTENT(INOUT)::Y(N_Y,N_RL)
REAL,INTENT(INOUT)::V_OUT(N_Z,N_R)
REAL,INTENT(OUT) :: REL_L1_NORM,ABS_TR_ERR,TR_ERR,NU_BAR,INNOV_2MOM
INTEGER I,ZSUM,J,N_Z2(1),SZMU0(2),SZCV(2),K,SEED,KGLOBAL,YI,YF,P,O
INTEGER,DIMENSION(:),ALLOCATABLE::IZ,zuse_local
REAL :: DIST(N_YP,N_YP), C,ZD,TR_CNUNU,TR_CZZ,&
  L1_NORM_DIFF,L1_NORM_EXP,COL_SUM
REAL,DIMENSION(:,:),ALLOCATABLE::Z,CYZ,CZZ,CV,CZZCV,KALM,CZZCVINV,V,MU0,&
  ZTILDA,LOCALSSEZZ,GLOBALSSEZZ,LOCALSSEYZ,GLOBALSSEYZ,YTILDA,V_TRANS,C0,&
  LOCALSSENUNU,GLOBALSSENUNU,CNUNU
REAL,DIMENSION(:),ALLOCATABLE::ZMEAS,LOCALSUMZ,GLOBALSUMZ,ZMEAN,YMEAN,&
  LOCALSUMY,GLOBALSUMY,NU,NUTCI,LOCALSUMNU,GLOBALSUMNU
!REAL NODATA,YMEAN(N_Y),LOCALSUMY(N_Y),GLOBALSUMY(N_Y),YTILDA(N_Y,N_RL)
REAL NODATA,mymax(n_y),mymin(n_y)
integer :: imax,imin
CHARACTER(32) N_YPC

ALLOCATE(YMEAN(N_Y),LOCALSUMY(N_Y),GLOBALSUMY(N_Y),YTILDA(N_Y,N_RL))

!DEFINE NO DATA
NODATA=-9999.

!  READ DISTANCE MATRIX FROM FILE
OPEN(UNIT=2,FILE='dist.in',STATUS='OLD')
!if(rank.eq.0) print *, 'n_yp=',n_yp
WRITE(N_YPC,*) N_YP !FORM A CHARACTER VARIABLE WITH VALUE OF N_YP
DO I=1,N_YP
  READ(2,'('// N_YPC // '(F10.4))') (DIST(I,J),J=1,N_YP)
END DO
CLOSE(2)

ALLOCATE(ZUSE_LOCAL(N_Z))
ZUSE_LOCAL=ZUSE

! DETERMINE WHICH MEASUREMENTS TO USE...
IF(MEAS_SWITCH.EQ.0)THEN !USE ALL MEASUREMENTS
  ZUSE_LOCAL=1
ELSEIF(MEAS_SWITCH.EQ.1)THEN !USE PM MEASUREMENTS ONLY
  ZUSE_LOCAL=0
  ZUSE_LOCAL(1:12)=1
ELSEIF(MEAS_SWITCH.EQ.2)THEN !USE PM+ALBEDO MEASUREMENTS ONLY
  ZUSE_LOCAL=0
  ZUSE_LOCAL(1:637)=1
ELSEIF(MEAS_SWITCH.EQ.3)THEN !USE ALBEDO MEASUREMENTS ONLY
  ZUSE_LOCAL=0
  ZUSE_LOCAL(13:637)=1
ELSEIF(MEAS_SWITCH.EQ.9)THEN !NO UPDATES
  ZUSE_LOCAL=0.
ELSE
  ! STOP IF THE SWITCH IS NOT EXPECTED VALUES
  PRINT *, 'ILLEGAL VALUE OF MEAS_SWITCH, ABORTING!!'
  CALL MPI_FINALIZE(ERR)
  STOP
END IF

ZSUM=0
DO I=1,N_Z
  IF (ZMEAS_IN(I).EQ.NODATA.OR.ZUSE_local(I).EQ.0) ZSUM=ZSUM+1
  IF (ZMEAS_IN(I).EQ.NODATA) THEN
    ZUSE_local(I)=0
  END IF
END DO

IF (ZSUM.EQ.N_Z) THEN
  !DON'T UPDATE Y AND GIVE A DUMMY VALUE TO V_OUT, OTHER OUTPUTS
  V_OUT=0.
  REL_L1_NORM=0.
  ABS_TR_ERR=0.
  TR_ERR=0.
  NU_BAR=0.
  INNOV_2MOM=0.
ELSE
  !1) RESHAPE ARRAYS BASED ON ZUSE INFO
  ALLOCATE(IZ(SUM(ZUSE_local)))
  J=0
  DO I=1,N_Z
    IF (ZUSE_local(I).EQ.1) THEN
      J=J+1
      IZ(J)=I
    END IF
  END DO
  N_Z2=SHAPE(IZ)  ! THIS NOW REPLACES N_Z ARRAY BASED ON NEW SIZE
                  ! AFTER RESHAPING ARRAY BASED ON ZUSE

  ! 2) ALLOCATION BASED ON RESHAPED DATA
  ALLOCATE(Z(N_Z2(1),N_RL),ZMEAN(N_Z2(1)),ZMEAS(N_Z2(1)),&
    CYZ(N_Y,N_Z2(1)),CZZ(N_Z2(1),N_Z2(1)),CV(N_Z2(1),N_Z2(1)),&
    CZZCV(N_Z2(1),N_Z2(1)),KALM(N_Y,N_Z2(1)),CZZCVINV(N_Z2(1),N_Z2(1)),&
    V(N_Z2(1),N_R),V_TRANS(N_R,N_Z2(1)),MU0(N_Z2(1),1),LOCALSUMZ(N_Z2(1)),&
    GLOBALSUMZ(N_Z2(1)),LOCALSSEZZ(N_Z2(1),N_Z2(1)),ZTILDA(N_Z2(1),N_RL),&
    GLOBALSSEZZ(N_Z2(1),N_Z2(1)),LOCALSSEYZ(N_Y,N_Z2(1)),&
    GLOBALSSEYZ(N_Y,N_Z2(1)),LOCALSSENUNU(N_Z2(1),N_Z2(1)),&
    GLOBALSSENUNU(N_Z2(1),N_Z2(1)),CNUNU(N_Z2(1),N_Z2(1)),LOCALSUMNU(N_Z2(1)),&
    GLOBALSUMNU(N_Z2(1)) )

  ! 3) RESHAPE Z,ZMEAS, AND CV BASED ON ZUSE
  CV=0
  DO I=1,N_Z2(1)
    J=IZ(I)
    Z(I,:)=Z_IN(J,:)
    ZMEAS(I)=ZMEAS_IN(J)
    CV(I,I)=CV_IN(J,J)
  END DO

  ! 4) COMPUTE LOCAL AND GLOBAL TILDA (DEVIATION FROM THE RESPECTIVE MEAN) 
  !      BY REDUCING ACROSS ALL PROCESSES

  ! COMPUTE LOCAL SUM OF Y AND Z VARIABLES
  LOCALSUMY=0.
  DO I=1,N_Y
    DO K=1,N_RL
      LOCALSUMY(I)=LOCALSUMY(I)+Y(I,K)
    END DO
  END DO
  LOCALSUMZ=0.
  DO I=1,N_Z2(1)
    DO K=1,N_RL
      LOCALSUMZ(I)=LOCALSUMZ(I)+Z(I,K)
    END DO
  END DO

  ! COMPUTE GLOBAL SUMS BY REDUCING LOCAL SUMS ACROSS ALL PROCESSES
  CALL MPI_ALLREDUCE(LOCALSUMY,GLOBALSUMY,N_Y,MPI_REAL,MPI_SUM,MY_COMM,ERR)
  CALL MPI_ALLREDUCE(LOCALSUMZ,GLOBALSUMZ,N_Z2(1),MPI_REAL,MPI_SUM,MY_COMM,ERR)

  YMEAN=GLOBALSUMY/N_R
  ZMEAN=GLOBALSUMZ/N_R

  ! COMPUTE VARIATION FROM THE RESPECTIVE MEANS
  DO I=1,N_Y
    DO K=1,N_RL
      YTILDA(I,K)=Y(I,K)-YMEAN(I)
    END DO
  END DO
  DO I=1,N_Z2(1)
    DO K=1,N_RL
      ZTILDA(I,K)=Z(I,K)-ZMEAN(I)
    END DO
  END DO

  ! 5) COMPUTE CZZ BY REDUCING A LOCAL SUM OF SQUARED DEVIANCE FROM THE MEAN
  !      ACROSS ALL PROCESSES
  LOCALSSEZZ=0.0
  DO I=1,N_Z2(1)
    DO J=1,N_Z2(1)
      DO K=1,N_RL
        LOCALSSEZZ(I,J)=LOCALSSEZZ(I,J)+ZTILDA(I,K)*ZTILDA(J,K)
      END DO
    END DO
  END DO
  CALL MPI_REDUCE(LOCALSSEZZ,GLOBALSSEZZ,N_Z2(1)*N_Z2(1),MPI_REAL,MPI_SUM,&
    0,MY_COMM,ERR)
  CZZ=GLOBALSSEZZ/(N_R-1)

  ! 6) COMPUTE CYZ BY REDUCING A LOCAL SUM OF THE PRODUCT OF Y AND Z 
  !      DEVIANCE FROM THEIR RESPECTIVE MEANS ACROSS ALL PROCESSES
  LOCALSSEYZ=0.0
  DO I=1,N_Y
    DO J=1,N_Z2(1)
      DO K=1,N_RL
        LOCALSSEYZ(I,J)=LOCALSSEYZ(I,J)+YTILDA(I,K)*ZTILDA(J,K)
      END DO
    END DO
  END DO
  CALL MPI_REDUCE(LOCALSSEYZ,GLOBALSSEYZ,N_Y*N_Z2(1),MPI_REAL,MPI_SUM,&
    0,MY_COMM,ERR)
  CYZ=GLOBALSSEYZ/(N_R-1)

  !IMPLEMENTATION OF SMOOTHING OF COVARIANCE MATRICES BY SHUR PRODUCT WITH
  ! COMPACTLY SUPPORTED COVARIANCE MATRICES. IT IS NOT NECESSARY TO SET UP A 
  ! SEPARATE SHUR PRODUCT FOR ALBEDO + PM TOGETHER, SINCE THEY WILL
  ! NEVER APPEAR IN THE SAME UPDATE AS LONG AS NO DAYTIME PM MEASUREMENTS ARE
  ! USED. 
   
  !FIRST, DEFINE COMPACTLY SUPPORTED COVARIANCE MATRIX, C0, AS GIVEN IN
  ! GASPARI AND COHN, QUART. JOURN. R. MET. SOCIETY, 1999, 125; EQ 4.10
  ALLOCATE(C0(N_YP,N_YP))
  C=RANGE_A/2.
  DO P=1,N_YP !STATE PIXEL LOOP
    DO J=1,N_YP !MEASUREMENT PIXEL LOOP
      ZD=DIST(P,J)
      IF( ZD.LE.C )THEN
        C0(P,J)=-0.25*(ZD/C)**5+0.5*(ZD/C)**4+0.625*(ZD/C)**3-&
          5./3.*(ZD/C)**2+1.
      ELSEIF( ZD.LE.2.*C ) THEN
        C0(P,J)=1./12.*(ZD/C)**5-0.5*(ZD/C)**4+0.625*(ZD/C)**3+&
          5./3.*(ZD/C)**2-5.*(ZD/C)+4.-2./3.*(C/ZD) 
      ELSE
        C0(P,J)=0.
      END IF
      !SOMETIMES, C0 ELEMENTS THAT SHOULD BE ZERO TURN OUT TO BE VERY SMALL
      !  NEGATIVE NUMBERS.  CHECK FOR THESE AND ZERO THEM.  IF THE NUMBER IS
      !  LARGE, REPORT AN ERROR.
      IF(C0(P,J).LT.0.)THEN
        IF(ABS(C0(P,J)).GT.1.E-6)THEN
          IF(RANK.EQ.0) PRINT *, 'ERROR IN KUPDATE: C0(P,J)<0!'
          C0(P,J)=0.
        ELSE
          C0(P,J)=0.
        END IF
      END IF
    END DO 
  END DO 

  !NEXT, MULTIPLY APPROPRIATE PARTS OF COVARIANCE MATRICES BY C0
  !NOTE: IF MEAS_SWITCH=1 (PM ONLY), WE TAKE NO ACTION, BECAUSE ALL STATES ARE
  !   CORRELATED WITH ALL MEAS.
  IF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.2)THEN
    !BECAUSE PM IS ONLY NIGHTTIME AND NIR IS ONLY DAYTIME, N_Z2(1) IS EITHER 12
    !  OR 625 IN THIS CASE
    IF(N_Z2(1).EQ.625)THEN
      !FOR PM + ALBEDO UPDATES, THIS MUST ONLY HAPPEN DURING THE ALBEDO UPDATE
      !  WHEN THE MEASUREMENT VECTOR HAS A LENGTH OF 625
      DO P=1,625 !ROW PIXEL LOOP
        DO J=1,625 !COLUMN PIXEL LOOP
          !MULTIPLY THE COVARIANCE MATRICES ELEMENTWISE BY C0. HERE, I 
          ! MULTIPLIED ALL OF THE STATE COVARIANCES AT EACH PIXEL BY THE SAME
          ! C0 ELEMENT.  I CAN'T SEE WHY THAT WOULDN'T WORK.
          YI=N_YR*(P-1)+1
          YF=N_YR*P
          DO K=YI,YF
            CYZ(K,J)=CYZ(K,J)*C0(P,J)
          END DO
          CZZ(P,J)=CZZ(P,J)*C0(P,J)
        END DO
      END DO
    END IF
  ELSEIF(MEAS_SWITCH.EQ.0)THEN
    !DURING THE DAYTIME, THERE ARE NIR+TIR, SO N_Z2(1)=1250. DURING THE NIGHTTIME,
    !  THERE ARE PM+TIR, N_Z2(1)=637
    IF (N_Z2(1).EQ.1250) THEN
      !DAYTIME: WE HAVE NIR+TIR
      !MY STRATEGY HERE WAS TO COPY THE CODE FOR MEAS_SWITCH 2 OR 3 TWICE FOR DIFFERENT
      !  MEASUREMENT NUMBERS.  IT COULD BE DONE SLIGHTLY MORE COMPACTLY. IT IS NECESSARY
      !  TO APPLY THE SCHUR PRODUCT SEPARATELY FOR THE CORRELATIONS BETWEEN THE DIFFERENT
      !  MEASUREMENT TYPES: FOR CZZ, IT MUST ESSENTIALLY BE APPLIED FOUR TIMES DURING
      !  THE DAYTIME MEASUREMNT.

      !NIR MEASUREMENTS
      DO P=1,625 !ROW PIXEL LOOP
        DO J=1,625 !COLUMN PIXEL LOOP
          !MULTIPLY THE COVARIANCE MATRICES ELEMENTWISE BY C0
          YI=N_YR*(P-1)+1
          YF=N_YR*P
          !CALCULATE THE INDECES TO ACCESS THE COVARIANCES
          I=J
          O=P
          DO K=YI,YF
            CYZ(K,I)=CYZ(K,I)*C0(P,J)
          END DO
          CZZ(O,I)=CZZ(O,I)*C0(P,J)
        END DO
      END DO
      !TIR MEASUREMENTS
      DO P=1,625 !ROW PIXEL LOOP
        DO J=1,625 !COLUMN PIXEL LOOP
          !MULTIPLY THE COVARIANCE MATRICES ELEMENTWISE BY C0
          YI=N_YR*(P-1)+1
          YF=N_YR*P
          !CALCULATE THE INDEX TO ACCESS THE COVARIANCES
          I=J+625
          O=P+625
          DO K=YI,YF
            CYZ(K,I)=CYZ(K,I)*C0(P,J)
          END DO
          CZZ(O,I)=CZZ(O,I)*C0(P,J)
        END DO
      END DO
      !NIR-TIR CORRELATIONS
      DO P=1,625 !ROW PIXEL LOOP
        DO J=1,625 !COLUMN PIXEL LOOP
          !CALCULATE THE INDEX TO ACCESS THE COVARIANCES
          I=J+625 !COLUMN
          O=P !ROW
          CZZ(O,I)=CZZ(O,I)*C0(P,J)
        END DO
      END DO
      !TIR-NIR CORRELATIONS
      DO P=1,625 !ROW PIXEL LOOP
        DO J=1,625 !COLUMN PIXEL LOOP
          !CALCULATE THE INDEX TO ACCESS THE COVARIANCES
          I=J !COLUMN
          O=P+625 !ROW
          CZZ(O,I)=CZZ(O,I)*C0(P,J)
        END DO
      END DO
    ELSEIF (N_Z2(1).EQ.637) THEN
      !NIGHTTIME: WE HAVE PM+TIR
      !TIR MEASUREMENTS
      DO P=1,625 !STATE PIXEL LOOP
        DO J=1,625 !MEASUREMENT PIXEL LOOP
          !MULTIPLY THE COVARIANCE MATRICES ELEMENTWISE BY C0
          YI=N_YR*(P-1)+1
          YF=N_YR*P
          !CALCULATE THE INDEX TO ACCESS THE COVARIANCES
          I=J+12
          O=P+12
          DO K=YI,YF
            CYZ(K,I)=CYZ(K,I)*C0(P,J)
          END DO
          CZZ(O,I)=CZZ(O,I)*C0(P,J)
        END DO
      END DO
    ELSEIF (N_Z2(1).EQ.14) THEN
      !NO CHNAGE NEEDS TO BE TAKEN
    ELSE
      !I WROTE THIS ASSUMING THAT N_Z2(1) WOULD BE 637 OR 1250... STOP
      !  HERE IF IT IS SOME OTHER NUMBER
      PRINT *, 'KUPDATE: INCORRECT CODE CONFIGURATION...'
      CALL MPI_FINALIZE(ERR)
      STOP
    END IF
  END IF

  ! 7) ON 0 NODE, COMPUTE GAIN THEN DISTRIBUTE TO INDIVIDUAL PROCESSES
  IF (RANK.EQ.0) THEN

!    print *, 'kupdate: starting step 7'

    CZZCV=CZZ+CV
!    DO I=1,N_Z2(1)
!      WRITE(16,'(625(E10.4,1X))') (CZZ(I,J),J=1,N_Z2(1))
!    END DO
    CALL INVERT_MATRIX(CZZCV,CZZCVINV,N_Z2(1))
    KALM=MATMUL(CYZ,CZZCVINV)
  END IF 

  CALL MPI_BCAST(KALM,N_Y*N_Z2(1),MPI_REAL,0,MY_COMM,ERR)

  ! 8) COMPUTE RANDOM VARIATES V FOR UPDATE
!    print *, 'kupdate: starting step 8'
  DO I=1,N_Z2(1)
    MU0(I,1)=0.
  END DO
  SZMU0=SHAPE(MU0)
  SZCV=SHAPE(CV)
  ! MVNRND RETURNS THE VARIATES WITH CASES (N_R, HERE) AS THE FIRST 
  ! DIMENSION.  FOR CONSISTENCY, THIS IS TRANSPOSED TO YIELD V
  ! FOR RANDOM NUMBER GENERATOR SEED SCHEME, SEE JOURNAL JUNE 23, 2006
  CALL MVNRND(MU0,CV,V_TRANS,N_R,SZMU0(1),SZMU0(2),SZCV(1),SZCV(2),70+M)
  V=TRANSPOSE(V_TRANS)

  ! ONLY THE FIRST N_Z2(1) VALUES ARE MAPPED TO V_OUT, THOUGH THE DIMENSION 
  !   OF V_OUT IS ACTUALLY N_Z
  V_OUT=0.0
  V_OUT(1:N_Z2(1),1:N_R)=V(1:N_Z2(1),1:N_R)

  ! 9) PERFORM UPDATE
!    print *, 'kupdate: starting step 9'
  DO K=1,N_RL
    ! CALCULATE GLOBAL REPLICATE NUMBER, KGLOBAL TO INDEX V
    KGLOBAL=RANK*N_RL+K
    DO I=1,N_Y
      DO J=1,N_Z2(1)
        Y(I,K)=Y(I,K)+KALM(I,J)*(ZMEAS(J)+V(J,KGLOBAL)-Z(J,K))
      END DO
    END DO
  END DO

  ! 10) CALCULATE FILTER INNOVATION (THIS COULD BE DONE AFTER 7), FOLLOWING
  !       CROW AND VAN LOON, 2006. 
!    print *, 'kupdate: starting step 10'
  ALLOCATE(NU(N_Z2(1)),NUTCI(N_Z2(1)))
  IF(RANK.EQ.0) THEN
    !A) COMPUTE MEAN INNOVATIONS, NU
    NU=ZMEAS-ZMEAN
    WRITE(15,'(625(E10.4,1X))') (NU(I),I=1,N_Z2(1))

    !B) COMPUTE PRODUCT OF NUT * CZZCVINV
    NUTCI=0.
    DO I=1,N_Z2(1)
      DO J=1,N_Z2(1)
        NUTCI(I)=NUTCI(I)+NU(J)*CZZCVINV(J,I)
      END DO
    END DO

    !C) COMPUTE PRODUCT OF NUTCI * NU
    INNOV_2MOM=0.
    DO I=1,N_Z2(1)
      INNOV_2MOM=INNOV_2MOM+NUTCI(I)*NU(I)/REAL(N_Z2(1))
    END DO
  END IF

  ! 11) COMPUTE ANOTHER METRIC FOR INNOVATION RIGHTNTESS: RELATIVE
  !    DIFFERENCE OF THE TRACE OF ACTUAL AND EXPECTED INNOVATION
  !    COVARIANCE
  
  !   A) COMPUTE INNOVATION COVARIANCE ACROSS THE ENSEMBLE BY REDUCING A 
  !      LOCAL SUM OF THE SQUARE OF INNOVATION DEVIANCE FROM ITS EXPECTED
  !      MEAN ACROSS ALL PROCESSES
!  print *, 'starting 11'
  LOCALSSENUNU=0.0
  DO I=1,N_Z2(1)
    DO J=1,N_Z2(1)
      DO K=1,N_RL
        LOCALSSENUNU(I,J)=LOCALSSENUNU(I,J)+(ZMEAS(I)-Z(I,K))*&
          (ZMEAS(J)-Z(J,K))
      END DO
    END DO
  END DO
  CALL MPI_REDUCE(LOCALSSENUNU,GLOBALSSENUNU,N_Z2(1)*N_Z2(1),MPI_REAL,MPI_SUM,&
    0,MY_COMM,ERR)
  CNUNU=GLOBALSSENUNU/(N_R-1)
  
  !   B) COMPUTE TRACE OF CNUNU AND CZZ+CV
  TR_CNUNU=0. !INITIALIZE TRACE OF CNUNU
  TR_CZZ=0. !INITIALIZE TRACE OF CZZ
  DO I=1,N_Z2(1)
    TR_CNUNU=TR_CNUNU+CNUNU(I,I)
    TR_CZZ=TR_CZZ+CZZ(I,I)+CV(I,I)
  END DO

  !   C) COMPUTE RELATIVE DIFFERENCE IN TRACES
  TR_ERR=(TR_CNUNU-TR_CZZ)/TR_CZZ

  ! 12) COMPUTE A SIMILAR METRIC, BUT ABSOLUTE VALUE
  
!    print *, 'kupdate: starting step 12'
  ! COMPUTE ABSOLUTE VALUE OF DIFFERENCES BETWEEN COVARIANCE
  !  MATRICES
  ABS_TR_ERR=0.
  DO I=1,N_Z2(1)
    ABS_TR_ERR=ABS_TR_ERR+ABS(CNUNU(I,I)-CZZ(I,I)-CV(I,I))/(CZZ(I,I)+CV(I,I))
  END DO

  ! 13) COMPUTE RELATIVE L1-NORM
  !    A) COMPUTE L1-NORM OF MATRIX DIFFERENCES
!    print *, 'kupdate: starting step 13'
  L1_NORM_DIFF=0.
  DO J=1,N_Z2(1)
    COL_SUM=0.
    DO I=1,N_Z2(1)
      COL_SUM=COL_SUM+ABS(CNUNU(I,J)-CZZ(I,J)-CV(I,J))
    END DO
    L1_NORM_DIFF=MAX(L1_NORM_DIFF,COL_SUM)
  END DO
  !    B) COMPUTE L1-NORM OF EXPECTED VALUE
  L1_NORM_EXP=0.
  DO J=1,N_Z2(1)
    COL_SUM=0.
    DO I=1,N_Z2(1)
      COL_SUM=COL_SUM+ABS(CZZ(I,J)+CV(I,J))
    END DO
    L1_NORM_EXP=MAX(L1_NORM_EXP,COL_SUM)
  END DO
  !    C) COMPUTE RELATIVE L1-NORM
  REL_L1_NORM=L1_NORM_DIFF/L1_NORM_EXP

  ! 14) COMPUTE AVERAGE OF INNOVATION MEANS
!    print *, 'kupdate: starting step 14'
  NU_BAR=0.
  DO I=1,N_Z2(1)
    NU_BAR=NU_BAR+NU(I)/REAL(N_Z2(1))
  END DO

!    print *, 'kupdate: after step 14'

  DEALLOCATE(Z,ZMEAN,ZMEAS,CYZ,CZZ,CV,KALM,CZZCVINV,V,MU0,LOCALSUMZ,&
    GLOBALSUMZ,LOCALSSEZZ,ZTILDA,GLOBALSSEZZ,LOCALSSEYZ,GLOBALSSEYZ,&
    zuse_local,V_TRANS,CZZCV,LOCALSSENUNU,GLOBALSSENUNU,CNUNU)
  DEALLOCATE(YMEAN,LOCALSUMY,GLOBALSUMY,YTILDA)
  DEALLOCATE(C0,NU,NUTCI)

!    print *, 'kupdate: after deallocation'

END IF

END SUBROUTINE KUPDATE