enkf.f90

!
 PROGRAM ENKF
!
! ----------------------------------------------------------------------

!  THIS IS AN IMPLEMENTATION OF THE ENSEMBLE KALMAN FILTER FOR MPI 
!  PARALLEL COMPUTING. MULTI-SCALE MEASUREMENTS ARE ACCOMMODATED BY 
!  STORING ALL THE MEASUREMENTS IN A SINGLE ARRAY UNTIL THE 
!  UPDATE SUBROUTINE. IN THIS IMPLEMENTATION, THE SAME NUMBER OF 
!  ENSEMBLE MEMBERS ARE STORED ON EACH OF THE PROCESSORS, AND THE UPDATE
!  IS PERFORMED USING 'REDUCE' COMMANDS IN THE KUPDATE SUBROUTINE.
!  THUS, THE NUMBER OF REPLICATES N_R MUST BE DIVISIBLE BY THE NUMBER 
!  OF PROCESSES NPROCS. FOR MEMORY CONSIDERATIONS, ONLY THE KEY ENSEMBLE
!  STATISTICS ARE STORED AT EACH TIMESTEP, AND ONLY ONE COPY OF THE 
!  ENSEMBLE IS CREATED IN MEMORY.  THE ENSEMBLE STATISTICS ARE WRITTEN 
!  EVERY FOUR HOURS TO A BINARY FILE.  COARSE SCALE  NOMINAL FORCING 
!  DATA ARE READ IN, AND RANDOMLY DISAGGREGATED WITHIN THE 
!  THE MEASUREMENT LOOP TO AVOID STORING THE ENTIRE RANDOMIZED FORCING 
!  ARRAY IN MEMORY, THOUGH THE PARAMETERS ARE RANDOMIZED AND STORED ON 
!  THE LOCAL PROCESS. IN ORDER TO AID IN DEBUGGING, ALL OF THE ENSEMBLE 
!  INFORMATION IS WRITTEN TO FILES ON EACH PROCESSOR AT THE MEASUREMENT  
!  TIME SPECIFIED IN THIS SOURCE CODE. MEASUREMENT INNOVATIONS ARE 
!  WRITTEN OUT TO A FILE.
!
!  NOTE: THIS CODE WAS LAST UPDATED TO THE RUN_PARAMS SUBROUTINE IN JULY 
!    2007, AND IS LIKELY OUT OF DATE, AT THAT POINT.
!
!  KEY VARIABLES
!  N_Y ~ NUMBER OF PIXELS TIMES THE NUMBER OF STATES PER PIXEL
!  N_R ~ TOTAL SIZE OF THE ENSEMBLE
!  N_RL ~ NUMBER OF ENSEMBLE MEMBERS STORED ON EACH LOCAL PROCESSOR
!  N_X ~ NUMBER OF AUXILIARY STATES AT EACH PIXEL
!  N_Z ~ TOTAL NUMBER OF MULTI-SCALE MEASUREMENTS
!  N_T ~ TOTAL NUMBER OF FORCING TIMES
!  N_YS ~ NUMBER OF STATISTICS CALCULATED FROM THE ENSEMBLE (MAX,MIN...)
!
!  Y ~ ARRAY OF STATES: N_Y x N_RL
!  X ~ ARRAY OF AUXILIARY STATES: N_X x N_RL. X IS NOT UPDATED, BUT 
!         NEEDED TO ALLOW STATE MODEL TO RESUME EXECUTION CORRECTLY
!  Z ~ ARRAY OF MEASUREMENTS N_Z x N_RL
!  YSTATS ~ ARCHIVE OF STATE VARIABLES OVER THE ENTIRE EXECUTION TIME,
!      INCLUDING BOTH PRIOR AND POSTERIOR STATES: N_Y x N_T+N_M+1 x N_YS
!
! CALLS: INTERFACEY, INTERFACEZ, KUPDATE, COMPILEYSTATS
!           SET UP FOR MPI - MD - AUG 2005
!           INITIALLY IMPLEMENTED FOR SNOW MODELS - MD - DEC 2005
!           ADJUSTED TO DISAGGREGATE FORCING - MD - JAN 2006
!           SET UP TO INTERFACE WITH NEW RUN_PARAMS - MD - JUL 2007
!
! OUTLINE:
!
! 0) DECLARATIONS
! 1) INITIAL SETUP (BEFORE MEASUREMENT LOOP BEGINS)
! 2) MEASUREMENT LOOP TO PERFORM FILTER OPERATIONS
! 3) COMPUTE ELAPSED TIME
! 4) ON HEAD NODE, WRITE YSTATS AND YTIME FROM FILTER RUN TO FILES

IMPLICIT NONE
INCLUDE 'mpif.h'
!INCLUDE 'mpef.h'



! 0.1) DECLARATIONS FOR PART 1
! 0.1.1) INITIALIZE MPI/MPE
INTEGER::RANK,NPROCS,STAT(MPI_STATUS_SIZE),IERR,N_RL
DOUBLE PRECISION STARTTIME,ENDTIME
! 0.1.2) READ SIZES OF ARRAYS FROM FILE
INTEGER::N_Y,N_YR,N_T,N_M,N_YP,N_YS,N_Z,N_ZR,N_R,N_U,N_A,N_X,N_C,I,N_PF,N_RP
INTEGER,DIMENSION(:),ALLOCATABLE::RTM_CTRL,N_ZP
CHARACTER(32) :: N_ZRC
! 0.1.4) GET RUN PARAMETERS
INTEGER :: MEAS_SWITCH,SPATIAL_PIC,ATM_SWITCH,GEN_SEEDS(5),ENS_NUM
INTEGER,DIMENSION(:),ALLOCATABLE :: MEASTIMES,ISTART,IEND,ZUSE,MSTART,MEND,ZUSE1
REAL :: PRECIP_SCALE,CVEG,LAM_VEG,DUMMY_INV_PRCP,RANGE_U,&
  RANGE_A,DUMMY_B_1,DUMMY_B_2!,F_RADT
 ! LAM_RADT,C_RADT,PSCALE,RANGE_B,RANGE_RADT,RANGE_CHI_P,C_CHI_P,&
 ! LAM_CHI_P 
