!!  Copyright (C) Stichting Deltares, 2005-2020.
!!
!!  This file is part of iMOD.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!  Contact: imod.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands.
!!
MODULE MOD_IPEST_IES

USE WINTERACTER
USE MOD_IDF_PAR
USE MOD_IDF, ONLY : IDFREAD,IDFREADSCALE,IDFALLOCATEX,IDFNULLIFY,IDFWRITE,IDFDEALLOCATEX,IDFGETLOC,IDFCOPY
USE IMODVAR, ONLY : DP_KIND,IDPROC
USE MOD_PMANAGER_PAR, ONLY : PEST,PBMAN,PRJIDF,PRJNLAY
USE MOD_UTL, ONLY  : UTL_GETUNIT,UTL_CAP,ITOS,RTOS,UTL_GETMED,UTL_CREATEDIR,UTL_GOODNESS_OF_FIT,UTL_NASH_SUTCLIFFE,UTL_STDEF, &
 UTL_GETGAMMA,UTL_MATMUL
USE MOD_IPEST_GLM_PAR
USE MOD_IPEST_GLM, ONLY : IPEST_GLM_GETJ,IPEST_GLM_ALLOCATEMSR,IPEST_GLM_RUNMODELS,IPEST_GLM_SETGROUPS
USE MOD_PMANAGER_UTL, ONLY : PMANAGER_SAVEMF2005_MOD_U2DREL,PMANAGER_SAVEMF2005_DEALLOCATE
USE MOD_OSD
USE MOD_IPEST_IES_PAR
USE MOD_LUDCMP
USE MOD_BATCH_UTL, ONLY : IMODBATCH_CREATEENSEMBLES_CHOLESKY

