!!  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,ILINE,KLINE,JLINE,SPACEDIM,TIMEDIM,OBSDIM
 REAL(KIND=DP_KIND) :: MINJ,L,L0,L1,M0,M1,MINL,F
 REAL(KIND=DP_KIND),DIMENSION(:),ALLOCATABLE :: MO
  
 !## 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')

 !## set dimensions
 SPACEDIM=PRJIDF%NROW*PRJIDF%NCOL
 TIMEDIM =PEST%NREALS
 ALLOCATE(MO(TIMEDIM)); MO=0.0D0

 CALL WMESSAGEENABLE(TIMEREXPIRED,1); MSR%PJ=HUGE(1.0D0); ITER=0
 DO WHILE (ITER.LE.PEST%PE_MXITER)
  
  !## next iteration
  WRITE(IUPESTOUT,'(/A/)') ' *** Next Optimization Cycle '//TRIM(ITOS(ITER+1))//' ***'
  WRITE(*,'(/A/)')         ' *** Next Optimization Cycle '//TRIM(ITOS(ITER+1))//' ***'

  !## no residual output
  IUPESTRESIDUAL=0
  !## initial run - no lambda testing
  IF(ITER.EQ.0)THEN
   WRITE(IUPESTPROGRESS,'(/A)') 'Initial run:'
   !## 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; JLINE=1; CALL IPEST_IES_SAVEREALS(DIR,ITER,JLAMBDA,JLINE,MNAME)
   !## run and process all realizations-testing
   IF(.NOT.IPEST_GLM_RUNMODELS(IBATCH,RNG,GPARAM,'R',DIR,MNAME,ITER,LAMBDA,JLAMBDA))EXIT
   !## add model-error (=0.0 here) to it
   MO=0.0D0; JE(:,ITER)=JE(:,ITER)+MO; CALL IPEST_IES_PROGRESS_OUT(ITER)
   WRITE(IUPESTEFFICIENCY,'(I5,A5,1X,F15.3,6(16X),1X,F15.3)') ITER,'-',SUM(JE(:,ITER)),LAMBDA
   !## dump dhg as backstepping might be possible
   CALL IPEST_IES_WRITE_DHG(DIR,ITER,JLAMBDA,JLINE,MO)    
  !## carry out the lambda tests
  ELSE
   !## simulate ensembles of different lambdas
   MINJ=10.0D20; JLAMBDA=0
   DO ILAMBDA=1,PBMAN%NLAMBDASEARCH
    MINL=10.0D20; KLINE=0
    DO ILINE=1,PBMAN%NLINESEARCH
     F=PBMAN%LINE_SEARCH(ILINE)
     !## read from idf file and save the realisations in lpf package - apply the correct IDF files and export them for the realizations
     CALL IPEST_IES_SAVEREALS(DIR,ITER,ILAMBDA,ILINE,MNAME)
     L=LAMBDA*PBMAN%LAMBDA_TEST(ILAMBDA)
     !## run and process all realizations-testing
     WRITE(IUPESTPROGRESS,'(/A)') 'Lambda Testing '//TRIM(ITOS(ILAMBDA))//' ('//TRIM(RTOS(L,'G',7))//'); Line-Search '//TRIM(ITOS(ILINE))//' (Factor '//TRIM(RTOS(F,'G',7))//')'
     IF(.NOT.IPEST_GLM_RUNMODELS(IBATCH,RNG,GPARAM,'R',DIR,MNAME,ITER,L,ILAMBDA))EXIT
     !## compute model error
     IF(PEST%PE_REGULARISATION.EQ.1)CALL IPEST_IES_COMPUTE_MODEL_ERROR(DIR,MO,SPACEDIM,TIMEDIM,OBSDIM,ILOG,ITER,ILAMBDA,ILINE)
     !## add model-error to it
     JE(:,ITER)=JE(:,ITER)+MO; CALL IPEST_IES_PROGRESS_OUT(ITER)
     WRITE(IUPESTEFFICIENCY,'(2I5,1X,F15.3,6(16X),2(1X,F15.3))') ITER,ILINE,SUM(JE(:,ITER)),L,F
     !## dump dhg as backstepping might be possible
     CALL IPEST_IES_WRITE_DHG(DIR,ITER,ILAMBDA,ILINE,MO)
     !## evaluate which lambda is best
     IF(SUM(JE(:,ITER)).LT.MINL)THEN
      MINL=SUM(JE(:,ITER)); KLINE=ILINE
     ENDIF
    ENDDO
    !## check whether this line-search is minimal compared to previous values, if so store it
    IF(MINL.LT.MINJ)THEN
     MINJ=MINL; JLAMBDA=ILAMBDA; JLINE=KLINE
    ENDIF
   ENDDO
   IF(JLAMBDA.EQ.0)THEN; STOP 'SOMETHING WENT WRONG: JLAMBDA=0'; ENDIF
   !## echo selected lambda
   L=LAMBDA*PBMAN%LAMBDA_TEST(JLAMBDA)
   F=PBMAN%LINE_SEARCH(JLINE)
   WRITE(IUPESTPROGRESS,'(/A/)') 'Selected lambda = '//TRIM(ITOS(JLAMBDA))//' ('//TRIM(RTOS(L,'F',3))// &
                                ') Selected linesearch = '//TRIM(ITOS(JLINE))//' ('//TRIM(RTOS(F,'F',3))// &
           ') has minimal Objective Function Value ('//TRIM(RTOS(MINJ,'F',3))//')'
   WRITE(IUPESTEFFICIENCY,'(/A/)') 'Selected lambda = '//TRIM(ITOS(JLAMBDA))//' ('//TRIM(RTOS(L,'F',3))// &
                                ') Selected linesearch = '//TRIM(ITOS(JLINE))//' ('//TRIM(RTOS(F,'F',3))// &
           ') has minimal Objective Function Value ('//TRIM(RTOS(MINJ,'F',3))//')'
   !## read residuals of all ensembles for selected jlambda and linesearch
   CALL IPEST_IES_READ_DHG(DIR,ITER,JLAMBDA,JLINE,MO)    
   
  ENDIF

  !## get update of realisations
  OBSDIM=MSR%NOBS
  IF(.NOT.IPEST_IES_UPDATE_ENSEMBLES(DIR,ITER,JLAMBDA,JLINE,IBATCH,MNAME,SPACEDIM,TIMEDIM,OBSDIM))EXIT  
  !## save mean/stdev heads
  CALL IPEST_IES_SAVE_HEAD_STATS(DIR,ITER,JLAMBDA)
  !## save mean/stdev parameters
  CALL IPEST_IES_SAVE_PARAMETERS_STATS(DIR,ITER,JLAMBDA,JLINE)
  
  !## 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
 WRITE(IUPESTOUT,'(/A10,999(1X,A15))') 'ITER','TOTAL',('VARIANCE_SIM'//TRIM(ITOS(IGRAD)),IGRAD=1,PEST%NREALS)
 DO ITER=0,PEST%PE_MXITER
  WRITE(IUPESTOUT,'(I10,999(1X,F15.7))') ITER,SUM(SIGMA(1:PEST%NREALS,ITER))/REAL(PEST%NREALS),(SIGMA(IGRAD,ITER),IGRAD=1,PEST%NREALS)
 ENDDO
 WRITE(IUPESTOUT,'(/A10,999(1X,A15))') 'ITER','TOTAL',('IMPROVEM._SIM'//TRIM(ITOS(IGRAD)),IGRAD=1,PEST%NREALS)
 DO ITER=0,PEST%PE_MXITER
  L0=SUM(SIGMA(1:PEST%NREALS,0))/REAL(PEST%NREALS)
  M0=SUM(MU(1:PEST%NREALS,0))/REAL(PEST%NREALS)
  L1=SUM(SIGMA(1:PEST%NREALS,ITER))/REAL(PEST%NREALS)
  M1=SUM(MU(1:PEST%NREALS,ITER))/REAL(PEST%NREALS)
  L =100.0D0*(1.0D0-(M0+L0)/(M1+L1))
  WRITE(IUPESTOUT,'(I10,999(1X,F15.7))') ITER,L, &
      (100.0D0*(1.0D0-((MU(IGRAD,ITER)+SIGMA(IGRAD,ITER))/(MU(IGRAD,0)+SIGMA(IGRAD,0)))),IGRAD=1,PEST%NREALS)
 ENDDO
 WRITE(IUPESTOUT,'(/A10,999(1X,A15))') 'ITER','TOTAL',('AVERAGE_SIM'//TRIM(ITOS(IGRAD)),IGRAD=1,PEST%NREALS)
 DO ITER=0,PEST%PE_MXITER
  WRITE(IUPESTOUT,'(I10,999(1X,F15.7))') ITER,SUM(MU(1:PEST%NREALS,ITER))/REAL(PEST%NREALS),(MU(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)
 !#####=================================================================
 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
  ENDIF
 ENDDO
 CALL IDFDEALLOCATEX(STDEV); CALL IDFDEALLOCATEX(MEAN)
 
 IPEST_IES_CREATEENSEMBLES=.TRUE.
 
 END FUNCTION IPEST_IES_CREATEENSEMBLES

 !#####=================================================================
 SUBROUTINE IPEST_IES_SAVE_HEAD_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_HEAD_STATS
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_SAVE_PARAMETERS_STATS(DIR,ITER,JLAMBDA,JLINE)
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,INTENT(IN) :: ITER,JLAMBDA,JLINE
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: TMP1
 TYPE(IDFOBJ),ALLOCATABLE,DIMENSION(:) :: PAR
 TYPE(IDFOBJ) :: MPR,VPR
 REAL(KIND=DP_KIND) :: AVG,STD
 INTEGER :: I,J,IROW,ICOL,N
 
 !## compute mean and variance heads
 ALLOCATE(PAR(PEST%NREALS)); DO I=1,PEST%NREALS; CALL IDFNULLIFY(PAR(I)); ENDDO; CALL IDFNULLIFY(MPR); CALL IDFNULLIFY(VPR)
 DO I=1,SIZE(PEST%PARAM)
  DO J=1,PEST%NREALS
   IF(ITER.EQ.0)THEN
    PAR(J)%FNAME=TRIM(DIR)//'\ENSEMBLES\PARAMETER'//TRIM(ITOS(I))//'_R'//TRIM(ITOS(J))//'.IDF'
   ELSE
    PAR(J)%FNAME=TRIM(DIR)//'\ENSEMBLES\PARAMETER'//TRIM(ITOS(I))//'_R'//TRIM(ITOS(J))//'_ITER'//TRIM(ITOS(ITER))// &
                 '_ILAMBDA'//TRIM(ITOS(JLAMBDA))//'_ILINE'//TRIM(ITOS(JLINE))//'.IDF'
   ENDIF
   IF(.NOT.IDFREAD(PAR(J),PAR(J)%FNAME,1))STOP
  ENDDO
  IF(I.EQ.1)THEN; CALL IDFCOPY(PAR(1),MPR); CALL IDFCOPY(PAR(1),VPR); ENDIF
  ALLOCATE(TMP1(PEST%NREALS)); DO IROW=1,MPR%NROW; DO ICOL=1,MPR%NCOL
   TMP1(1)=0.0D0; DO J=1,PEST%NREALS; TMP1(J)=PAR(J)%X(ICOL,IROW); ENDDO
   CALL UTL_STDEF(TMP1,PEST%NREALS,PAR(1)%NODATA,STD,AVG,N)
   MPR%X(ICOL,IROW)=AVG; VPR%X(ICOL,IROW)=STD
  ENDDO; ENDDO
  DEALLOCATE(TMP1)
  CALL UTL_CREATEDIR(TRIM(DIR)//'\STATPARS')
  MPR%FNAME=TRIM(DIR)//'\STATPARS\MEAN_PARAMETER'//TRIM(ITOS(I))//'_ITER_'//TRIM(ITOS(ITER))//'.IDF'
  IF(.NOT.IDFWRITE(MPR,MPR%FNAME,1))STOP
  VPR%FNAME=TRIM(DIR)//'\STATPARS\STDEV_PARAMETER'//TRIM(ITOS(I))//'_ITER_'//TRIM(ITOS(ITER))//'.IDF'
  IF(.NOT.IDFWRITE(VPR,VPR%FNAME,1))STOP
 ENDDO
 
 END SUBROUTINE IPEST_IES_SAVE_PARAMETERS_STATS
 
 !#####=================================================================
 LOGICAL FUNCTION IPEST_IES_UPDATE_ENSEMBLES(DIR,ITER,JLAMBDA,JLINE,IBATCH,MNAME,SPACEDIM,TIMEDIM,OBSDIM) 
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR,MNAME
 INTEGER,INTENT(IN) :: IBATCH,SPACEDIM,TIMEDIM,OBSDIM
 INTEGER,INTENT(INOUT) :: ITER
 INTEGER,INTENT(INOUT) :: JLAMBDA,JLINE
 INTEGER :: I,II,JJ,J,K,N,IROW,ICOL,IX,ILOG,ILAMBDA,IPARAM,NPOP,ILINE
 REAL(KIND=DP_KIND) :: X,VAR,A,L,TL,MAXTL,GAMMA,F1,F2,F3,J0,MLT
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: DM,X1,X2,X3
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: SO,MO,O,TMP1,PCLASS,RCLASS
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: REALS
 INTEGER,ALLOCATABLE,DIMENSION(:,:) :: PDF 
 CHARACTER(LEN=256) :: FNAME
 LOGICAL :: LMODELERROR
 INTEGER,SAVE :: KLAMBDA
 
 WRITE(*,'(/A/)') 'Update ensembles ...'
 
 LMODELERROR=.FALSE.; IF(PEST%PE_REGULARISATION.EQ.1)LMODELERROR=.TRUE.
 
 !## apply log transformation
 ILOG=1
 GAMMA=4.0D0
 MAXTL=100.0D0
 IPEST_IES_UPDATE_ENSEMBLES=.TRUE.

 !## 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) 
   SO(J)=SO(J)+MSR%DHG(J,K)**2.0D0*MSR%W(K)
   O(J)=O(J)+MSR%DHG(J,K)
  ENDDO
  CALL UTL_STDEF(TMP1,OBSDIM,99999.0D0,SIGMA(J,ITER),MU(J,ITER),NPOP)
  !## get pdf of error distribution
  CALL IPEST_IES_GETPDF(TMP1,PDF(:,J),RCLASS) 
 ENDDO

 !## compute model-error per ensemble
 ALLOCATE(MO(TIMEDIM)); MO=0.0D0
 IF(LMODELERROR.AND.ITER.GT.0)CALL IPEST_IES_COMPUTE_MODEL_ERROR(DIR,MO,SPACEDIM,TIMEDIM,OBSDIM,ILOG,ITER,JLAMBDA,JLINE)

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

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

 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,'(I5,A5,99(1X,F15.3))') ITER,'-',MSR%TJ,F3,SUM(SO),SUM(MO),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
    !## stop optimization, no progression anymore
    IPEST_IES_UPDATE_ENSEMBLES=.FALSE.; RETURN
   ELSE

    !## get update for lambda
    CALL IPEST_IES_UPDATE_LAMBDA(DIR,ITER,GAMMA,MO)

    !## store lambda that "won-the-prize" for back-stepping possibility within the next iteration
    KLAMBDA=JLAMBDA
    WRITE(IUPESTOUT,'(/A/)') 'Awesome, significant reduction of the objective function reduction: proceed'
    !## save previous objective function value
    MSR%PJ=MSR%TJ

   ENDIF
  ELSE
  
   WRITE(IUPESTEFFICIENCY,'(/A,F15.3)') 'OLD LAMBDA',LAMBDA

   LAMBDA=LAMBDA*PBMAN%LAMBDA_TEST(JLAMBDA)
   N=PBMAN%NLAMBDASEARCH/2; IF(MOD(PBMAN%NLAMBDASEARCH,2).NE.0)N=N+1
   IF(JLAMBDA.LT.N)THEN
    LAMBDA=1.0D0/PBMAN%LAMBDA_TEST(PBMAN%NLAMBDASEARCH)*LAMBDA
   ELSEIF(JLAMBDA.GT.N)THEN
    LAMBDA=1.0D0/PBMAN%LAMBDA_TEST(1)*LAMBDA
   ENDIF

!   IF(PBMAN%LAMBDA_TEST(PBMAN%NLAMBDASEARCH).NE.1.0D0)THEN
!    LAMBDA=LAMBDA*PBMAN%LAMBDA_TEST(PBMAN%NLAMBDASEARCH)
!   ELSE
!    LAMBDA=LAMBDA*GAMMA
!   ENDIF
   !## reset ensembles to previous iteration
   ITER=ITER-1
   !## store lambda that "won-the-prize" for back-stepping possibility within the next iteration
   JLAMBDA=MAX(1,KLAMBDA)
   WRITE(IUPESTOUT,'(/A/)') 'No objective function reduction: reset to previous iteration'
   
   WRITE(IUPESTEFFICIENCY,'(/A,2I10)') 'Read ensembles for lambda/linesearch',JLAMBDA,JLINE
   WRITE(IUPESTEFFICIENCY,'(/A,F15.3)') 'NEW LAMBDA',LAMBDA
   
   !## read residuals of all ensembles for selected jlambda
   CALL IPEST_IES_READ_DHG(DIR,ITER,JLAMBDA,JLINE,MO)
   WRITE(IUPESTEFFICIENCY,'(/A,F15.3)') 'SUM J',JE(:,ITER)
   
  ENDIF
 ELSE
  MSR%PJ=MSR%TJ
 ENDIF
 
 !## try to find acceptable ensemble updates for all lambda to be tested
 WRITE(IUPESTOUT,'(/A)') 'Process Lambdas ('//TRIM(ITOS(PBMAN%NLAMBDASEARCH))//')'
 DO ILAMBDA=1,PBMAN%NLAMBDASEARCH
  MLT=PBMAN%LAMBDA_TEST(ILAMBDA)
  DO
   WRITE(IUPESTOUT,'(/A)') 'Process Lambda '//TRIM(ITOS(ILAMBDA))//' ('//TRIM(RTOS(MLT*LAMBDA,'G',7))//')' 
 
   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,JLINE) 

     !## 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
    
    ENDIF
    
    !## save the updates for the next iteration
    ALLOCATE(TMP1(SPACEDIM)) 
    
    DO ILINE=1,PBMAN%NLINESEARCH

     !## read realisatons again from previous cycle
     WRITE(IUPESTOUT,'(/A)') 'Get parameters from previous cycle'
     CALL IPEST_IES_COMPUTE_GETREALS(DIR,I,TIMEDIM,SPACEDIM,REALS,ITER,ILOG,JLAMBDA,JLINE) 

     !## processing the inverse of the prior parameter covariance matrix
     CALL IPEST_IES_COMPUTE_PARTIAL_M(ITER,I,TIMEDIM,SPACEDIM,OBSDIM,DM,A,REALS,L,LMODELERROR,ILINE); TL=TL+L
     WRITE(IUPESTOUT,'(/A,F15.7/)') 'LENGTH OF INDIVIDUAL ENSEMBLE UPDATE=',L

     IF(.NOT.IDFALLOCATEX(PRJIDF))STOP
 
     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,ILINE,'Write')
      !## 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
    
     CALL IDFDEALLOCATEX(PRJIDF)  

    ENDDO
    
    DEALLOCATE(TMP1)

   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,MO,DM,REALS)
 
 FLUSH(IUPESTOUT); FLUSH(IUPESTPROGRESS); FLUSH(IUPESTEFFICIENCY)
 
 END FUNCTION IPEST_IES_UPDATE_ENSEMBLES
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_COMPUTE_MODEL_ERROR(DIR,MO,SPACEDIM,TIMEDIM,OBSDIM,ILOG,ITER,JLAMBDA,JLINE)
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 REAL(KIND=DP_KIND),DIMENSION(:),INTENT(OUT) :: MO
 INTEGER,INTENT(IN) :: SPACEDIM,TIMEDIM,OBSDIM,ILOG,ITER,JLAMBDA,JLINE
 INTEGER :: IPARAM,I,J
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: X1,X2,X3
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: REALS
 
 ALLOCATE(X1(SPACEDIM,TIMEDIM),X2(SPACEDIM,SPACEDIM),X3(SPACEDIM,TIMEDIM))
 IF(ALLOCATED(REALS))DEALLOCATE(REALS); ALLOCATE(REALS(TIMEDIM,SPACEDIM))
 
 DO IPARAM=1,SIZE(PEST%PARAM)
  !## read realisatons for model-error
  WRITE(IUPESTOUT,'(/A)') 'Read previous set of parameter to compute the MODEL-ERROR'
  CALL IPEST_IES_COMPUTE_GETREALS(DIR,IPARAM,TIMEDIM,SPACEDIM,REALS,ITER,ILOG,JLAMBDA,JLINE) 
  X2=PEST%PARAM(IPARAM)%ICOV
  DO I=1,TIMEDIM
   DO J=1,SPACEDIM; X1(J,I)=REALS(I,J)-PEST%PARAM(IPARAM)%MPR(I,J); ENDDO
  ENDDO
  CALL UTL_MATMUL(X1,X2,X3); DEALLOCATE(X1); ALLOCATE(X1(TIMEDIM,SPACEDIM))
  X1=TRANSPOSE(X3)
  DEALLOCATE(X2); ALLOCATE(X2(TIMEDIM,TIMEDIM))
  CALL UTL_MATMUL(X3,X1,X2)
  DO I=1,TIMEDIM; MO(I)=X2(I,I); ENDDO    
  DEALLOCATE(X1,X2,X3)
 ENDDO

 MO=MO/(REAL(SPACEDIM,8)/REAL(OBSDIM,8))

 END SUBROUTINE IPEST_IES_COMPUTE_MODEL_ERROR
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_UPDATE_LAMBDA(DIR,ITER,GAMMA,MO)
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 REAL(KIND=DP_KIND),DIMENSION(:),INTENT(OUT) :: MO
 INTEGER,INTENT(IN) :: ITER
 REAL(KIND=DP_KIND),INTENT(IN) :: GAMMA
 INTEGER :: I,J
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: XLAMBDA,YLAMBDA,CPARABOL
 REAL(KIND=DP_KIND) :: X,MINJ,JJ

 WRITE(IUPESTEFFICIENCY,'(/A,F15.3)') 'OLD LAMBDA ',LAMBDA

 IF(PBMAN%NLAMBDASEARCH.GT.1)THEN 

  I=PBMAN%NLAMBDASEARCH
  ALLOCATE(XLAMBDA(I),YLAMBDA(I),CPARABOL(3))
  DO I=1,PBMAN%NLAMBDASEARCH
   MINJ=1.0D21; DO J=1,PBMAN%NLINESEARCH
    CALL IPEST_IES_READ_DHG(DIR,ITER,I,J,MO)
    JJ=SUM(JE(:,ITER))+SUM(MO)
    IF(JJ.LT.MINJ)MINJ=JJ
   ENDDO
   YLAMBDA(I)=MINJ 
   XLAMBDA(I)=PBMAN%LAMBDA_TEST(I)*LAMBDA
  ENDDO
  
  WRITE(IUPESTEFFICIENCY,'(/2A10,1X,2A15)') 'ITER','ILAMBDA','LAMBDA','J'
  DO I=1,PBMAN%NLAMBDASEARCH
   WRITE(IUPESTEFFICIENCY,'(2I10,1X,2F15.3)') ITER,I,YLAMBDA(I),XLAMBDA(I)
  ENDDO
 
  CALL UTL_FIT_PARABOLA(CPARABOL,XLAMBDA,YLAMBDA)
  !## valley at this lambda value
  X=-CPARABOL(2)/(2.0D0*CPARABOL(3))

  !## bot-parabol
  IF(CPARABOL(3).GT.0.0D0)THEN
   LAMBDA=MAX(LAMBDA/GAMMA,X)
   WRITE(IUPESTEFFICIENCY,'(/A,F15.3)') 'BOT LAMBDA ', X 
  ELSE
   LAMBDA=LAMBDA/GAMMA
   WRITE(IUPESTEFFICIENCY,'(/A,F15.3)') 'TOP LAMBDA ', X 
  ENDIF
  WRITE(IUPESTEFFICIENCY,'(A,F15.3/)') 'NEW LAMBDA ',LAMBDA
  DEALLOCATE(XLAMBDA,YLAMBDA,CPARABOL)
 ELSE
  X=LAMBDA; LAMBDA=LAMBDA/GAMMA
  WRITE(IUPESTEFFICIENCY,'(/A,F15.3/)') 'NEW LAMBDA ',LAMBDA
 ENDIF

 END SUBROUTINE IPEST_IES_UPDATE_LAMBDA

 !#####=================================================================
 SUBROUTINE IPEST_IES_WRITE_DHG(DIR,ITER,ILAMBDA,ILINE,MO)
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 REAL(KIND=DP_KIND),DIMENSION(:),INTENT(IN) :: MO
 INTEGER,INTENT(IN) :: ITER,ILAMBDA,ILINE
 INTEGER :: IU,I,J
 REAL(KIND=DP_KIND) :: W
 
 IU=UTL_GETUNIT(); OPEN(IU,FILE=TRIM(DIR)//'\IIES\DHG_ITER'//TRIM(ITOS(ITER))//'_ILAMBDA'//TRIM(ITOS(ILAMBDA))// &
              '_ILINE'//TRIM(ITOS(ILINE))//'.TXT',STATUS='UNKNOWN',ACTION='WRITE')
 WRITE(IU,'(/A/)') 'TOTAL-ERROR'
 WRITE(IU,'(2A10)') 'REAL.','ERROR'
 DO I=1,SIZE(MSR%DHG,1); WRITE(IU,'(I10,F15.3)') I,JE(I,ITER); ENDDO
 WRITE(IU,'(20A1)') ('-',I=1,20)
 WRITE(IU,'(A10,F15.3)') 'TOTAL',SUM(JE(:,ITER))
 WRITE(IU,'(/A/)') 'DATA-ERROR'
 WRITE(IU,'(2A10,99A15)') 'MEASURE','STDEV',('DELTA_H_'//TRIM(ITOS(I)),I=1,SIZE(MSR%DHG,1))
 DO J=1,MSR%NOBS 
  W=SQRT(1.0D0/MSR%W(J)); WRITE(IU,'(I10,99F15.7)') J,W,(MSR%DHG(I,J),I=1,SIZE(MSR%DHG,1))
 ENDDO
 WRITE(IU,'(/A/)') 'MODEL-ERROR'
 WRITE(IU,'(2A10)') 'REAL.','ERROR'
 DO I=1,SIZE(MO)
  WRITE(IU,'(I10,F15.3)') I,MO(I)
 ENDDO
 CLOSE(IU)
 
 END SUBROUTINE IPEST_IES_WRITE_DHG
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_READ_DHG(DIR,ITER,ILAMBDA,ILINE,MO)
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,INTENT(IN) :: ITER,ILAMBDA,ILINE
 REAL(KIND=DP_KIND),DIMENSION(:),INTENT(OUT) :: MO
 INTEGER :: IU,I,J
 REAL(KIND=DP_KIND) :: W
 
 IU=UTL_GETUNIT()
 IF(ITER.EQ.0)THEN
  OPEN(IU,FILE=TRIM(DIR)//'\IIES\DHG_ITER'//TRIM(ITOS(ITER))//'_ILAMBDA'// &
       TRIM(ITOS(1))//'_ILINE'//TRIM(ITOS(1))//'.TXT',STATUS='OLD',ACTION='READ')
 ELSE
  OPEN(IU,FILE=TRIM(DIR)//'\IIES\DHG_ITER'//TRIM(ITOS(ITER))//'_ILAMBDA'// &
       TRIM(ITOS(ILAMBDA))//'_ILINE'//TRIM(ITOS(ILINE))//'.TXT',STATUS='OLD',ACTION='READ')
 ENDIF
 DO I=1,4; READ(IU,*); ENDDO
 DO I=1,SIZE(MSR%DHG,1)
  READ(IU,'(10X,F15.3)') JE(I,ITER)
 ENDDO
 DO I=1,6; READ(IU,*); ENDDO
 DO J=1,MSR%NOBS
  READ(IU,'(10X,99F15.7)') W,(MSR%DHG(I,J),I=1,SIZE(MSR%DHG,1))
  MSR%W(J)=1.0D0/W**2.0D0
!  DO I=1,SIZE(MSR%DHG,1)
!   JE(I,ITER)=JE(I,ITER)+MSR%W(J)*MSR%DHG(I,J)**2.0D0  
!  ENDDO
 ENDDO
 DO I=1,4; READ(IU,*); ENDDO
 DO I=1,SIZE(MO)
  READ(IU,'(10X,F15.3)') MO(I)
 ENDDO
 CLOSE(IU)
 
 END SUBROUTINE IPEST_IES_READ_DHG

 !#####=================================================================
 SUBROUTINE IPEST_IES_PROGRESS_OUT(ITER)
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ITER
 REAL(KIND=DP_KIND) :: F1,F2
 INTEGER :: I
 
 WRITE(IUPESTPROGRESS,'(125A1)') ('-',I=1,125)
 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,L3
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: XT1,XT2,XT3,XT4
 
 !## carried out only once
 IF(ASSOCIATED(PEST%PARAM(IPARAM)%ICOV))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) - use double precision
 CALL IDFNULLIFY(IDF); IDF%ITYPE=8
 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)%SQRTCOV(N,N),PEST%PARAM(IPARAM)%ICOV(N,N))
 IF(.NOT.IDFALLOCATEX(IDF))RETURN
 
 INQUIRE(FILE=TRIM(DIR)//'\ENSEMBLES\ICOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF'    ,EXIST=L1)
 INQUIRE(FILE=TRIM(DIR)//'\ENSEMBLES\SQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF' ,EXIST=L2)
 INQUIRE(FILE=TRIM(DIR)//'\ENSEMBLES\ISQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF',EXIST=L3)
 IF(L1.AND.L2.AND.L3)THEN
  WRITE(*,'(/A)') 'Found existing:'
  WRITE(*,'(A)') '- ICOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF'
  WRITE(*,'(A)') '- SQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF'
  WRITE(*,'(A)') '- ISQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF'
  WRITE(*,'(/A$)') 'Do you want to use them (Y/y) ? '
  READ(*,'(A1)') YN
  IF(YN.EQ.'Y'.OR.YN.EQ.'y')THEN
   WRITE(*,'(A)') 'Reading '//TRIM(DIR)//'\ENSEMBLES\ICOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF ...'
   IF(.NOT.IDFREAD(IDF,TRIM(DIR)//'\ENSEMBLES\ICOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF',1))STOP
   PEST%PARAM(IPARAM)%ICOV=IDF%X
   WRITE(*,'(A)') 'Reading '//TRIM(DIR)//'\ENSEMBLES\SQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF ...'
   IF(.NOT.IDFREAD(IDF,TRIM(DIR)//'\ENSEMBLES\SQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF',1))STOP
   PEST%PARAM(IPARAM)%SQRTCOV=IDF%X
   WRITE(*,'(A)') 'Reading '//TRIM(DIR)//'\ENSEMBLES\ISQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF ...'
   IF(.NOT.IDFREAD(IDF,TRIM(DIR)//'\ENSEMBLES\ISQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF',1))STOP
   PEST%PARAM(IPARAM)%ISQRTCOV=IDF%X
  ELSE
   L1=.FALSE.; L2=.FALSE.; L3=.FALSE.
  ENDIF
 ENDIF
 IF(.NOT.L1.OR..NOT.L2)THEN
  ALLOCATE(XT1(N,N),XT2(N,N),XT3(N,N),XT4(N,N))
  XT1=PEST%PARAM(IPARAM)%COV
  XT2=0.0D0; XT3=0.0D0; XT4=0.0D0
  call IPEST_IES_DUMP(Xt1,'cov',0)
  CALL LUDCMP_CALC_SQRTROOTINVERSE(N,XT1,XT2,XT3,XT4)
  !## use variances only (not the covariates)
!  XT2=0.0D0; DO J=1,N; XT2(J,J)=1.0D0/PEST%PARAM(IPARAM)%COV(J,J); ENDDO
  PEST%PARAM(IPARAM)%ICOV    =XT2
  PEST%PARAM(IPARAM)%SQRTCOV =XT3
  PEST%PARAM(IPARAM)%ISQRTCOV=XT4
!## checkt as it need to be identity matrix
!  xt3=PEST%PARAM(IPARAM)%COV
!  CALL UTL_MATMUL(xt3,xt2,xt4)
  DEALLOCATE(XT1,XT2,XT3,XT4)
  DO J=1,N; DO K=1,N; IDF%X(J,K)=PEST%PARAM(IPARAM)%ICOV(J,K); ENDDO; ENDDO
  IDF%FNAME=TRIM(DIR)//'\ENSEMBLES\ICOV_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)%SQRTCOV(J,K); ENDDO; ENDDO
  IDF%FNAME=TRIM(DIR)//'\ENSEMBLES\SQRTCOV_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)%ISQRTCOV(J,K); ENDDO; ENDDO
  IDF%FNAME=TRIM(DIR)//'\ENSEMBLES\ISQRTCOV_IPARAM'//TRIM(ITOS(IPARAM))//'.IDF'; IF(.NOT.IDFWRITE(IDF,IDF%FNAME,1))STOP
 ENDIF
 CALL IDFDEALLOCATEX(IDF)
 !## save rest in single precision
 IDF%ITYPE=4
 DEALLOCATE(PEST%PARAM(IPARAM)%COV)
 
 END SUBROUTINE IPEST_IES_INITIATE_INVERSE_MATRICES
  
 !#####=================================================================
 SUBROUTINE IPEST_IES_GETPDF(X,PDF,XCLASS) 
 !#####=================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),DIMENSION(:) :: XCLASS
 REAL(KIND=DP_KIND),DIMENSION(:) :: X
 INTEGER,DIMENSION(:) :: PDF
 REAL(KIND=DP_KIND) :: LX 
 INTEGER :: I,J 
 
 PDF=0; DO I=1,SIZE(X) 
  DO J=1,SIZE(XCLASS)-1
   LX=X(I)
   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,LMODELERROR,ILINE)
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ITER,IPARAM,TIMEDIM,SPACEDIM,OBSDIM,ILINE
 LOGICAL,INTENT(IN) :: LMODELERROR
 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,F
 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))
 ALLOCATE(X2(SPACEDIM,SPACEDIM)); X2=-PEST%PARAM(IPARAM)%SQRTCOV
 CALL UTL_MATMUL(X2,DM,X1)
 DEALLOCATE(X2)
 CALL IPEST_IES_DUMP(PEST%PARAM(IPARAM)%SQRTCOV,'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
  ALLOCATE(X2(SPACEDIM,AM_TRUNCATED))
  ALLOCATE(X3(SPACEDIM,SPACEDIM))
  X3=PEST%PARAM(IPARAM)%ISQRTCOV
  CALL UTL_MATMUL(X1,X3,X2)
  DEALLOCATE(X1,X3)
  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
  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)
  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))
  ALLOCATE(X2(SPACEDIM,SPACEDIM)); X2=-PEST%PARAM(IPARAM)%SQRTCOV
  CALL UTL_MATMUL(X2,DM,X1)
  ALLOCATE(PM2(TIMEDIM,SPACEDIM))
  CALL UTL_MATMUL(X1,X7,PM2) 
  DEALLOCATE(X1,X2,X7)

 ENDIF
 
 ALLOCATE(X1(TIMEDIM,SPACEDIM))
 
 CALL IPEST_IES_DUMP(PM1,'PM1',ITER)    
 IF(ITER.GT.0.AND.LMODELERROR)CALL IPEST_IES_DUMP(PM2,'PM2',ITER)
 
 F=PBMAN%LINE_SEARCH(ILINE)
 X1=REALS+F*PM1; IF(ITER.GT.0)X1=X1+F*PM2
! 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

 CALL UTL_CREATEDIR('D:\IMOD-MODELS\IES\MODELS\TMP')
 
 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,ILAMBDA,ILINE) 
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,INTENT(IN) :: IPARAM,TIMEDIM,SPACEDIM,ITER,ILOG,ILAMBDA,ILINE
 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,ILAMBDA,ILINE,'Read')
  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,ILAMBDA,ILINE) 
 !#####=================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,INTENT(INOUT) :: ITER
 INTEGER,INTENT(IN) :: IPARAM,TIMEDIM,SPACEDIM,ILOG,ILAMBDA,ILINE
 REAL(KIND=DP_KIND),DIMENSION(TIMEDIM,SPACEDIM) :: DM,REALS
 INTEGER :: J,K
 REAL(KIND=DP_KIND) :: RN,XSUM
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: X1
 
 !## read realisatons
 WRITE(IUPESTOUT,'(/A)') 'Read previous set of parameter to compute DM'
 CALL IPEST_IES_COMPUTE_GETREALS(DIR,IPARAM,TIMEDIM,SPACEDIM,REALS,ITER,ILOG,ILAMBDA,ILINE) 
 
 !## 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
 ALLOCATE(X1(SPACEDIM,SPACEDIM)); X1=PEST%PARAM(IPARAM)%ISQRTCOV
 CALL UTL_MATMUL(X1,REALS,DM)
 DEALLOCATE(X1)
 
 !## read realisatons again
 WRITE(IUPESTOUT,'(/A)') 'Read previous set of parameter to refresh for AM'
 CALL IPEST_IES_COMPUTE_GETREALS(DIR,IPARAM,TIMEDIM,SPACEDIM,REALS,ITER,ILOG,ILAMBDA,ILINE) 

 END SUBROUTINE IPEST_IES_COMPUTE_DM
 
 !#####=================================================================
 SUBROUTINE IPEST_IES_SAVEREALS(DIR,ITER,ILAMBDA,ILINE,MNAME) 
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ITER
 INTEGER,INTENT(IN) :: ILAMBDA,ILINE
 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
 
 WRITE(IUPESTOUT,'(/A/)') 'Compute Realizations'

 !## 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,ILINE,'Read')
    
    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)//'\IIES_R#'//TRIM(ITOS(I))//'_L#'//TRIM(ITOS(ILAMBDA))//'"'
     ELSE
      LINE(K)=FNAME
     ENDIF
    ENDIF
   ENDDO
   CLOSE(IU)
   IF(J.EQ.1)ALLOCATE(LINE(K))
  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)
  DEALLOCATE(LINE)
! # 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,ILAMBDA,ILINE,TXT)
 !#####=================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ITER,ILAMBDA,ILINE
 CHARACTER(LEN=*),INTENT(IN) :: FNAME_IN,DIR,TXT
 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))// &
         '_ILAMBDA'//TRIM(ITOS(ILAMBDA))//'_ILINE'//TRIM(ITOS(ILINE))//'.IDF'
  FNAME_OUT=TRIM(DIR)//'\ENSEMBLES\'//TRIM(FNAME_OUT)
 ENDIF
 
 WRITE(IUPESTOUT,'(A)') TRIM(TXT)//': '//TRIM(FNAME_OUT)
 
 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.
 
 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
 ALLOCATE(SIGMA(N,0:PEST%PE_MXITER)); SIGMA=0.0D0
 ALLOCATE(MU   (N,0:PEST%PE_MXITER)); MU   =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,'(2A5,10(1X,A15))') 'ITER','LTEST','TOTAL-ERROR_NEW','TOTAL-ERROR_OLD','DATA-ERROR','MODEL-ERROR','AVG-RESIDUAL','RELATIVE','ABSOLUTE','LAMBDA','MEAN','STDEV'
 WRITE(IUPESTEFFICIENCY,'(10X,10(1X,A15))') ('[-]',I=1,4),'[M]','[%]','[%]','[-]','[-]','[-]'

 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(ALLOCATED(SIGMA))  DEALLOCATE(SIGMA)
 IF(ALLOCATED(MU))     DEALLOCATE(MU)
 
 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