REAL,DIMENSION(:),ALLOCATABLE:: ALPHA_BAR, FREQ, THETA,TB_OUTPUT,TB_OUTPUTALL,CV1
REAL,DIMENSION(:,:),ALLOCATABLE::CV,Y_IN,X_IN,TBRAW,LAM_U,LAM_ALPHA!
REAL,DIMENSION(:,:,:),ALLOCATABLE::ALBEDO,CU,C_ALPHA
LOGICAL :: GEN_SWITCH
! 0.1.5) READ FROM FILE
INTEGER :: M,J,N_BDRF,rp
INTEGER,DIMENSION(:),ALLOCATABLE::VEG_MAP
REAL,DIMENSION(:),ALLOCATABLE::ELEVATION,VCOVER
REAL,DIMENSION(:,:),ALLOCATABLE::ZMEAS,ZS,ZF,BDRF
CHARACTER(32) :: N_TC,N_ZC
! 0.1.7) INITIAL CONDITIONS
INTEGER :: R,XI,XF,SEED0
REAL,DIMENSION(:,:), ALLOCATABLE :: Y0_IN,X0_IN
REAL :: T0(5),TF(5)
! 0.1.8) CREATE ENSEMBLE OF PARAMETERS
REAL,DIMENSION(:,:,:), ALLOCATABLE :: ALPHA,FALPHA
INTEGER :: KGLOBAL
! 0.1.9) CREATE STATIC ENSEMBLE OF FORCING PERTURBATIONS
INTEGER :: SEED
INTEGER,DIMENSION(:),ALLOCATABLE :: MAP 
REAL,DIMENSION(:,:,:),ALLOCATABLE :: F_U
! 0.1.10) CREATE STATIC ENSEMBLE OF VEGETATION PERTURBATIONS
REAL,DIMENSION(:,:),ALLOCATABLE :: F_VEG ! MATCHES DIMENSION IN MVNRND
! 0.2 MEASUREMENT LOOP
! 0.2.1) INDICES
INTEGER YINDEX1,YINDEX2,UINDEX1,UINDEX2,N_UT,N_S,FIRST
! 0.2.2) DEFINE INITIAL CONDITIONS
REAL,DIMENSION(:,:), ALLOCATABLE :: Y0,X0
! 0.2.4) STATE MODEL INTEGRATION
INTEGER :: YI,YF,P,K,MONTH,IFLAG
LOGICAL :: LIMIT
REAL,DIMENSION(:), ALLOCATABLE :: Y0_RAW,ZMEAS1
REAL,DIMENSION(:,:), ALLOCATABLE :: Y,Z,V_OUT,ALBEDO_IN,YPRIOR,F,X,FT_1,Z1,V_OUT1
REAL,DIMENSION(:,:,:),ALLOCATABLE::YT,U_NOM,XT,U_RTM,U_DISAG,U
! 0.2.6) DEFINE TIME VECTOR
DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: YTIME
! 0.2.8) UPDATE 
REAL,DIMENSION(:),ALLOCATABLE :: INNOV_2MOM,TR_ERR,REL_L1_NORM,ABS_TR_ERR,&
  NU_BAR
! 0.2.9) RECORD STATE STATISTICS
REAL,DIMENSION(:,:,:),ALLOCATABLE::YSTATS, DUMMY
! 0.4) WIRTE OUTPUT FILES
INTEGER :: T
LOGICAL :: OUTPUT
DOUBLE PRECISION :: TOUTPUT(6), DIFF
!0.)  DECLARATIONS FOR FILE NAMES...
CHARACTER(9):: FNAMES
CHARACTER(22):: FNAMES2

!INTEGER :: UPDATE,MINUPDATE

! 1) INITIAL SETUP (BEFORE MEASUREMENT LOOP BEGINS)
! 1.1) INITIALIZE MPI AND MPE VARIABLES
CALL MPI_INIT(IERR)
CALL MPI_COMM_RANK(MPI_COMM_WORLD,RANK,IERR)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NPROCS,IERR)

!1.2) CODE USED TO OPEN FILES FOR EACH PROCESS
!FNAMES= "123456789"
!FNAMES2="1011121314151617181920"
!I=RANK+1
!IF(I.LT.10)THEN
!  OPEN(FILE=FNAMES(I:I),unit=100-I,status='unknown')!ELSE
!  J=I-9
!  OPEN(FILE=FNAMES2(2*J-1:2*J),unit=100-I,status='UNKNOWN')
!END IF

!OPEN FILE FOR INNOVATIONS
OPEN(UNIT=15,FILE='innovs.out',STATUS='UNKNOWN')

! INITIALIZE MPE
!CALL MPE_INIT_LOG()

! START TIMER
IF(RANK.EQ.0) STARTTIME=MPI_WTIME()

! ENTER MPE PROPERTIES
IF(RANK.EQ.0)THEN
!  CALL MPE_DESCRIBE_STATE(1,2,'PROPAGATION','white:gray')
!  CALL MPE_DESCRIBE_STATE(11,12,'CALCULATION OF MEASUREMENTS','green:gray')
!  CALL MPE_DESCRIBE_STATE(13,14,'UPDATE','cyan:gray')
!  CALL MPE_DESCRIBE_STATE(15,16,'COMPILE STATS','yellow:gray')
END IF

! 1.2) READ SIZES OF THINGS FROM FILE, DEFINE SECONDARY SIZES
OPEN(UNIT=1,FILE='sizes.in',STATUS='OLD')