CONTAINS

 !#####=================================================================
 SUBROUTINE IPEST_IES_MAIN(DIR,MNAME,IBATCH)
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IBATCH
 CHARACTER(LEN=*),INTENT(IN) :: DIR,MNAME
 INTEGER :: ITER,IGRAD,ILOG,ILAMBDA,JLAMBDA
 REAL(KIND=DP_KIND) :: MINJ,L
 
 !## log10 transform
 ILOG=1
 
 IF(.NOT.IPEST_IES_ALLOCATE(DIR))RETURN 
 IF(.NOT.IPEST_IES_CREATEENSEMBLES(DIR,ILOG))RETURN
 
 !## organise groups
 CALL IPEST_GLM_SETGROUPS()

 CALL UTL_CREATEDIR(TRIM(DIR)//'\STATHEADS')
 CALL UTL_CREATEDIR(TRIM(DIR)//'\ENSEMBLES')

 CALL WMESSAGEENABLE(TIMEREXPIRED,1); MSR%PJ=HUGE(1.0D0); ITER=0
 DO WHILE (ITER.LE.PEST%PE_MXITER)
  
  !## no residual output
  IUPESTRESIDUAL=0
  !## initial run - no lambda testing
  IF(ITER.EQ.0)THEN
   !## read from idf file and save the realisations in lpf package - apply the correct IDF files and export them in for the realizations
   JLAMBDA=1; CALL IPEST_IES_SAVEREALS(DIR,ITER,JLAMBDA,MNAME)
   !## run and process all realizations-testing
   CALL IPEST_GLM_RUNMODELS(IBATCH,RNG,GPARAM,'R',DIR,MNAME,ITER,LAMBDA)
  !## carry out the lambda tests
  ELSE
   !## simulate different lambdas
   MINJ=10.0D20; JLAMBDA=0; DO ILAMBDA=1,PBMAN%NLINESEARCH
    !## read from idf file and save the realisations in lpf package - apply the correct IDF files and export them in for the realizations
    CALL IPEST_IES_SAVEREALS(DIR,ITER,ILAMBDA,MNAME)
    L=LAMBDA*PBMAN%LAMBDA_TEST(ILAMBDA)
    !## run and process all realizations-testing
    WRITE(IUPESTPROGRESS,'(/A5,4A15,15X,A15)') 'RS','REALIZATION','LAMBDA','TOT_J','RED_J','SIM_TIME(SEC)' 
    CALL IPEST_GLM_RUNMODELS(IBATCH,RNG,GPARAM,'R',DIR,MNAME,ITER,L)
    !## evaluate which lambda is best
    IF(SUM(JE(:,1:PEST%NREALS)).LT.MINJ)THEN
     MINJ=SUM(JE(:,1:PEST%NREALS)); JLAMBDA=ILAMBDA
    ENDIF
    CALL IPEST_IES_PROGRESS_OUT(ITER)
   ENDDO
   IF(JLAMBDA.EQ.0)THEN; STOP 'SOMETHING WENT WRONG: JLAMBDA=0'; ENDIF
   !## get appropriate msr%dhg() results
  ENDIF
  
  !## get update of realisations
  IF(.NOT.IPEST_IES_UPDATE_ENSEMBLES(DIR,ITER,JLAMBDA))EXIT
!  !## save mean/stdev heads
!  CALL IPEST_IES_SAVE_RESULT_STATS(DIR,ITER,JLAMBDA)

! !## save residuals
! IF(IUPESTRESIDUAL.GT.0)CLOSE(IUPESTRESIDUAL)
! IUPESTRESIDUAL=UTL_GETUNIT(); OPEN(IUPESTRESIDUAL,FILE=TRIM(DIR)//'\IPEST\LOG_PEST_RESIDUAL_'// &
!   TRIM(ITOS(ITER))//'.TXT',STATUS='UNKNOWN',ACTION='WRITE')
! IF(.NOT.IPEST_GLM_GETJ(DIR,IJ,LPARAM(IJ),'L',IBATCH))RETURN
! CLOSE(IUPESTRESIDUAL)
  
  !## next cycle
  ITER=ITER+1

 ENDDO
 
 !## finished write summary of objective function values
 WRITE(IUPESTOUT,'(A10,999(1X,A15))') 'ITER','TOTAL',('OBJF_SIM'//TRIM(ITOS(IGRAD)),IGRAD=1,PEST%NREALS)
 DO ITER=0,PEST%PE_MXITER
  WRITE(IUPESTOUT,'(I10,999(1X,F15.7))') ITER,SUM(JE(1:PEST%NREALS,ITER)),(JE(IGRAD,ITER),IGRAD=1,PEST%NREALS)
 ENDDO

 CALL WMESSAGETIMER(0); CALL WMESSAGEENABLE(TIMEREXPIRED,0)
 PBMAN%IIES=0; CALL PMANAGER_SAVEMF2005_DEALLOCATE(); CALL IPEST_IES_DEALLOCATE
 
 END SUBROUTINE IPEST_IES_MAIN

 !#####=================================================================
 LOGICAL FUNCTION IPEST_IES_CREATEENSEMBLES(DIR,ILOG)!,IBATCH)
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 REAL(KIND=DP_KIND) :: K1,K2
 INTEGER,INTENT(IN) :: ILOG
 INTEGER :: I,J
 TYPE(IDFOBJ) :: STDEV,MEAN
 
 IPEST_IES_CREATEENSEMBLES=.FALSE.
 CALL IDFCOPY(PRJIDF,STDEV); CALL IDFCOPY(PRJIDF,MEAN)
 
 DO I=1,SIZE(PEST%PARAM)
  IF(ILOG.EQ.1)THEN
   K1=LOG10(ABS(PEST%PARAM(I)%PMIN))
   K2=LOG10(ABS(PEST%PARAM(I)%PMAX))
  ELSE
   K1=PEST%PARAM(I)%PMIN
   K2=PEST%PARAM(I)%PMAX
  ENDIF
  PEST%PARAM(I)%PARMEAN=(K1+K2)/2.0D0
  PEST%PARAM(I)%PARSTD=(K2-K1)/4.0D0
  WRITE(*,'(A,4F15.7)')  'LOG10: MIN,MAX,MEAN,STDEV',K1,K2,PEST%PARAM(I)%PARMEAN,PEST%PARAM(I)%PARSTD
  WRITE(*,'(A,4F15.7)')  '       MIN,MAX,MEAN,STDEV',10.0D0**K1,10.0D0**K2,10.0D0**PEST%PARAM(I)%PARMEAN,10.0D0**PEST%PARAM(I)%PARSTD
  IF(PEST%ICREATEENSEMBLES.EQ.1)THEN
   STDEV%X=PEST%PARAM(I)%PARSTD; STDEV%NODATA=-999.99D0
   MEAN%X=PEST%PARAM(I)%PARMEAN; MEAN%NODATA=-999.99D0
   CALL IMODBATCH_CREATEENSEMBLES_CHOLESKY(STDEV,MEAN,TRIM(DIR)//'\ENSEMBLES', &
       PEST%PARAM(I)%PARRANGE,PEST%NREALS,'PARAMETER'//TRIM(ITOS(I)),ILOG)
   DO J=1,PEST%NREALS
    PEST%PARAM(I)%REALSFNAME(J)=TRIM(DIR)//'\ENSEMBLES\PARAMETER'//TRIM(ITOS(I))//'_R'//TRIM(ITOS(J))//'.IDF'
   ENDDO
!  PEST%PARAM(I)%ICREATEENSEMBLES
  ENDIF
 ENDDO
 CALL IDFDEALLOCATEX(STDEV); CALL IDFDEALLOCATEX(MEAN)
 
 IPEST_IES_CREATEENSEMBLES=.TRUE.
 
 END FUNCTION IPEST_IES_CREATEENSEMBLES

 !#####=================================================================
 SUBROUTINE IPEST_IES_SAVE_RESULT_STATS(DIR,ITER,JLAMBDA)
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,INTENT(IN) :: ITER,JLAMBDA
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: TMP1
 TYPE(IDFOBJ),ALLOCATABLE,DIMENSION(:) :: HDS
 TYPE(IDFOBJ) :: MHD,VHD
 REAL(KIND=DP_KIND) :: AVG,STD
 INTEGER :: I,IROW,ICOL,N,ILAY
 
 !## compute mean and variance heads
 ALLOCATE(HDS(PEST%NREALS)); DO I=1,PEST%NREALS; CALL IDFNULLIFY(HDS(I)); ENDDO; CALL IDFNULLIFY(MHD); CALL IDFNULLIFY(VHD)
 DO ILAY=1,PRJNLAY
  DO I=1,PEST%NREALS
   HDS(I)%FNAME=TRIM(DIR)//'\IIES_R#'//TRIM(ITOS(I))//'_L#'//TRIM(ITOS(JLAMBDA))//'\HEAD\HEAD_STEADY-STATE_L'//TRIM(ITOS(ILAY))//'.IDF'
   IF(.NOT.IDFREAD(HDS(I),HDS(I)%FNAME,1))STOP
  ENDDO
  IF(ILAY.EQ.1)THEN; CALL IDFCOPY(HDS(1),MHD); CALL IDFCOPY(HDS(1),VHD); ENDIF
  ALLOCATE(TMP1(PEST%NREALS)); DO IROW=1,MHD%NROW; DO ICOL=1,MHD%NCOL
   TMP1(1)=0.0D0; DO I=1,PEST%NREALS; TMP1(I)=HDS(I)%X(ICOL,IROW); ENDDO
   CALL UTL_STDEF(TMP1,PEST%NREALS,HDS(1)%NODATA,STD,AVG,N)
   MHD%X(ICOL,IROW)=AVG; VHD%X(ICOL,IROW)=STD
  ENDDO; ENDDO
  DEALLOCATE(TMP1)
  CALL UTL_CREATEDIR(TRIM(DIR)//'\STATHEADS')
  MHD%FNAME=TRIM(DIR)//'\STATHEADS\MEAN_HEAD_STEADY-STATE_L'//TRIM(ITOS(ILAY))//'_ITER_'//TRIM(ITOS(ITER))//'.IDF'
  IF(.NOT.IDFWRITE(MHD,MHD%FNAME,1))STOP
  VHD%FNAME=TRIM(DIR)//'\STATHEADS\STDEV_HEAD_STEADY-STATE_L'//TRIM(ITOS(ILAY))//'_ITER_'//TRIM(ITOS(ITER))//'.IDF'
  IF(.NOT.IDFWRITE(VHD,VHD%FNAME,1))STOP
 ENDDO
 
 END SUBROUTINE IPEST_IES_SAVE_RESULT_STATS
 
 !#####=================================================================
 LOGICAL FUNCTION IPEST_IES_UPDATE_ENSEMBLES(DIR,ITER,JLAMBDA) 
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,INTENT(INOUT) :: ITER
 INTEGER,INTENT(IN) :: JLAMBDA
 INTEGER :: I,II,J,JJ,K,N,SPACEDIM,TIMEDIM,OBSDIM,IROW,ICOL,IX,ILOG,ILAMBDA
 REAL(KIND=DP_KIND) :: X,VAR,A,L,TL,MAXTL,GAMMA,F1,F2,F3,J0,MLT
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: DM 
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: SO,O,TMP1,PCLASS,RCLASS
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: REALS
 INTEGER,ALLOCATABLE,DIMENSION(:,:) :: PDF 
 CHARACTER(LEN=256) :: FNAME
  
 WRITE(*,'(/A/)') 'Update ensembles ...'
 
 !## apply log transformation
 ILOG=1
 GAMMA=4.0D0
 MAXTL=100.0D0
 IPEST_IES_UPDATE_ENSEMBLES=.TRUE.

 SPACEDIM=PRJIDF%NROW*PRJIDF%NCOL
 TIMEDIM=PEST%NREALS
 OBSDIM=MSR%NOBS

 !## log10 parameter classes for pdf
 X=-3.0D0; N=1; DO; N=N+1; X=X+0.25D0; IF(X.GT.3.0D0)EXIT; ENDDO
 ALLOCATE(PCLASS(N)); PCLASS=0.0D0
 PCLASS(1)=-3.0D0; N=1; DO; N=N+1; PCLASS(N)=PCLASS(N-1)+0.25D0; IF(PCLASS(N).GT.3.0D0)EXIT; ENDDO
 !## residual classes for pdf
 X=-10.5D0; N=1; DO; N=N+1; X=X+0.5D0; IF(X.GT.10.0D0)EXIT; ENDDO
 ALLOCATE(RCLASS(N)); RCLASS=0.0D0
 RCLASS(1)=-10.5D0; N=1; DO; N=N+1; RCLASS(N)=RCLASS(N-1)+0.5D0; IF(RCLASS(N).GT.10.0D0)EXIT; ENDDO
 IF(ALLOCATED(PDF))DEALLOCATE(PDF); ALLOCATE(PDF(SIZE(RCLASS),TIMEDIM))

 ALLOCATE(TMP1(OBSDIM)) 

 !## get residual objective function value
 ALLOCATE(SO(TIMEDIM),O(TIMEDIM)); SO=0.0D0; O=0.0D0
 !## use weigth (get variance from weights) - get matrix of objective function values
 DO J=1,TIMEDIM
  TMP1=0.0D0
  DO K=1,OBSDIM
   TMP1(K)=MSR%DHG(J,K) 
   VAR=1.0D0/MSR%W(K); SO(J)=SO(J)+MSR%DHG(J,K)**2.0D0*(1.0D0/VAR)
   O(J)=O(J)+MSR%DHG(J,K)**2.0D0
   
   !## volgens mij moet error parameters er nog bij
   
  ENDDO
  !## get pdf of error distribution
  CALL IPEST_IES_GETPDF(TMP1,PDF(:,J),RCLASS) 
 ENDDO

 !## get current sum of weighted objective function values
 MSR%TJ=SUM(SO); MSR%J=SUM(O)/DBLE(TIMEDIM)/DBLE(OBSDIM)

 !## write total objective function value
 WRITE(IUPESTOUT,'(27A1)') ('-',I=1,27)
 WRITE(IUPESTOUT,'(I5,7X,F15.7)') ITER,MSR%TJ

 CALL IPEST_IES_PROGRESS_OUT(ITER)
 
 WRITE(IUPESTOUT,'(/A/)') 'PDF OBJECTIVE FUNCTION VALUES'
 WRITE(IUPESTOUT,'(A15,999(1X,A5))') 'CLASS','ALL',('R'//TRIM(ITOS(II)),II=1,TIMEDIM)
 DO JJ=1,SIZE(PDF,1)
  WRITE(IUPESTOUT,'(F15.7,999(1X,I5))') RCLASS(JJ),SUM(PDF(JJ,1:TIMEDIM)),(PDF(JJ,II),II=1,TIMEDIM) 
 ENDDO
 DEALLOCATE(TMP1)

 !## determine initial lambda
 IF(ITER.EQ.0)THEN; X=MSR%TJ/DBLE(2.0D0*OBSDIM); X=LOG10(X); IX=FLOOR(X); LAMBDA=10.0D0**IX; ENDIF
 
 !## evaluate progression
 IF(ITER.GT.0)THEN

  J0=SUM(JE(1:PEST%NREALS,0)); F3=J0
  
  !## do something with statistics, mean and stdev
  IF(ITER.GT.1)F3=MSR%PJ
  F1=100.0D0*(1.0D0-MSR%TJ/F3)
  F2=100.0D0*(1.0D0-MSR%TJ/J0)
  WRITE(IUPESTEFFICIENCY,'(I10,99(1X,F15.3))') ITER,MSR%TJ,F3,MSR%J,F1,F2,LAMBDA, &
     (PEST%PARAM(J)%MEAN(TIMEDIM+1),PEST%PARAM(J)%STD(TIMEDIM+1),J=1,SIZE(PEST%PARAM))

  WRITE(IUPESTOUT,'(/2A10,4A15)') 'ITER','ENSEMBLE','MEAN','STDEV','10^MEAN','10^STDEV'
  DO J=1,SIZE(PEST%PARAM)
   DO I=1,TIMEDIM
    WRITE(IUPESTOUT,'(2I10,4F15.7)') I,J,PEST%PARAM(J)%MEAN(I),PEST%PARAM(J)%STD(I), &
     10.0D0**PEST%PARAM(J)%MEAN(I),10.0D0**PEST%PARAM(J)%STD(I)
   ENDDO
   WRITE(IUPESTOUT,'(80A1)') ('-',K=1,80) 
   I=TIMEDIM+1
   WRITE(IUPESTOUT,'(20X,4F15.7/)') PEST%PARAM(J)%MEAN(I),PEST%PARAM(J)%STD(I), &
     10.0D0**PEST%PARAM(J)%MEAN(I),10.0D0**PEST%PARAM(J)%STD(I)
  ENDDO

  IF(MSR%TJ.LE.F3)THEN 
   !## too less improvement
   IF(1.0D0-MSR%TJ/F3.LT.PEST%PE_STOP)THEN
    IPEST_IES_UPDATE_ENSEMBLES=.FALSE.; RETURN
   ELSE
    !## YES!!! Continue to next iteration
    LAMBDA=LAMBDA/GAMMA
    WRITE(IUPESTOUT,'(/A/)') 'Awesome, significant reduction of the objective function reduction: proceed'
    !## save previous objective function value
    MSR%PJ=MSR%TJ
   ENDIF
  ELSE
   LAMBDA=LAMBDA*GAMMA
   !## reset ensembles to previous iteration
   ITER=ITER-1
   WRITE(IUPESTOUT,'(/A/)') 'No objective function reduction: reset to previous iteration'
  ENDIF
 ENDIF
 
 !## try to find acceptable ensemble updates for all lambda to be tested
 DO ILAMBDA=1,PBMAN%NLINESEARCH
  MLT=PBMAN%LAMBDA_TEST(ILAMBDA)
  DO
 
!   A=LAMBDA+1.0D0
   TL=0.0D0; A=MLT*LAMBDA+1.0D0

   !## processing the inverse of the prior parameter covariance matrix
   DO I=1,SIZE(PEST%PARAM)

    IF(ILAMBDA.EQ.1)THEN
     !## initialisation
     IF(ITER.EQ.0)CALL IPEST_IES_INITIATE_INVERSE_MATRICES(DIR,I,TIMEDIM,SPACEDIM)
  
     !## compute dm-matrix
     IF(ALLOCATED(REALS))DEALLOCATE(REALS); IF(ALLOCATED(DM))DEALLOCATE(DM)
     ALLOCATE(REALS(TIMEDIM,SPACEDIM),DM(TIMEDIM,SPACEDIM))
     !## compute dm based upon last succesfull set of parameters
     CALL IPEST_IES_COMPUTE_DM(DIR,I,TIMEDIM,SPACEDIM,REALS,DM,ITER,ILOG,JLAMBDA) 

     !## store prior-estimated parameters - need to be done per parameter
     IF(ITER.EQ.0)THEN
      !## might be coming here twice as improvement of objective function is not enough
      IF(.NOT.ASSOCIATED(PEST%PARAM(I)%MPR))THEN
       ALLOCATE(PEST%PARAM(I)%MPR(TIMEDIM,SPACEDIM)); PEST%PARAM(I)%MPR=REALS
       !## compute am-matrix initially through svd - this is specific per parameter
       CALL IPEST_IES_COMPUTE_AM(I,TIMEDIM,SPACEDIM,DM)
      ENDIF
     ENDIF
    
    ELSE
     !## read realisatons again
     CALL IPEST_IES_COMPUTE_GETREALS(DIR,I,TIMEDIM,SPACEDIM,REALS,ITER,ILOG,JLAMBDA) 
    ENDIF
    
    !## processing the inverse of the prior parameter covariance matrix
    CALL IPEST_IES_COMPUTE_PARTIAL_M(ITER,I,TIMEDIM,SPACEDIM,OBSDIM,DM,A,REALS,L); TL=TL+L
    WRITE(IUPESTOUT,'(/A,F15.7/)') 'LENGTH OF INDIVIDUAL ENSEMBLE UPDATE=',L

    !## save the updates for the next iteration
    IF(.NOT.IDFALLOCATEX(PRJIDF))STOP
    ALLOCATE(TMP1(SPACEDIM)) 

    IF(ALLOCATED(PDF))DEALLOCATE(PDF); ALLOCATE(PDF(SIZE(PCLASS),TIMEDIM))
    DO K=1,TIMEDIM
     II=0; DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL; II=II+1
      TMP1(II)=REALS(K,II)
      IF(ILOG.EQ.1)THEN; PRJIDF%X(ICOL,IROW)=10.0D0**TMP1(II) 
      ELSE; PRJIDF%X(ICOL,IROW)=TMP1(II); ENDIF
     ENDDO; ENDDO
     !## get statistics (variance is outcome)
     CALL UTL_STDEF(TMP1,SPACEDIM,PRJIDF%NODATA,PEST%PARAM(I)%STD(K),PEST%PARAM(I)%MEAN(K),N)
     !## get pdf of parameters
     CALL IPEST_IES_GETPDF(TMP1,PDF(:,K),PCLASS) 
     !## get name of ensemble
     CALL IPEST_IES_GETENSEMBLENAME(ITER+1,DIR,PEST%PARAM(I)%REALSFNAME(K),FNAME,ILAMBDA)
     !## save ensemble for current iteration and lambda test
     IF(.NOT.IDFWRITE(PRJIDF,FNAME,1))STOP
    ENDDO
    PEST%PARAM(I)%MEAN(TIMEDIM+1)=SUM(PEST%PARAM(I)%MEAN(1:TIMEDIM))/REAL(TIMEDIM,8)
    PEST%PARAM(I)%STD(TIMEDIM+1) =SUM(PEST%PARAM(I)%STD(1:TIMEDIM)) /REAL(TIMEDIM,8)
    WRITE(IUPESTOUT,'(/A/)') 'PDF PARAMETERS '//TRIM(ITOS(I))
    WRITE(IUPESTOUT,'(A15,999(1X,A5))') 'CLASS','ALL',('R'//TRIM(ITOS(II)),II=1,TIMEDIM)
    DO JJ=1,SIZE(PDF,1)
     WRITE(IUPESTOUT,'(F15.7,999(1X,I5))') 10.0D0**PCLASS(JJ),SUM(PDF(JJ,1:TIMEDIM)),(PDF(JJ,II),II=1,TIMEDIM) 
    ENDDO
  
    DEALLOCATE(TMP1); CALL IDFDEALLOCATEX(PRJIDF)  

   ENDDO
   
   !## do not allow enormous updates
   IF(TL.LT.MAXTL)EXIT
  
   LAMBDA=LAMBDA*GAMMA
   WRITE(IUPESTOUT,'(/A,F15.7/)') 'No acceptable ensemble update (too large), increasing lambda ',LAMBDA

  ENDDO
 ENDDO

! !## possible update vector is less that criterion
! WRITE(IUPESTOUT,'(/A,F15.7/)') 'TOTAL LENGTH OF ENSEMBLE UPDATE=',TL
! IF(TL.LT.PEST%PE_PADJ)IPEST_IES_UPDATE_ENSEMBLES=.FALSE.
  
 DEALLOCATE(SO,DM,REALS)
 
 FLUSH(IUPESTOUT); FLUSH(IUPESTPROGRESS); FLUSH(IUPESTEFFICIENCY)
 
 END FUNCTION IPEST_IES_UPDATE_ENSEMBLES
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_PROGRESS_OUT(ITER)
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ITER
 REAL(KIND=DP_KIND) :: F1,F2
 INTEGER :: I
 
 WRITE(IUPESTPROGRESS,'(95A1)') ('-',I=1,95)
 F1=SUM(JE(1:PEST%NREALS,ITER))
 IF(ITER.EQ.0)THEN
  WRITE(IUPESTPROGRESS,'(5X,A15,15X,F15.3)') 'TOTAL',F1
 ELSE
  F2=SUM(JE(1:PEST%NREALS,ITER-1))
  WRITE(IUPESTPROGRESS,'(5X,A15,15X,2F15.3)') 'TOTAL',F1,F1-F2
 ENDIF
 FLUSH(IUPESTPROGRESS)
 
 END SUBROUTINE IPEST_IES_PROGRESS_OUT

 !#####=================================================================
 SUBROUTINE IPEST_IES_INITIATE_INVERSE_MATRICES(DIR,IPARAM,TIMEDIM,SPACEDIM)
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,INTENT(IN) :: IPARAM,TIMEDIM,SPACEDIM
 INTEGER :: IU,N,J,K,IROW,ICOL,IROW2,ICOL2
 TYPE(IDFOBJ) :: IDF
 CHARACTER(LEN=1) :: YN
 REAL(KIND=DP_KIND) :: X1,Y1,X2,Y2,GAMMA,XCOR
 LOGICAL :: L1,L2
 
 !## carried out only once
 IF(ASSOCIATED(PEST%PARAM(IPARAM)%INVCOV))RETURN

 ALLOCATE(PEST%PARAM(IPARAM)%MEAN(TIMEDIM+1)); PEST%PARAM(IPARAM)%MEAN=0.0D0
 ALLOCATE(PEST%PARAM(IPARAM)%STD(TIMEDIM+1));  PEST%PARAM(IPARAM)%STD=0.0D0
 IF(TRIM(PEST%PARAM(IPARAM)%ICOVFNAME).NE.'')THEN
  IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=PEST%PARAM(IPARAM)%ICOVFNAME,STATUS='OLD',ACTION='READ')
  READ(IU,*); READ(IU,*) N; IF(N.NE.SPACEDIM)STOP 'N.NE.SPACEDIM IN COVARIANCE MATRIX'; ALLOCATE(PEST%PARAM(IPARAM)%COV(SPACEDIM,SPACEDIM))
  DO J=1,SPACEDIM; READ(IU,*) (PEST%PARAM(IPARAM)%COV(J,K),K=1,SPACEDIM); ENDDO; CLOSE(IU)
 ELSE
  !## fill in the covariance matrix with covariates using the exponential variogram with range ... sill=1.0
  ALLOCATE(PEST%PARAM(IPARAM)%COV(SPACEDIM,SPACEDIM))
  PEST%PARAM(IPARAM)%COV=0.0D0
  DO J=1,SPACEDIM
   PEST%PARAM(IPARAM)%COV(J,J)=PEST%PARAM(IPARAM)%PARSTD**2.0D0
  ENDDO
  IROW=1; ICOL=0; DO J=1,SPACEDIM
   ICOL=ICOL+1; IF(ICOL.GT.PRJIDF%NCOL)THEN; ICOL=1; IROW=IROW+1; ENDIF
   CALL IDFGETLOC(PRJIDF,IROW,ICOL,X1,Y1)
   IROW2=IROW; ICOL2=ICOL; DO K=J+1,SPACEDIM
    ICOL2=ICOL2+1; IF(ICOL2.GT.PRJIDF%NCOL)THEN; ICOL2=1; IROW2=IROW2+1; ENDIF
    CALL IDFGETLOC(PRJIDF,IROW2,ICOL2,X2,Y2)
    !## gets the variogram values
    GAMMA=UTL_GETGAMMA(X1,Y1,X2,Y2,PEST%PARAM(IPARAM)%PARRANGE,1.0D0,0.0D0,3)
    XCOR=1.0D0-GAMMA
    PEST%PARAM(IPARAM)%COV(J,K)=SQRT(PEST%PARAM(IPARAM)%COV(J,J))*SQRT(PEST%PARAM(IPARAM)%COV(K,K))*XCOR
    PEST%PARAM(IPARAM)%COV(K,J)=PEST%PARAM(IPARAM)%COV(J,K)
   ENDDO
  ENDDO
  N=SPACEDIM
 ENDIF
 !## save as IDF (for fun)
 CALL IDFNULLIFY(IDF)
 IDF%DX=1.0D0; IDF%DY=1.0D0; IDF%NCOL=N; IDF%NROW=N
 IDF%XMIN=0.0D0; IDF%XMAX=IDF%XMIN+IDF%NCOL*IDF%DX
 IDF%YMIN=0.0D0; IDF%YMAX=IDF%YMIN+IDF%NROW*IDF%DY
 IF(.NOT.IDFALLOCATEX(IDF))RETURN
 DO J=1,N; DO K=1,N; IDF%X(J,K)=PEST%PARAM(IPARAM)%COV(J,K); ENDDO; ENDDO
 IDF%FNAME=TRIM(DIR)//'\ENSEMBLES\COV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF'; IF(.NOT.IDFWRITE(IDF,IDF%FNAME,1))STOP
 CALL IDFDEALLOCATEX(IDF)
  
 ALLOCATE(PEST%PARAM(IPARAM)%ISQRTCOV(N,N),PEST%PARAM(IPARAM)%INVCOV(N,N))
 IF(.NOT.IDFALLOCATEX(IDF))RETURN
 
 INQUIRE(FILE=TRIM(DIR)//'\ENSEMBLES\ISQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF',EXIST=L1)
 INQUIRE(FILE=TRIM(DIR)//'\ENSEMBLES\ICOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF',EXIST=L2)
 IF(L1.AND.L2)THEN
  WRITE(*,'(/A$)') 'Found existing ISQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF and ICOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF, do you want to use them ?'
  READ(*,'(A1)') YN
  IF(YN.EQ.'Y'.OR.YN.EQ.'y')THEN
   IF(.NOT.IDFREAD(IDF,TRIM(DIR)//'\ENSEMBLES\ISQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF',1))STOP
   PEST%PARAM(IPARAM)%ISQRTCOV=IDF%X
   IF(.NOT.IDFREAD(IDF,TRIM(DIR)//'\ENSEMBLES\ICOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF',1))STOP
   PEST%PARAM(IPARAM)%INVCOV=IDF%X
  ELSE
   L1=.FALSE.; L2=.FALSE.
  ENDIF
 ENDIF
 IF(.NOT.L1.OR..NOT.L2)THEN
  CALL LUDCMP_CALC_SQRTROOTINVERSE(N,PEST%PARAM(IPARAM)%COV,PEST%PARAM(IPARAM)%ISQRTCOV,PEST%PARAM(IPARAM)%INVCOV)
  DO J=1,N; DO K=1,N; IDF%X(J,K)=PEST%PARAM(IPARAM)%ISQRTCOV(J,K); ENDDO; ENDDO
  IDF%FNAME=TRIM(DIR)//'\ENSEMBLES\ISQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF'; IF(.NOT.IDFWRITE(IDF,IDF%FNAME,1))STOP
  DO J=1,N; DO K=1,N; IDF%X(J,K)=PEST%PARAM(IPARAM)%INVCOV(J,K); ENDDO; ENDDO
  IDF%FNAME=TRIM(DIR)//'\ENSEMBLES\ICOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF'; IF(.NOT.IDFWRITE(IDF,IDF%FNAME,1))STOP
 ENDIF
 CALL IDFDEALLOCATEX(IDF)
 DEALLOCATE(PEST%PARAM(IPARAM)%COV)
 
 END SUBROUTINE IPEST_IES_INITIATE_INVERSE_MATRICES
  
 !#####=================================================================
 SUBROUTINE IPEST_IES_GETPDF(X,PDF,XCLASS) !,CT)
 !#####=================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),DIMENSION(:) :: XCLASS
! CHARACTER(LEN=1),INTENT(IN) :: CT
 REAL(KIND=DP_KIND),DIMENSION(:) :: X
 INTEGER,DIMENSION(:) :: PDF
 REAL(KIND=DP_KIND) :: LX 
 INTEGER :: I,J 
 
 PDF=0; DO I=1,SIZE(X) !SPACEDIM
  DO J=1,SIZE(XCLASS)-1
!   IF(CT.EQ.'L')THEN
!    LX=LOG10(X(I))
!   ELSE
   LX=X(I)
!   ENDIF
   IF(J.EQ.1)THEN
    IF(LX.LT.XCLASS(J))PDF(J)=PDF(J)+1
   ELSEIF(J.EQ.SIZE(XCLASS))THEN
    IF(LX.GE.XCLASS(J))PDF(SIZE(PDF))=PDF(SIZE(PDF))+1
   ELSE
    IF(LX.GE.XCLASS(J).AND.LX.LT.XCLASS(J+1))THEN
     PDF(J)=PDF(J)+1; EXIT
    ENDIF
   ENDIF
  ENDDO
 ENDDO
 
 END SUBROUTINE IPEST_IES_GETPDF
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_COMPUTE_PARTIAL_M(ITER,IPARAM,TIMEDIM,SPACEDIM,OBSDIM,DM,A,REALS,L)
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ITER,IPARAM,TIMEDIM,SPACEDIM,OBSDIM
 REAL(KIND=DP_KIND),INTENT(OUT) :: L
 REAL(KIND=DP_KIND),INTENT(IN) :: A
 REAL(KIND=DP_KIND),DIMENSION(TIMEDIM,SPACEDIM),INTENT(IN) :: DM
 REAL(KIND=DP_KIND),DIMENSION(TIMEDIM,SPACEDIM),INTENT(INOUT) :: REALS
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: X1,X2,X3,X4,X5,X6,X7,PM1,PM2,DD
 REAL(KIND=DP_KIND) :: RN,XSUM,VAR
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: W
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: V,VV,U,TMP1
 INTEGER :: I,J,MAXDIM
 
 MAXDIM=MAX(TIMEDIM,OBSDIM)
 ALLOCATE(X1(TIMEDIM,OBSDIM),DD(TIMEDIM,OBSDIM),X2(OBSDIM,OBSDIM)); X1=0.0D0; DD=0.0D0; X2=0.0D0

 RN=SQRT(DBLE(TIMEDIM)-1.0D0)

 !## compute ensemble mean
 DO I=1,OBSDIM
  !## compute average observational error
  XSUM=0.0D0; DO J=1,TIMEDIM; XSUM=XSUM+MSR%DHG(J,I); ENDDO; XSUM=XSUM/DBLE(TIMEDIM)
  DO J=1,TIMEDIM; X1(J,I)=(MSR%DHG(J,I)-XSUM)/RN; ENDDO
 ENDDO

 X2=0.0D0; DO I=1,OBSDIM; VAR=1.0D0/MSR%W(I); X2(I,I)=SQRT(1.0D0/VAR); ENDDO
 !## dd is [timedim,obsdim]
 CALL IPEST_IES_DUMP(X1,'X1',ITER) !EQUAL
 CALL IPEST_IES_DUMP(X2,'X2',ITER) !EQUAL
 CALL UTL_MATMUL(X2,X1,DD); DEALLOCATE(X1)
 CALL IPEST_IES_DUMP(DD,'X3',ITER) !EQUAL
 
 ALLOCATE(W(MAXDIM),U(OBSDIM,OBSDIM),VV(TIMEDIM,OBSDIM),TMP1(OBSDIM,OBSDIM))
 W=0.0D0; U=0.0D0; VV=0.0D0; TMP1=0.0D0
 !## compute svd
 TMP1=DD; CALL LUDCMP_SVD_MAIN(OBSDIM,TIMEDIM,MAXDIM,TMP1,U,W,VV)
 CALL IPEST_IES_DUMP(U,'X4',ITER)  !MINUS DIFF
 CALL IPEST_IES_DUMP(VV,'X5',ITER) !MINUS DIFF

 !## truncate for non-zeros
 DD_TRUNCATED=0; DO I=1,MAXDIM; IF(W(I).LE.0.0001D0)EXIT; DD_TRUNCATED=DD_TRUNCATED+1; ENDDO
 
 !## store actual vector v
 ALLOCATE(V(DD_TRUNCATED,TIMEDIM))
 DO I=1,DD_TRUNCATED; DO J=1,TIMEDIM
  V(I,J)=VV(J,I)
 ENDDO; ENDDO
 CALL IPEST_IES_DUMP(V,'X6',ITER) !SMALL DIFF

 !## compute x1 from chen et al.
 ALLOCATE(X1(OBSDIM,DD_TRUNCATED),X3(OBSDIM,DD_TRUNCATED)); X1=0.0D0
 DO I=1,OBSDIM; DO J=1,DD_TRUNCATED; X1(I,J)=U(J,I); ENDDO; ENDDO 
 CALL IPEST_IES_DUMP(X1,'X7',ITER) !minus
 CALL IPEST_IES_DUMP(X2,'X8',ITER) !SMALL DIFF
 CALL UTL_MATMUL(X1,X2,X3)
 DEALLOCATE(X1); ALLOCATE(X1(TIMEDIM,DD_TRUNCATED))

 !## msr%dhg() can have too many rows, clip it to obsdim
 ALLOCATE(X4(TIMEDIM,OBSDIM))
 DO I=1,TIMEDIM; DO J=1,OBSDIM; X4(I,J)=MSR%DHG(I,J); ENDDO; ENDDO
 CALL IPEST_IES_DUMP(X3,'X9',ITER)  !SMALL DIFF
 CALL IPEST_IES_DUMP(X4,'X10',ITER) !EQUAL
 CALL UTL_MATMUL(X3,X4,X1)
 CALL IPEST_IES_DUMP(X1,'X12',ITER) !DIFFER
 DEALLOCATE(X4)
! STOP
 
 !## compute x2 from chen et al.
 DEALLOCATE(X3); ALLOCATE(X3(DD_TRUNCATED,DD_TRUNCATED)); X3=0.0D0
 DO I=1,DD_TRUNCATED; X3(I,I)=1.0D0/(A+W(I)**2.0D0); ENDDO
 DEALLOCATE(X2); ALLOCATE(X2(TIMEDIM,DD_TRUNCATED))
 CALL IPEST_IES_DUMP(X3,'X11',ITER) !EQUAL
 CALL IPEST_IES_DUMP(X1,'X12',ITER) !DIFFER
 CALL UTL_MATMUL(X3,X1,X2)
 
 !## compute x3 from chen et al.
 DEALLOCATE(X3); ALLOCATE(X3(DD_TRUNCATED,DD_TRUNCATED)); X3=0.0D0
 DO I=1,DD_TRUNCATED; X3(I,I)=W(I); ENDDO
 ALLOCATE(X4(DD_TRUNCATED,TIMEDIM)); X4=0.0D0
 DO I=1,DD_TRUNCATED; DO J=1,TIMEDIM; X4(I,J)=V(I,J); ENDDO; ENDDO
 ALLOCATE(X5(DD_TRUNCATED,TIMEDIM))
 CALL IPEST_IES_DUMP(X4,'X13',ITER)
 CALL IPEST_IES_DUMP(X3,'X14',ITER)
 CALL UTL_MATMUL(X4,X3,X5)
 DEALLOCATE(X3,X4)
 ALLOCATE(X3(TIMEDIM,TIMEDIM))
 CALL IPEST_IES_DUMP(X5,'X15',ITER)
 CALL IPEST_IES_DUMP(X2,'X16',ITER)
 CALL UTL_MATMUL(X5,X2,X3)
 DEALLOCATE(X1,X2,X5)

 ALLOCATE(X1(TIMEDIM,SPACEDIM))
 CALL UTL_MATMUL(-PEST%PARAM(IPARAM)%INVCOV,DM,X1)
 CALL IPEST_IES_DUMP(PEST%PARAM(IPARAM)%INVCOV,'X17',ITER)
 CALL IPEST_IES_DUMP(DM,'X18',ITER)
 ALLOCATE(PM1(TIMEDIM,SPACEDIM))
 !## compute pm1
 CALL IPEST_IES_DUMP(X1,'X19',ITER) 
 CALL IPEST_IES_DUMP(X3,'X20',ITER) 
 CALL UTL_MATMUL(X1,X3,PM1); DEALLOCATE(X1,X3)

 IF(ITER.GT.0)THEN
 
  !## calculate x4 from paper of chen et al.
  ALLOCATE(X1(SPACEDIM,AM_TRUNCATED)); X1=0.0D0
  DO I=1,SPACEDIM; DO J=1,AM_TRUNCATED; X1(I,J)=PEST%PARAM(IPARAM)%AM(J,I); ENDDO; ENDDO
!  DO I=1,SPACEDIM; DO J=1,AM_TRUNCATED; X1(I,J)=AM(J,I); ENDDO; ENDDO
  ALLOCATE(X2(SPACEDIM,AM_TRUNCATED))
  CALL UTL_MATMUL(X1,PEST%PARAM(IPARAM)%ISQRTCOV,X2)
  DEALLOCATE(X1)
  ALLOCATE(X4(TIMEDIM,AM_TRUNCATED),X1(TIMEDIM,SPACEDIM))
  DO I=1,TIMEDIM; DO J=1,SPACEDIM; X1(I,J)=REALS(I,J)-PEST%PARAM(IPARAM)%MPR(I,J); ENDDO; ENDDO
!  DO I=1,TIMEDIM; DO J=1,SPACEDIM; X1(I,J)=REALS(I,J)-MPR(I,J); ENDDO; ENDDO
  CALL UTL_MATMUL(X2,X1,X4)
  DEALLOCATE(X1,X2)

  !## calculate x5 from paper of chen et al.
  ALLOCATE(X5(TIMEDIM,SPACEDIM))
  CALL UTL_MATMUL(PEST%PARAM(IPARAM)%AM,X4,X5)
!  CALL UTL_MATMUL(AM,X4,X5)
  DEALLOCATE(X4)

  !## calculate x6 from paper of chen et al.
  ALLOCATE(X4(SPACEDIM,TIMEDIM))
  DO I=1,SPACEDIM; DO J=1,TIMEDIM; X4(I,J)=DM(J,I); ENDDO; ENDDO
  ALLOCATE(X6(TIMEDIM,TIMEDIM))
  CALL UTL_MATMUL(X4,X5,X6)
  DEALLOCATE(X4,X5)

  !## calculate x7 from paper of chen et al.
  !## part 1
  ALLOCATE(X1(DD_TRUNCATED,DD_TRUNCATED)); X1=0.0D0 
  DO I=1,DD_TRUNCATED; IF(W(I).GT.0.0D0)X1(I,I)=1.0D0/(A+W(I)**2.0D0); ENDDO
  ALLOCATE(X2(DD_TRUNCATED,TIMEDIM),X3(DD_TRUNCATED,TIMEDIM)); X2=0.0D0; X3=0.0D0
  DO I=1,DD_TRUNCATED; DO J=1,TIMEDIM; X3(I,J)=V(I,J); ENDDO; ENDDO
  CALL UTL_MATMUL(X3,X1,X2)
  DEALLOCATE(X1,X3)

  ALLOCATE(X1(TIMEDIM,DD_TRUNCATED),X3(TIMEDIM,TIMEDIM))
  DO I=1,TIMEDIM; DO J=1,DD_TRUNCATED; X1(I,J)=V(J,I); ENDDO; ENDDO
  CALL UTL_MATMUL(X2,X1,X3)
  DEALLOCATE(X2,X1)

  ALLOCATE(X7(TIMEDIM,TIMEDIM))
  CALL UTL_MATMUL(X3,X6,X7)
  DEALLOCATE(X3,X6)

  ALLOCATE(X1(TIMEDIM,SPACEDIM))
  CALL UTL_MATMUL(-PEST%PARAM(IPARAM)%INVCOV,DM,X1)
  ALLOCATE(PM2(TIMEDIM,SPACEDIM))
  CALL UTL_MATMUL(X1,X7,PM2) 
  DEALLOCATE(X1,X7)

 ENDIF
 
 ALLOCATE(X1(TIMEDIM,SPACEDIM))
 
 CALL IPEST_IES_DUMP(PM1,'PM1',ITER)    
 IF(ITER.GT.0)CALL IPEST_IES_DUMP(PM2,'PM2',ITER)
 
 X1=REALS+PM1; IF(ITER.GT.0)X1=X1+PM2
 DEALLOCATE(PM1); IF(ALLOCATED(PM2))DEALLOCATE(PM2)

 !## check update
 L=0.0D0; DO I=1,TIMEDIM; DO J=1,SPACEDIM; L=L+(X1(I,J)-REALS(I,J))**2.0D0; ENDDO; ENDDO; L=SQRT(L)
 
 !## copy update to reals
 REALS=X1
 
 DEALLOCATE(X1,TMP1,U,V,W)
 
 END SUBROUTINE IPEST_IES_COMPUTE_PARTIAL_M
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_DUMP(X,TXT,ITER) 
 !#####=================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(:,:) :: X
 CHARACTER(LEN=*),INTENT(IN) :: TXT
 INTEGER,INTENT(IN) :: ITER
 INTEGER :: I,J,IU
 
 RETURN
 
 IU=UTL_GETUNIT()
#if (defined(DEBUG))
 OPEN(IU,FILE='D:\IMOD-MODELS\IES\MODELS\TMP\DEBUG_'//  TRIM(TXT)//'_ITER'//TRIM(ITOS(ITER))//'.M',STATUS='UNKNOWN',ACTION='WRITE')
#elif (defined(RELEASE))
 OPEN(IU,FILE='D:\IMOD-MODELS\IES\MODELS\TMP\RELEASE_'//TRIM(TXT)//'_ITER'//TRIM(ITOS(ITER))//'.M',STATUS='UNKNOWN',ACTION='WRITE')
#endif
! WRITE(IU,*) 'COL,ROW,VALUE'
 WRITE(IU,'(A)') TRIM(TXT)//'=['
 DO J=1,SIZE(X,2)
  IF(J.LT.SIZE(X,2))THEN
   WRITE(IU,'(99(F15.7,A1))') (X(I,J),',',I=1,SIZE(X,1)-1),X(SIZE(X,1),J),';'
  ELSE
   WRITE(IU,'(99(F15.7,A1))') (X(I,J),',',I=1,SIZE(X,1)-1),X(SIZE(X,1),J) 
  ENDIF
 ENDDO
 WRITE(IU,'(A)') ']'
 CLOSE(IU)
 
 END SUBROUTINE IPEST_IES_DUMP
 
 !#####================================================================= 
 SUBROUTINE IPEST_IES_COMPUTE_AM(IPARAM,TIMEDIM,SPACEDIM,DM) 
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: TIMEDIM,SPACEDIM,IPARAM
 REAL(KIND=DP_KIND),DIMENSION(TIMEDIM,SPACEDIM),INTENT(IN) :: DM
 INTEGER :: I,J,NT,NS
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: W
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: A,V,U
 LOGICAL :: LTEST
 
 LTEST=.FALSE.

 IF(LTEST)THEN
  NT=2; NS=4
 ELSE
  NT=TIMEDIM; NS=SPACEDIM
 ENDIF
 
 ALLOCATE(W(MAX(NT,NS)),A(NT,NS),U(NS,NS),V(NT,NS))
 W=0.0D0; A=0.0D0; U=0.0D0; V=0.0D0
 
 IF(.NOT.LTEST)THEN
  A=DM
 ELSE
  A=0.0D0
!  A(1,1)=2.0D0
!  A(2,1)=4.0D0
!  A(1,2)=1.0D0
!  A(2,2)=3.0D0
 ENDIF

! !## see whether matlab can solve it
! IU=UTL_GETUNIT(); OPEN(IU,FILE='D:\TMP.M',STATUS='OLD',ACTION='WRITE')
! DO I=1,NS
!  WRITE(IU,'(9999F10.2)') (A(J,I),J=1,NT)
! ENDDO
! CLOSE(IU)
 
 !## svd computation
 CALL LUDCMP_SVD_MAIN(NS,NT,MAX(NS,NT),A,U,W,V)
 
 !## truncate all zero and get inverse of singular values
 IF(.NOT.LTEST)THEN
  !## truncate for non-zeros
  AM_TRUNCATED=0; DO I=1,SPACEDIM; IF(W(I).LE.0.0001D0)EXIT; AM_TRUNCATED=AM_TRUNCATED+1; ENDDO
  
!  IF(ALLOCATED(AM))DEALLOCATE(AM); ALLOCATE(AM(AM_TRUNCATED,SPACEDIM)); AM=0.0D0

  ALLOCATE(PEST%PARAM(IPARAM)%AM(AM_TRUNCATED,SPACEDIM)); PEST%PARAM(IPARAM)%AM=0.0D0

  DO I=1,SPACEDIM; DO J=1,AM_TRUNCATED
   IF(W(J).GT.0.0D0)PEST%PARAM(IPARAM)%AM(J,I)=U(J,I)*1.0D0/W(J)
!   IF(W(J).GT.0.0D0)AM(J,I)=U(J,I)*1.0D0/W(J)
  ENDDO; ENDDO

 ENDIF

 CALL IPEST_IES_DUMP(PEST%PARAM(IPARAM)%AM,'AM_X1',0) 
 CALL IPEST_IES_DUMP(U,'AM_X2',0) 
 CALL IPEST_IES_DUMP(DM,'AM_X3',0) 
 
 IF(ALLOCATED(U))DEALLOCATE(U); IF(ALLOCATED(W))DEALLOCATE(W)
 IF(ALLOCATED(V))DEALLOCATE(V); IF(ALLOCATED(A))DEALLOCATE(A)

 END SUBROUTINE IPEST_IES_COMPUTE_AM

 !#####=================================================================
 SUBROUTINE IPEST_IES_COMPUTE_GETREALS(DIR,IPARAM,TIMEDIM,SPACEDIM,REALS,ITER,ILOG,JLAMBDA) 
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,INTENT(IN) :: IPARAM,TIMEDIM,SPACEDIM,ITER,ILOG,JLAMBDA
 REAL(KIND=DP_KIND),DIMENSION(TIMEDIM,SPACEDIM) :: REALS
 INTEGER :: J,K,IROW,ICOL
 CHARACTER(LEN=256) :: FNAME
 
 !## compute delta-m
 REALS=0.0D0
 DO J=1,TIMEDIM
  CALL IPEST_IES_GETENSEMBLENAME(ITER,DIR,PEST%PARAM(IPARAM)%REALSFNAME(J),FNAME,JLAMBDA)
  IF(.NOT.IDFREADSCALE(FNAME,PRJIDF,10,0,1.0D0,0))RETURN
  K=0; DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL
   K=K+1
   REALS(J,K)=PRJIDF%X(ICOL,IROW)
   IF(ILOG)REALS(J,K)=LOG10(REALS(J,K))
  ENDDO; ENDDO
  CALL IDFDEALLOCATEX(PRJIDF)
 ENDDO
 
 END SUBROUTINE IPEST_IES_COMPUTE_GETREALS
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_COMPUTE_DM(DIR,IPARAM,TIMEDIM,SPACEDIM,REALS,DM,ITER,ILOG,JLAMBDA) 
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,INTENT(INOUT) :: ITER
 INTEGER,INTENT(IN) :: IPARAM,TIMEDIM,SPACEDIM,ILOG,JLAMBDA
 REAL(KIND=DP_KIND),DIMENSION(TIMEDIM,SPACEDIM) :: DM,REALS
 INTEGER :: J,K
 REAL(KIND=DP_KIND) :: RN,XSUM
 
 !## read realisatons
 CALL IPEST_IES_COMPUTE_GETREALS(DIR,IPARAM,TIMEDIM,SPACEDIM,REALS,ITER,ILOG,JLAMBDA) 
 
 !## compute delta-m
 RN=SQRT(DBLE(TIMEDIM)-1.0D0)

 DO K=1,SPACEDIM
  XSUM=0.0D0; DO J=1,TIMEDIM; XSUM=XSUM+REALS(J,K); ENDDO; XSUM=XSUM/DBLE(TIMEDIM)
  DO J=1,TIMEDIM
   REALS(J,K)=(REALS(J,K)-XSUM)/RN
  ENDDO
 ENDDO
 
 CALL UTL_MATMUL(PEST%PARAM(IPARAM)%ISQRTCOV,REALS,DM)

 !## read realisatons again
 CALL IPEST_IES_COMPUTE_GETREALS(DIR,IPARAM,TIMEDIM,SPACEDIM,REALS,ITER,ILOG,JLAMBDA) 

 END SUBROUTINE IPEST_IES_COMPUTE_DM
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_SAVEREALS(DIR,ITER,ILAMBDA,MNAME) 
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(INOUT) :: ITER
 INTEGER,INTENT(IN) :: ILAMBDA
 CHARACTER(LEN=*),INTENT(IN) :: DIR,MNAME
 INTEGER :: I,J,K,ILAY,SCL_UP,SCL_D,IU,IOS
 CHARACTER(LEN=256) :: FNAME
 CHARACTER(LEN=256),ALLOCATABLE,DIMENSION(:) :: LINE
 
 !## read idf
 DO I=1,PEST%NREALS
 
  DO J=1,SIZE(PEST%PARAM)
   DO K=1,SIZE(PEST%PARAM(J)%ILS)
    ILAY=PEST%PARAM(J)%ILS(K)
    IF(ILAY.LE.0.OR.ILAY.GT.PRJNLAY)CYCLE
    
    CALL IPEST_IES_GETENSEMBLENAME(ITER,DIR,PEST%PARAM(J)%REALSFNAME(I),FNAME,ILAMBDA)
    
    SELECT CASE (PEST%PARAM(J)%PPARAM)
     CASE ('KH')
 
      !## it is not allowed to scale here as the covariance matrix does not scale with it
      SCL_UP=10; SCL_D=0      
      
      !## read/clip/scale idf file
      IF(.NOT.IDFREADSCALE(FNAME,PRJIDF,SCL_UP,SCL_D,1.0D0,0))RETURN
      !## save array, do not correct for boundary condition as we not yet know for what layer the zone will apply
      IF(.NOT.PMANAGER_SAVEMF2005_MOD_U2DREL(TRIM(DIR)//'\MODELINPUT\LPF7\HK_L'//TRIM(ITOS(ILAY))//'_R#'//TRIM(ITOS(I))//'.ARR',PRJIDF,0,0,1,0))RETURN

    END SELECT
   ENDDO
  ENDDO
 
  !## adjust the MET-file, make sure the results are written in the right folder
  DO J=1,2
   FNAME=TRIM(DIR)//'\MODELINPUT\'//TRIM(MNAME)//'_R#'//TRIM(ITOS(I))//'.MET7'
   IU=UTL_GETUNIT(); OPEN(IU,FILE=FNAME,STATUS='OLD',ACTION='READ'); K=0
   DO
    READ(IU,'(A256)',IOSTAT=IOS) FNAME; IF(IOS.NE.0)EXIT; FNAME=UTL_CAP(FNAME,'U')
    K=K+1; IF(J.EQ.2)THEN
     IF(INDEX(FNAME,'RESULTDIR').GT.0)THEN
      LINE(K)='RESULTDIR "'//TRIM(DIR)//'IES_R#'//TRIM(ITOS(I))//'_L#'//TRIM(ITOS(ILAMBDA))//'"'
     ELSE
      LINE(K)=FNAME
     ENDIF
    ENDIF
   ENDDO
   CLOSE(IU)
  ENDDO
  FNAME=TRIM(DIR)//'\MODELINPUT\'//TRIM(MNAME)//'_R#'//TRIM(ITOS(I))//'.MET7'
  IU=UTL_GETUNIT(); OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',ACTION='WRITE'); K=0
  DO J=1,SIZE(LINE); WRITE(IU,'(A)') TRIM(LINE(J)); ENDDO; CLOSE(IU)
! # MET7 File Generated by Beta Build Intel v2019.1.144 [13-03-2020 14:57]-iMOD [V5_1 X64 Debug] !!! Expiring date: 30/4/2020 !!!
!COORD_XLL 183950.000
!COORD_YLL 377150.000
!COORD_XLL_NB 183950.000
!COORD_YLL_NB 377150.000
!COORD_XUR_NB 184850.000
!COORD_YUR_NB 377850.000
!RESULTDIR "D:\IMOD-MODELS\IES\MODELS\RESULTS\MODEL_V2_LAMBDAS\IIES_R#1"
!SAVEDOUBLE 0
 ENDDO
 
 END SUBROUTINE IPEST_IES_SAVEREALS
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_GETENSEMBLENAME(ITER,DIR,FNAME_IN,FNAME_OUT,JLAMBDA)
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ITER,JLAMBDA
 CHARACTER(LEN=*),INTENT(IN) :: FNAME_IN,DIR
 CHARACTER(LEN=*),INTENT(OUT) :: FNAME_OUT

 IF(ITER.EQ.0)THEN
  FNAME_OUT=FNAME_IN
 ELSE
  FNAME_OUT=FNAME_IN(INDEX(FNAME_IN,'\',.TRUE.)+1:)
!  FNAME_OUT=FNAME_OUT(:INDEX(FNAME_OUT,'.',.TRUE.)-1)//'_ITER'//TRIM(ITOS(ITER))//'.IDF'
  FNAME_OUT=FNAME_OUT(:INDEX(FNAME_OUT,'.',.TRUE.)-1)//'_ITER'//TRIM(ITOS(ITER))// &
         '_ILAMBDA'//TRIM(ITOS(JLAMBDA))//'.IDF'
  FNAME_OUT=TRIM(DIR)//'\ENSEMBLES\'//TRIM(FNAME_OUT)
 ENDIF

 END SUBROUTINE IPEST_IES_GETENSEMBLENAME
 
 !#####=================================================================
 LOGICAL FUNCTION IPEST_IES_ALLOCATE(DIR)
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER :: I,N
 
 IPEST_IES_ALLOCATE=.FALSE.
 
 PEST%PE_REGULARISATION=0
 
 CALL IPEST_IES_DEALLOCATE()
 N=PEST%NREALS; ALLOCATE(RNG(N),IPROC(2,N),GPARAM(N),ISTATUS(N),STIME(N),ETIME(N))

 ALLOCATE(JE(N,0:PEST%PE_MXITER)); JE=0.0D0
 
 !## set realisation runbatch-files
 DO I=1,PEST%NREALS; RNG(I)=TRIM(DIR)//'\RUN_R#'//TRIM(ITOS(I))//'.BAT'; GPARAM(I)=I; ENDDO

 !## allocate memory
 IF(.NOT.IPEST_GLM_ALLOCATEMSR())RETURN
 
 !## open output files
 CALL UTL_CREATEDIR(TRIM(DIR)//'\IIES')
 IUPESTOUT=UTL_GETUNIT();        OPEN(IUPESTOUT,       FILE=TRIM(DIR)//'\IIES\LOG_PEST.TXT'           ,STATUS='UNKNOWN',ACTION='WRITE')
 IUPESTPROGRESS=UTL_GETUNIT();   OPEN(IUPESTPROGRESS,  FILE=TRIM(DIR)//'\IIES\LOG_PEST_PROGRESS.TXT'  ,STATUS='UNKNOWN',ACTION='WRITE')
 IUPESTEFFICIENCY=UTL_GETUNIT(); OPEN(IUPESTEFFICIENCY,FILE=TRIM(DIR)//'\IIES\LOG_PEST_EFFICIENCY.TXT',STATUS='UNKNOWN',ACTION='WRITE')
 WRITE(IUPESTEFFICIENCY,'(/60X,A30)') '--------IMPROVEMENT (%)-------'
 WRITE(IUPESTEFFICIENCY,'(A10,6(1X,A15))') 'ITER','CURRENT_TJ','PREVIOUS_PJ','AVG-RESIDUAL','RELATIVE','ABSOLUTE','LAMBDA'

 IPEST_IES_ALLOCATE=.TRUE.
 
 END FUNCTION IPEST_IES_ALLOCATE

 !#####=================================================================
 SUBROUTINE IPEST_IES_DEALLOCATE()
 !#####=================================================================
 IMPLICIT NONE

 IF(ALLOCATED(RNG))    DEALLOCATE(RNG)
 IF(ALLOCATED(IPROC))  DEALLOCATE(IPROC)
 IF(ALLOCATED(GPARAM)) DEALLOCATE(GPARAM)
 IF(ALLOCATED(ISTATUS))DEALLOCATE(ISTATUS)
 IF(ALLOCATED(DM))     DEALLOCATE(DM)
 IF(ALLOCATED(JE))     DEALLOCATE(JE)
 
 IF(IUPESTOUT.GT.0)        CLOSE(IUPESTOUT);         IUPESTOUT=0
 IF(IUPESTPROGRESS.GT.0)   CLOSE(IUPESTPROGRESS);    IUPESTPROGRESS=0
 IF(IUPESTEFFICIENCY.GT.0) CLOSE(IUPESTEFFICIENCY);  IUPESTEFFICIENCY=0

 END SUBROUTINE IPEST_IES_DEALLOCATE
 
END MODULE MOD_IPEST_IES