!NUMBER OF STATES AT EACH PIXEL, N_YR
READ(1,*) 
READ(1,10) N_YR
!NUMBER OF TIMESTEPS, N_T
READ(1,*) 
READ(1,10) N_T
!NUMBER OF MEASUREMENT TIMES, N_M
READ(1,*) 
READ(1,10) N_M
!TOTAL NUMBER OF Y PIXELS, N_YP
READ(1,*)
READ(1,10) N_YP
!NUMBER OF Y STATISTICS TO BE RECORDED, N_YS
READ(1,*) 
READ(1,10) N_YS
!NUMBER OF FORCING DATA REQ'D FOR EACH PIXEL AND TIMESTEP, N_U
READ(1,*) 
READ(1,10) N_U
!NUMBER OF REPLICATES: DIVISIBLE BY NPROCS IN THIS VERSION, N_R
READ(1,*) 
READ(1,10) N_R
!NUMBER OF MEASUREMENT TYPES AT EACH PIXEL, N_ZR
READ(1,*) 
READ(1,10) N_ZR
!TOTAL NUMBER OF MEAS. FOR EACH RESOLUTION (MUST BE N_ZR VALUES HERE),N_ZP
ALLOCATE(N_ZP(N_ZR))
WRITE(N_ZRC,*) N_ZR
READ(1,*) 
READ(1,'('//N_ZRC//'(I4,1X))') (N_ZP(I),I=1,N_ZR) 
!NUMBER OF PARAMETERS, N_A
READ(1,*) 
READ(1,10) N_A
!NUMBER OF AUXILIARY VARIABLES FOR STATE MODEL, N_X
READ(1,*)
READ(1,10) N_X
!NUMBER OF RTM CONTROL VARIABLES
READ(1,*)
READ(1,10) N_C
ALLOCATE(RTM_CTRL(N_C))
RTM_CTRL=0 !INITIALIZE ARRAY
!NUMBER OF AUXILIARY INPUTS FOR SNOW RTM
READ(1,*)
READ(1,10) RTM_CTRL(5)
!NUMBER OF SNOW INPUTS FOR SNOW RTM 
READ(1,*)
READ(1,10) RTM_CTRL(6)
!NUMBER OF INPUTS FOR CANOPY RTM 
READ(1,*)
READ(1,10) RTM_CTRL(7)
!NUMBER OF INPUTS FOR ATMOSPHERIC RTM 
READ(1,*)
READ(1,10) RTM_CTRL(8)
!NUMBER OF  PASSIVE MICROWAVE CHANNELS
READ(1,*)
READ(1,10) RTM_CTRL(9)
!NUMBER OF NLDAS FORCING PIXELS
READ(1,*)
READ(1,10) N_PF
10 FORMAT(I5)
CLOSE(1)
N_Y=N_YR*N_yp  !COMPUTE TOTAL NUMBER OF STATES
N_Z=SUM(N_ZP)  !COMPUTE TOTAL NUMBER OF MEASUREMENTS


! 1.3) CHECK TO BE SURE THAT N_R IS DIVISIBLE BY NPROCS, COMPUTE DIMENSION OF
! LOCAL STATE VARIABLE, THEN ALLOCATE Y,Y0,Z AND V_OUT ARRAYS
IF (MOD(N_R,NPROCS).NE.0) THEN
  PRINT *, 'ERROR IN ENKF!  NUMBER OF REPLICATES MUST BE DIVISIBLE BY ',&
           'NUMBER OF PROCESSES.  ABORTING'
  STOP
ELSEIF(N_R.EQ.1) THEN
  PRINT *, 'FOR THIS CODE, 1 REPLICATE IS NOT A VALID OPTION.  SPECIFY',&
    'TWO OR MORE REPLICATES'
ELSE
  N_RL=N_R/NPROCS   !NUMBER OF REPLICATES ON LOCAL PROCESS
END IF

! 1.4) GET RUN PARAMETERS 

!ALLOCATIONS
ALLOCATE(MEASTIMES(N_M),CV(N_Z,N_Z),ZUSE(N_Z),CU(N_YP,N_YP,N_U),&
  LAM_U(N_YP,N_U),ISTART(NPROCS),IEND(NPROCS),ALPHA_BAR(N_A),&
  FREQ(RTM_CTRL(9)),THETA(RTM_CTRL(9)),MSTART(NPROCS),MEND(NPROCS),&
  C_ALPHA(N_YP,N_YP,N_A),LAM_ALPHA(N_YP,N_A))

GEN_SWITCH=.FALSE.

! THE DUMMY_INV_PRCP ARGUMENT IS TO HANDLE THE RUN_PARAMS INVERT_PRECIP 
! ARGUMENT THAT IS USED TO IVNERT THE LAPSE RATE IN THE GENERATOR.  IT 
! IS PASSED TO THE DISAGGREGATE SUBROUTINE, BUT IT IS NOT USED

CALL RUN_PARAMS(N_YR,N_T,N_M,N_YP,N_Y,N_YS,N_U,N_R,N_ZR,N_ZP,N_Z,&
  MEASTIMES,CU,CV,LAM_U,ISTART,IEND,NPROCS,N_A,ALPHA_BAR,ZUSE,RTM_CTRL(9),&
  FREQ,THETA,MSTART,MEND,C_ALPHA,LAM_ALPHA,PRECIP_SCALE,MEAS_SWITCH,&
  GEN_SWITCH,SPATIAL_PIC,CVEG,LAM_VEG,DUMMY_INV_PRCP,RANK,&
  ATM_SWITCH,RANGE_U,RANGE_A,GEN_SEEDS,ENS_NUM,DUMMY_B_1,DUMMY_B_2)

RTM_CTRL(3)=ATM_SWITCH

! 1.5) READ MEASUREMENTS, LATITUDE/LONGITUDE DATA, VEGETATION AND FORCING DATA
! 1.5.1) READ MEASUREMENTS
ALLOCATE(ZMEAS(N_Z,N_M))
IF(RANK.EQ.0)THEN
  WRITE(N_ZC,*) N_Z
  OPEN(UNIT=2,FILE='zmeas.in',STATUS='OLD')
  DO M=1,N_M
!    DO I=1,N_Z
      READ(2,'('//N_ZC//'F12.4)') (ZMEAS(I,M),I=1,N_Z)
!    END DO
  END DO
  CLOSE(2)
!  25 FORMAT(F12.4)
END IF
CALL MPI_BCAST(ZMEAS,N_M*N_Z,MPI_REAL,0,MPI_COMM_WORLD,IERR)

! 1.5.2) READ STATE PIXEL LAT LON PAIRS
ALLOCATE(ZS(N_YP,2))
OPEN(UNIT=3,FILE='state_pixels.in',STATUS='OLD')
DO I=1,N_YP
  READ(3,30) ZS(I,1),ZS(I,2)
END DO
CLOSE(3)

! 1.5.3) READ NLDAS (COARSE FORCING) PIXEL LAT LON PAIRS
ALLOCATE(ZF(N_PF,2))
OPEN(UNIT=4,FILE='nldas_pixels.in',STATUS='OLD')
DO I=1,N_PF
    READ(4,30) ZF(I,1),ZF(I,2)
END DO
30 FORMAT(2F12.4)
CLOSE(4)

! 1.5.4) READ ELEVATIONS AT STATE PIXELS
!        THIS FILE IS ARRANGED BY PRINTING THE 25 VALUES OF THE LOWERMOST ROW 
!        ON THE FIRST LINE, ETC.
ALLOCATE(ELEVATION(N_YP))
OPEN(UNIT=5,FILE='state_elev.in',STATUS='OLD')
DO I=1,N_YP
  READ(5,40) ELEVATION(I)
END DO
CLOSE(5)
40 FORMAT(F12.5) 

! 1.5.5) READ BDRF LOOKUP TABLE FOR INTERFACEZ
OPEN(UNIT=60,FILE='bdrf_lookup.in',STATUS='OLD')
READ(60,*) N_BDRF
ALLOCATE(BDRF(N_BDRF,2)) 
DO I=1,N_BDRF
  READ(60,47) BDRF(I,1),BDRF(I,2)
END DO
CLOSE(60)
45 FORMAT(I10)
47 FORMAT(2(F12.6))

! 1.5.6) READ VEGETATION DATA FROM FILE FOR SSIB
ALLOCATE(VEG_MAP(N_YP))
CALL READ_VEG(N_YP,VEG_MAP)

! 1.5.7) READ VEGETATION COVER DATA
ALLOCATE(VCOVER(N_YP))
OPEN(UNIT=7,FILE='veg_density.in',STATUS='OLD')
DO I=1,N_YP
  READ(7,50) VCOVER(I)
END DO
50 FORMAT(F10.4) 
CLOSE(7)

! 1.5.8) READ FORCING DATA FROM FILE ON 0 NODE, READ NLDAS DATA, BROADCAST
ALLOCATE(U_NOM(N_U,N_T,N_PF))
IF(RANK.EQ.0)THEN
  WRITE(N_TC,*) N_T
  OPEN(FILE='nldas_forcing.in',UNIT=8,STATUS='OLD')
  DO I=1,N_PF
    DO J=1,N_U
      READ(8,'('//N_TC//'F12.4)') (U_NOM(J,K,I),K=1,N_T)
    END DO
    READ(8,*) !PICK UP THE BLANK LINE HERE...
  END DO
  CLOSE(8)

  !1.6) SCALE AND ADJUST PRECIP FOR SYNTHETIC TEST / BROADCAST
  DO I=1,N_T
    U_NOM(3,I,1:N_PF)=U_NOM(3,I,1:N_PF)*PRECIP_SCALE
  END DO
END IF

!if (rank.eq.0) then
! print *, (u_nom(3,k,1) ,k=240:312)
!end if
!stop

!SEND TO OTHER PROCESSES FROM 0
CALL MPI_BCAST(U_NOM,N_T*N_U*N_PF,MPI_REAL,0,MPI_COMM_WORLD,IERR)

!IF SPATIAL_PIC PARAMETER IS SET TO 0, ALL NLDAS FORCING WILL BE USED
!OTHERWISE, NLDAS FORCING WILL BE SET TO THE FORCING FROM SPATIAL_PIC PIXEL
IF(SPATIAL_PIC.GT.0)THEN
  DO I=1,N_PF
    U_NOM(1:N_U,1:N_T,I)= U_NOM(1:N_U,1:N_T,SPATIAL_PIC) 
  END DO
END IF

! 1.7) ASSIGN INITIAL CONDITION OF STATES, AUXILIARY
ALLOCATE(Y0_IN(N_Y,N_RL),X0_IN(N_X*N_YP,N_RL),Y0_RAW(N_YR))
Y0_IN=0. !INITIALIZE ALL SNOW VARIABLES AT 0.
DO I=14,N_Y,N_YR
  DO R=1,N_RL
    Y0_IN(I,R)=265.
  END DO
END DO
!INITIALIZE X0 AS IN THE GENERATOR
X0_IN=0.
DO K=1,N_RL
  DO P=1,N_YP
    XI=1+N_X*(P-1)
    DO I=XI,XI+3
      X0_IN(I,K)=270.
    END DO
    DO I=XI+4,XI+6
      X0_IN(I,K)=0.15
    END DO
    X0_IN(XI+9,K)=0.5E-3
    X0_IN(XI+10,K)=4.3
  END DO
END DO

! TIME DATA BASED ON MIDNIGHT UTC DATA...
T0(1)=18
T0(2)=9
T0(3)=30
T0(4)=273
T0(5)=2002
!T0(1)=18   !HOUR
!T0(2)=10   !MONTH
!T0(3)=31   !DAY OF THE MONTH
!T0(4)=304  !JULIAN DAY

! 1.8) CREATE ENSEMBLE OF PARAMETERS
ALLOCATE(ALPHA(N_RL,N_YP,N_A),FALPHA(N_R,N_YP,N_A))
! FOR RANDOM NUMBER GENERATOR SEED SCHEME, SEE Reports/Models/EnKF/seeds.odt
DO I=1,N_A
  IF(C_ALPHA(1,1,I).EQ.0.)THEN
    ! MVNRND CANNOT BE CALLED WHEN THE COVARIANCE HAS ZERO DIAGONAL VALUES,
    ! BECAUSE IT IS SET UP TO CALL CHOLESKY, WHICH CANNOT WORK FOR ZERO
    ! DIAGONAL VALUES IN ITS' CURRENT SETUP.  SET F_ALPHA TO 0.
    FALPHA(1:N_R,1:N_YP,I)=0.
  ELSE
    CALL MVNRND(LAM_ALPHA(1:N_YP,I),C_ALPHA(1:N_YP,1:N_YP,I),&
      FALPHA(1:N_R,1:N_YP,I),N_R,N_YP,1,N_YP,N_YP,40+I)
  END IF
END DO

! 1.8.3) FORM PARAMETER ARRAY
DO I=1,N_A
  DO J=1,N_YP
    DO K=1,N_RL
      KGLOBAL=RANK*N_RL+K
      ALPHA(K,J,I)=ALPHA_BAR(I)*EXP(FALPHA(KGLOBAL,J,I))
    END DO
  END DO
END DO
!LIMIT ALPHA_5
DO I=5,5
  DO J=1,N_YP
    DO K=1,N_RL
      ALPHA(K,J,I)=MIN(ALPHA(K,J,I),1.0E-3)
    END DO
  END DO
END DO

!1.9) CREATE STATIC ENSEMBLE OF FORCING PERTURBATIONS
ALLOCATE(F_U(N_R,N_YP,N_U),MAP(N_YP))
! FOR RANDOM NUMBER GENERATOR SEED SCHEME, SEE Reports/Models/EnKF/seeds.odt
DO I=1,N_U
  IF(CU(1,1,I).EQ.0)THEN
    ! MVNRND CANNOT BE CALLED WHEN THE COVARIANCE HAS ZERO DIAGONAL VALUES,
    ! BECAUSE IT IS SET UP TO CALL CHOLESKY, WHICH CANNOT WORK FOR ZERO
    ! DIAGONAL VALUES IN ITS' CURRENT SETUP.  SET F_U TO 0.
    F_U(1:N_R,1:N_YP,I)=0.
  ELSE
    !FOR MODIFICATION OF RANDOM NUMBER SCHEME, SEE JOURNAL ON OCT 09, 2006
    SEED=50+I+1000*ENS_NUM
    CALL MVNRND(LAM_U(1:N_YP,I),CU(1:N_YP,1:N_YP,I),F_U(1:N_R,1:N_YP,I),&
      N_R,N_YP,1,N_YP,N_YP,SEED)
  END IF
END DO

!1.10) CREATE STATIC ENSEMBLE OF VEGETATION PERTURBATIONS
! FOR RANDOM NUMBER GENERATOR SEED SCHEME, Reports/Models/EnKF/seeds.odt 
ALLOCATE(F_VEG(N_R,1))
CALL MVNRND(LAM_VEG,CVEG,F_VEG,N_R,1,1,1,1,60)

!2.) MEASUREMENT LOOP TO PERFORM FILTER OPERATIONS

! 2.0) INITIALIZATIONS AND ALLOCATIONS FOR MEASUREMENT LOOP
ALLOCATE(Y(N_Y,N_RL),Z(N_Z,N_RL),V_OUT(N_Z,N_R),X(N_X*N_YP,N_RL),&
  YPRIOR(N_Y,N_RL),Y0(N_Y,N_RL),X0(N_X*N_YP,N_RL),&
  YTIME(N_T+N_M+1), YSTATS(N_Y,N_T+N_M+1,N_YS),ALBEDO_IN(N_YP,N_RL),&
  INNOV_2MOM(N_M),TR_ERR(N_M),REL_L1_NORM(N_M),ABS_TR_ERR(N_M),&
  NU_BAR(N_M))
ALLOCATE(Z1(1,N_RL)) ! RK 11/25/2014

! INITIALIZE STATE ARRAY AND 'FIRST' INDEX
Y=Y0_IN
X=X0_IN
FIRST=0


!OPEN FILE FOR SPATIALLY-AVERAGED TB OUTPUT DIAGNOSTIC
IF(RANK.EQ.0) THEN
  OPEN(file='tbraw_t.out',unit=9,status='unknown')
END IF

!ALLOCATE VARIABLES FOR SPATIALLY-AVERAGED TB OUTPUT DIAGNOSTIC
ALLOCATE(TB_OUTPUT(N_ZP(1)),TB_OUTPUTALL(N_ZP(1)))

!MEASUREMENT LOOP
DO M=1,N_M+1
  IF(RANK.EQ.0) PRINT *,'BEGINNING MEASUREMENT INTERVAL ',M,'/',N_M+1

  !2.1) COMPUTE INDECES
  !2.1.1)DETERMINE NUMBER OF FORCING STEPS FOR THIS INTERVAL (ALL NODES)
  IF (M.EQ.N_M+1) THEN
    N_S=N_T-MEASTIMES(M-1)
  ELSE
    N_S=MEASTIMES(M)-FIRST
  END IF

  !DETERMINE INDECES FOR STATE ARRAY Y AND FORCING ARRAY U (ALL NODES)
  IF (M.EQ.1) THEN
    YINDEX1=1
  ELSEIF (M.EQ.2) THEN
    YINDEX1=YINDEX1+(MEASTIMES(M-1))+1
  ELSE
    YINDEX1=YINDEX1+(MEASTIMES(M-1)-MEASTIMES(M-2))+1
  END IF
  IF (M.EQ.N_M+1) THEN
    YINDEX2=YINDEX1+N_S
  ELSE
    YINDEX2=YINDEX1+N_S+1
  END IF

  IF (M.EQ.N_M+1) THEN
    UINDEX1=MEASTIMES(M-1)+1
    UINDEX2=N_T
  ELSE
    UINDEX1=FIRST+1
    UINDEX2=MEASTIMES(M)
  END IF
  N_UT=UINDEX2-UINDEX1+1 !NUMBER OF FORCING STEPS

  !2.2) DEFINE INITIAL CONDITIONS
  Y0=Y
  X0=X
  IF(M>1) T0=TF !USE PREVIOUS UPDATE TIME TO INITIA

! CALL MPE_LOG_EVENT(9,M,'BEGIN PROP.')

  !2.3) ALLOCATIONS FOR MODEL INTEGRATION
  IF(M<N_M+1)THEN
    ALLOCATE(YT(N_Y,N_S+2,N_RL)) ! N_S PLUS I.C. AND POSTERIOR
  ELSE
    ALLOCATE(YT(N_Y,N_S+1,N_RL)) ! N_S PLUS I.C.; NO UPDATE
  END IF
  ALLOCATE(XT(N_X*N_YP,N_S,N_RL),ALBEDO(N_S,N_YP,N_RL),U_DISAG(N_U,N_UT,N_YP),&
  U(N_U,N_UT,N_YP),U_RTM(N_U,N_RL,N_YP))

  !2.4) REPLICATE LOOP FOR INTEGRATING MODEL
  DO K=1,N_RL
    !2.4.1) DISAGGREGATE FORCING
    ALLOCATE(DUMMY(N_U,N_UT,N_PF))
    DUMMY=U_NOM(1:N_U,UINDEX1:UINDEX2,1:N_PF) !UINDEX1:UINDEX2 2014/9/9 by R.Kim

    CALL DISAGGREGATE(N_YP,N_PF,N_U,N_UT,DUMMY,ELEVATION,ZF,ZS,U_DISAG,&
      ALPHA(K,1,7),ALPHA(K,1,8),MAP,GEN_SWITCH,DUMMY_INV_PRCP)

    DEALLOCATE(DUMMY)
   
 !2.4.2) PERTURB PRECIPITATION AT THE DOMAIN SCALE
    
    ! SET FORCING EQUAL TO NOMINAL DISAGGREGATED FORCING
    U=U_DISAG

    ! CALCULATE GLOBAL REPLICATE NUMBER, KGLOBAL TO INDEX F_U.  THIS VARIABLE
    !   USED IN SECTION 2.4.2 AND IN SECTION 2.4.3.4, CALL TO SSIB
    KGLOBAL=RANK*N_RL+K
    ! MULTIPLY BY F_U AND FU_MEAN
    DO T=1,N_UT
      DO P=1,N_YP
        do i=1,n_u
          U(I,T,P)=U_DISAG(I,T,P)*EXP(F_U(KGLOBAL,P,I))
        END DO
        !APPLY LIMIT TO PRECIPITATION
        IF(U_DISAG(3,T,P).GT.0.0.AND.U(3,T,P).LT.1E-9)THEN
           U(3,T,P)=1E-9
        END IF
      END DO
    END DO
   
    !2.4.3) PIXEL LOOP FOR INTEGRATING MODEL

     DO P=1,N_YP!PIXEL LOOP

      !2.4.3.1) COMPUTE INDICES
      !COMPUTE YI and YF, INITIAL AND FINAL STATE INDECES AT THIS PIXEL
      YI=N_YR*(P-1)+1  
      YF=N_YR*P
      XI=N_X*(P-1)+1
      XF=N_X*P

      !2.4.3.2) APPLY LIMITS TO INITIAL QUANTITIES, ONLY IF THERE IS 
      !          ACTUALLY SNOW (G.T. 0.01 MM) IN THIS PIXEL, REPLICATE.
      !          THE SECOND CONDITION MAKES SURE LIMITS ARE APPLIED IF
      !          SNOW DEPTH IS NEGATIVE.
      LIMIT=.TRUE.
      IF(Y0(YI,K).LT.1E-5) LIMIT=.FALSE. !TESTING FOR ZERO SNOW DEPTH
      DO I=YI,YF
        IF(Y0(I,K).LT.0.) LIMIT=.TRUE.    !IF STATE < 0, DO LIMITS
      END DO
      IF(Y0(YF,K).LT.200) LIMIT=.TRUE.   !APPLY IF TGS IS TOO LOW
      IF(Y0(YF,K).GT.285) LIMIT=.TRUE.   !APPLY IF TGS IS TOO HIGH
      IF(M.EQ.1) LIMIT=.FALSE.           !DON'T APPLY ON FIRST INTERVAL

      IF( LIMIT.EQ..TRUE. ) THEN
        Y0_RAW=Y0(YI:YF,K)
        !SNOW DEPTH
        IFLAG=0
        IF(Y0(YI,K)<0.05) THEN
          IF(Y0(YI,K)<0.0) THEN
            Y0(YI,K)=0.0
            IFLAG=IFLAG+1
          END IF
        END IF
        IF(Y0(YI,K)>4.) THEN
          Y0(YI,K)=4.
          IFLAG=IFLAG+1
        END IF 
        !SNOW DENSITY
        DO J=YI+1,YI+3
          IF(Y0(J,K)>400.)THEN
            Y0(J,K)=400.
            IFLAG=IFLAG+1
          END IF
          IF(Y0(J,K)<100.0)THEN
            Y0(J,K)=100.0
            IFLAG=IFLAG+1
          END IF
        END DO
        !SNOW TEMPERATURE
        DO J=YI+4,YI+6
          IF(Y0(J,K)>273.16)THEN !TEMP MAX = MELTING TEMP
            Y0(J,K)=273.16
            IFLAG=IFLAG+1
          END IF
          IF(Y0(J,K)<240.00)THEN
            Y0(J,K)=240.00
            IFLAG=IFLAG+1
          END IF
        END DO
        !SNOW LIQUID WATER CONTENT
        DO J=YI+7,YI+9
          IF(Y0(J,K)>0.02)THEN !
            Y0(J,K)=0.02
            IFLAG=IFLAG+1
          END IF
          IF(Y0(J,K)<0.0)THEN
            Y0(J,K)=0.0
            IFLAG=IFLAG+1
          END IF
        END DO
        !SNOW GRAIN SIZE
        DO J=YI+10,YI+12
          IF(Y0(J,K)>3.0E-3)THEN
            Y0(J,K)=3.0E-3
            IFLAG=IFLAG+1
          END IF
          IF(Y0(J,K)<0.1E-4)THEN
            Y0(J,K)=0.1E-4
            IFLAG=IFLAG+1
          END IF
        END DO
        !GROUND TEMPERATURE
        IF(Y0(YI+13,K)>285.00)THEN
          Y0(YI+13,K)=285.00
          IFLAG=IFLAG+1
        END IF
        IF(Y0(YI+13,K)<240.00)THEN
          Y0(YI+13,K)=240.00
          IFLAG=IFLAG+1
        END IF
!       IF(IFLAG.GT.3.AND.M.NE.1)THEN
!         PRINT *, 'FOUND MORE THAN 3 ILLEGAL INITIAL STATES, RESETTING I.C.',&
!           'M=',M,'RANK=',RANK,'P=',P,'K=',K,'Y0_RAW=',Y0_RAW,'YPRIOR=',&
!           YPRIOR(YI:YF,K) 
!         Y0(YI:YF,K)=YPRIOR(YI:YF,K)
!       END IF
      END IF

      !2.4.3.3) SET INITIAL CONDITION IN YT
      YT(YI:YF,1,K)=Y0(YI:YF,K)

      !2.4.3.4) CALL TO SSIB TO PROPAGATE STATES.
      CALL SSIB3(VEG_MAP(P),N_U,N_S,U(1:N_U,1:N_UT,P),ALPHA(K,P,1:6),&
        T0,X0(XI:XF,K),Y0(YI:YF,K),ZS(P,1:2),YT(YI:YF,2:N_S+1,K),P,N_YR,N_X,&
        XT(XI:XF,1:N_S,K),ALBEDO(1:N_S,P,K),TF,K,RANK,M,6,F_VEG(KGLOBAL,1),&
        VCOVER(P))


      !2.4.3.6) ISOLATE FINAL VALUES FROM TIME SERIES

      Y(YI:YF,K)=YT(YI:YF,N_S+1,K) !ISOLATE FINAL STATE VALUE
      X(XI:XF,K)=XT(XI:XF,N_S,K) !ISOLATE FINAL AUXILIARY VALUE 
      U_RTM(1:N_U,K,P)=U(1:N_U,N_UT,P) !ISOLATE FINAL FORCING VALUE
    END DO !END PIXEL LOOP
  END DO ! END REPLICATE LOOP 

  !2.5) DEALLOCATE VARIABLES USED ONLY FOR INTEGRATING MODEL
  DEALLOCATE(U,U_DISAG)
    
  !2.6) SET TIME VECTOR
  !NOTE THAT THE T VECTOR IS INDEXED AS HOUR, MONTH, DAY OF MONTH, JULIAN, YEAR

  !FIRST COMPUTE SERIAL DATE OF INITIAL TIME
  CALL SERIAL_DATE(INT(T0(5)),INT(T0(2)),INT(T0(3)),INT(T0(1)),YTIME(YINDEX1))

  !COMPUTE REST OF SERIAL DATES FROM SUM OF PREVIOUS DATE AND ONE HOUR
  DO I=2,N_UT+1
    YTIME(YINDEX1+I-1)=YTIME(YINDEX1+I-2)+1./24.
  END DO

!  CALL MPE_LOG_EVENT(10,M,'END PROP.')

  !2.7) UPDATE STATES
  IF(M<N_M+1) THEN !NO UPDATE ON LAST MEASUREMENT INTERVAL
!    CALL MPE_LOG_EVENT(11,M,'START MEAS CALC.')
    !2.7.1) COMPUTE PREDICTED OBSERVATIONS, Z
    ALLOCATE(TBRAW(N_ZP(1)*N_YP,N_RL))!,TBRAW_T(N_ZP(1)*N_YP,N_T,N_RL))    
    MONTH=12
    ALBEDO_IN=ALBEDO(N_S,1:N_YP,1:N_RL)

    CALL INTERFACEZ(Y,Z,N_RL,N_R,N_Y,N_YP,N_YR,N_Z,N_ZP,N_ZR,N_C,N_X,RTM_CTRL,&
      FREQ,THETA,X,N_A,N_U,U_RTM,MONTH,TBRAW,ALBEDO_IN,M,RANK,VEG_MAP,IERR,&
      F_VEG,BDRF,N_BDRF,MEAS_SWITCH,VCOVER)

!open(file='tbraw.out',unit=29,status='unknown')
!write(29,'(12(f10.4))') tbraw


    !2.7.1B OUTPUT SPATIALLY-AVERAGED TB FOR DIAGNOSTIC
    DO I=1,N_ZP(1)
     TB_OUTPUT(I)=SUM(Z(I,:)) 
    END DO

    CALL MPI_REDUCE(TB_OUTPUT,TB_OUTPUTALL,N_ZP(1),MPI_REAL,MPI_SUM,&
      0,MPI_COMM_WORLD,IERR)

    IF(RANK.EQ.0)THEN
      WRITE(9,'(12(F10.4))') TB_OUTPUTALL/N_R
    END IF

    DEALLOCATE(TBRAW)

   !2.7.1C PREDICTED OBSERVATIONS, Z1, FOR SINGLE FREQUENCY(37V) ! RK 11/25/2014
    DO I=1,N_RL
    Z1(1,I)=Z(11,I)
    END DO  

    ALLOCATE(ZMEAS1(1),ZUSE1(1),CV1(1),V_OUT1(1,N_R)) 
    ZMEAS1=ZMEAS(11,M)
    CV1=4.0  

    !2.7.2) SAVE PRIOR STATES
    YPRIOR=Y

    !2.7.3) COMPUTE KALMAN GAIN AND UPDATE STATE AND 'FIRST' INDEX
!    CALL MPE_LOG_EVENT(13,M,'START UPDATE')

     CALL KUPDATE(Y,Z1,ZUSE1,CV1,ZMEAS1,V_OUT1,N_RL,N_Y,1,N_R,&
       MPI_COMM_WORLD,RANK,IERR,M,MEAS_SWITCH,RANGE_U,RANGE_A,1,1,&
       N_YP,N_YR,INNOV_2MOM(M),TR_ERR(M),ABS_TR_ERR(M),REL_L1_NORM(M),&
       NU_BAR(M))

!     CALL KUPDATE(Y,Z,ZUSE,CV,ZMEAS(1:N_Z,M),V_OUT,N_RL,N_Y,N_Z,N_R,&
!       MPI_COMM_WORLD,RANK,IERR,M,MEAS_SWITCH,RANGE_U,RANGE_A,N_ZR,N_ZP,&
!       N_YP,N_YR,INNOV_2MOM(M),TR_ERR(M),ABS_TR_ERR(M),REL_L1_NORM(M),&
!       NU_BAR(M))

     !2.7.4 WRITE OUT UPDATE DATA TO FILES ON THE LOCAL NODES
!     IF(M.EQ.1)THEN
!        I=RANK+1
!        WRITE(100-I,*) 'PRIOR STATE...'
!        DO K=1,N_RL
!          WRITE(100-I,'(8750(e12.6,1x))') (YPRIOR(J,K),J=1,N_Y)
!        END DO
!        WRITE(100-I,*) 'PREDICTED MEASUREMENT...'
!        DO K=1,N_RL
!          WRITE(100-I,'(1262(e12.6,1x))') (Z(J,K),J=1,N_Z)
!        END DO
!        WRITE(100-I,*) 'RANDOM ERROR ...'
!        DO K=1,N_RL
!          KGLOBAL=RANK*N_RL+K
!          WRITE(100-I,'(1262(e12.6,1x))') (V_OUT(J,KGLOBAL),J=1,N_Z)
!        END DO
!        WRITE(100-I,*) 'POSTERIOR STATE...'
!        DO K=1,N_RL
!          WRITE(100-I,'(8750(e12.6,1x))') (Y(J,K),J=1,N_Y)
!        END DO

!       THESE LINES CAN BE UNCOMMENTED TO STOP EXECUTION AT THIS POINT
!        CALL MPI_FINALIZE(IERR)
!        PRINT *, 'STOPPING...,RANK=',RANK
!        STOP
!     END IF 

!UPDATE 'FIRST' INDEX - modified RK 2014/3/4
   IF (M.LT.N_M+1) THEN
    FIRST=MEASTIMES(M)
   END IF
!    CALL MPE_LOG_EVENT(14,M,'END UPDATE')
 END IF

  !2.8) COMPILE STATISTICS AND DEALLOCATE YT,XT
  !CALL MPE_LOG_EVENT(15,M,'START STATS')

  !NOTE THAT ONLY 0 PROCESS HAS COMPLETE STATISTICS, BUT THEN THIS IS THE ONLY
  !  ONE THAT WILL WRITE THEM OUT.
  IF(M<N_M+1)THEN
    YT(1:N_Y,N_S+2,1:N_RL)=Y !ADD POSTERIOR STATES TO YT
    ! USE A DUMMY ARG TO PASS ARRAY... NOT SURE WHY IT DOESN'T WORK OTHERWISE
    ALLOCATE(DUMMY(N_Y,N_S+2,N_YS))
    CALL COMPILEYSTATS(YT,DUMMY,N_R,N_S+2,&
      N_Y,N_YS,N_RL,RANK,MPI_COMM_WORLD,IERR)
    YSTATS(1:N_Y,YINDEX1:YINDEX2,1:N_YS)=DUMMY
    DEALLOCATE(DUMMY) 
  ELSE !FOR THE LAST MEASUREMENT INTERVAL WITH NO UPDATE
    ALLOCATE(DUMMY(N_Y,N_S+1,N_YS))
    CALL COMPILEYSTATS(YT,DUMMY,N_R,N_S+1,&
         N_Y,N_YS,N_RL,RANK,MPI_COMM_WORLD,IERR)
    YSTATS(1:N_Y,YINDEX1:YINDEX2,1:N_YS)=DUMMY
    DEALLOCATE(DUMMY) 
  END IF

!  CALL MPE_LOG_EVENT(16,M,'END STATS')

  !2.10) DEALLOCATIONS
  DEALLOCATE(YT,XT,ALBEDO,U_RTM)

END DO  !End measurement loop

!3) COMPUTE ELAPSED TIME
IF(RANK.EQ.0) THEN
  ENDTIME=MPI_WTIME()
END IF

! 4) ON O NODE, WRITE YSTATS AND YTIME FROM FILTER RUN TO FILES
IF(RANK.EQ.0)THEN
  PRINT *, 'WRITING YSTATS FILE'
  !4.1 DEFINE THE TOUTPUT VECTOR, WHICH CONTAINS THE HOURS AT WHICH YSTATS
  !      VECTOR WILL BE WRITTEN TO FILE
  J=0
  DO I=1,24,4 ! I=1,5,9,13,17,21
    J=J+1
    TOUTPUT(J)=DBLE(I)/24.
  END DO
  
  !4.2 OPEN OUTPUT FILES
  OPEN(FILE='ystats.out',UNIT=13,STATUS='UNKNOWN',FORM='UNFORMATTED')
  OPEN(FILE='ytime.out',UNIT=10,STATUS='UNKNOWN')
  OPEN(FILE='ytime_all.out',UNIT=11,STATUS='UNKNOWN')
  OPEN(FILE='innov_stats.out',UNIT=12,STATUS='UNKNOWN')
!  OPEN(FILE='tbraw_t.out',UNIT=14,STATUS='UNKNOWN')
 
 !4.3 LOOP OVER DATA, WRITING OUTPUT CORRESPONDING TO TOUTPUT
  DO T=1,N_M+N_T+1
    !ASSUME THAT WE ARE NOT AT AN OUTPUT TIME, THEN IF THE FRACTION PART OF 
    !  THE SERIAL DATE IN YTIME MATCHES AN OUTPUT TIME, SET OUTPUT = .TRUE.
    OUTPUT=.FALSE. 
    DO I=1,6
      ! THIS REPRESENTS THE DIFFERENCE BETWEEN THE FRACTION PORTION OF YTIME
      !   AND THE TOUTPUT TIME.  TO BE USED IN THE TEST OF WHETHER IT IS AN
      !   OUTPUT TIME.
      DIFF=ABS(YTIME(T)-INT(YTIME(T))-TOUTPUT(I))
      IF(DIFF.LT.1E-5) OUTPUT=.TRUE.
    END DO
    IF(OUTPUT.EQ..TRUE.) THEN
      DO I=1,4
        WRITE(13)  (YSTATS(J,T,I),J=1,N_Y)
      END DO
      WRITE(10,*) YTIME(T)
    END IF
    WRITE(11,*) YTIME(T) !THIS CONTAINS ALL OF THE TIME INDECES
  END DO
  !4.4 WRITE OUT SOME INNOVATION STATISTICS
  DO I=1,N_M
    WRITE(12,'(5(E14.6))') INNOV_2MOM(I),TR_ERR(I),ABS_TR_ERR(I),&
      REL_L1_NORM(I),NU_BAR(I)
  END DO
  !4.5 WRITE OUT ZOUT AND TBRAW
!  DO I=1,N_RL
!    WRITE(9) (Z(J,I),J=1,N_Z)
!  END DO
!  DO I=1,N_YP*N_ZP(1)
!    WRITE(14,'(F10.3,1X,F10.3,1X,F10.3,1X,F10.3)') (TBRAW(I,J),J=1,N_RL)
!  END DO 
  CLOSE(9) !CLOSE z output file
  CLOSE(10)!CLOSE YTIME FILE
  CLOSE(11)!CLOSE FILE FOR ALL YTIME INDECES
  CLOSE(12)!CLOSE INNOVATIONS FILE
  CLOSE(13)!CLOSE YSTATISTICS FILE
END IF

CLOSE(15) !FILE FOR MEASUREMENT INNOVATIONS
close(14)
close(29)

!CALL MPE_FINISH_LOG('enkf_mpe')
CALL MPI_FINALIZE(IERR)

END PROGRAM ENKF