!! Copyright (C) Stichting Deltares, 2005-2023. !! !! 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 . !! !! Contact: imod.support@deltares.nl !! Stichting Deltares !! P.O. Box 177 !! 2600 MH Delft, The Netherlands. !! MODULE MOD_IPEST_GLM #if(defined(DEFPARALLEL)) USE OMP_LIB #endif USE WINTERACTER USE MOD_IDF_PAR USE MOD_IDF, ONLY : IDFNULLIFY,IDFALLOCATEX,IDFWRITE,IDFREAD,IDFDEALLOCATEX,IDFWRITEFREE_ROW,IDFIROWICOL USE IMODVAR, ONLY : DP_KIND USE MOD_OSD, ONLY : OSD_OPEN USE MOD_PMANAGER_PAR, ONLY : PEST,SIM,PARAM,PRJNLAY,PBMAN,PRJNPER,BND,HNOFLOW,PCG,PRJIDF,PSTMEASURE,TOPICS,TKDW,TKHV,TKVV,TKVA,TSTO,TSPY USE MOD_PMANAGER_UTL, ONLY : PMANAGER_SAVEMF2005_MOD_U2DREL USE MOD_UTL, ONLY : UTL_GETUNIT,UTL_CAP,VTOS,UTL_GETMED,UTL_CREATEDIR,UTL_GOODNESS_OF_FIT,UTL_NASH_SUTCLIFFE,& VTOS,UTL_MF2005_MAXNO,UTL_DEL1TREE,YMDHMSTOITIME,UTL_COMPUTE_GXG,UTL_CONTINUE USE MOD_IPEST_GLM_PAR USE MOD_KRIGING_PAR USE MOD_KRIGING USE MOD_LUDCMP !INCLUDE 'for_iosdef.for' CONTAINS !#####================================================================= SUBROUTINE IPEST_GLM_MAIN(DIR,MNAME,IBATCH,ITER) !#####================================================================= IMPLICIT NONE REAL(KIND=DP_KIND) :: GAMMA !## lambda increase/decrease INTEGER,INTENT(IN) :: IBATCH INTEGER,INTENT(OUT) :: ITER CHARACTER(LEN=*),INTENT(IN) :: DIR,MNAME CHARACTER(LEN=256) :: LINE,STRING,DIRNAME,DIRO,DIRP REAL(KIND=DP_KIND) :: LAMBDA,DUMMY INTEGER :: I,J,N,IX,ILOG,NITER,IOS,IP,IACT LOGICAL :: LEX !## natural log transformation ILOG=1 !## initialize gamma GAMMA=4.0D0 !## organise groups IF(.NOT.IPEST_GLM_SETGROUPS(IBATCH))RETURN !## allocate memory for simulations CALL IPEST_GLM_ALLOCATE(DIR) !## set initial values for alpha() DO I=1,SIZE(PEST%PARAM); IF(.NOT.IPEST_GLM_CHK(I,IBATCH))RETURN; ENDDO !## remove all nam-files unneeded DO I=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(I)%PACT).EQ.1)CYCLE INQUIRE(FILE=TRIM(DIR)//'\'//TRIM(MNAME)//'_P#'//TRIM(VTOS(I))//'.NAM',EXIST=LEX) IF(LEX)CALL IOSDELETEFILE(TRIM(DIR)//'\'//TRIM(MNAME)//'_P#'//TRIM(VTOS(I))//'.NAM') ENDDO !## output map DIRP=TRIM(DIR)//'\IPEST'; IF(TRIM(PBMAN%IPESTPOUTPUT).NE.'')DIRP=PBMAN%IPESTPOUTPUT DIRO=DIR; IF(PBMAN%OUTPUT.NE.'')DIRO=PBMAN%OUTPUT CALL IOSDIRNAME(DIRNAME) !## make sure diro and pbman%output are all absolute pathnames STRING=DIRO; CALL UTL_RELPATHNAME(DIRNAME,STRING,DIRO) STRING=DIRP; CALL UTL_RELPATHNAME(DIRNAME,STRING,DIRP) STRING=PBMAN%OUTPUT; CALL UTL_RELPATHNAME(DIRNAME,STRING,PBMAN%OUTPUT) !## restart, if so read previous files IF(PBMAN%RESTART.GT.0)THEN INQUIRE(FILE=TRIM(DIRP)//'\LOG_PEST_EFFICIENCY.TXT',EXIST=LEX) !## read initial objective function value IF(LEX)THEN IUPESTEFFICIENCY=UTL_GETUNIT(); OPEN(IUPESTEFFICIENCY, FILE=TRIM(DIRP)//'\LOG_PEST_EFFICIENCY.TXT',STATUS='OLD',ACTION='READ') DO I=1,2; READ(IUPESTEFFICIENCY,*); ENDDO; READ(IUPESTEFFICIENCY,*) MSR%TJRESTART,DUMMY,MSR%RJRESTART NITER=0; DO READ(IUPESTEFFICIENCY,*,IOSTAT=IOS) IF(IOS.NE.0)EXIT; NITER=NITER+1 ENDDO LINE='COPY '//TRIM(DIRP)//'\LOG_PEST_EFFICIENCY.TXT '//TRIM(DIRP)//'\LOG_PEST_EFFICIENCY_COPY.TXT' CLOSE(IUPESTEFFICIENCY); CALL SYSTEM(LINE) ENDIF INQUIRE(FILE=TRIM(DIRP)//'\LOG_PEST_RUNFILE.TXT',EXIST=LEX) !## read initial objective function value IF(LEX)THEN IUPESTRUNFILE=UTL_GETUNIT(); OPEN(IUPESTRUNFILE, FILE=TRIM(DIRP)//'\LOG_PEST_RUNFILE.TXT' ,STATUS='OLD',ACTION='READ') DO READ(IUPESTRUNFILE,'(A256)',IOSTAT=IOS) LINE IF(IOS.NE.0)THEN WRITE(*,*) 'Reading previous parameters did not read in correctly' STOP ENDIF LINE=UTL_CAP(LINE,'U') WRITE(STRING,'(A,1X,I10)') 'Copy in the runfile, iteration',NITER; STRING=UTL_CAP(STRING,'U') IF(INDEX(LINE,STRING).GT.0)THEN WRITE(*,*) 'New values:' !## reset parameters DO IP=1,SIZE(PEST%PARAM) READ(IUPESTRUNFILE,*) IACT,PEST%PARAM(IP)%PPARAM,PEST%PARAM(IP)%PILS,PEST%PARAM(IP)%PIZONE,PEST%PARAM(IP)%PINI !## make upper case PEST%PARAM(IP)%PPARAM=UTL_CAP(PEST%PARAM(IP)%PPARAM,'U') WRITE(*,'(A,1X,2(I5,1X),F15.3)') TRIM(PEST%PARAM(IP)%PPARAM),PEST%PARAM(IP)%PILS,PEST%PARAM(IP)%PIZONE,PEST%PARAM(IP)%PINI IF(PEST%PARAM(IP)%PLOG.EQ.1)THEN PEST%PARAM(IP)%PINI=LOG(PEST%PARAM(IP)%PINI) ELSEIF(PEST%PARAM(IP)%PLOG.EQ.2)THEN PEST%PARAM(IP)%PINI=LOG10(PEST%PARAM(IP)%PINI) ENDIF PEST%PARAM(IP)%ALPHA(1)=PEST%PARAM(IP)%PINI !## current alpha PEST%PARAM(IP)%ALPHA(2)=PEST%PARAM(IP)%PINI !## previous alpha PEST%PARAM(IP)%LALPHA=PEST%PARAM(IP)%ALPHA(1) ENDDO EXIT ENDIF ENDDO LINE='COPY '//TRIM(DIRP)//'\LOG_PEST_RUNFILE.TXT '//TRIM(DIRP)//'\LOG_PEST_RUNFILE_COPY.TXT' CLOSE(IUPESTRUNFILE); CALL SYSTEM(LINE) ENDIF ENDIF !## open output files CALL UTL_CREATEDIR(TRIM(DIRP)) IF(PBMAN%NCYCLE.GT.0)THEN IUPESTOUT=UTL_GETUNIT(); OPEN(IUPESTOUT, FILE=TRIM(DIRP)//'\LOG_PEST_ICYCLE'//TRIM(VTOS(PBMAN%ICYCLE))//'.TXT' ,STATUS='REPLACE',ACTION='WRITE') IUPESTPROGRESS=UTL_GETUNIT(); OPEN(IUPESTPROGRESS, FILE=TRIM(DIRP)//'\LOG_PEST_PROGRESS_ICYCLE'//TRIM(VTOS(PBMAN%ICYCLE))//'.TXT' ,STATUS='REPLACE',ACTION='WRITE') IUPESTEFFICIENCY=UTL_GETUNIT(); OPEN(IUPESTEFFICIENCY, FILE=TRIM(DIRP)//'\LOG_PEST_EFFICIENCY_ICYCLE'//TRIM(VTOS(PBMAN%ICYCLE))//'.TXT' ,STATUS='REPLACE',ACTION='WRITE') IUPESTSENSITIVITY=UTL_GETUNIT(); OPEN(IUPESTSENSITIVITY,FILE=TRIM(DIRP)//'\LOG_PEST_SENSITIVITY_ICYCLE'//TRIM(VTOS(PBMAN%ICYCLE))//'.TXT',STATUS='REPLACE',ACTION='WRITE') IUPESTRUNFILE=UTL_GETUNIT(); OPEN(IUPESTRUNFILE, FILE=TRIM(DIRP)//'\LOG_PEST_RUNFILE_ICYCLE'//TRIM(VTOS(PBMAN%ICYCLE))//'.TXT' ,STATUS='REPLACE',ACTION='WRITE') IUPESTJACOBIAN=UTL_GETUNIT(); OPEN(IUPESTJACOBIAN, FILE=TRIM(DIRP)//'\LOG_PEST_JACOBIAN_ICYCLE'//TRIM(VTOS(PBMAN%ICYCLE))//'.TXT' ,STATUS='REPLACE',ACTION='WRITE') ELSE IUPESTOUT=UTL_GETUNIT(); OPEN(IUPESTOUT, FILE=TRIM(DIRP)//'\LOG_PEST.TXT' ,STATUS='REPLACE',ACTION='WRITE') IUPESTPROGRESS=UTL_GETUNIT(); OPEN(IUPESTPROGRESS, FILE=TRIM(DIRP)//'\LOG_PEST_PROGRESS.TXT' ,STATUS='REPLACE',ACTION='WRITE') IUPESTEFFICIENCY=UTL_GETUNIT(); OPEN(IUPESTEFFICIENCY, FILE=TRIM(DIRP)//'\LOG_PEST_EFFICIENCY.TXT' ,STATUS='REPLACE',ACTION='WRITE') IUPESTSENSITIVITY=UTL_GETUNIT(); OPEN(IUPESTSENSITIVITY,FILE=TRIM(DIRP)//'\LOG_PEST_SENSITIVITY.TXT',STATUS='REPLACE',ACTION='WRITE') IUPESTRUNFILE=UTL_GETUNIT(); OPEN(IUPESTRUNFILE, FILE=TRIM(DIRP)//'\LOG_PEST_RUNFILE.TXT' ,STATUS='REPLACE',ACTION='WRITE') IUPESTJACOBIAN=UTL_GETUNIT(); OPEN(IUPESTJACOBIAN, FILE=TRIM(DIRP)//'\LOG_PEST_JACOBIAN.TXT' ,STATUS='REPLACE',ACTION='WRITE') ENDIF WRITE(IUPESTOUT,'(A)') '===============' WRITE(IUPESTOUT,'(A)') 'iPESTP-Settings' WRITE(IUPESTOUT,'(A)') '===============' WRITE(IUPESTOUT,'(A)') 'Number of CPU used: '//TRIM(VTOS(PBMAN%NCPU)) WRITE(IUPESTOUT,'(A)') 'Number of Lambda-searches: '//TRIM(VTOS(PBMAN%NLAMBDASEARCH)) WRITE(IUPESTOUT,'(A,99F10.3)') '- Lambdas: ',(PBMAN%LAMBDA_TEST(I),I=1,PBMAN%NLAMBDASEARCH) WRITE(IUPESTOUT,'(A)') 'Number of Line-searches: '//TRIM(VTOS(PBMAN%NLINESEARCH)) WRITE(IUPESTOUT,'(A,99F10.3)') '- Line-searches: ',(PBMAN%LINE_SEARCH(I),I=1,PBMAN%NLINESEARCH) WRITE(IUPESTOUT,'(A)') 'Number of Pest Iterations: '//TRIM(VTOS(PEST%PE_MXITER)) IF(PEST%PE_MXITER.EQ.0)WRITE(IUPESTOUT,'(A)') ' Sensitivity Analyses Started' WRITE(IUPESTOUT,'(A)') 'Stop Criterium Objective Function Value: '//TRIM(VTOS(PEST%PE_STOP,'F',3)) WRITE(IUPESTOUT,'(A)') 'Sensitivity to Exclude Parameter (temporarily): '//TRIM(VTOS(ABS(PEST%PE_SENS),'F',3)) IF(PEST%PE_SENS.GE.0.0D0) WRITE(IUPESTOUT,'(A)') '- Parameters are turned off TEMPORARILY whenever they have less sensitivity' IF(PEST%PE_SENS.LT.0.0D0) WRITE(IUPESTOUT,'(A)') '- Parameters are turned off PERMENANTLY whenever they have less sensitivity' WRITE(IUPESTOUT,'(A)') 'Acceptable residual value: '//TRIM(VTOS(PEST%PE_DRES,'F',3)) WRITE(IUPESTOUT,'(A)') 'Absolute Residual Target: '//TRIM(VTOS(PEST%PE_TARGET(1),'F',3)) WRITE(IUPESTOUT,'(A)') 'Relative Residual Target: '//TRIM(VTOS(PEST%PE_TARGET(2),'F',3)) SELECT CASE (PEST%PE_SCALING) CASE (0); WRITE(IUPESTOUT,'(A)') 'No Scaling, No SVD' CASE (1); WRITE(IUPESTOUT,'(A)') 'Yes Scaling, No SVD' CASE (2); WRITE(IUPESTOUT,'(A)') 'Yes Scaling, Yes SVD' CASE (3); WRITE(IUPESTOUT,'(A)') 'No Scaling, Yes SVD' END SELECT SELECT CASE (ABS(PEST%PE_KTYPE)) CASE (1); WRITE(IUPESTOUT,'(A)') 'Linear Semivariogram is used (if neccessary)' CASE (2); WRITE(IUPESTOUT,'(A)') 'Spherical Semivariogram is used (if neccessary)' CASE (3); WRITE(IUPESTOUT,'(A)') 'Exponential Semivariogram is used (if neccessary)' CASE (4); WRITE(IUPESTOUT,'(A)') 'Gaussian Semivariogram is used (if neccessary)' CASE (5); WRITE(IUPESTOUT,'(A)') 'Power Semivariogram is used (if neccessary)' CASE (6); WRITE(IUPESTOUT,'(A)') 'Thiessen Interpolation is used (if neccessary)' END SELECT WRITE(IUPESTOUT,'(A)') 'Kriging Range: '//TRIM(VTOS(PEST%PE_KRANGE,'F',3)) WRITE(IUPESTOUT,'(A)') 'Termination Criterion for Parameter Adjustments (vectorlength): '//TRIM(VTOS(PEST%PE_PADJ,'F',3)) WRITE(IUPESTOUT,'(/A)') 'Parameters' WRITE(IUPESTOUT,'(A2,1X,A5,1X,A5,1X,A7,5(1X,A15),3A10,5A15)') 'AC','PTYPE','ILS','IZN','INITIAL','DELTA','MINIMUM','MAXIMUM','FADJ','IGROUP','LTRANS', & 'NODES','ACRONYM','PPRIOR','STDEV','SDATE','EDATE' DO I=1,SIZE(PEST%PARAM); CALL IPEST_GLM_ECHO_PARAMETERS(I); ENDDO; WRITE(IUPESTOUT,*) WRITE(IUPESTOUT,'(/1X,2(A,I10),A/)') 'Optimizing ',SIZE(RNG),' Parameters/Groups out of ',SIZE(PEST%PARAM),' parameters' WRITE(IUPESTEFFICIENCY,'(8(A15,1X))') 'TOTAL_J','MEAS._J','PARAM._J','RMSE_TJ','ADJUSTMENTS','CUR_IMPROVEMENT','TOT_IMPROVEMENT' WRITE(IUPESTEFFICIENCY,'(8(A15,1X))') '(L2)', '(L2)', '(L2)', '(L)', '(-)', '(%)', '(%)' IF(PBMAN%RESTART.EQ.1)THEN WRITE(IUPESTEFFICIENCY,'(3(F15.3,1X),A)') MSR%TJRESTART,MSR%TJRESTART-MSR%RJRESTART,MSR%RJRESTART,' >>> restart org value <<< ' ENDIF WRITE(IUPESTSENSITIVITY,'(A)') 'Sensitivity (%):' N=0; DO J=1,SIZE(PEST%PARAM); IF(ABS(PEST%PARAM(J)%PACT).EQ.1)N=N+1; ENDDO !.AND.PEST%PARAM(J)%PIGROUP.GT.0)N=N+1; ENDDO CALL IPEST_GLM_WRITEHEADER('Iteration ',N,IUPESTSENSITIVITY) !## get and remember actual iMOD run location + switch to temporal simulation directory CALL IOSDIRNAME(DIRNAME); CALL IOSDIRCHANGE(TRIM(DIR)//'\') CALL WMESSAGEENABLE(TIMEREXPIRED,1); MSR%PJ=HUGE(1.0D0) !## start optimization cycle ITER=0 MAINLOOP: DO !## run and process all lambda-testing !## save initial residuals IF(ITER.EQ.0)THEN IF(IUPESTRESIDUAL.GT.0)CLOSE(IUPESTRESIDUAL); IUPESTRESIDUAL=UTL_GETUNIT() IF(PBMAN%NCYCLE.GT.0)THEN OPEN(IUPESTRESIDUAL,FILE=TRIM(DIRP)//'\LOG_PEST_RESIDUAL_ICYCLE'//TRIM(VTOS(PBMAN%ICYCLE))//'.TXT',STATUS='REPLACE',ACTION='WRITE') ELSE OPEN(IUPESTRESIDUAL,FILE=TRIM(DIRP)//'\LOG_PEST_RESIDUAL_0.TXT',STATUS='REPLACE',ACTION='WRITE') ENDIF IF(IUPESTPRESIDUAL.GT.0)CLOSE(IUPESTPRESIDUAL); IUPESTPRESIDUAL=UTL_GETUNIT() IF(PBMAN%NCYCLE.GT.0)THEN OPEN(IUPESTPRESIDUAL,FILE=TRIM(DIRP)//'\LOG_PESTP_RESIDUAL_ICYCLE'//TRIM(VTOS(PBMAN%ICYCLE))//'.TXT',STATUS='REPLACE',ACTION='WRITE') ELSE OPEN(IUPESTPRESIDUAL,FILE=TRIM(DIRP)//'\LOG_PESTP_RESIDUAL_0.TXT',STATUS='REPLACE',ACTION='WRITE') ENDIF ELSE IUPESTRESIDUAL=0; IUPESTPRESIDUAL=0 ENDIF !## run lambda-testing IF(.NOT.IPEST_GLM_RUNMODELS(DIRNAME,IBATCH,RNL,LPARAM,'L',DIR,MNAME,ITER,LAMBDA,0))EXIT !## initial lambda IF(ITER.GT.0)THEN !## determine best lambda and continue, otherwise determine another set of lambda's IF(.NOT.IPEST_GLM_NEXT(IBATCH,ITER,DIR,DIRO,DIRP,LAMBDA,GAMMA,MNAME))THEN IF(GAMMA.GT.0.0D0)THEN !## determine new gradient I=IPEST_GLM_GRADIENT(IBATCH,ITER,LAMBDA,ABS(GAMMA)) ELSE !## overrule stopping IF(GAMMA.LT.0.0D0)I=-1 ENDIF SELECT CASE (I) !## quit CASE (-1); ITER=ITER-1; EXIT MAINLOOP !## start next cycle CASE (0 ); EXIT !## repeat current cycle CASE (1 ); CYCLE END SELECT ENDIF ELSE IF(MSR%NOBS.GT.0)THEN !## set current h on correct position DO I=1,MSR%NOBS MSR%HL(0,I) =MSR%HL(1,I) MSR%DHL(0,I)=MSR%DHL(1,I) ENDDO LAMBDA=MSR%TJ/DBLE(2.0D0*MSR%NOBS) IF(ILOG.EQ.1)THEN LAMBDA=LOG(LAMBDA); IX=FLOOR(LAMBDA); LAMBDA=EXP(REAL(IX,8)) ELSEIF(ILOG.EQ.2)THEN LAMBDA=LOG10(LAMBDA); IX=FLOOR(LAMBDA); LAMBDA=10.0D0**IX ENDIF !## write initial lambda WRITE(IUPESTOUT,'(/A/)') 'Initial Lambda_0 computed as '//TRIM(VTOS(LAMBDA,'G',3)) MSR%TJ_H(ITER)=MSR%TJ; MSR%RJ_H(ITER)=MSR%RJ; MSR%GOF_H(ITER)=MSR%GOF(1); MSR%NSC_H(ITER)=MSR%NSC(1) !## save initial heads with these (only for modflow2005/modflow6) and if isteady=1 icm transient IF(PBMAN%IFORMAT.EQ.2.OR.PBMAN%IFORMAT.EQ.3)CALL IPEST_GLM_UPDATEHEADS(DIR,DIRO,1) ENDIF ENDIF IF(.NOT.IPEST_GLM_ECHOPARAMETERS(IBATCH,ITER))EXIT !## if mxiter <0 stop IF(PEST%PE_MXITER.LT.0)EXIT !## next cycle ITER=ITER+1 !## copy current objective function value to previous objective function value IF(MSR%NOBS.GT.0)MSR%PJ=MSR%TJ !## "melt" all parameters for next cycle DO I=1,SIZE(PEST%PARAM); IF(PEST%PARAM(I)%PACT.EQ.-1)PEST%PARAM(I)%PACT=1; ENDDO !## run and process all parameters-sensitivities IUPESTRESIDUAL=0; IUPESTPRESIDUAL=0; IF(.NOT.IPEST_GLM_RUNMODELS(DIRNAME,IBATCH,RNG,GPARAM,'P',DIR,MNAME,ITER,LAMBDA,0))EXIT !## determine new gradient I=IPEST_GLM_GRADIENT(IBATCH,ITER,LAMBDA,GAMMA) !## sensitivities finished IF(PEST%PE_MXITER.EQ.0)EXIT SELECT CASE (I) !## quit CASE (-1); ITER=ITER-1; EXIT MAINLOOP END SELECT ENDDO MAINLOOP CALL WMESSAGETIMER(0); CALL WMESSAGEENABLE(TIMEREXPIRED,0) !## bring the iMOD run location back to he origional directory CALL IOSDIRCHANGE(DIRNAME) END SUBROUTINE IPEST_GLM_MAIN !#####================================================================= LOGICAL FUNCTION IPEST_GLM_RUNMODELS(DIRNAME,IBATCH,RN,RPARAM,RT,DIR,MNAME,ITER,LAMBDA,ILAMBDA) !#####================================================================= IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN),DIMENSION(:) :: RN REAL(KIND=DP_KIND),INTENT(IN) :: LAMBDA INTEGER,INTENT(IN),DIMENSION(:) :: RPARAM INTEGER,INTENT(IN) :: ITER,IBATCH,ILAMBDA CHARACTER(LEN=1),INTENT(IN) :: RT CHARACTER(LEN=*),INTENT(IN) :: DIR,MNAME,DIRNAME INTEGER :: IGRAD,NDONE,NCPU,IFLAGS,I,T,ET,ST,CLOCK_RATE INTEGER,DIMENSION(2) :: IDPROC LOGICAL :: LRESTART IPEST_GLM_RUNMODELS=.FALSE. IF(RT.EQ.'L')THEN WRITE(IUPESTPROGRESS,'(/A5,4A15,15X,4A15)') 'LS','PARAMETER','LAMBDA','TOT_J','RED_J','GOODNESS FIT','NASH SUTCLIFFE','SIM_TIME(SEC)','PP_TIME(SEC)' WRITE(IUPESTOUT,'(/A/)') ' *** Lambda Cycle ***'; WRITE(*,'(/A/)') ' *** Lambda Cycle ***' ELSEIF(RT.EQ.'P')THEN WRITE(IUPESTPROGRESS,'(/A5,5A15,30X,2A15)') 'GD','PARAMETER','FACTOR','TOT_J','RED_J','CUR_PARAMETER','SIM_TIME(SEC)','PP_TIME(SEC)' WRITE(IUPESTOUT,'(/A/)') ' *** Sensitivity Cycle ***'; WRITE(*,'(/A/)') ' *** Sensitivity Cycle ***' ELSEIF(RT.EQ.'R')THEN WRITE(IUPESTPROGRESS,'(/A5,4A15,15X,4A15)') 'RS','REALIZATION','LAMBDA','TOT_J','RED_J','GOODNESS FIT','NASH SUTCLIFFE','SIM_TIME(SEC)','PP_TIME(SEC)' WRITE(IUPESTOUT,'(/A/)') ' *** Realization Cycle '//TRIM(VTOS(ILAMBDA))//' ***' WRITE(*, '(/A/)') ' *** Realization Cycle '//TRIM(VTOS(ILAMBDA))//' ***' ENDIF !## see whether to reuse sensitivities LRESTART=.FALSE. IF(RT.EQ.'P'.AND.ITER.EQ.1)THEN !## reuse sensitivities if not equal to mod(iter,usesens) !LRESTART=.FALSE.; IF(MOD(ITER-1,PBMAN%USESENS).NE.0) IF(PBMAN%USESENS.EQ.1)LRESTART=.TRUE. IF(LRESTART)WRITE(*,'(/A/)') '>>> BE AWARE SENSITIVIES ARE RE-USED <<<' ENDIF CALL WMESSAGETIMER(NCSECS,IREPEAT=1); IGRAD=0; NCPU=0; IPROC=0; ISTATUS=-1; INSENS=0; IF(LRESTART)INSENS=-1 !## executes lamdba-testing on commandtool such that commands alike 'dir' etc. works IFLAGS=PROCCMDPROC; IF(PBMAN%CMDHIDE.EQ.1)IFLAGS=IFLAGS+PROCSILENT !## all done initially except the first IF(RT.EQ.'L'.AND.ITER.EQ.0)THEN; ISTATUS=1; ISTATUS(1)=-1; ENDIF !## exclude parameters that are turned off IF(RT.EQ.'P'.AND.PEST%PE_SENS.LT.0.0D0)THEN !## determine what model to be ran IGRAD=0; DO I=1,SIZE(PEST%PARAM) !## associated parameters to existing groups inactive for gradient computation IF(ABS(PEST%PARAM(I)%PACT).EQ.1)THEN !## flag as done IGRAD=IGRAD+1; IF(PEST%PARAM(I)%INSENS.EQ.1)INSENS(IGRAD)=1 ENDIF ENDDO ENDIF !## start the lambda-analyses DO !## start processes DO !## find gradient simulation to be (re)carried out DO IGRAD=1,SIZE(RN); IF(ISTATUS(IGRAD).EQ.-1)EXIT; ENDDO; IF(IGRAD.GT.SIZE(RN))EXIT !## number of cpu full NCPU=NCPU+1 !## parameter insensitive, skip simulating this parameter again IF(ABS(INSENS(IGRAD)).EQ.1)THEN IPROC(:,IGRAD)=1; LRESTART=.TRUE. ELSE SELECT CASE (RT) CASE ('L','P') !## adjust alpha for current igrad CALL IPEST_GLM_NEXTGRAD(RPARAM(IGRAD),RT,IGRAD) !## update arr-files for modflow6/seawat IF(PBMAN%IFORMAT.EQ.3.OR.PBMAN%IFORMAT.EQ.6)THEN !## update files for this simulation CALL IPEST_GLM_SAVE_PARAMETERS(IBATCH,DIR,RPARAM(IGRAD),RT,IGRAD) ELSE !## define update in pst file IF(.NOT.IPEST_GLM_PST(DIRNAME,DIR,MNAME,IGRAD,RPARAM(IGRAD),RT))THEN CALL IPEST_GLM_ERROR(IBATCH,'ERROR CREATED PST1 FILE FOR '//RT//'#'//TRIM(VTOS(RPARAM(IGRAD)))); RETURN ENDIF ENDIF END SELECT !## wait before starting a new process IF(PBMAN%NSWAIT.GT.0)THEN WRITE(*,'(/A)') '>>> Wait '//TRIM(VTOS(PBMAN%NSWAIT/100))//' seconds ...' CALL SYSTEM_CLOCK(ST) !## check residual, if you have to wait anyhow NDONE=IPEST_GLM_EVALUATE(IBATCH,ITER,NCPU,DIR,RN,RT,RPARAM,LAMBDA,ILAMBDA,MNAME) CALL SYSTEM_CLOCK(ET,CLOCK_RATE) T=100*((ET-ST)/CLOCK_RATE) WRITE(*,*) ST,ET,T,PBMAN%NSWAIT-T T=PBMAN%NSWAIT-T !## wait remaining time IF(T.GT.0)CALL IOSWAIT(PBMAN%NSWAIT) WRITE(*,'(A/)') ' Finished waiting <<<' ENDIF !## clear error I=WINFOERROR(1); IDPROC=0 CALL IOSCOMMAND(TRIM(RN(IGRAD)),IFLAGS=IFLAGS,IDPROC=IDPROC); IPROC(:,IGRAD)=IDPROC IF(WINFOERROR(1).EQ.ERROSCOMMAND)THEN CALL IPEST_GLM_ERROR(IBATCH,'FAILED TO START MODEL '//RT//'#'//TRIM(VTOS(RPARAM(IGRAD)))); RETURN ENDIF CALL IPEST_GLM_PROGRESS_NUMBERS(RN,NDONE) ENDIF !## save start time CALL SYSTEM_CLOCK(STIME(IGRAD)) !## started ISTATUS(IGRAD)=1 !## all started DO IGRAD=1,SIZE(RN); IF(ISTATUS(IGRAD).EQ.-1)EXIT; ENDDO; IF(IGRAD.GT.SIZE(RN))EXIT !## maximum number of cpu reached IF(NCPU.GE.PBMAN%NCPU)EXIT ENDDO !## evaluate processes that are finished NDONE=IPEST_GLM_EVALUATE(IBATCH,ITER,NCPU,DIR,RN,RT,RPARAM,LAMBDA,ILAMBDA,MNAME) !## finished if all succesfully completed IF(NDONE.EQ.SIZE(RN))EXIT !## first iteration only need a single run IF(RT.EQ.'L'.AND.ITER.EQ.0.AND.ISTATUS(1).EQ.0)EXIT ENDDO IPEST_GLM_RUNMODELS=.TRUE. END FUNCTION IPEST_GLM_RUNMODELS !###==================================================================== SUBROUTINE IPEST_GLM_PROGRESS_NUMBERS(RN,NDONE) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN),DIMENSION(:) :: RN INTEGER,INTENT(OUT) :: NDONE INTEGER :: IGRAD,NBUSY REAL(KIND=DP_KIND) :: F !## how much done? NDONE=0; NBUSY=0 DO IGRAD=1,SIZE(RN) IF(ISTATUS(IGRAD).EQ.0)NDONE=NDONE+1 IF(ISTATUS(IGRAD).EQ.1)NBUSY=NBUSY+1 ENDDO F=DBLE(NDONE)*100.0D0/DBLE(SIZE(RN)) WRITE(6,'(A)') '+Still running '//TRIM(VTOS(NBUSY))//' models; completed: '//TRIM(VTOS(F,'F',2))//'% (total '// & TRIM(VTOS(NDONE))//' out of '//TRIM(VTOS(SIZE(RN)))//' simulations completed)' END SUBROUTINE IPEST_GLM_PROGRESS_NUMBERS !###==================================================================== INTEGER FUNCTION IPEST_GLM_EVALUATE(IBATCH,ITER,NCPU,DIR,RN,RT,RPARAM,LAMBDA,ILAMBDA,MNAME) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: IBATCH,ITER,ILAMBDA INTEGER,INTENT(IN),DIMENSION(:) :: RPARAM INTEGER,INTENT(INOUT) :: NCPU REAL(KIND=DP_KIND),INTENT(IN) :: LAMBDA CHARACTER(LEN=*),INTENT(IN),DIMENSION(:) :: RN CHARACTER(LEN=*),INTENT(IN) :: DIR,RT,MNAME INTEGER :: ITYPE,N,JGRAD,IEXCOD !,IW TYPE(WIN_MESSAGE) :: MESSAGE INTEGER,DIMENSION(2) :: IDPROC INTEGER :: ILIN,NDONE IPEST_GLM_EVALUATE=0 CALL IPEST_GLM_PROGRESS_NUMBERS(RN,NDONE) ILIN=0 DO CALL WMESSAGE(ITYPE,MESSAGE) !## timer expired IF(ITYPE.EQ.TIMEREXPIRED)THEN N=0; DO JGRAD=1,SIZE(RN) !## all handled process IF(IPROC(1,JGRAD)+IPROC(2,JGRAD).EQ.0)CYCLE !## insensitive IF(ABS(INSENS(JGRAD)).EQ.1)THEN ISTATUS(JGRAD)=0; IEXCOD=0 !## sensitive ELSE !## check running status IDPROC=IPROC(:,JGRAD) CALL IOSCOMMANDCHECK(IDPROC,ISTATUS(JGRAD),IEXCOD=IEXCOD) ENDIF !## stopped running IF(ISTATUS(JGRAD).EQ.0)THEN !## save end time CALL SYSTEM_CLOCK(ETIME(JGRAD)) !## error occured IF(IEXCOD.NE.0)THEN CALL IPEST_GLM_ERROR(IBATCH,'ERROR OCCURED RUNNING MODEL '//RT//'#'//TRIM(VTOS(RPARAM(JGRAD)))//'; Error code '//TRIM(VTOS(IEXCOD))) !## question to continue? IF(.NOT.UTL_CONTINUE('Inspect model '//TRIM(RN(JGRAD))//'; something went wrong in the SIMULATION; Error code '//TRIM(VTOS(IEXCOD))))STOP !## try again - need to run again ISTATUS(JGRAD)=-1 ENDIF !## set part of objective function IF(ISTATUS(JGRAD).EQ.0)THEN !## check run IF(.NOT.IPEST_GLM_CHECKRUN(DIR,MNAME,RT,RPARAM(JGRAD)))THEN CALL IPEST_GLM_ERROR(IBATCH,'WARNING, AN ERROR OCCURED EXAMINING THE LIST FILE FOR MODEL '//RT//'#'//TRIM(VTOS(RPARAM(JGRAD)))//' IMOD KEEPS TRYING AGAIN') IF(.NOT.UTL_CONTINUE('Inspect model '//TRIM(RN(JGRAD))//'; something went wrong in the LIST EXAMINING'))STOP !## try again - need to run again ISTATUS(JGRAD)=-1 ENDIF IF(ISTATUS(JGRAD).EQ.0)THEN !## save start time CALL SYSTEM_CLOCK(STIMEJ(JGRAD)) IF(INSENS(JGRAD).EQ.0)THEN IF(.NOT.IPEST_GLM_GETJ_AVG(DIR,JGRAD,RPARAM(JGRAD),RT,IBATCH,ILAMBDA,MNAME))THEN !## question to continue? --- causes by MF6 that cannot write obs.txt file CALL IPEST_GLM_ERROR(IBATCH,'WARNING, AN ERROR OCCURED COLLECTING RESIDUALS FOR MODEL '//RT//'#'//TRIM(VTOS(RPARAM(JGRAD)))//' IMOD KEEPS TRYING AGAIN') IF(.NOT.UTL_CONTINUE('Inspect model '//TRIM(RN(JGRAD))//'; something went wrong in the RESIDUAL COMPUTATION'))STOP !## try again - need to run again ISTATUS(JGRAD)=-1 ENDIF IF(MSR%NOBS.GT.0)THEN IF(RT.EQ.'L'.OR.RT.EQ.'R')THEN MSR%DHL_J(JGRAD)=MSR%TJ MSR%GOF(JGRAD)=UTL_GOODNESS_OF_FIT(GF_H,GF_O,MSR%NOBS) MSR%NSC(JGRAD)=UTL_NASH_SUTCLIFFE(GF_H,GF_O,MSR%NOBS) ELSEIF(RT.EQ.'P')THEN MSR%DHG_J(JGRAD)=MSR%TJ ENDIF ENDIF!### insensitive set equal to reference head ELSE IF(MSR%NOBS.GT.0)MSR%HG(JGRAD,:)=MSR%HL(0,:) ENDIF ENDIF !## save end time CALL SYSTEM_CLOCK(ETIMEJ(JGRAD)) !## write echo CALL IPEST_GLM_PROGRESS(ITER,JGRAD,RPARAM(JGRAD),RT,LAMBDA); IF(IUPESTPROGRESS.GT.0)FLUSH(IUPESTPROGRESS) ENDIF !## reset running handle proces IPROC(:,JGRAD)=0 !## release cpu so another can be started NCPU=NCPU-1 !## finished, quit this loop, and try to restart a new new EXIT ENDIF ENDDO CALL IPEST_GLM_PROGRESS_NUMBERS(RN,NDONE) !## start another one as a proces has been stopped and there might be still one waiting in the que IF(NCPU.LT.PBMAN%NCPU)EXIT ENDIF ENDDO ! write(*,*) 'see whether to start another one' IPEST_GLM_EVALUATE=NDONE END FUNCTION IPEST_GLM_EVALUATE !###====================================================================== LOGICAL FUNCTION IPEST_GLM_CHECKRUN(DIR,MNAME,CTYPE,IPARAM) !###====================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: DIR,CTYPE,MNAME INTEGER,INTENT(IN) :: IPARAM INTEGER :: I,J,IU,IOS,N CHARACTER(LEN=256) :: FNAME,DIRO,LINE INTEGER,ALLOCATABLE,DIMENSION(:) :: ITER REAL(KIND=DP_KIND) :: QS REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: QPERC IPEST_GLM_CHECKRUN=.FALSE. ALLOCATE(ITER(PRJNPER),QPERC(PRJNPER)) IF(PBMAN%IFORMAT.EQ.2)THEN DIRO=DIR; IF(PBMAN%OUTPUT.NE.'')DIRO=PBMAN%OUTPUT FNAME=TRIM(DIRO)//'\'//TRIM(MNAME)//'_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IPARAM))//'.LIST' WRITE(*,'(A)') 'Checking LIST-file='//TRIM(FNAME) !## check whether converged IU=UTL_GETUNIT(); OPEN(IU,FILE=FNAME,STATUS='OLD',FORM='FORMATTED',ACTION='READ') N=0; DO READ(IU,'(A256)',IOSTAT=IOS) LINE IF(IOS.NE.0)EXIT IF(INDEX(LINE,'TOTAL ITERATIONS').GT.0)THEN N=N+1 READ(LINE,*,IOSTAT=IOS) ITER(N) IF(IOS.NE.0)THEN WRITE(*,'(/A)') '>>> ERROR READING TOTAL ITERATIONS' WRITE(*,'(A/)') TRIM(LINE)//' <<<' RETURN ENDIF ENDIF IF(INDEX(LINE,'PERCENT DISCREPANCY').GT.0.AND.INDEX(LINE,'PERCENT DISCREPANCY IS').EQ.0)THEN READ(LINE,'(63X,F15.2)',IOSTAT=IOS) QPERC(N); QPERC(N)=ABS(QPERC(N)) IF(IOS.NE.0)THEN WRITE(*,'(/A)') '>>> ERROR READING PERCENT DISCREPANCY' WRITE(*,'(A/)') TRIM(LINE)//' <<<' RETURN ENDIF ENDIF IF(INDEX(LINE,'FAILURE TO MEET SOLVER CONVERGENCE CRITERIA').GT.0)THEN ITER(N)=-1*ITER(N) ENDIF ENDDO CLOSE(IU) !## something went wrong IF(N.NE.PRJNPER)THEN WRITE(*,'(/A)') '>>> ERROR READING ALL TIME STEPS IN' WRITE(*,'(A/)') TRIM(FNAME)//' <<<' RETURN ENDIF !## check failure DO I=1,N; IF(ITER(I).LE.0)EXIT; ENDDO IF(I.LE.N)THEN QS=0.0D0; J=0; DO I=1,N WRITE(*,'(2I10,F15.2)') I,ITER(I),QPERC(I) IF(ITER(I).LE.0)THEN J=J+1 ELSE QS=QS+QPERC(I) ENDIF ENDDO WRITE(*,'(/A)') '>>> FAILED TO CONVERGE FOR '//TRIM(VTOS(J))//' TIME STEPS IN' WRITE(*,'(A/)') TRIM(FNAME)//' <<<' QS=QS/REAL(N-J,8) !## thresshold 10 times average budget error QS=10.0D0*QS J=0; DO I=1,N IF(ITER(I).LE.0)THEN IF(QPERC(I).GT.QS)J=J+1 ENDIF ENDDO IF(J.EQ.0)THEN WRITE(*,'(/A)') '>>> CONTINUING AS BUDGET ERROR OF FAILED MODEL IS LESS/EQUAL THE AVERAGE BUDGET ERROR OF '//TRIM(VTOS(QS,'F',2))//' % <<<' ELSE WRITE(*,'(/A)') '>>> PROCESS PAUSED AS BUDGET ERROR OF FAILED MODEL IS GREATER THAN 10 TIMES AVERAGE BUDGET ERROR OF '//TRIM(VTOS(QS,'F',2))//' %' WRITE(*,'(A/)') ' CLICK KEYBOARD TO CONTINUE OR CRTL-DEL TO STOP' PAUSE !; STOP ENDIF ENDIF ENDIF DEALLOCATE(ITER,QPERC) IPEST_GLM_CHECKRUN=.TRUE. END FUNCTION IPEST_GLM_CHECKRUN !###==================================================================== SUBROUTINE IPEST_GLM_ERROR(IBATCH,TXT) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: IBATCH CHARACTER(LEN=*),INTENT(IN) :: TXT IF(IBATCH.EQ.1)WRITE(*,'(/A/)') TRIM(TXT) IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,TRIM(TXT),'Error') WRITE(IUPESTOUT,'(/A/)') TRIM(TXT) END SUBROUTINE IPEST_GLM_ERROR !###==================================================================== LOGICAL FUNCTION IPEST_GLM_CHK(IP,IBATCH) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: IP,IBATCH INTEGER :: I,N REAL(KIND=DP_KIND) :: K1,K2 IPEST_GLM_CHK=.FALSE. DO I=1,SIZE(PARAM) PEST%PARAM(IP)%PPARAM=UTL_CAP(PEST%PARAM(IP)%PPARAM,'U') IF(TRIM(PEST%PARAM(IP)%PPARAM).EQ.TRIM(PARAM(I)))EXIT ENDDO IF(I.GT.SIZE(PARAM))THEN IF(IBATCH.EQ.1)THEN WRITE(*,'(/A)') 'Error can not recognize parameter type:'//TRIM(PEST%PARAM(IP)%PPARAM) WRITE(*,'(/A)') ' Choose from:' DO I=1,SIZE(PARAM); WRITE(*,'(A)') ' - '//TRIM(PARAM(I)); ENDDO WRITE(*,'(A)') ELSEIF(IBATCH.EQ.0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Error can not recognize parameter type:'//TRIM(PEST%PARAM(IP)%PPARAM),'Error') ENDIF RETURN ENDIF IF(PEST%PARAM(IP)%PMIN.GE.PEST%PARAM(IP)%PMAX)THEN IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'No proper parameter width defined for parameter '//TRIM(VTOS(IP)) IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'No proper parameter width defined for parameter '//TRIM(VTOS(IP)),'Error') RETURN ENDIF !## estimate variance of parameter IF(PEST%PARAM(IP)%PARSTD.EQ.0.0D0)THEN IF(PEST%PARAM(IP)%PLOG.EQ.1)THEN K1=LOG(ABS(PEST%PARAM(IP)%PMIN)) K2=LOG(ABS(PEST%PARAM(IP)%PMAX)) ELSEIF(PEST%PARAM(IP)%PLOG.EQ.2)THEN K1=LOG10(ABS(PEST%PARAM(IP)%PMIN)) K2=LOG10(ABS(PEST%PARAM(IP)%PMAX)) ELSE K1=PEST%PARAM(IP)%PMIN K2=PEST%PARAM(IP)%PMAX ENDIF PEST%PARAM(IP)%PARSTD=(K2-K1)/4.0D0 ENDIF IF(PEST%PARAM(IP)%PINI.LT.PEST%PARAM(IP)%PMIN.OR.PEST%PARAM(IP)%PINI.GT.PEST%PARAM(IP)%PMAX)THEN IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Parameter '//TRIM(VTOS(IP))//' outside parameter width' IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Parameter '//TRIM(VTOS(IP))//' outside parameter width','Error') RETURN ENDIF SELECT CASE (TRIM(PEST%PARAM(IP)%PPARAM)) CASE ('KD','KH','SC','AF','VA','SY') IF(PEST%PARAM(IP)%PILS.LE.0.OR.PEST%PARAM(IP)%PILS.GT.PRJNLAY)THEN IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Parameter '//TRIM(VTOS(IP))//': ILS exceeds NLAY' IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Parameter '//TRIM(VTOS(IP))//': ILS exceeds NLAY','Error') RETURN ENDIF CASE ('VC','KV') IF(PEST%PARAM(IP)%PILS.LE.0.OR.PEST%PARAM(IP)%PILS.GT.PRJNLAY-1)THEN IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Parameter '//TRIM(VTOS(IP))//': ILS exceeds NLAY-1' IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Parameter '//TRIM(VTOS(IP))//': ILS exceeds NLAY-1','Error') RETURN ENDIF CASE ('EP','RE') IF(PEST%PARAM(IP)%PILS.NE.1)THEN IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Parameter '//TRIM(VTOS(IP))//': ILS need to be equal to 1' IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Parameter '//TRIM(VTOS(IP))//': ILS need to be equal to 1','Error') RETURN ENDIF END SELECT !## scaling IF(PEST%PARAM(IP)%PLOG.GT.0)THEN IF(PEST%PARAM(IP)%PDELTA.EQ.1.0D0)WRITE(*,'(/A/)') 'You can not specify delta alpha eq 1.0' IF(PEST%PARAM(IP)%PMIN .EQ.0.0D0)WRITE(*,'(/A/)') 'You can not specify minimal value eq 0.0 for log-transformed parameters' IF(PEST%PARAM(IP)%PLOG.EQ.1)THEN PEST%PARAM(IP)%PINI =LOG(PEST%PARAM(IP)%PINI) PEST%PARAM(IP)%PMIN =LOG(PEST%PARAM(IP)%PMIN) PEST%PARAM(IP)%PMAX =LOG(PEST%PARAM(IP)%PMAX) PEST%PARAM(IP)%PDELTA=LOG(PEST%PARAM(IP)%PDELTA) PEST%PARAM(IP)%PPRIOR=LOG(PEST%PARAM(IP)%PPRIOR) ELSEIF(PEST%PARAM(IP)%PLOG.EQ.2)THEN PEST%PARAM(IP)%PINI =LOG10(PEST%PARAM(IP)%PINI) PEST%PARAM(IP)%PMIN =LOG10(PEST%PARAM(IP)%PMIN) PEST%PARAM(IP)%PMAX =LOG10(PEST%PARAM(IP)%PMAX) PEST%PARAM(IP)%PDELTA=LOG10(PEST%PARAM(IP)%PDELTA) PEST%PARAM(IP)%PPRIOR=LOG10(PEST%PARAM(IP)%PPRIOR) ENDIF ENDIF PEST%PARAM(IP)%ALPHA(1)=PEST%PARAM(IP)%PINI !## current alpha PEST%PARAM(IP)%ALPHA(2)=PEST%PARAM(IP)%PINI !## previous alpha ALLOCATE(PEST%PARAM(IP)%ALPHA_HISTORY(0:ABS(PEST%PE_MXITER))); PEST%PARAM(IP)%ALPHA_HISTORY=0.0D0 !## allocate memory for running the models N=0 ; DO I=1,SIZE(PEST%PARAM) !## associated parameters to existing groups inactive for gradient computation IF(PEST%PARAM(I)%PACT.EQ.1)N=N+1 ENDDO ALLOCATE(PEST%PARAM(IP)%GALPHA(N)); PEST%PARAM(IP)%GALPHA=0.0D0 !## set initial lalpha to initial value N=PBMAN%NLAMBDASEARCH; ALLOCATE(PEST%PARAM(IP)%LALPHA(N)); PEST%PARAM(IP)%LALPHA=PEST%PARAM(IP)%ALPHA(1) !## fill in default acronym IF(PEST%PARAM(IP)%ACRONYM.EQ.'')THEN WRITE(PEST%PARAM(IP)%ACRONYM,'(A2,2I5.5,I3.3)') PEST%PARAM(IP)%PPARAM,PEST%PARAM(IP)%PILS,PEST%PARAM(IP)%PIZONE,ABS(PEST%PARAM(IP)%PIGROUP) ENDIF PEST%PARAM(IP)%ACRONYM=ADJUSTR(PEST%PARAM(IP)%ACRONYM) !## final check, check whether the optimized parameter is part of the model SELECT CASE (PEST%PARAM(IP)%PPARAM) CASE ('KD'); IF(TOPICS(TKDW)%IACT_MODEL.EQ.0)STOP 'PARAMETER KD CANNOT BE OPTIMIZED AS KDW IS NOT PRESENT' CASE ('KH'); IF(TOPICS(TKHV)%IACT_MODEL.EQ.0)STOP 'PARAMETER KH CANNOT BE OPTIMIZED AS KHV IS NOT PRESENT' CASE ('KV'); IF(TOPICS(TKVV)%IACT_MODEL.EQ.0)STOP 'PARAMETER KV CANNOT BE OPTIMIZED AS KVV IS NOT PRESENT' CASE ('SC'); IF(TOPICS(TSTO)%IACT_MODEL.EQ.0)STOP 'PARAMETER SC CANNOT BE OPTIMIZED AS STO IS NOT PRESENT' CASE ('SY'); IF(TOPICS(TKHV)%IACT_MODEL.EQ.0)STOP 'PARAMETER SY CANNOT BE OPTIMIZED AS SPY IS NOT PRESENT' CASE ('VA'); IF(TOPICS(TKVA)%IACT_MODEL.EQ.0)STOP 'PARAMETER VA CANNOT BE OPTIMIZED AS KVA IS NOT PRESENT' END SELECT IPEST_GLM_CHK=.TRUE. END FUNCTION IPEST_GLM_CHK !###==================================================================== LOGICAL FUNCTION IPEST_GLM_READ_ZONES_OPENFILE(DIR,JU,NCOL,NROW) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: DIR INTEGER,INTENT(OUT) :: JU,NROW,NCOL INTEGER :: N,IOS IPEST_GLM_READ_ZONES_OPENFILE=.FALSE. JU=UTL_GETUNIT(); OPEN(JU,FILE=TRIM(DIR)//'\PARAM_DUMP_IPEST_IMOD.DAT',STATUS='OLD',ACTION='READ',FORM='FORMATTED') READ(JU,*) READ(JU,*) READ(JU,'(22X,I10)',IOSTAT=IOS) NCOL READ(JU,'(22X,I10)',IOSTAT=IOS) NROW READ(JU,'(22X,I10)',IOSTAT=IOS) N IF(IOS.NE.0)THEN CLOSE(JU); WRITE(*,'(/A)') 'OOPS, AN OLD VERSION OF A PARAM_DUMP_IPEST_IMOD FOUND, NEED TO BE RECREATED'; RETURN ENDIF IF(N.NE.SIZE(PEST%PARAM))THEN CLOSE(JU); WRITE(*,'(/A)') 'NUMBER OF ZONES IN THE PARAM_DUMP_IPEST_IMOD DOES NOT COINCIDE WITH THE ACTUAL NUMBER OF PARAMETERS'; RETURN ENDIF IPEST_GLM_READ_ZONES_OPENFILE=.TRUE. END FUNCTION IPEST_GLM_READ_ZONES_OPENFILE !###==================================================================== SUBROUTINE IPEST_GLM_READ_ZONES(DIR,NCOL,NROW) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: DIR INTEGER,INTENT(OUT) :: NCOL,NROW INTEGER :: JU,I,J IF(.NOT.IPEST_GLM_READ_ZONES_OPENFILE(DIR,JU,NCOL,NROW))STOP DO I=1,SIZE(PEST%PARAM) READ(JU,*) READ(JU,*) READ(JU,'(2I10)') PEST%PARAM(I)%NODES,PEST%PARAM(I)%ZTYPE IF(PEST%PARAM(I)%ZTYPE.EQ.0)THEN ALLOCATE(PEST%PARAM(I)%IROW(PEST%PARAM(I)%NODES), & PEST%PARAM(I)%ICOL(PEST%PARAM(I)%NODES), & PEST%PARAM(I)%F( PEST%PARAM(I)%NODES)) READ(JU,*) DO J=1,PEST%PARAM(I)%NODES READ(JU,'(2I10,F10.4)') PEST%PARAM(I)%IROW(J),PEST%PARAM(I)%ICOL(J),PEST%PARAM(I)%F(J) ENDDO ELSE ALLOCATE(PEST%PARAM(I)%XY(PEST%PARAM(I)%NODES,2)) READ(JU,*) DO J=1,PEST%PARAM(I)%NODES READ(JU,*) PEST%PARAM(I)%XY(J,1),PEST%PARAM(I)%XY(J,2) ENDDO ENDIF ! IF(PEST%PARAM(I)%PPARAM.EQ.'HF')THEN ! PEST%PARAM(I)%NODES=0 !## one single cell used as zone for horizontal barrier module ! ELSE IF(PEST%PARAM(I)%NODES.EQ.0)PEST%PARAM(I)%PACT=0 ! ENDIF ENDDO CLOSE(JU) END SUBROUTINE IPEST_GLM_READ_ZONES !###==================================================================== SUBROUTINE IPEST_GLM_RESET_PARAMETER(DIR) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: DIR INTEGER :: I LOGICAL :: LEX !## scaling WRITE(*,*) '>>> Cleaning for '//TRIM(VTOS(SIZE(PEST%PARAM)))//' parameters ...' DO I=1,SIZE(PEST%PARAM) PEST%PARAM(I)%PACT=ABS(PEST%PARAM(I)%PACT) IF(PEST%PARAM(I)%PLOG.EQ.1)THEN PEST%PARAM(I)%PINI =EXP(PEST%PARAM(I)%PINI) PEST%PARAM(I)%PMIN =EXP(PEST%PARAM(I)%PMIN) PEST%PARAM(I)%PMAX =EXP(PEST%PARAM(I)%PMAX) PEST%PARAM(I)%PDELTA=EXP(PEST%PARAM(I)%PDELTA) PEST%PARAM(I)%PPRIOR=EXP(PEST%PARAM(I)%PPRIOR) ELSEIF(PEST%PARAM(I)%PLOG.EQ.2)THEN PEST%PARAM(I)%PINI =10.0D0**(PEST%PARAM(I)%PINI) PEST%PARAM(I)%PMIN =10.0D0**(PEST%PARAM(I)%PMIN) PEST%PARAM(I)%PMAX =10.0D0**(PEST%PARAM(I)%PMAX) PEST%PARAM(I)%PDELTA=10.0D0**(PEST%PARAM(I)%PDELTA) PEST%PARAM(I)%PPRIOR=10.0D0**(PEST%PARAM(I)%PPRIOR) ENDIF IF(PBMAN%ICLEAN.EQ.1)THEN IF(ABS(PEST%PARAM(I)%PACT).EQ.1.AND.PEST%PARAM(I)%PIGROUP.GT.0)THEN INQUIRE(DIRECTORY=TRIM(DIR)//'\IPEST_P#'//TRIM(VTOS(I)),EXIST=LEX) WRITE(*,'(A)') ' - '//TRIM(DIR)//'\IPEST_P#'//TRIM(VTOS(I)) IF(LEX)THEN IF(.NOT.UTL_DEL1TREE(TRIM(DIR)//'\IPEST_P#'//TRIM(VTOS(I)),IQUESTION=0))THEN WRITE(*,'(/A)') '>>> Could not delete folder:' WRITE(*,'(A/)') ' '//TRIM(DIR)//'\IPEST_P#'//TRIM(VTOS(I)) ENDIF ENDIF ENDIF ENDIF ENDDO WRITE(*,*) ' Finished <<<' WRITE(*,*) ' Finished <<<' END SUBROUTINE IPEST_GLM_RESET_PARAMETER !###==================================================================== SUBROUTINE IPEST_GLM_CLOSE_FILES() !###==================================================================== IMPLICIT NONE ! return WRITE(*,*) 'Closing files' !return ! WRITE(*,*) IUPESTOUT IF(IUPESTOUT.GT.0) CLOSE(IUPESTOUT); IUPESTOUT=0 ! WRITE(*,*) IUPESTPROGRESS IF(IUPESTPROGRESS.GT.0) CLOSE(IUPESTPROGRESS); IUPESTPROGRESS=0 ! WRITE(*,*) IUPESTEFFICIENCY IF(IUPESTEFFICIENCY.GT.0) CLOSE(IUPESTEFFICIENCY); IUPESTEFFICIENCY=0 ! WRITE(*,*) IUPESTSENSITIVITY IF(IUPESTSENSITIVITY.GT.0)CLOSE(IUPESTSENSITIVITY); IUPESTSENSITIVITY=0 ! WRITE(*,*) IUPESTRUNFILE IF(IUPESTRUNFILE.GT.0) CLOSE(IUPESTRUNFILE); IUPESTRUNFILE=0 ! WRITE(*,*) IUPESTJACOBIAN IF(IUPESTJACOBIAN.GT.0) CLOSE(IUPESTJACOBIAN); IUPESTJACOBIAN=0 END SUBROUTINE IPEST_GLM_CLOSE_FILES !###==================================================================== SUBROUTINE IPEST_GLM_ECHO_PARAMETERS(IP) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: IP INTEGER :: I CHARACTER(LEN=12) :: STRING IF(PBMAN%IPESTP.EQ.1)THEN STRING=TRIM(VTOS(PEST%PARAM(IP)%NODES)) ELSE STRING='REPLACE' ENDIF I=PEST%PARAM(IP)%ORG_PACT !## parameter off in resulting run-file IF(PEST%PE_SENS.LT.0.0D0.AND.PEST%PARAM(IP)%INSENS.EQ.1)I=0 IF(PEST%PARAM(IP)%PLOG.EQ.0)THEN WRITE(IUPESTOUT,'(I2,1X,A5,1X,I5,1X,I7,5(1X,F15.3),I10,I10,A10,A15,2F15.3,2A15)') I,PEST%PARAM(IP)%PPARAM,PEST%PARAM(IP)%PILS, & PEST%PARAM(IP)%PIZONE,PEST%PARAM(IP)%PINI,PEST%PARAM(IP)%PDELTA,PEST%PARAM(IP)%PMIN,PEST%PARAM(IP)%PMAX,PEST%PARAM(IP)%PINCREASE, & ABS(PEST%PARAM(IP)%PIGROUP),PEST%PARAM(IP)%PLOG,TRIM(STRING),PEST%PARAM(IP)%ACRONYM,PEST%PARAM(IP)%PPRIOR,PEST%PARAM(IP)%PARSTD, & PEST%PARAM(IP)%SDATE,PEST%PARAM(IP)%EDATE PEST%PARAM(IP)%ALPHA_HISTORY(0)=PEST%PARAM(IP)%PINI ELSEIF(PEST%PARAM(IP)%PLOG.EQ.1)THEN WRITE(IUPESTOUT,'(I2,1X,A5,1X,I5,1X,I7,5(1X,F15.3),I10,I10,A10,A15,2F15.3,2A15)') I,PEST%PARAM(IP)%PPARAM,PEST%PARAM(IP)%PILS, & PEST%PARAM(IP)%PIZONE,EXP(PEST%PARAM(IP)%PINI),EXP(PEST%PARAM(IP)%PDELTA),EXP(PEST%PARAM(IP)%PMIN),EXP(PEST%PARAM(IP)%PMAX),PEST%PARAM(IP)%PINCREASE, & ABS(PEST%PARAM(IP)%PIGROUP),PEST%PARAM(IP)%PLOG,TRIM(STRING),PEST%PARAM(IP)%ACRONYM,EXP(PEST%PARAM(IP)%PPRIOR),EXP(PEST%PARAM(IP)%PARSTD), & PEST%PARAM(IP)%SDATE,PEST%PARAM(IP)%EDATE PEST%PARAM(IP)%ALPHA_HISTORY(0)=EXP(PEST%PARAM(IP)%PINI) ELSEIF(PEST%PARAM(IP)%PLOG.EQ.2)THEN WRITE(IUPESTOUT,'(I2,1X,A5,1X,I5,1X,I7,5(1X,F15.3),I10,I10,A10,A15,2F15.3,2A15)') I,PEST%PARAM(IP)%PPARAM,PEST%PARAM(IP)%PILS, & PEST%PARAM(IP)%PIZONE,10.0D0**(PEST%PARAM(IP)%PINI),10.0D0**PEST%PARAM(IP)%PDELTA,10.0D0**(PEST%PARAM(IP)%PMIN),10.0D0**(PEST%PARAM(IP)%PMAX),PEST%PARAM(IP)%PINCREASE, & ABS(PEST%PARAM(IP)%PIGROUP),PEST%PARAM(IP)%PLOG,TRIM(STRING),PEST%PARAM(IP)%ACRONYM,10.0D0**(PEST%PARAM(IP)%PPRIOR),10.0D0**PEST%PARAM(IP)%PARSTD, & PEST%PARAM(IP)%SDATE,PEST%PARAM(IP)%EDATE PEST%PARAM(IP)%ALPHA_HISTORY(0)=10.0D0**(PEST%PARAM(IP)%PINI) ENDIF END SUBROUTINE IPEST_GLM_ECHO_PARAMETERS !#####================================================================= LOGICAL FUNCTION IPEST_GLM_SETGROUPS(IBATCH) !#####================================================================= IMPLICIT NONE INTEGER,INTENT(IN) :: IBATCH INTEGER :: I,J,NG,NP IPEST_GLM_SETGROUPS=.TRUE.; IF(PBMAN%IPESTP+PBMAN%IIES.EQ.0)RETURN IPEST_GLM_SETGROUPS=.FALSE. ! !## reset groups ! DO I=1,SIZE(PEST%PARAM) ! PEST%PARAM(I)%PACT=PEST%PARAM(I)%ORG_PACT; PEST%PARAM(I)%PIGROUP=ABS(PEST%PARAM(I)%PIGROUP) ! ENDDO !## get number of active parameters (number of groups need to be zero) NP=0; NG=0; DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PACT.EQ.1)NP=NP+1 IF(PEST%PARAM(I)%PACT.EQ.2)NG=NG+1 ENDDO WRITE(*,'(/A)') 'Got Initially '//TRIM(VTOS(NP+NG))//' defined parameters; number of active parameters/groups '//TRIM(VTOS(NP)) !## set igroup lt 0 for followers in group - check whether factors within group are equal --- need to be !## make sure group with active nodes is positive rest is negative IF(IUPESTOUT.NE.0)WRITE(IUPESTOUT,'(/A10,1X,A2,1X,A15,1X,A10)') 'NO','P','ACRONYM','GROUP_SIZE' IF(IUPESTOUT.EQ.0)WRITE(* ,'(/A10,1X,A2,1X,A15,1X,A10)') 'NO','P','ACRONYM','GROUP_SIZE' NP=0; DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PACT.NE.1)CYCLE !## get all parameters in same group for current parameter NG=1; NP=NP+1; DO J=1,SIZE(PEST%PARAM) IF(I.EQ.J)CYCLE IF(PEST%PARAM(J)%PACT.NE.1)CYCLE IF(PEST%PARAM(J)%PIGROUP.EQ.PEST%PARAM(I)%PIGROUP)THEN !## check factor IF(PEST%PARAM(J)%PINI.NE.PEST%PARAM(I)%PINI)THEN IF(IBATCH.EQ.1)THEN WRITE(*,'(/A)') 'Initial factor in an group need to be identicial' IF(PEST%PARAM(J)%PACT.EQ.0)WRITE(*,'(A)') '>>> Parameter is inactive, but it needs to be active if a parameter in that group is active too <<<' WRITE(*,'(A/)') 'Check initial factors for group '//TRIM(VTOS(PEST%PARAM(J)%PIGROUP)) STOP ELSE CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Parameters within the same group need to have idential initial values','Error') RETURN ENDIF ENDIF PEST%PARAM(J)%PIGROUP=-1*ABS(PEST%PARAM(J)%PIGROUP) !## make sure groupsnames are equally PEST%PARAM(J)%ACRONYM= PEST%PARAM(I)%ACRONYM PEST%PARAM(J)%PACT=2 NG=NG+1 ENDIF ENDDO IF(IUPESTOUT.NE.0)WRITE(IUPESTOUT,'(I10,1X,A2,1X,A15,1X,I10)') NP,PEST%PARAM(I)%PPARAM,PEST%PARAM(I)%ACRONYM,NG IF(IUPESTOUT.EQ.0)WRITE(* ,'(I10,1X,A2,1X,A15,1X,I10)') NP,PEST%PARAM(I)%PPARAM,PEST%PARAM(I)%ACRONYM,NG ENDDO !## check labels equal and unique DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PACT.NE.1)CYCLE !## get all parameters in same group for current parameter DO J=1,SIZE(PEST%PARAM) IF(I.EQ.J)CYCLE !## skip inactive parameters IF(PEST%PARAM(J)%PACT.EQ.0)CYCLE !## skip group parameters IF(PEST%PARAM(J)%PACT.EQ.2)CYCLE IF(PEST%PARAM(I)%ACRONYM.EQ.PEST%PARAM(J)%ACRONYM)THEN WRITE(*,'(3I10,A15)') I,PEST%PARAM(I)%PACT,PEST%PARAM(I)%PIGROUP,PEST%PARAM(I)%ACRONYM WRITE(*,'(3I10,A15)') J,PEST%PARAM(J)%PACT,PEST%PARAM(J)%PIGROUP,PEST%PARAM(J)%ACRONYM WRITE(*,'(/A/)') '>>> Duplicate groupslabels assigned which is not allowed <<<'; PAUSE ENDIF ENDDO ENDDO !## check whether pact=0 in combination with igroup DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PIGROUP.LT.0.AND.PEST%PARAM(I)%PACT.NE.2)THEN WRITE(*,'(/2A10)') 'ACTIVE','GROUP' WRITE(*,'(2I10)') PEST%PARAM(I)%PACT,PEST%PARAM(I)%PIGROUP WRITE(*,'(/A/)') '>>> Groups are not assigned correctly <<<'; PAUSE ENDIF ENDDO !## get number of active parameters NP=0; NG=0; DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PACT.EQ.1)NP=NP+1 IF(PEST%PARAM(I)%PACT.EQ.2)NG=NG+1 ENDDO WRITE(*,'(/A)') 'Got Initially '//TRIM(VTOS(NP+NG))//' defined parameters; number of active parameters/groups '//TRIM(VTOS(NP)) IPEST_GLM_SETGROUPS=.TRUE. END FUNCTION IPEST_GLM_SETGROUPS !#####================================================================= LOGICAL FUNCTION IPEST_GLM_PST(DIRNAME,DIR,MNAME,IGRAD,IPARAM,CTYPE) !#####================================================================= IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: DIR,CTYPE,MNAME,DIRNAME INTEGER,INTENT(IN) :: IGRAD,IPARAM CHARACTER(LEN=256) :: FNAME,FN1,FN2 INTEGER :: I,N,IU,JU,IOS,NCOL,NROW,ICOL,IROW,IEQ REAL(KIND=DP_KIND) :: X1,Y1,X2,Y2,DX !,A REAL(KIND=DP_KIND),DIMENSION(:),ALLOCATABLE :: SX,SY IPEST_GLM_PST=.FALSE. FNAME=TRIM(DIR)//'\MODELINPUT\'//TRIM(MNAME)//'_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IPARAM))//'.PST1' IU=UTL_GETUNIT(); OPEN(IU,FILE=FNAME,STATUS='OLD',ACTION='READ',IOSTAT=IOS) IF(IOS.NE.0)THEN WRITE(*,'(/A)') 'Cannot read the appropriate PST1-file for parameter '//TRIM(VTOS(IGRAD)) WRITE(*,'(A/)') TRIM(FNAME) STOP ENDIF FNAME=TRIM(DIR)//'\MODELINPUT\'//TRIM(MNAME)//'_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IPARAM))//'.PST1_' JU=UTL_GETUNIT(); OPEN(JU,FILE=FNAME,STATUS='REPLACE',ACTION='WRITE',IOSTAT=IOS) IF(IOS.NE.0)THEN WRITE(*,'(/A)') 'Cannot update the appropriate PST-file for parameter '//TRIM(VTOS(IGRAD)) WRITE(*,'(A/)') TRIM(FNAME) STOP ENDIF !## copy header READ(IU,'(A256)',IOSTAT=IOS) LINE; IF(IOS.NE.0)RETURN; WRITE(JU,'(A)') TRIM(LINE) READ(IU,'(A256)',IOSTAT=IOS) LINE; IF(IOS.NE.0)RETURN; WRITE(JU,'(A)') TRIM(LINE) READ(LINE,*,IOSTAT=IOS) NCOL,NROW; IF(IOS.NE.0)RETURN READ(IU,*,IOSTAT=IOS) X1,Y1,X2,Y2,IEQ; IF(IOS.NE.0)RETURN; WRITE(JU,*) X1,Y1,X2,Y2,IEQ IF(IEQ.EQ.0)THEN READ(IU,*,IOSTAT=IOS) DX; IF(IOS.NE.0)RETURN; WRITE(JU,*) DX ELSE ALLOCATE(SX(0:NCOL)); READ(IU,*,IOSTAT=IOS) (SX(ICOL),ICOL=0,NCOL); IF(IOS.NE.0)RETURN; WRITE(JU,*) (SX(ICOL),ICOL=0,NCOL); DEALLOCATE(SX) ALLOCATE(SY(0:NROW)); READ(IU,*,IOSTAT=IOS) (SY(IROW),IROW=0,NROW); IF(IOS.NE.0)RETURN; WRITE(JU,*) (SY(IROW),IROW=0,NROW); DEALLOCATE(SY) ENDIF !## copy measurements READ(IU,*,IOSTAT=IOS) N; IF(IOS.NE.0)RETURN; WRITE(JU,*) N DO I=1,ABS(N) READ(IU,'(A256)',IOSTAT=IOS) LINE IF(IOS.NE.0)RETURN !## ensure a global path here IF(INDEX(LINE,'.\').GT.0)THEN CALL UTL_SUBST(LINE,'.\',TRIM(DIRNAME)//'\') CALL UTL_SUBST(LINE,'\\','\') ENDIF WRITE(JU,'(A)') TRIM(LINE) ENDDO !## copy parameters READ(IU,*,IOSTAT=IOS) N; IF(IOS.NE.0)RETURN; WRITE(JU,*) N READ(IU,'(A256)',IOSTAT=IOS) LINE; IF(IOS.NE.0)RETURN; WRITE(JU,'(A)') TRIM(LINE) !## write blankout idf IF(PEST%PE_KTYPE.LT.0)THEN READ(IU,'(A256)',IOSTAT=IOS) LINE; IF(IOS.NE.0)RETURN; WRITE(JU,'(A)') TRIM(LINE) ENDIF N=0; IF(ASSOCIATED(PEST%S_PERIOD))N=SIZE(PEST%S_PERIOD) IF(N.GT.0)THEN DO I=1,SIZE(PEST%S_PERIOD) READ(IU,'(A256)',IOSTAT=IOS) LINE; IF(IOS.NE.0)RETURN; WRITE(JU,'(A)') TRIM(LINE) ENDDO ENDIF !## no batchfiles no matter what in pst-files, iMOD takes care of it N=0 !; IF(ASSOCIATED(PEST%B_FRACTION))N=SIZE(PEST%B_FRACTION) IF(N.GT.0)THEN DO I=1,SIZE(PEST%B_FRACTION) READ(IU,'(A256)',IOSTAT=IOS) LINE; IF(IOS.NE.0)RETURN; WRITE(JU,'(A)') TRIM(LINE) ENDDO ENDIF DO I=1,SIZE(PEST%PARAM) !## read old settings READ(IU,'(A256)',IOSTAT=IOS) LINE; IF(IOS.NE.0)RETURN IF(PEST%PARAM(I)%PLOG.EQ.0)THEN LINE=TRIM(VTOS(MIN(1,ABS(PEST%PARAM(I)%PACT))))//','// & TRIM(PEST%PARAM(I)%PPARAM) //','// & TRIM(VTOS(PEST%PARAM(I)%PILS)) //','// & TRIM(VTOS(PEST%PARAM(I)%PIZONE)) //','// & TRIM(VTOS(PEST%PARAM(I)%ALPHA(1),'G',7)) //','// & TRIM(VTOS(PEST%PARAM(I)%PDELTA,'G',7)) //','// & TRIM(VTOS(PEST%PARAM(I)%PMIN/10.0D0,'G',7)) //','// & TRIM(VTOS(PEST%PARAM(I)%PMAX*10.0D0,'G',7)) //','// & TRIM(VTOS(PEST%PARAM(I)%PINCREASE,'G',7))//','// & TRIM(VTOS(ABS(PEST%PARAM(I)%PIGROUP))) //','// & TRIM(VTOS(PEST%PARAM(I)%PLOG))//','// & TRIM(PEST%PARAM(I)%ACRONYM)//','// & TRIM(VTOS(PEST%PARAM(I)%PPRIOR,'G',7))//','// & TRIM(VTOS(PEST%PARAM(I)%PARSTD,'G',7))//','// & TRIM(PEST%PARAM(I)%SDATE)//','// & TRIM(PEST%PARAM(I)%EDATE) ELSEIF(PEST%PARAM(I)%PLOG.EQ.1)THEN LINE=TRIM(VTOS(MIN(1,ABS(PEST%PARAM(I)%PACT))))//','// & TRIM(PEST%PARAM(I)%PPARAM) //','// & TRIM(VTOS(PEST%PARAM(I)%PILS)) //','// & TRIM(VTOS(PEST%PARAM(I)%PIZONE)) //','// & TRIM(VTOS(EXP(PEST%PARAM(I)%ALPHA(1)),'G',7))//','// & TRIM(VTOS(EXP(PEST%PARAM(I)%PDELTA),'G',7))//','// & TRIM(VTOS(EXP(PEST%PARAM(I)%PMIN)/10.0D0,'G',7)) //','// & TRIM(VTOS(EXP(PEST%PARAM(I)%PMAX)*10.0D0,'G',7)) //','// & TRIM(VTOS(PEST%PARAM(I)%PINCREASE,'G',7)) //','// & TRIM(VTOS(ABS(PEST%PARAM(I)%PIGROUP))) //','// & TRIM(VTOS(PEST%PARAM(I)%PLOG))//','// & TRIM(PEST%PARAM(I)%ACRONYM)//','// & TRIM(VTOS(EXP(PEST%PARAM(I)%PPRIOR),'G',7))//','// & TRIM(VTOS(EXP(PEST%PARAM(I)%PARSTD),'G',7))//','// & TRIM(PEST%PARAM(I)%SDATE)//','// & TRIM(PEST%PARAM(I)%EDATE) ELSEIF(PEST%PARAM(I)%PLOG.EQ.2)THEN LINE=TRIM(VTOS(MIN(1,ABS(PEST%PARAM(I)%PACT))))//','// & TRIM(PEST%PARAM(I)%PPARAM) //','// & TRIM(VTOS(PEST%PARAM(I)%PILS)) //','// & TRIM(VTOS(PEST%PARAM(I)%PIZONE)) //','// & TRIM(VTOS(10.0D0**(PEST%PARAM(I)%ALPHA(1)),'G',7))//','// & TRIM(VTOS(10.0D0**PEST%PARAM(I)%PDELTA,'G',7))//','// & TRIM(VTOS(10.0D0**(PEST%PARAM(I)%PMIN)/10.0D0,'G',7)) //','// & TRIM(VTOS(10.0D0**(PEST%PARAM(I)%PMAX)*10.0D0,'G',7)) //','// & TRIM(VTOS(PEST%PARAM(I)%PINCREASE,'G',7)) //','// & TRIM(VTOS(ABS(PEST%PARAM(I)%PIGROUP))) //','// & TRIM(VTOS(PEST%PARAM(I)%PLOG))//','// & TRIM(PEST%PARAM(I)%ACRONYM)//','// & TRIM(VTOS(10.0D0**PEST%PARAM(I)%PPRIOR,'G',7))//','// & TRIM(VTOS(10.0D0**PEST%PARAM(I)%PARSTD,'G',7))//','// & TRIM(PEST%PARAM(I)%SDATE)//','// & TRIM(PEST%PARAM(I)%EDATE) ENDIF WRITE(JU,'(A)') TRIM(LINE) IF(TRIM(CTYPE).EQ.'L')THEN PEST%PARAM(I)%LALPHA(IGRAD)=PEST%PARAM(I)%ALPHA(1) ELSEIF(TRIM(CTYPE).EQ.'P')THEN PEST%PARAM(I)%GALPHA(IGRAD)=PEST%PARAM(I)%ALPHA(1) ENDIF ENDDO !## copy zones READ(IU,*,IOSTAT=IOS) N; IF(IOS.NE.0)RETURN; WRITE(JU,*) N DO I=1,N; READ(IU,'(A256)',IOSTAT=IOS) LINE; IF(IOS.NE.0)RETURN; WRITE(JU,'(A)') TRIM(LINE); ENDDO CLOSE(IU,STATUS='DELETE'); CLOSE(JU) FN1=TRIM(DIR)//'\MODELINPUT\'//TRIM(MNAME)//'_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IPARAM))//'.PST1_' FN2=TRIM(DIR)//'\MODELINPUT\'//TRIM(MNAME)//'_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IPARAM))//'.PST1' CALL IOSRENAMEFILE(FN1,FN2) IPEST_GLM_PST=.TRUE. END FUNCTION IPEST_GLM_PST !#####================================================================= SUBROUTINE IPEST_GLM_ALLOCATE(DIR) !#####================================================================= IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: DIR INTEGER :: I,N,M LOGICAL :: LEX CALL IPEST_GLM_DEALLOCATE() !## allocate memory for running the models N=0 ; DO I=1,SIZE(PEST%PARAM) !## associated parameters to existing groups inactive for gradient computation IF(PEST%PARAM(I)%PACT.EQ.1)N=N+1 ENDDO M=PBMAN%NLAMBDASEARCH WRITE(*,'(/1X,A,I10,A/)') 'Optimizing (nett groups included): ',N,' Parameters/Groups' ALLOCATE(RNG(N),RNL(M),GPARAM(N),LPARAM(M)) !## allocate memory for process-status memory M=MAX(M,N); ALLOCATE(ISTATUS(M),INSENS(M),IPROC(2,M),STIME(M),ETIME(M),STIMEJ(M),ETIMEJ(M)) !## set linesearch-runbatch-files N=PBMAN%NLAMBDASEARCH; IF(PEST%PE_MXITER.LT.0)N=1 DO I=1,N !## modflow6 IF(PBMAN%IFORMAT.EQ.3)THEN RNL(I)=TRIM(DIR)//'\IPEST_L#'//TRIM(VTOS(I))//'\RUN_L#'//TRIM(VTOS(I))//'.BAT' ELSE RNL(I)=TRIM(DIR)//'\RUN_L#'//TRIM(VTOS(I))//'.BAT' ENDIF LPARAM(I)=I INQUIRE(FILE=RNL(I),EXIST=LEX) IF(.NOT.LEX)THEN WRITE(*,'(/A)') 'Cannot read the appropriate BAT-file for parameter '//TRIM(VTOS(I)) WRITE(*,'(A/)') TRIM(RNL(I)) STOP ENDIF RNL(I)='"'//TRIM(RNL(I))//'"' ENDDO !## set gradient-runbatch-files - equal to number of parameter to be estimated IF(PEST%PE_MXITER.GE.0)THEN N=0; DO I=1,SIZE(PEST%PARAM) !## parameter IF(PEST%PARAM(I)%PACT.NE.1)CYCLE; N=N+1 !.OR.PEST%PARAM(I)%PIGROUP.LT.0)CYCLE; N=N+1 IF(PBMAN%IFORMAT.EQ.3)THEN RNG(N)=TRIM(DIR)//'\IPEST_P#'//TRIM(VTOS(I))//'\RUN_P#'//TRIM(VTOS(I))//'.BAT' ELSE RNG(N)=TRIM(DIR)//'\RUN_P#'//TRIM(VTOS(I))//'.BAT' ENDIF INQUIRE(FILE=RNG(N),EXIST=LEX) IF(.NOT.LEX)THEN WRITE(*,'(/A)') 'Cannot read the appropriate BAT-file for parameter '//TRIM(VTOS(I)) WRITE(*,'(A/)') TRIM(RNG(N)) STOP ENDIF GPARAM(N)=I; RNG(N)='"'//TRIM(RNG(N))//'"' ENDDO ENDIF END SUBROUTINE IPEST_GLM_ALLOCATE !#####================================================================= SUBROUTINE IPEST_GLM_DEALLOCATE() !#####================================================================= IMPLICIT NONE IF(ALLOCATED(RNG)) DEALLOCATE(RNG) IF(ALLOCATED(RNL)) DEALLOCATE(RNL) IF(ALLOCATED(IPROC)) DEALLOCATE(IPROC) IF(ALLOCATED(GPARAM)) DEALLOCATE(GPARAM) IF(ALLOCATED(LPARAM)) DEALLOCATE(LPARAM) IF(ALLOCATED(ISTATUS))DEALLOCATE(ISTATUS) IF(ALLOCATED(INSENS)) DEALLOCATE(INSENS) IF(ALLOCATED(STIME)) DEALLOCATE(STIME) IF(ALLOCATED(ETIME)) DEALLOCATE(ETIME) IF(ALLOCATED(STIMEJ)) DEALLOCATE(STIMEJ) IF(ALLOCATED(ETIMEJ)) DEALLOCATE(ETIMEJ) END SUBROUTINE IPEST_GLM_DEALLOCATE !#####================================================================= LOGICAL FUNCTION IPEST_GLM_NEXT(IBATCH,ITER,DIR,DIRO,DIRP,LAMBDA,GAMMA,MNAME) !#####================================================================= IMPLICIT NONE INTEGER,INTENT(IN) :: ITER,IBATCH CHARACTER(LEN=*),INTENT(IN) :: DIR,DIRO,DIRP,MNAME REAL(KIND=DP_KIND),INTENT(INOUT) :: GAMMA REAL(KIND=DP_KIND),INTENT(INOUT) :: LAMBDA REAL(KIND=DP_KIND) :: MJ,X INTEGER :: I,IJ REAL(KIND=DP_KIND),DIMENSION(3) :: XLAMBDA,YLAMBDA,CPARABOL IPEST_GLM_NEXT=.FALSE. !## get optimal model with lowest objective function value IJ=0; MJ=HUGE(1.0D0); DO I=1,PBMAN%NLAMBDASEARCH XLAMBDA(I)=PBMAN%LAMBDA_TEST(I)*LAMBDA YLAMBDA(I)=MSR%DHL_J(I) IF(MSR%DHL_J(I).LT.MJ)THEN MJ=MSR%DHL_J(I) IJ=I ENDIF ENDDO !## store winning lambda LAMBDAS(ITER)=XLAMBDA(IJ) !## determine lambda CALL UTL_FIT_PARABOLA(CPARABOL,XLAMBDA,YLAMBDA) ! !## valley at this lambda value ! X=-CPARABOL(2)/(2.0D0*CPARABOL(3)) IF(MJ.LT.MSR%PJ)THEN !## bot-parabola and objective function reduced IF(CPARABOL(3).GT.0.0D0)THEN !## valley at this lambda value X=-CPARABOL(2)/(2.0D0*CPARABOL(3)) !## update lambda LAMBDA=MAX(LAMBDA/GAMMA,X) WRITE(IUPESTPROGRESS,*) WRITE(IUPESTPROGRESS,'(2(A,F15.3))') 'BOT PARABOLA @ LAMBDA= ', X,' NEW LAMBDA= ',LAMBDA WRITE(IUPESTPROGRESS,'(1(A,F15.3))') 'ESTIMATED OBJECTIVE FUNCTION VALUE (NEXT TIME) = ', CPARABOL(1)+CPARABOL(2)*X+CPARABOL(3)*X**2.0D0 ELSE LAMBDA=LAMBDA*PBMAN%LAMBDA_TEST(IJ) ENDIF !## no reduction start with largest lambda to test another set of lambda's ELSE LAMBDA=PBMAN%LAMBDA_TEST(PBMAN%NLAMBDASEARCH)*LAMBDA ENDIF WRITE(IUPESTPROGRESS,'(/A)') 'New Lambda_0 '//TRIM(VTOS(LAMBDA,'F',7))//' (for upper boundary of lambda-tests); Objective Function Value '//TRIM(VTOS(MJ,'F',7)) !## is there no reduction achieved - return and try another set of lambdas IF(MJ.GE.MSR%PJ)THEN !## stop, no progression found anymore GAMMA=GAMMA*2.0D0 IF(GAMMA.GT.10.0D0)THEN GAMMA=-1.0D0*GAMMA !## save alphas for history-logging DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PACT.EQ.0)CYCLE IF(PEST%PARAM(I)%PLOG.EQ.1)THEN PEST%PARAM(I)%ALPHA_HISTORY(ITER)=EXP(PEST%PARAM(I)%ALPHA(1)) ELSEIF(PEST%PARAM(I)%PLOG.EQ.2)THEN PEST%PARAM(I)%ALPHA_HISTORY(ITER)=10.0D0**(PEST%PARAM(I)%ALPHA(1)) ELSE PEST%PARAM(I)%ALPHA_HISTORY(ITER)=PEST%PARAM(I)%ALPHA(1) ENDIF ENDDO WRITE(IUPESTPROGRESS,'(/A)') 'No reduction of objective function found in lambda search, process finished.' ENDIF !## enlarge gamma to be steepest descent method more than Gauss-Newton LAMBDA=LAMBDA*GAMMA RETURN ENDIF !## update new objective function value to be minimized MSR%TJ=MJ !## save residuals IF(IUPESTRESIDUAL.GT.0)CLOSE(IUPESTRESIDUAL); IUPESTRESIDUAL=0 IUPESTRESIDUAL=UTL_GETUNIT(); OPEN(IUPESTRESIDUAL,FILE=TRIM(DIRP)//'\LOG_PEST_RESIDUAL_'// & TRIM(VTOS(ITER))//'.TXT',STATUS='REPLACE',ACTION='WRITE') IF(IUPESTPRESIDUAL.GT.0)CLOSE(IUPESTPRESIDUAL); IUPESTPRESIDUAL=0 IUPESTPRESIDUAL=UTL_GETUNIT(); OPEN(IUPESTPRESIDUAL,FILE=TRIM(DIRP)//'\LOG_PESTP_RESIDUAL_'// & TRIM(VTOS(ITER))//'.TXT',STATUS='REPLACE',ACTION='WRITE') IF(.NOT.IPEST_GLM_GETJ_AVG(DIRO,IJ,LPARAM(IJ),'L',IBATCH,0,MNAME))RETURN !## set current h on correct position DO I=1,MSR%NOBS MSR%HL(0,I) =MSR%HL(IJ,I) MSR%DHL(0,I)=MSR%DHL(IJ,I) ENDDO MSR%GOF(1)=MSR%GOF(IJ); MSR%NSC(1)=MSR%NSC(IJ) MSR%GOF_H(ITER)=MSR%GOF(1); MSR%NSC_H(ITER)=MSR%NSC(1) MSR%TJ_H(ITER)=MSR%TJ; MSR%RJ_H(ITER)=MSR%RJ !## set correct set of parameters DO I=1,SIZE(PEST%PARAM) PEST%PARAM(I)%ALPHA(1)=PEST%PARAM(I)%LALPHA(IJ) ENDDO IF(PBMAN%IFORMAT.EQ.2.OR.PBMAN%IFORMAT.EQ.3)CALL IPEST_GLM_UPDATEHEADS(DIR,DIRO,IJ) IPEST_GLM_NEXT=.TRUE. END FUNCTION IPEST_GLM_NEXT !#####================================================================= SUBROUTINE IPEST_GLM_UPDATEHEADS(DIR,DIRO,IJ) !#####================================================================= IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: DIR,DIRO INTEGER,INTENT(IN) :: IJ INTEGER :: ILAY,IROW,ICOL TYPE(IDFOBJ) :: SHD TYPE(IDFOBJ),SAVE :: HED LOGICAL :: LEX !## only reset starting head conditions for steady-state solutions IF(PBMAN%ISS.EQ.1)THEN IF(PBMAN%ISTEADY.EQ.0)RETURN ENDIF !## copy results from "winner" as starting heads IF(ALLOCATED(BND))THEN CALL IDFDEALLOCATE(BND,SIZE(BND)); DEALLOCATE(BND) ENDIF ALLOCATE(BND(1)) !## process from nlay to top DO ILAY=PRJNLAY,1,-1 !## only if it exists INQUIRE(FILE=TRIM(DIRO)//'\IPEST_L#'//TRIM(VTOS(IJ))//'\HEAD\HEAD_STEADY-STATE_L'//TRIM(VTOS(ILAY))//'.IDF',EXIST=LEX) IF(LEX)THEN IF(IDFREAD(SHD,TRIM(DIRO)//'\IPEST_L#'//TRIM(VTOS(IJ))//'\HEAD\HEAD_STEADY-STATE_L'//TRIM(VTOS(ILAY))//'.IDF',1))THEN IF(.NOT.ASSOCIATED(HED%X))THEN CALL IDFCOPY(SHD,HED); HED%X=SHD%X ENDIF IF(.NOT.ASSOCIATED(BND(1)%X))CALL IDFCOPY(SHD,BND(1)) !## fill in nodata if needed DO IROW=1,SHD%NROW; DO ICOL=1,SHD%NCOL IF(SHD%X(ICOL,IROW).EQ.SHD%NODATA)SHD%X(ICOL,IROW)=HED%X(ICOL,IROW) ENDDO; ENDDO IF(.NOT.PMANAGER_SAVEMF2005_MOD_U2DREL(TRIM(DIR)//'\MODELINPUT\BAS6\STRT_L'//TRIM(VTOS(ILAY))//'.ARR', & SHD,0,0,1,-1,'F10.3'))THEN ENDIF HED=SHD ENDIF CALL IDFDEALLOCATEX(SHD) ENDIF ENDDO IF(ALLOCATED(BND))THEN CALL IDFDEALLOCATE(BND,SIZE(BND)); DEALLOCATE(BND) ENDIF CALL IDFDEALLOCATEX(SHD) ! CALL IDFDEALLOCATEX(HED) END SUBROUTINE IPEST_GLM_UPDATEHEADS !#####================================================================= SUBROUTINE IPEST_GLM_NEXTGRAD(IPARAM,RTYPE,ISIM) !#####================================================================= IMPLICIT NONE INTEGER,INTENT(IN) :: IPARAM,ISIM CHARACTER(LEN=1),INTENT(IN) :: RTYPE INTEGER :: I REAL(KIND=DP_KIND) :: FCT IF(RTYPE.EQ.'L')THEN DO I=1,SIZE(PEST%PARAM) PEST%PARAM(I)%ALPHA(1)=PEST%PARAM(I)%LALPHA(ISIM) ENDDO ELSEIF(RTYPE.EQ.'P')THEN !## reset all alpha's PEST%PARAM%ALPHA(1)=PEST%PARAM%ALPHA(2) !## adjust all parameters within the same group DO I=1,SIZE(PEST%PARAM) !## skip inactive parameters IF(PEST%PARAM(I)%PACT.EQ.0)CYCLE IF(ABS(PEST%PARAM(I)%PIGROUP).EQ.ABS(PEST%PARAM(IPARAM)%PIGROUP))THEN IF(PEST%PARAM(I)%PLOG.EQ.1)THEN PEST%PARAM(I)%ALPHA(1)=PEST%PARAM(I)%ALPHA(2)+PEST%PARAM(I)%PDELTA FCT=EXP(PEST%PARAM(I)%ALPHA(1)) ELSEIF(PEST%PARAM(I)%PLOG.EQ.2)THEN PEST%PARAM(I)%ALPHA(1)=PEST%PARAM(I)%ALPHA(2)+PEST%PARAM(I)%PDELTA FCT=10.0D0**(PEST%PARAM(I)%ALPHA(1)) ELSE PEST%PARAM(I)%ALPHA(1)=PEST%PARAM(I)%ALPHA(2)+PEST%PARAM(I)%PDELTA FCT=PEST%PARAM(I)%ALPHA(1) ENDIF IF(PEST%PARAM(I)%PIGROUP.GT.0)THEN WRITE(IUPESTOUT,'(A)') 'Adjusting Parameter '//TRIM(PEST%PARAM(I)%PPARAM)// & '['//TRIM(PEST%PARAM(I)%ACRONYM)//']'// & ';ils='//TRIM(VTOS(PEST%PARAM(I)%PILS))// & ';izone='//TRIM(VTOS(PEST%PARAM(I)%PIZONE))// & ';igroup='//TRIM(VTOS(PEST%PARAM(I)%PIGROUP))// & ';factor='//TRIM(VTOS(FCT,'*',1)) ENDIF ENDIF ENDDO ENDIF END SUBROUTINE IPEST_GLM_NEXTGRAD !#####================================================================= SUBROUTINE IPEST_GLM_SAVE_PARAMETERS(IBATCH,DIR,IPARAM,RT,IGRAD) !#####================================================================= IMPLICIT NONE INTEGER,INTENT(IN) :: IGRAD,IPARAM,IBATCH CHARACTER(LEN=*),INTENT(IN) :: DIR,RT INTEGER :: I,J,K,ILAY,ISUB,IROW,ICOL,NCOL,NROW,IPER,ISYS INTEGER,DIMENSION(4) :: IOS INTEGER(KIND=DP_KIND) :: SDATE_PAR,EDATE_PAR,SDATE_MOD,EDATE_MOD REAL(KIND=DP_KIND) :: FCT,F REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:) :: XORG,XADJ,FRAC,XCNT REAL(KIND=DP_KIND),POINTER,DIMENSION(:) :: FADJ INTEGER,POINTER,DIMENSION(:,:) :: IDX CHARACTER(LEN=2),DIMENSION(12) :: IPTYPE CHARACTER(LEN=5),DIMENSION(12) :: FPTYPE CHARACTER(LEN=4),DIMENSION(12) :: SVFILE INTEGER,DIMENSION(12) :: NCLIST,ACLIST,NCILST CHARACTER(LEN=256) :: FNAME,CFOLDER LOGICAL :: LEX IPTYPE= ['KH ' ,'VA' ,'SC' ,'SY' ,'DC' ,'RC' ,'GC', 'RE' , 'QR', 'KD', 'VC', 'HF'] !# modflow6 IF(PBMAN%IFORMAT.EQ.3)THEN FPTYPE=['K' ,'K33' ,'SS' ,'SY' ,'DRN' ,'RIV' ,'GHB' ,'RCH' , 'WEL','NaN','NaN', 'HFB'] SVFILE=['NPF6','NPF6','STO6','STO6','DRN6','RIV6','GHB6','RCH6','WEL6','NaN','NaN','HFB6'] ELSE IF(TOPICS(TKDW)%IACT_MODEL.EQ.1)THEN FPTYPE=['HK' ,'VKA' ,'SF1' ,'SF2' ,'DRN' ,'RIV' ,'GHB' ,'RCH' , 'WEL','TRAN','VCONT','NaN'] SVFILE=['LPF7','LPF7','BCF6','BCF6','DRN7','RIV7','GHB7','RCH7','WEL7','BCF6','BCF6' ,'NaN'] ELSE FPTYPE=['HK' ,'VKA' ,'SF1' ,'SF2' ,'DRN' ,'RIV' ,'GHB' ,'RCH' , 'WEL','TRAN','VCONT','NaN'] SVFILE=['LPF7','LPF7','LPF7','LPF7','DRN7','RIV7','GHB7','RCH7','WEL7','BCF6','BCF6' ,'NaN'] ENDIF ENDIF !## number of columns to be read as reals modflow6 IF(PBMAN%IFORMAT.EQ.3)THEN !## number of reals to be read per package NCLIST=[ 0, 0, 0, 0, 2, 3, 2 , 1 , 1 , 0, 0, 1] !## number of integer to be read per package NCILST=[ 0, 0, 0, 0, 3, 3, 3 , 3 , 3 , 0, 0, 6] IF(PBMAN%DDRN.NE.0.0D0)NCLIST(5)=3 !## number of columns to be read as reals imod-wq ELSE !## number of reals to be read per package NCLIST=[ 0, 0, 0, 0, 2, 4, 3 , 1 , 1 , 0, 0, 0] !## number of integer to be read per package NCILST=[ 0, 0, 0, 0, 3, 3, 3 , 3 , 3 , 0, 0, 6] ENDIF !## parameter to be adjusted (from the reals in clist) ACLIST= [ 0, 0, 0, 0, 2, 2, 2 , 1 , 1 , 0, 0, 1] !## modify parameters per submodel DO ISUB=1,PBMAN%NSUBMODEL IF(PBMAN%IFORMAT.EQ.3)CFOLDER='\GWF_'//TRIM(VTOS(ISUB))//'\' IF(PBMAN%IFORMAT.EQ.6)CFOLDER='\' !## read the zones per submodel CALL IPEST_GLM_READ_ZONES(TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT',NCOL,NROW) ALLOCATE(FRAC(NCOL,NROW)) DO I=1,SIZE(IPTYPE) !## read from list of number SELECT CASE (IPTYPE(I)) CASE ('DC','RC','GC','QR','RE','HF') !## see whether this parameter is optimized DO J=1,SIZE(PEST%PARAM) !## check whether the parameter is available IF(PEST%PARAM(J)%PPARAM.EQ.IPTYPE(I))EXIT ENDDO; IF(J.GT.SIZE(PEST%PARAM))CYCLE ISYS=0 DO ISYS=ISYS+1 !## mf6 IF(PBMAN%IFORMAT.EQ.3)THEN IF(IPTYPE(I).EQ.'HF')THEN FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\'//RT//'#'//TRIM(VTOS(IPARAM)) ELSE FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\SYS'//TRIM(VTOS(ISYS)); IF(.NOT.IOSDIREXISTS(FNAME))EXIT FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\SYS'//TRIM(VTOS(ISYS))//'\'//RT//'#'//TRIM(VTOS(IPARAM)) ENDIF CALL UTL_CREATEDIR(FNAME) !## imod-wq ELSE FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\'//RT//'#'//TRIM(VTOS(IPARAM)) CALL UTL_CREATEDIR(FNAME) ENDIF DO IPER=1,PRJNPER !## mf6 IF(PBMAN%IFORMAT.EQ.3)THEN !## read original parameter values IF(IPTYPE(I).EQ.'HF')THEN FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\'//TRIM(FPTYPE(I))//'_T'//TRIM(VTOS(IPER))//'.ARR' ELSE FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\SYS'//TRIM(VTOS(ISYS))//'\'//TRIM(FPTYPE(I))//'_T'//TRIM(VTOS(IPER))//'.ARR' ENDIF !## imod-wq ELSE !## read original parameter values FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\'//TRIM(FPTYPE(I))//'_T'//TRIM(VTOS(IPER))//'.ARR' ENDIF !## skip it for this period INQUIRE(FILE=FNAME,EXIST=LEX); IF(.NOT.LEX)CYCLE IF(PBMAN%IFORMAT.EQ.6.AND.IPTYPE(I).EQ.'RE')THEN IF(.NOT.IPEST_GLM_READ_ARRFILE(FNAME,XORG))RETURN IF(.NOT.ASSOCIATED(XADJ))ALLOCATE(XADJ(SIZE(XORG,1),SIZE(XORG,2))) IF(.NOT.ASSOCIATED(XCNT))ALLOCATE(XCNT(SIZE(XORG,1),SIZE(XORG,2))) XADJ=XORG; XCNT=0.0D0 ELSE IF(.NOT.IPEST_GLM_READ_LSTFILE(FNAME,NCLIST(I),XORG,NCILST(I),IDX))RETURN IF(.NOT.ASSOCIATED(XADJ))ALLOCATE(XADJ(SIZE(XORG,1),SIZE(XORG,2))) XADJ=XORG; ALLOCATE(FADJ(SIZE(IDX,2))); FADJ=0.0D0 ENDIF !## process parameter values DO J=1,SIZE(PEST%PARAM) !## check whether the parameter is available - leave this out otherwise parameters that are kept constant are excluded IF(PEST%PARAM(J)%PPARAM.NE.IPTYPE(I))CYCLE !## check date if needed IF(PEST%PARAM(J)%SDATE.NE.''.AND.PEST%PARAM(J)%EDATE.NE.'')THEN IF(TRIM(SIM(IPER)%CDATE).EQ.'STEADY-STATE')THEN IF(TRIM(PEST%PARAM(J)%SDATE).EQ.'STEADY-STATE'.AND.TRIM(PEST%PARAM(J)%EDATE).EQ.'STEADY-STATE')THEN IF(TRIM(SIM(IPER)%CDATE).NE.'STEADY-STATE')CYCLE ELSE !## skip this parameter CYCLE ENDIF ELSE READ(PEST%PARAM(J)%SDATE,*,IOSTAT=IOS(1)) SDATE_PAR READ(PEST%PARAM(J)%EDATE,*,IOSTAT=IOS(2)) EDATE_PAR SDATE_MOD=YMDHMSTOITIME(SIM(IPER )%IYR,SIM(IPER )%IMH,SIM(IPER )%IDY,SIM(IPER )%IHR,SIM(IPER )%IMT,SIM(IPER )%ISC) EDATE_MOD=YMDHMSTOITIME(SIM(IPER+1)%IYR,SIM(IPER+1)%IMH,SIM(IPER+1)%IDY,SIM(IPER+1)%IHR,SIM(IPER+1)%IMT,SIM(IPER+1)%ISC) !## ready to evaluate IF(SUM(IOS).EQ.0)THEN !## make sure all are in same dimensions SDATE_MOD=UTL_COMPLETEDATE(SDATE_MOD); EDATE_MOD=UTL_COMPLETEDATE(EDATE_MOD) SDATE_PAR=UTL_COMPLETEDATE(SDATE_PAR); EDATE_PAR=UTL_COMPLETEDATE(EDATE_PAR) IF(EDATE_PAR.LE.SDATE_MOD)CYCLE; IF(SDATE_PAR.GT.EDATE_MOD)CYCLE ENDIF ENDIF ENDIF !## found correct parameter in time : apply factor FCT=PEST%PARAM(J)%ALPHA(1) IF(PEST%PARAM(J)%PLOG.EQ.1)FCT=EXP(FCT) IF(PEST%PARAM(J)%PLOG.EQ.2)FCT=10.0D0**FCT IF(PEST%PARAM(J)%ZTYPE.EQ.0)THEN !## fill array with location of a zone FRAC=0.0D0 DO K=1,PEST%PARAM(J)%NODES IROW=PEST%PARAM(J)%IROW(K) ICOL=PEST%PARAM(J)%ICOL(K) FRAC(ICOL,IROW)=PEST%PARAM(J)%F(K) ENDDO IF(PBMAN%IFORMAT.EQ.6.AND.IPTYPE(I).EQ.'RE')THEN DO K=1,PEST%PARAM(J)%NODES IROW=PEST%PARAM(J)%IROW(K); ICOL=PEST%PARAM(J)%ICOL(K) F =FRAC(ICOL,IROW) IF(XCNT(ICOL,IROW).EQ.0.0D0)THEN XADJ(ICOL,IROW)=XORG(ICOL,IROW)*F*FCT ELSE XADJ(ICOL,IROW)=XADJ(ICOL,IROW)+XORG(ICOL,IROW)*F*FCT ENDIF XCNT(ICOL,IROW)=XCNT(ICOL,IROW)+F ENDDO ELSE IF(IPTYPE(I).EQ.'HF')THEN DO K=1,SIZE(IDX,2) !## skip wrong system number IF(IDX(4,K).NE.PEST%PARAM(J)%PILS)CYCLE IROW=IDX(2,K); ICOL=IDX(3,K); F=FRAC(ICOL,IROW); IF(F.LE.0.0D0)CYCLE IROW=IDX(5,K); ICOL=IDX(6,K); F=FRAC(ICOL,IROW); IF(F.LE.0.0D0)CYCLE IF(FADJ(K).EQ.0.0D0)THEN XADJ(ACLIST(I),K)=XORG(ACLIST(I),K)*F*FCT ELSE !## modify current package for this zone XADJ(ACLIST(I),K)=XADJ(ACLIST(I),K)+XORG(ACLIST(I),K)*F*FCT ENDIF FADJ(K)=FADJ(K)+F ENDDO ELSE DO K=1,SIZE(IDX,2) !## skip wrong system number IF(IDX(4,K).NE.PEST%PARAM(J)%PILS)CYCLE IROW=IDX(2,K); ICOL=IDX(3,K); F=FRAC(ICOL,IROW); IF(F.LE.0.0D0)CYCLE IF(FADJ(K).EQ.0.0D0)THEN XADJ(ACLIST(I),K)=XORG(ACLIST(I),K)*F*FCT ELSE !## modify current package for this zone XADJ(ACLIST(I),K)=XADJ(ACLIST(I),K)+XORG(ACLIST(I),K)*F*FCT ENDIF FADJ(K)=FADJ(K)+F ENDDO ENDIF ENDIF ENDIF ENDDO IF(PBMAN%IFORMAT.EQ.6.AND.IPTYPE(I).EQ.'RE')THEN !## check whether adjustments are okay DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL IF(XCNT(ICOL,IROW).GT.1.01D0.AND.XORG(ICOL,IROW).NE.XADJ(ICOL,IROW))THEN WRITE(*,'(A)') '>>> Parameter: '//IPTYPE(I)//' <<<' WRITE(*,'(A,F15.7,A,2I5,A)') 'Too much altered (',XCNT(ICOL,IROW),') in a single cell [rc]: (',IROW,ICOL,')' WRITE(*,'(2(A,F15.7))') 'Original Value= ',XORG(ICOL,IROW),' ; Modified Value= ',XADJ(ICOL,IROW) PAUSE ENDIF ENDDO; ENDDO ELSE DO K=1,SIZE(IDX,2) IROW=IDX(2,K); ICOL=IDX(3,K) IF(FADJ(K).GT.1.01D0.AND.XORG(ACLIST(I),K).NE.XADJ(ACLIST(I),K))THEN WRITE(*,'(A)') '>>> Parameter :'//IPTYPE(I)//' for system '//TRIM(VTOS(K))//' <<<' WRITE(*,'(A,F15.7,A,I10,A,3I5,A)') 'Too much altered (',FADJ(K),') in a single entry (',K,') for cell [lrc]: ',IDX(4,K),IROW,ICOL,')' WRITE(*,'(2(A,F15.7))') 'Original Value= ',XORG(ACLIST(I),K),' ; Modified Value= ',XADJ(ACLIST(I),K) PAUSE ENDIF ENDDO ENDIF !## add if needed pilot point interpolation IF(IPEST_GLM_COMPUTE_PILOTPOINTS(IPTYPE(I),ILAY,IBATCH))THEN IF(PBMAN%IFORMAT.EQ.6.AND.IPTYPE(I).EQ.'RE')THEN DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL !## copy interpolated data into current adjusted dataset IF(PRJIDF%X(ICOL,IROW).NE.PRJIDF%NODATA)XADJ(ICOL,IROW)=XADJ(ICOL,IROW)*PRJIDF%X(ICOL,IROW) ENDDO; ENDDO ELSE IF(IPTYPE(I).EQ.'HF')THEN DO K=1,SIZE(IDX,2) !## skip wrong system number IF(IDX(7,K).NE.PEST%PARAM(J)%PILS)CYCLE !## inside from node IROW=IDX(2,K); ICOL=IDX(3,K); IF(PRJIDF%X(ICOL,IROW).EQ.PRJIDF%NODATA)CYCLE !## inside to node IROW=IDX(5,K); ICOL=IDX(6,K); IF(PRJIDF%X(ICOL,IROW).EQ.PRJIDF%NODATA)CYCLE !## copy interpolated data into current adjusted dataset XADJ(ACLIST(I),K)=XADJ(ACLIST(I),K)*PRJIDF%X(ICOL,IROW) ENDDO ELSE DO K=1,SIZE(IDX,2) !## skip wrong system number IF(IDX(4,K).NE.PEST%PARAM(J)%PILS)CYCLE IROW=IDX(2,K); ICOL=IDX(3,K) !## copy interpolated data into current adjusted dataset IF(PRJIDF%X(ICOL,IROW).NE.PRJIDF%NODATA)XADJ(ACLIST(I),K)=XADJ(ACLIST(I),K)*PRJIDF%X(ICOL,IROW) ENDDO ENDIF ENDIF ENDIF !## mf6 IF(PBMAN%IFORMAT.EQ.3)THEN !## save adjusted package-values IF(IPTYPE(I).EQ.'HF')THEN FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\'//RT//'#'//TRIM(VTOS(IPARAM))//'\'// & TRIM(FPTYPE(I))//'_T'//TRIM(VTOS(IPER))//'_'//RT//'#'//TRIM(VTOS(IPARAM))//'.ARR' ELSE FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\SYS'//TRIM(VTOS(ISYS))//'\'//RT//'#'//TRIM(VTOS(IPARAM))//'\'// & TRIM(FPTYPE(I))//'_T'//TRIM(VTOS(IPER))//'_'//RT//'#'//TRIM(VTOS(IPARAM))//'.ARR' ENDIF !## imod-wq ELSE !## save adjusted package-values FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\'//RT//'#'//TRIM(VTOS(IPARAM))//'\'// & TRIM(FPTYPE(I))//'_T'//TRIM(VTOS(IPER))//'_'//RT//'#'//TRIM(VTOS(IPARAM))//'.ARR' ENDIF IF(PBMAN%IFORMAT.EQ.6.AND.IPTYPE(I).EQ.'RE')THEN IF(.NOT.IPEST_GLM_WRITE_ARRFILE(FNAME,XADJ))RETURN ELSE IF(.NOT.IPEST_GLM_WRITELSTFILE(FNAME,XADJ,IDX))RETURN ENDIF !## deallocate memory DEALLOCATE(XORG,XADJ) IF(ASSOCIATED(IDX))DEALLOCATE(IDX) IF(ASSOCIATED(FADJ))DEALLOCATE(FADJ) IF(ASSOCIATED(XCNT))DEALLOCATE(XCNT) ENDDO IF(PBMAN%IFORMAT.EQ.6)EXIT IF(PBMAN%IFORMAT.EQ.3.AND.IPTYPE(I).EQ.'HF')EXIT ENDDO CASE ('KH','VA','SC','SY','KD','VC') !## see whether this parameter is optimized DO J=1,SIZE(PEST%PARAM) !## check whether the parameter is available - leave this out otherwise parameters that are kept constant are excluded IF(PEST%PARAM(J)%PPARAM.EQ.IPTYPE(I))EXIT ENDDO; IF(J.GT.SIZE(PEST%PARAM))CYCLE FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\'//RT//'#'//TRIM(VTOS(IPARAM)) CALL UTL_CREATEDIR(FNAME) DO ILAY=1,PRJNLAY DO J=1,SIZE(PEST%PARAM) !## check whether the parameter is available - leave this out otherwise parameters that are kept constant are excluded IF(PEST%PARAM(J)%PPARAM.EQ.IPTYPE(I))THEN IF(PEST%PARAM(J)%PILS.EQ.ILAY)EXIT ENDIF ENDDO !## skip this layer as it is not in the list optimized parameters IF(J.GT.SIZE(PEST%PARAM))CYCLE !## read original parameter values FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\'//TRIM(FPTYPE(I))//'_L'//TRIM(VTOS(ILAY))//'.ARR' INQUIRE(FILE=FNAME,EXIST=LEX) IF(.NOT.LEX)THEN WRITE(*,'(/A)') 'Cannot find '//TRIM(FNAME) WRITE(*,'(A/)') 'Probably you are trying to calibrate a parameter not-existing in the current model' STOP ENDIF IF(.NOT.IPEST_GLM_READ_ARRFILE(FNAME,XORG))RETURN IF(.NOT.ASSOCIATED(XADJ))ALLOCATE(XADJ(SIZE(XORG,1),SIZE(XORG,2))) IF(.NOT.ASSOCIATED(XCNT))ALLOCATE(XCNT(SIZE(XORG,1),SIZE(XORG,2))) IF(IPTYPE(I).EQ.'VC')THEN DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL IF(XORG(ICOL,IROW).NE.0.0D0)XORG(ICOL,IROW)=1.0D0/XORG(ICOL,IROW) ENDDO; ENDDO ENDIF XADJ=XORG; XCNT=0.0D0 !## process zonation DO J=1,SIZE(PEST%PARAM) !## check whether the parameter is available - leave this out otherwise parameters that are kept constant are excluded IF(PEST%PARAM(J)%PPARAM.NE.IPTYPE(I))CYCLE IF(PEST%PARAM(J)%PILS.NE.ILAY)CYCLE !## found correct parameter - apply factor FCT=PEST%PARAM(J)%ALPHA(1) IF(PEST%PARAM(J)%PLOG.EQ.1)FCT=EXP(FCT) IF(PEST%PARAM(J)%PLOG.EQ.2)FCT=10.0D0**FCT IF(PEST%PARAM(J)%ZTYPE.EQ.0)THEN DO K=1,PEST%PARAM(J)%NODES IROW=PEST%PARAM(J)%IROW(K) ICOL=PEST%PARAM(J)%ICOL(K) F =PEST%PARAM(J)%F(K) IF(XCNT(ICOL,IROW).EQ.0.0D0)THEN XADJ(ICOL,IROW)=XORG(ICOL,IROW)*F*FCT ELSE XADJ(ICOL,IROW)=XADJ(ICOL,IROW)+XORG(ICOL,IROW)*F*FCT ENDIF XCNT(ICOL,IROW)=XCNT(ICOL,IROW)+F ENDDO ENDIF ENDDO !## check whether adjustments are okay DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL IF(XCNT(ICOL,IROW).GT.1.01D0.AND.XORG(ICOL,IROW).NE.XADJ(ICOL,IROW))THEN WRITE(*,'(A)') '>>> Parameter :'//IPTYPE(I)//' model layer '//TRIM(VTOS(ILAY))//' <<<' WRITE(*,'(A,F15.7,A,3I5,A)') 'Too much altered (',XCNT(ICOL,IROW),') in a single cell [lrc]: (',ILAY,IROW,ICOL,')' WRITE(*,'(2(A,F15.7))') 'Original Value= ',XORG(ICOL,IROW),' ; Modified Value= ',XADJ(ICOL,IROW) PAUSE ENDIF ENDDO; ENDDO !## add if needed pilot point interpolation IF(IPEST_GLM_COMPUTE_PILOTPOINTS(IPTYPE(I),ILAY,IBATCH))THEN DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL !## copy interpolated data into current adjusted dataset IF(PRJIDF%X(ICOL,IROW).NE.PRJIDF%NODATA)XADJ(ICOL,IROW)=XADJ(ICOL,IROW)*PRJIDF%X(ICOL,IROW) ENDDO; ENDDO ENDIF IF(IPTYPE(I).EQ.'VC')THEN DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL IF(XORG(ICOL,IROW).NE.0.0D0)XORG(ICOL,IROW)=1.0D0/XORG(ICOL,IROW) ENDDO; ENDDO ENDIF !## save adjusted values FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//SVFILE(I)//'\'//RT//'#'//TRIM(VTOS(IPARAM))//'\'// & TRIM(FPTYPE(I))//'_L'//TRIM(VTOS(ILAY))//'_'//RT//'#'//TRIM(VTOS(IPARAM))//'.ARR' IF(.NOT.IPEST_GLM_WRITE_ARRFILE(FNAME,XADJ))RETURN !## deallocate memory DEALLOCATE(XORG,XADJ,XCNT) ENDDO END SELECT ENDDO !## clean memory DO I=1,SIZE(PEST%PARAM) IF(ASSOCIATED(PEST%PARAM(I)%IROW))DEALLOCATE(PEST%PARAM(I)%IROW) IF(ASSOCIATED(PEST%PARAM(I)%ICOL))DEALLOCATE(PEST%PARAM(I)%ICOL) IF(ASSOCIATED(PEST%PARAM(I)%F)) DEALLOCATE(PEST%PARAM(I)%F) IF(ASSOCIATED(PEST%PARAM(I)%XY)) DEALLOCATE(PEST%PARAM(I)%XY) ENDDO ENDDO DEALLOCATE(FRAC) DO I=1,SIZE(PEST%PARAM) IF(RT.EQ.'L')THEN PEST%PARAM(I)%LALPHA(IGRAD)=PEST%PARAM(I)%ALPHA(1) ELSEIF(RT.EQ.'P')THEN PEST%PARAM(I)%GALPHA(IGRAD)=PEST%PARAM(I)%ALPHA(1) ENDIF ENDDO END SUBROUTINE IPEST_GLM_SAVE_PARAMETERS !###==================================================================== LOGICAL FUNCTION IPEST_GLM_COMPUTE_PILOTPOINTS(IPTYPE,ILAY,IBATCH) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: IPTYPE !,RT INTEGER,INTENT(IN) :: ILAY,IBATCH REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: XD,YD,ZD,WD,PD REAL(KIND=DP_KIND) :: FCT INTEGER :: J,K,ND,ILOG,IROW,ICOL !,PE_KTYPE_COPY REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:) :: XBLNK IPEST_GLM_COMPUTE_PILOTPOINTS=.FALSE. !## process pilot points if needed ND=0; DO J=1,SIZE(PEST%PARAM) IF(PEST%PARAM(J)%PPARAM.NE.IPTYPE)CYCLE IF(PEST%PARAM(J)%PILS.NE.ILAY)CYCLE IF(PEST%PARAM(J)%ZTYPE.EQ.1)ND=ND+PEST%PARAM(J)%NODES ENDDO !## no pilotpoints IF(ND.EQ.0)RETURN ALLOCATE(XD(ND),YD(ND),ZD(ND),WD(ND),PD(ND)) ND=0; DO J=1,SIZE(PEST%PARAM) IF(PEST%PARAM(J)%PPARAM.NE.IPTYPE)CYCLE IF(PEST%PARAM(J)%PILS.NE.ILAY)CYCLE IF(PEST%PARAM(J)%ZTYPE.EQ.0)CYCLE FCT=PEST%PARAM(J)%ALPHA(1) DO K=1,PEST%PARAM(J)%NODES ND=ND+1 !## kriging on log-values XD(ND)=PEST%PARAM(J)%XY(K,1) YD(ND)=PEST%PARAM(J)%XY(K,2) ZD(ND)=0.0D0 !## set log-type (be aware not possible to use mix of them) ILOG=PEST%PARAM(J)%PLOG PD(ND)=FCT WD(ND)=1.0D0 ENDDO ENDDO ! PE_KTYPE_COPY=PEST%PE_KTYPE ! IF(TRIM(RT).EQ.'P')PEST%PE_KTYPE=6 !## number of pilotpoint to be used in interpolation MAXPNT=10; IF(ABS(PEST%PE_KTYPE).EQ.6)MAXPNT=1 IBLANKOUT=0; IF(PEST%PE_KTYPE.LE.0)IBLANKOUT=1 KTYPE=-1*ABS(PEST%PE_KTYPE) ! !## ordinary kriging - spherical ! IF(KTYPE.EQ.3)KTYPE=-2 !3 ! IF(KTYPE.EQ.2)KTYPE=-2 !3 !2 ! !## simple kriging ! IF(KTYPE.EQ.1)KTYPE= 2 !3 !2 RANGE=PEST%PE_KRANGE PNTSEARCH=1 SILL =100.0D0 NUGGET=0.0D0 IQUADRANT=0 !## copy what has been adjusted so far IF(.NOT.IDFALLOCATEX(PRJIDF))RETURN !## read blank-out region IF(IBLANKOUT.EQ.1)THEN IF(.NOT.IPEST_GLM_READ_ARRFILE(PEST%PPBNDIDF,XBLNK))RETURN PRJIDF%X=0.0D0 !## interpolate for areas equal to nodata only DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL IF(XBLNK(ICOL,IROW).NE.0.0)PRJIDF%X(ICOL,IROW)=PRJIDF%NODATA ENDDO; ENDDO ELSE !## interpolate all PRJIDF%X=PRJIDF%NODATA; ALLOCATE(XBLNK(PRJIDF%NCOL,PRJIDF%NROW)); XBLNK=1.0D0 ENDIF CALL KRIGING_MAIN(SIZE(XD),XD,YD,ZD,PD,WD,PRJIDF,IBATCH,0.0D0) !## back-transformation DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL !## copy interpolated data into current adjusted dataset IF(XBLNK(ICOL,IROW).NE.0.0)THEN IF(ILOG.EQ.1)PRJIDF%X(ICOL,IROW)=EXP(PRJIDF%X(ICOL,IROW)) IF(ILOG.EQ.2)PRJIDF%X(ICOL,IROW)=10.0D0**PRJIDF%X(ICOL,IROW) ELSE PRJIDF%X(ICOL,IROW)=PRJIDF%NODATA ENDIF ENDDO; ENDDO IF(ASSOCIATED(XBLNK))DEALLOCATE(XBLNK) ! PEST%PE_KTYPE=PE_KTYPE_COPY IPEST_GLM_COMPUTE_PILOTPOINTS=.TRUE. END FUNCTION IPEST_GLM_COMPUTE_PILOTPOINTS !###==================================================================== LOGICAL FUNCTION IPEST_GLM_READ_LSTFILE(FNAME,NC,X,NCI,IDX) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: FNAME INTEGER,INTENT(IN) :: NC,NCI REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:) :: X INTEGER,POINTER,DIMENSION(:,:) :: IDX INTEGER :: IU,IOS,I,J,N CHARACTER(LEN=256) :: LINE IPEST_GLM_READ_LSTFILE=.FALSE. IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='OLD',FORM='FORMATTED',ACTION='READ',IOSTAT=IOS) IF(IOS.NE.0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot open file: '//CHAR(13)// & '['//TRIM(FNAME)//']'//CHAR(13)//'for reading','Error') RETURN ENDIF DO I=1,2 N=0; DO READ(IU,'(A256)',IOSTAT=IOS) LINE IF(IOS.NE.0)EXIT; IF(INDEX(LINE,'DIMENSIONS').GT.0)EXIT; N=N+1 IF(I.EQ.2)THEN READ(LINE,*) (IDX(J,N),J=1,NCI),(X(J,N),J=1,NC),IDX(NCI+1,N) ENDIF ENDDO IF(I.EQ.1)THEN ALLOCATE(IDX(NCI+1,N),X(NC,N)); IDX=0; X=0.0D0 ENDIF REWIND(IU) ENDDO CLOSE(IU) IPEST_GLM_READ_LSTFILE=.TRUE. END FUNCTION IPEST_GLM_READ_LSTFILE !###==================================================================== LOGICAL FUNCTION IPEST_GLM_WRITELSTFILE(FNAME,X,IDX) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: FNAME REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:) :: X INTEGER,POINTER,DIMENSION(:,:) :: IDX INTEGER :: IU,IOS,I,J,NCI,NCR CHARACTER(LEN=52) :: FRM IPEST_GLM_WRITELSTFILE=.FALSE. NCI=SIZE(IDX,1) NCR=SIZE(X ,1) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='REPLACE',FORM='FORMATTED',ACTION='WRITE',IOSTAT=IOS) IF(IOS.NE.0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot open file: '//CHAR(13)// & '['//TRIM(FNAME)//']'//CHAR(13)//'for writing','Error') RETURN ENDIF WRITE(FRM,'(A10,I2.2,A14)') '('//TRIM(VTOS(NCI-1))//'(I5,1X),',NCR,'(F15.3,1X),I5)' DO I=1,SIZE(IDX,2) WRITE(IU,FRM) (IDX(J,I),J=1,NCI-1),(X(J,I),J=1,NCR),IDX(NCI,I) ENDDO CLOSE(IU) IPEST_GLM_WRITELSTFILE=.TRUE. END FUNCTION IPEST_GLM_WRITELSTFILE !###==================================================================== LOGICAL FUNCTION IPEST_GLM_READ_ARRFILE(FNAME,X,IDBL) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: FNAME INTEGER,INTENT(OUT),OPTIONAL :: IDBL REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:) :: X INTEGER :: IU,NLINE,IOS,NCOL,NROW,IROW,ICOL CHARACTER(LEN=52) :: LINE IPEST_GLM_READ_ARRFILE=.FALSE. IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='OLD',FORM='FORMATTED',ACTION='READ',IOSTAT=IOS) IF(IOS.NE.0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot open file: '//CHAR(13)// & '['//TRIM(FNAME)//']'//CHAR(13)//'for reading','Error') RETURN ENDIF !## look for network dimensions NLINE=0; DO READ(IU,'(A52)') LINE IF(INDEX(LINE,'DIMENSIONS').GT.0)EXIT NLINE=NLINE+1 ENDDO READ(IU,'(A52)') LINE; READ(LINE(2:),*) NCOL READ(IU,'(A52)') LINE; READ(LINE(2:),*) NROW ALLOCATE(X(NCOL,NROW)) REWIND(IU) DO IROW=1,NROW READ(IU,*) (X(ICOL,IROW),ICOL=1,NCOL) ENDDO !## try to read double/single precision IF(PRESENT(IDBL))THEN !## default single precision READ(IU,*,IOSTAT=IOS) IDBL; IF(IOS.NE.0)IDBL=4 ENDIF CLOSE(IU) IPEST_GLM_READ_ARRFILE=.TRUE. END FUNCTION IPEST_GLM_READ_ARRFILE !###==================================================================== LOGICAL FUNCTION IPEST_GLM_WRITE_ARRFILE(FNAME,X) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: FNAME REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:) :: X INTEGER :: IU,IOS,NCOL,NROW,IROW REAL(KIND=DP_KIND) :: NODATA IPEST_GLM_WRITE_ARRFILE=.FALSE. IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='REPLACE',FORM='FORMATTED',ACTION='WRITE',IOSTAT=IOS) IF(IOS.NE.0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot open file: '//CHAR(13)// & '['//TRIM(FNAME)//']'//CHAR(13)//'for writing','Error') RETURN ENDIF NCOL=SIZE(X,1) NROW=SIZE(X,2) NODATA=HUGE(1.0) DO IROW=1,NROW CALL IDFWRITEFREE_ROW(IU,X(:,IROW),NCOL,NODATA,0,'G15.7') ENDDO CALL IDFWRITEFREE_HEADER(IU,PRJIDF) CLOSE(IU) IPEST_GLM_WRITE_ARRFILE=.TRUE. END FUNCTION IPEST_GLM_WRITE_ARRFILE !###==================================================================== INTEGER FUNCTION IPEST_GLM_GRADIENT(IBATCH,ITER,LAMBDA,GAMMA) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: ITER,IBATCH REAL(KIND=DP_KIND),INTENT(INOUT) :: LAMBDA REAL(KIND=DP_KIND),INTENT(IN) :: GAMMA REAL(KIND=DP_KIND) :: DJ1,DJ2,DP,WW,TS,DF1,EIGWTHRESHOLD,W,H1,H2,MARQUARDT,P1,P2,PMIN,PMAX,WF,SS,F INTEGER :: I,J,K,NP,IP1,NE,ISING,ILAMBDA,IBND,IPARAM,NN INTEGER,ALLOCATABLE,DIMENSION(:) :: INDX LOGICAL :: LJQJ REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: BJQJ,JQJ REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: U REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: JS,P,PT REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: N,RU REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: S REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: EIGV,B,M REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: EIGW,JQR LOGICAL :: LSVD IPEST_GLM_GRADIENT=0 SELECT CASE (PEST%PE_SCALING) CASE (0,1); LSVD=.FALSE. CASE (2,3); LSVD=.TRUE. END SELECT !## sensitivity NP=0; DO I=1,SIZE(PEST%PARAM); IF(ABS(PEST%PARAM(I)%PACT).NE.1)CYCLE; NP=NP+1; ENDDO IF(.NOT.ALLOCATED(S)) ALLOCATE(S(NP)) WRITE(*,'(A)') ' >>> Computing Sensitivities on '//TRIM(VTOS(UTL_OMP_GET_MAX_THREADS()))//' threads' S=0.0D0; IPARAM=0 DO I=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(I)%PACT).NE.1)CYCLE DF1=PEST%PARAM(I)%PDELTA IPARAM=IPARAM+1 !$OMP PARALLEL PRIVATE(J,H1,H2,W,TS) SHARED (MSR,DF1) TS=0.0D0 SS=0.0D0 !$OMP DO DO J=1,MSR%NOBS W=MSR%W(J); H1=MSR%HG(IPARAM,J); H2=MSR%HL(0,J) TS=TS+(W*((H1-H2)/DF1)) ENDDO !$OMP END DO !$OMP CRITICAL SS=SS+TS !$OMP END CRITICAL !$OMP END PARALLEL S(IPARAM)=SS F=REAL(IPARAM,8)/REAL(NP,8); WRITE(6,'(A,F10.2,A)') '+ >>> Progress ',F*100.0D0,'%' ENDDO WRITE(*,'(A)') ' >>> Computing Sensitivities on '//TRIM(VTOS(UTL_OMP_GET_MAX_THREADS()))//' threads <<<' DO I=1,NP; S(I)=S(I)/DBLE(MSR%NOBS); ENDDO ! !## remove small values 0.1% of biggest sensitivity ! MAXS=MAXVAL(ABS(S)); MAXS=MAXS/1000.0D0 ! DO I=1,NP; IF(ABS(S(I)).LE.MAXS)S(I)=0.0D0; ENDDO TS=SUM(ABS(S)); IPARAM=0; DO I=1,SIZE(PEST%PARAM) !## skip inactive parameters IF(ABS(PEST%PARAM(I)%PACT).NE.1)CYCLE; IPARAM=IPARAM+1; IF(TS.NE.0.0D0)S(IPARAM)=S(IPARAM)/TS ENDDO; S=ABS(S)*100.0D0 WRITE(IUPESTSENSITIVITY,'(I9,6X,99999F15.3)') ITER,(S(I),I=1,NP) IF(PEST%PE_MXITER.EQ.0)RETURN IPARAM=0; DO I=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(I)%PACT).NE.1)CYCLE; IPARAM=IPARAM+1 !## remove sensitivities of zero IF(S(IPARAM).LE.0.0D0)PEST%PARAM(I)%PACT=-1 ENDDO !##=================== !## write statistics of ALL parameters prior to the update !##=================== NP=0; DO I=1,SIZE(PEST%PARAM); IF(PEST%PARAM(I)%PACT.EQ.1)NP=NP+1; ENDDO IF(ALLOCATED(JQJ))DEALLOCATE(JQJ); ALLOCATE(JQJ (NP,NP)) IF(ALLOCATED(BJQJ))DEALLOCATE(BJQJ); ALLOCATE(BJQJ(NP,NP)) IF(ALLOCATED(EIGW))DEALLOCATE(EIGW); ALLOCATE(EIGW(NP )) IF(ALLOCATED(EIGV))DEALLOCATE(EIGV); ALLOCATE(EIGV(NP,NP)) !## write statistics (covariance/correlation) BJQJ=0.0D0; LJQJ=.TRUE.; CALL IPEST_GLM_JQJ(IBATCH,MARQUARDT,BJQJ,JQJ,NP,.TRUE.,LJQJ) !## reset insensitive parameters to log them in a list IPARAM=0; DO I=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(I)%PACT).NE.1)CYCLE; IPARAM=IPARAM+1 !## remove sensitivities of zero IF(S(IPARAM).LE.0.0D0)PEST%PARAM(I)%PACT=1 ENDDO !## reset parameters - alpha(2)=previous alpha DO I=1,SIZE(PEST%PARAM); PEST%PARAM(I)%ALPHA(1)=PEST%PARAM(I)%ALPHA(2); ENDDO !## "freeze"-insensitive parameters according to a specified sens.-criterion K=0; IPARAM=0; DO I=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(I)%PACT).NE.1)CYCLE; IPARAM=IPARAM+1 IF(S(IPARAM).LE.ABS(PEST%PE_SENS))THEN PEST%PARAM(I)%PACT=-1; K=K+1 IF(K.EQ.1)THEN WRITE(IUPESTOUT,'(/A)') 'List of Insensitive Parameter (Sensitivity <= '//TRIM(VTOS(ABS(PEST%PE_SENS),'F',7))//' %)' IF(PEST%PE_SENS.GE.0.0D0)WRITE(IUPESTOUT,'(A/)') ' turned off TEMPORARILY' IF(PEST%PE_SENS.LT.0.0D0)WRITE(IUPESTOUT,'(A/)') ' turned off PERMENANTLY' ENDIF IF(PEST%PE_SENS.LT.0.0D0)PEST%PARAM(I)%INSENS=1 WRITE(IUPESTOUT,'(A15,F15.7,A2)') PEST%PARAM(I)%ACRONYM,S(IPARAM),' %' ENDIF ENDDO WRITE(IUPESTOUT,'(/70A1)') ('=',I=1,70); WRITE(IUPESTOUT,'(A)') 'Start Lambda Testing'; WRITE(IUPESTOUT,'(70A1/)') ('=',I=1,70) !## set boundaries for parameters DO I=1,SIZE(PEST%PARAM); CALL IPEST_GLM_GETBOUNDARY(I,IBND,P1,P2,PMIN,PMAX); PEST%PARAM(I)%IBND=IBND; ENDDO !## compute update vector for lambdas ILAMBDA=0 DO ILAMBDA=ILAMBDA+1 !## finished IF(ILAMBDA.GT.PBMAN%NLAMBDASEARCH)EXIT !## store current size of matrices NN=0; IF(ALLOCATED(JQJ))NN=SIZE(JQR) !## initiate number of parameters NP=0; DO I=1,SIZE(PEST%PARAM); IF(PEST%PARAM(I)%PACT.EQ.1)NP=NP+1; ENDDO IF(NP.EQ.0)THEN; CALL IPEST_GLM_ERROR(IBATCH,'All (remaining) parameters are insensitive, process stopped!'); IPEST_GLM_GRADIENT=-1; RETURN; ENDIF !## reset matrices IF(NP.NE.NN)THEN LJQJ=.TRUE. IF(ALLOCATED(JQJ))DEALLOCATE(JQJ); ALLOCATE(JQJ (NP,NP )) IF(ALLOCATED(BJQJ))DEALLOCATE(BJQJ); ALLOCATE(BJQJ (NP,NP)) IF(ALLOCATED(JQR))DEALLOCATE(JQR); ALLOCATE(JQR (NP )) IF(ALLOCATED(U ))DEALLOCATE(U); ALLOCATE(U (NP )) IF(ALLOCATED(EIGW))DEALLOCATE(EIGW); ALLOCATE(EIGW(NP )) IF(ALLOCATED(EIGV))DEALLOCATE(EIGV); ALLOCATE(EIGV(NP,NP )) WRITE(*,'(A)') ' >>> Computing JQr on '//TRIM(VTOS(UTL_OMP_GET_MAX_THREADS()))//' threads <<<' !## construct jTqr (<--- r is residual for current parameter set) JQR=0.0; I=0; IPARAM=0 DO IP1=1,SIZE(PEST%PARAM) !## row IF(ABS(PEST%PARAM(IP1)%PACT).NE.1)CYCLE IPARAM=IPARAM+1; IF(PEST%PARAM(IP1)%PACT.NE.1)CYCLE DF1=PEST%PARAM(IP1)%PDELTA I=I+1 !$OMP PARALLEL PRIVATE(J,H1,H2,DJ1,DJ2,W,TS) SHARED (MSR,DF1) TS=0.0D0 SS=0.0D0 !$OMP DO DO J=1,MSR%NOBS H1=MSR%HG(IPARAM,J); H2=MSR%HL(0,J) DJ1=(H1-H2)/DF1 ; DJ2=MSR%DHL(0,J) W =MSR%W(J) TS=TS+(DJ1*W*DJ2) ! JQR(I)=JQR(I)+(DJ1*W*DJ2) ENDDO !$OMP END DO !$OMP CRITICAL SS=SS+TS !$OMP END CRITICAL !$OMP END PARALLEL JQR(I)=JQR(I)+SS ENDDO WRITE(*,'(A)') ' End Computing JQr <<<' !## add parameter regularisation IF(PEST%PE_REGULARISATION.EQ.1)THEN !## add weight to balance number of observations and parameters WF=IPEST_GLM_BALANCEFACTOR() I=0 DO IP1=1,SIZE(PEST%PARAM) !## row IF(PEST%PARAM(IP1)%PACT.NE.1)CYCLE I=I+1 WW=1.0D0/PEST%PARAM(IP1)%PARSTD**2.0D0 !## add balance weighting WW=WF*WW DP=(PEST%PARAM(IP1)%ALPHA(2)-PEST%PARAM(IP1)%PPRIOR)*WW JQR(I)=JQR(I)+DP ENDDO ENDIF ELSE !## reuse jqj LJQJ=.FALSE. ENDIF IF(LAMBDA.LT.0.0D0)THEN WRITE(*,*) 'LAMBDA < 0 IS NOT POSSIBLE'; PAUSE ENDIF WRITE(IUPESTOUT,'(/A/)') 'Lambda_0 (lambda number '//TRIM(VTOS(ILAMBDA))//') prior to the update vector '//TRIM(VTOS(LAMBDA,'G',3)) !## try to come with some suitable lambdas DO !## set marquardt-lambda (cannot become negative) MARQUARDT=MAX(LAMBDA*PBMAN%LAMBDA_TEST(ILAMBDA),0.0D0) WRITE(*,'(A)') 'Computing lambda '//TRIM(VTOS(LAMBDA,'F',7))//'; MARQUARDT is '//TRIM(VTOS(MARQUARDT,'F',7)) !## construct jqj - normal matrix/hessian CALL IPEST_GLM_JQJ(IBATCH,MARQUARDT,BJQJ,JQJ,NP,.FALSE.,LJQJ); LJQJ=.FALSE. !## project on important singular values IF(LSVD)THEN !## add eigenvalue decomposition WRITE(*,'(A)') 'Computing eigenvalues '//TRIM(VTOS(NP))//'x'//TRIM(VTOS(NP))//' matrix' CALL IPEST_GLM_EIGDECOM(IBATCH,JQJ,EIGW,EIGV,NP,.TRUE.) !.FALSE.) !## make percentage of eigenvalues EIGW=(EIGW*100.0D0)/SUM(EIGW) EIGWTHRESHOLD=0.0D0 !% explained variance DO NE=1,NP; EIGWTHRESHOLD=EIGWTHRESHOLD+EIGW(NE); IF(EIGWTHRESHOLD.GE.PBMAN%EIGV)EXIT; ENDDO; NE=MIN(NE,NP) WRITE(IUPESTOUT,'(A)') 'Problem projected on '//TRIM(VTOS(NE))//' eigenvalues (out of '//TRIM(VTOS(NP))//')' ALLOCATE(P(NP,NE)); P(:,1:NE)=EIGV(:,1:NE); ALLOCATE(M(NE,NE),N(NE),RU(NE),PT(NE,NP)) !## compute pp=pt(jqj) on eigen-space PT=0.0; DO I=1,NE; DO J=1,NP DO K=1,NP PT(I,J)=PT(I,J)+P(K,I)*JQJ(K,J) ENDDO ENDDO; ENDDO !## project jqj on eigen-space M=0.0; DO I=1,NE; DO J=1,NE DO K=1,NP M(I,J)=M(I,J)+PT(I,K)*P(K,J) ENDDO ENDDO; ENDDO !## project right hand side on eigenspace N=0.0; DO I=1,NE DO K=1,NP N(I)=N(I)+P(K,I)*JQR(K) END DO ENDDO !## compute inverse of (Pt(JQJ)P)-1 -> B IF(ALLOCATED(INDX))DEALLOCATE(INDX); ALLOCATE(INDX(NE)) IF(ALLOCATED(B ))DEALLOCATE(B); ALLOCATE(B (NE,NE)) CALL IPEST_LUDECOMP_DBL(M,INDX,NE,ISING) IF(ISING.EQ.1)THEN; CALL IPEST_GLM_ERROR(IBATCH,'Singular matrix after projection on eigenvectors which is rather odd, stopped.'); IPEST_GLM_GRADIENT=-1; RETURN; ENDIF B=0.0D0; DO I=1,NE; B(I,I)=1.0D0; ENDDO DO I=1,NE; CALL IPEST_LUBACKSUB_DBL(M,INDX,B(1,I),NE); ENDDO !## compute U=(M)-1*N RU=0.0D0; DO I=1,NE; DO J=1,NE RU(I)=RU(I)+(B(J,I)*N(J)) ENDDO; ENDDO !## reproject reduced gradient on original space !## compute U=(M)-1*N U=0.0D0; DO I=1,NP; DO J=1,NE U(I)=U(I)+(P(I,J)*RU(J)) ENDDO; ENDDO DEALLOCATE(P,PT,M,N,RU,INDX,B) ELSE !## compute inverse of (JQJ)-1 -> B IF(ALLOCATED(INDX))DEALLOCATE(INDX); ALLOCATE(INDX(NP)) IF(ALLOCATED(B))DEALLOCATE(B); ALLOCATE(B(NP,NP)) IF(NP.EQ.1)THEN B(1,1)=1.0D0/JQJ(1,1) ELSE CALL IPEST_LUDECOMP_DBL(JQJ,INDX,NP,ISING) IF(ISING.EQ.1)THEN; CALL IPEST_GLM_ERROR(IBATCH,'Singular matrix,try activating the SVD option to avoid this, stopped.'); IPEST_GLM_GRADIENT=-1; RETURN; ENDIF B=0.0D0; DO I=1,NP; B(I,I)=1.0D0; ENDDO DO I=1,NP; CALL IPEST_LUBACKSUB_DBL(JQJ,INDX,B(1,I),NP); ENDDO ENDIF !## compute (JQJ)-1*JQR U=0.0D0 DO I=1,NP; DO J=1,NP U(I)=U(I)+(B(J,I)*JQR(J)) ENDDO; ENDDO DEALLOCATE(INDX,B) ENDIF !## pointing downhill U=-1.0D0*U !## store gradient update vector in list I=IPEST_GLM_UPGRADEVECTOR_LAMBDA(ILAMBDA,U) SELECT CASE (I) CASE (0) !## step inappropriate (too large) try another lambda LAMBDA=LAMBDA*GAMMA CASE (1) EXIT CASE (2) !## try the same lambda without the excluded parameter due to a bondary hit ILAMBDA=ILAMBDA-1 EXIT END SELECT ENDDO ENDDO !## lambda-test-loop WRITE(IUPESTOUT,'(A)') 'LAMBDA_0 POSTERIOR TO THE UPDATE VECTOR '//TRIM(VTOS(LAMBDA,'G',3)) !## print parameters for lambda testing WRITE(IUPESTOUT,'(15X,99A15)') ('LAMBDA'//TRIM(VTOS(J)),J=1,PBMAN%NLAMBDASEARCH) WRITE(IUPESTOUT,'(A15,99G15.7)') 'Parameters',(LAMBDA*PBMAN%LAMBDA_TEST(J),J=1,PBMAN%NLAMBDASEARCH) WRITE(IUPESTOUT,'(999A1)') ('-',I=1,15*(PBMAN%NLAMBDASEARCH+1)) DO I=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(I)%PACT).NE.1)CYCLE IF(PEST%PARAM(I)%PLOG.EQ.0)THEN WRITE(IUPESTOUT,'(A15,99F15.7)') PEST%PARAM(I)%ACRONYM,(PEST%PARAM(I)%LALPHA(J),J=1,PBMAN%NLAMBDASEARCH) ELSEIF(PEST%PARAM(I)%PLOG.EQ.1)THEN WRITE(IUPESTOUT,'(A15,99F15.7)') PEST%PARAM(I)%ACRONYM,(EXP(PEST%PARAM(I)%LALPHA(J)),J=1,PBMAN%NLAMBDASEARCH) ELSEIF(PEST%PARAM(I)%PLOG.EQ.2)THEN WRITE(IUPESTOUT,'(A15,99F15.7)') PEST%PARAM(I)%ACRONYM,(10.0D0**(PEST%PARAM(I)%LALPHA(J)),J=1,PBMAN%NLAMBDASEARCH) ENDIF ENDDO WRITE(IUPESTOUT,'(/70A1)') ('=',I=1,70); WRITE(IUPESTOUT,'(A)') 'End Lambda Testing'; WRITE(IUPESTOUT,'(70A1/)') ('=',I=1,70) IF(IUPESTOUT.GT.0)FLUSH(IUPESTOUT) IF(ALLOCATED(JQJ)) DEALLOCATE(JQJ ) IF(ALLOCATED(BJQJ))DEALLOCATE(BJQJ) IF(ALLOCATED(U )) DEALLOCATE(U ) IF(ALLOCATED(EIGW))DEALLOCATE(EIGW) IF(ALLOCATED(EIGV))DEALLOCATE(EIGV) IF(ALLOCATED(JQR ))DEALLOCATE(JQR ) IF(ALLOCATED(S ))DEALLOCATE(S ) IF(ALLOCATED(JS ))DEALLOCATE(JS ) IPEST_GLM_GRADIENT=1 END FUNCTION IPEST_GLM_GRADIENT !###==================================================================== SUBROUTINE IPEST_GLM_GETBOUNDARY(IP1,IBND,P1,P2,PMIN,PMAX) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: IP1 INTEGER,INTENT(OUT) :: IBND REAL(KIND=DP_KIND),INTENT(OUT) :: P1,P2,PMIN,PMAX !## parameter adjustment hit the parameter boundary P1 =PEST%PARAM(IP1)%ALPHA(1) P2 =PEST%PARAM(IP1)%ALPHA(2) PMIN=PEST%PARAM(IP1)%PMIN PMAX=PEST%PARAM(IP1)%PMAX IBND=0 !## shoot over IF(P1.LE.PMIN)IBND=-1 IF(P1.GE.PMAX)IBND= 1 ! !## too close ! IF(ABS(P1-PMIN).LE.XPBND)IBND=-1; IF(ABS(PMAX-P1).LE.XPBND)IBND= 1 END SUBROUTINE IPEST_GLM_GETBOUNDARY !###==================================================================== INTEGER FUNCTION IPEST_GLM_UPGRADEVECTOR_LAMBDA(ILAMBDA,U) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: ILAMBDA REAL(KIND=DP_KIND),DIMENSION(:),INTENT(IN) :: U INTEGER :: I,IP1,IP2,IBND REAL(KIND=DP_KIND) :: F,PMAX,PMIN,P1,P2 !## exit code IPEST_GLM_UPGRADEVECTOR_LAMBDA=0 !## fill in by default DO IP1=1,SIZE(PEST%PARAM); PEST%PARAM(IP1)%ALPHA(1)=PEST%PARAM(IP1)%ALPHA(2); ENDDO !## update parameters I=0; DO IP1=1,SIZE(PEST%PARAM) IF(PEST%PARAM(IP1)%PACT.NE.1)CYCLE I=I+1; DO IP2=1,SIZE(PEST%PARAM) IF(PEST%PARAM(IP1)%PIGROUP.EQ.ABS(PEST%PARAM(IP2)%PIGROUP))THEN PEST%PARAM(IP2)%ALPHA(1)=PEST%PARAM(IP2)%ALPHA(2)+U(I) ENDIF ENDDO ENDDO !## check for size of adjustment DO IP1=1,SIZE(PEST%PARAM) !## inactive parameter IF(PEST%PARAM(IP1)%PACT.NE.1)CYCLE !## check size of adjustment IF(PEST%PARAM(IP1)%PLOG.EQ.1)THEN F=EXP(PEST%PARAM(IP1)%ALPHA(1))/EXP(PEST%PARAM(IP1)%ALPHA(2)) ELSEIF(PEST%PARAM(IP1)%PLOG.EQ.2)THEN F=10.0D0**(PEST%PARAM(IP1)%ALPHA(1))/10.0D0**(PEST%PARAM(IP1)%ALPHA(2)) ELSE F=PEST%PARAM(IP1)%ALPHA(1)/PEST%PARAM(IP1)%ALPHA(2) ENDIF !## adjustment too large try another lambda IF(F.LT.1.0D0/PEST%PARAM(IP1)%PINCREASE.OR. & F.GT. PEST%PARAM(IP1)%PINCREASE)RETURN CALL IPEST_GLM_GETBOUNDARY(IP1,IBND,P1,P2,PMIN,PMAX) !## parameter hits the boundary IF(IBND.NE.0)THEN !## hits the same boundary as before - skip it IF(IBND.EQ.PEST%PARAM(IP1)%IBND)THEN !## ignore this parameter (group) for now - reuse current lambda and search another update vector ignoring this parameters again PEST%PARAM(IP1)%PACT=-1; IPEST_GLM_UPGRADEVECTOR_LAMBDA=2 !; RETURN ENDIF ENDIF ENDDO !## if hit to boundary, return IF(IPEST_GLM_UPGRADEVECTOR_LAMBDA.EQ.2)RETURN !## corrects all gradients in case limits are overwritten !## adjust all parameters DO IP1=1,SIZE(PEST%PARAM) IF(PEST%PARAM(IP1)%PACT.NE.1)CYCLE PMIN=PEST%PARAM(IP1)%PMIN PMAX=PEST%PARAM(IP1)%PMAX !## adjust in case parameter exceeds limits IF(PEST%PARAM(IP1)%ALPHA(1).LT.PMIN)PEST%PARAM(IP1)%ALPHA(1)=PMIN IF(PEST%PARAM(IP1)%ALPHA(1).GT.PMAX)PEST%PARAM(IP1)%ALPHA(1)=PMAX ENDDO !## copy gradients to all groups DO IP1=1,SIZE(PEST%PARAM) IF(PEST%PARAM(IP1)%PACT.NE.1)CYCLE DO IP2=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(IP2)%PIGROUP).EQ.PEST%PARAM(IP1)%PIGROUP)PEST%PARAM(IP2)%ALPHA(1)=PEST%PARAM(IP1)%ALPHA(1) ENDDO ENDDO !## copy alphas to correct vector with updates per lambda DO IP1=1,SIZE(PEST%PARAM) PEST%PARAM(IP1)%LALPHA(ILAMBDA)=PEST%PARAM(IP1)%ALPHA(1) ENDDO !## correct update gradient found IPEST_GLM_UPGRADEVECTOR_LAMBDA=1 END FUNCTION IPEST_GLM_UPGRADEVECTOR_LAMBDA !###==================================================================== LOGICAL FUNCTION IPEST_GLM_ECHOPARAMETERS(IBATCH,ITER) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: ITER,IBATCH INTEGER :: IP1,N,I,J REAL(KIND=DP_KIND) :: C1,C2,C3 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: GRADUPDATE IPEST_GLM_ECHOPARAMETERS=.FALSE. N=0; C1=0.0D0 DO I=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(I)%PACT).EQ.1)THEN N=N+1 C1=C1+(PEST%PARAM(I)%ALPHA(2)-PEST%PARAM(I)%ALPHA(1))**2.0D0 ENDIF ENDDO C1=SQRT(C1) IF(ITER.EQ.0)THEN IF(PBMAN%RESTART.NE.0)THEN MSR%TJ=MSR%TJRESTART; MSR%RJ=MSR%RJRESTART ELSE WRITE(IUPESTEFFICIENCY,'(4(F15.3,1X))') MSR%TJ,MSR%TJ-MSR%RJ,MSR%RJ,REAL(SQRT(MSR%TJ))/REAL(MSR%NOBS,8) ENDIF ELSE C2=(1.0D0-MSR%TJ/MSR%PJ)*100.0D0 C3=(1.0D0-MSR%TJ/MSR%TJ_H(0))*100.0D0 WRITE(IUPESTEFFICIENCY,'(8(F15.3,1X))') MSR%TJ,MSR%TJ-MSR%RJ,MSR%RJ,REAL(SQRT(MSR%TJ))/REAL(MSR%NOBS,8),C1,C2,C3 ENDIF !## save alphas for history-logging DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PACT.EQ.0)CYCLE IF(PEST%PARAM(I)%PLOG.EQ.1)THEN PEST%PARAM(I)%ALPHA_HISTORY(ITER)=EXP(PEST%PARAM(I)%ALPHA(1)) ELSEIF(PEST%PARAM(I)%PLOG.EQ.2)THEN PEST%PARAM(I)%ALPHA_HISTORY(ITER)=10.0D0**(PEST%PARAM(I)%ALPHA(1)) ELSE PEST%PARAM(I)%ALPHA_HISTORY(ITER)=PEST%PARAM(I)%ALPHA(1) ENDIF ENDDO !## replace old by new parameter values PEST%PARAM%ALPHA(2)=PEST%PARAM%ALPHA(1) WRITE(IUPESTOUT,'(A/)') 'Optimization History:' WRITE(IUPESTOUT,'(A15,999(5X,A7,I3.3))') 'Statistics',('ITER',I,I=ITER,0,-1) WRITE(IUPESTOUT,'(999A1)') ('-',I=1,19+(ITER+1)*15) IF(MSR%NOBS.GT.0)THEN WRITE(IUPESTOUT,'(A15,99F15.3)') 'Total Obj. Val.',(MSR%TJ_H(I),I=ITER,0,-1) WRITE(IUPESTOUT,'(A15,99F15.3)') 'Meas. Obj. Val.',(MSR%TJ_H(I)-MSR%RJ_H(I),I=ITER,0,-1) WRITE(IUPESTOUT,'(A15,99F15.3)') 'Param Obj. Val.',(MSR%RJ_H(I),I=ITER,0,-1) WRITE(IUPESTOUT,'(A15,99I15)') 'Number of Obs. ',(MSR%NOBS,I=ITER,0,-1) WRITE(IUPESTOUT,'(A15,99F15.3)') 'Goodness of Fit',(MSR%GOF_H(I),I=ITER,0,-1) WRITE(IUPESTOUT,'(A15,99F15.3)') 'Nash Sutcliffe ',(MSR%NSC_H(I),I=ITER,0,-1) ELSE WRITE(IUPESTOUT,'(A15,99I15)') 'Number of Obs. ',(MSR%NOBS,I=ITER,0,-1) ENDIF WRITE(IUPESTOUT,'(/A15,999(5X,A7,I3.3))') 'Parameter',('ITER',I,I=ITER,0,-1) WRITE(IUPESTOUT,'(999A1)') ('-',I=1,15+(ITER+1)*15) ALLOCATE(GRADUPDATE(ITER)); GRADUPDATE=0.0D0 N=0 DO IP1=1,SIZE(PEST%PARAM) WRITE(BLINE,'(99(F15.7))') (PEST%PARAM(IP1)%ALPHA_HISTORY(I),I=ITER,0,-1) IF(ABS(PEST%PARAM(IP1)%PACT).EQ.1)THEN WRITE(IUPESTOUT,'(A15,A)') PEST%PARAM(IP1)%ACRONYM,TRIM(BLINE) N=N+1 DO I=1,ITER GRADUPDATE(I)=GRADUPDATE(I)+(PEST%PARAM(IP1)%ALPHA_HISTORY(I)-PEST%PARAM(IP1)%ALPHA_HISTORY(I-1))**2.0D0 ENDDO ENDIF ENDDO GRADUPDATE=SQRT(GRADUPDATE) WRITE(IUPESTOUT,'(999A1)') ('-',I=1,15+(ITER+1)*15) WRITE(IUPESTOUT,'(A15,99F15.7/)') 'Adjustment',(GRADUPDATE(I),I=ITER,1,-1) DEALLOCATE(GRADUPDATE) WRITE(IUPESTRUNFILE,'(/A,I10/)') 'Copy in the runfile, iteration ',ITER DO I=1,SIZE(PEST%PARAM) J=MAX(0,MIN(1,ABS(PEST%PARAM(I)%PACT))) IF(PEST%PARAM(I)%PLOG.EQ.1)THEN WRITE(IUPESTRUNFILE,'(I2,1X,A,1X,I5,1X,I7,1X,5(F10.3,1X),I7,1X,I2,1X,A15,2F10.3,2A15)') J, & !## iact PEST%PARAM(I)%PPARAM, & !## ptype PEST%PARAM(I)%PILS, & !## ilayer/system PEST%PARAM(I)%PIZONE, & !## zone number EXP(PEST%PARAM(I)%ALPHA(1)), & !## initial value EXP(PEST%PARAM(I)%PDELTA), & !## finite difference step EXP(PEST%PARAM(I)%PMIN), & !## minimal value EXP(PEST%PARAM(I)%PMAX),& !## maximal value PEST%PARAM(I)%PINCREASE,& !## maximal adjust factor ABS(PEST%PARAM(I)%PIGROUP),& !## group number PEST%PARAM(I)%PLOG,& !## log transformed TRIM(PEST%PARAM(I)%ACRONYM), & EXP(PEST%PARAM(I)%PPRIOR), & EXP(PEST%PARAM(I)%PARSTD), & TRIM(PEST%PARAM(I)%SDATE), & TRIM(PEST%PARAM(I)%EDATE) ELSEIF(PEST%PARAM(I)%PLOG.EQ.2)THEN WRITE(IUPESTRUNFILE,'(I2,1X,A,1X,I5,1X,I7,1X,5(F10.3,1X),I7,1X,I2,1X,A15,2F10.3,2A15)') J, & !## iact PEST%PARAM(I)%PPARAM, & !## ptype PEST%PARAM(I)%PILS, & !## ilayer/system PEST%PARAM(I)%PIZONE, & !## zone number 10.0D0**(PEST%PARAM(I)%ALPHA(1)), & !## initial value 10.0D0**(PEST%PARAM(I)%PDELTA), & !## finite difference step 10.0D0**(PEST%PARAM(I)%PMIN), & !## minimal value 10.0D0**(PEST%PARAM(I)%PMAX),& !## maximal value PEST%PARAM(I)%PINCREASE,& !## maximal adjust factor ABS(PEST%PARAM(I)%PIGROUP),& !## group number PEST%PARAM(I)%PLOG,& !## log transformed TRIM(PEST%PARAM(I)%ACRONYM), & 10.0D0**(PEST%PARAM(I)%PPRIOR), & 10.0D0**(PEST%PARAM(I)%PARSTD), & TRIM(PEST%PARAM(I)%SDATE), & TRIM(PEST%PARAM(I)%EDATE) ELSE WRITE(IUPESTRUNFILE,'(I2,1X,A,1X,I5,1X,I7,1X,5(F10.3,1X),I7,1X,I2,1X,A15,2F10.3,2A15)') J, & !## iact PEST%PARAM(I)%PPARAM, & !## ptype PEST%PARAM(I)%PILS, & !## ilayer/system PEST%PARAM(I)%PIZONE, & !## zone number PEST%PARAM(I)%ALPHA(1), & !## initial value PEST%PARAM(I)%PDELTA, & !## finite difference step PEST%PARAM(I)%PMIN, & !## minimal value PEST%PARAM(I)%PMAX,& !## maximal value PEST%PARAM(I)%PINCREASE,& !## maximal adjust factor ABS(PEST%PARAM(I)%PIGROUP),& !## group number PEST%PARAM(I)%PLOG, & !## log transformed TRIM(PEST%PARAM(I)%ACRONYM), & PEST%PARAM(I)%PPRIOR, & PEST%PARAM(I)%PARSTD, & TRIM(PEST%PARAM(I)%SDATE), & TRIM(PEST%PARAM(I)%EDATE) ENDIF ENDDO IF(IUPESTOUT.GT.0)FLUSH(IUPESTOUT); IF(IUPESTPROGRESS.GT.0)FLUSH(IUPESTPROGRESS); IF(IUPESTEFFICIENCY.GT.0)FLUSH(IUPESTEFFICIENCY) IF(IUPESTSENSITIVITY.GT.0)FLUSH(IUPESTSENSITIVITY); IF(IUPESTRUNFILE.GT.0)FLUSH(IUPESTRUNFILE); IF(IUPESTJACOBIAN.GT.0)FLUSH(IUPESTJACOBIAN) IF(ITER.GT.0)THEN WRITE(*,'(/A/)') 'Current Obj.Func. '//TRIM(VTOS(MSR%TJ,'F',7))//'; current/total improvement '//TRIM(VTOS(C2,'F',7))//';'//TRIM(VTOS(C3,'F',7))//'%' !## continue ? IF(ITER+1.GT.PEST%PE_MXITER)THEN CALL IPEST_GLM_ERROR(IBATCH,'Pest iteration terminated: PEST_ITER (='//TRIM(VTOS(PEST%PE_MXITER))//') = PEST_NITER (='// & TRIM(VTOS(PEST%PE_MXITER))//')'); RETURN ENDIF C2=(1.0D0-MSR%TJ_H(ITER-1)/MSR%TJ_H(0))*100.0D0 C3=(1.0D0-MSR%TJ_H(ITER )/MSR%TJ_H(0))*100.0D0 C1=C3-C2 IF(C1.LT.PEST%PE_STOP)THEN CALL IPEST_GLM_ERROR(IBATCH,'Pest iteration terminated as decrease objective function ('//TRIM(VTOS(C1,'G',7))// & '%) < PEST_JSTOP ('//TRIM(VTOS(PEST%PE_STOP,'G',7))//'%)'); RETURN ENDIF ENDIF IF(MSR%NOBS.EQ.0)MSR%TJ=0.0D0 IF(MSR%TJ.LE.0.0D0)THEN CALL IPEST_GLM_ERROR(IBATCH,'Objective Function <= 0.0 ('//TRIM(VTOS(MSR%TJ,'G',7))//')'); RETURN ENDIF !## next iteration WRITE(IUPESTOUT,'(/A/)') ' *** Next Optimization Cycle '//TRIM(VTOS(ITER+1))//' ***' WRITE(*,'(/A/)') ' *** Next Optimization Cycle '//TRIM(VTOS(ITER+1))//' ***' IPEST_GLM_ECHOPARAMETERS=.TRUE. END FUNCTION IPEST_GLM_ECHOPARAMETERS !###==================================================================== SUBROUTINE IPEST_GLM_JQJ(IBATCH,MARQUARDT,BJQJ,JQJ,NP,LCOV,LJQJ) !###==================================================================== IMPLICIT NONE INTEGER,PARAMETER :: PRINTLIMIT=50 REAL(KIND=DP_KIND),INTENT(IN) :: MARQUARDT REAL(KIND=DP_KIND),INTENT(INOUT),DIMENSION(NP,NP) :: BJQJ REAL(KIND=DP_KIND),INTENT(OUT),DIMENSION(NP,NP) :: JQJ INTEGER,INTENT(IN) :: NP,IBATCH LOGICAL,INTENT(IN) :: LCOV,LJQJ INTEGER :: I,J,IP1,IP2,II,JJ,N,IPARAM,ISING,IERROR,IOS,I1,I2 REAL(KIND=DP_KIND) :: DF1,DJ1,DJ2,B1,CB,W,ZW,Z1,Z2,H1,H2,F,Z REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: B REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: COR,COV REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: DIAG,DJ CHARACTER(LEN=15),ALLOCATABLE,DIMENSION(:) :: CPARAM INTEGER,ALLOCATABLE,DIMENSION(:) :: INDX INTEGER :: CLOCK_RATE LOGICAL :: LEX !## save jacobian for all IF(LCOV)THEN LEX=.FALSE. IF(PBMAN%BOSS.EQ.1.AND.PBMAN%NWORKERS.EQ.0)LEX=.TRUE. IF(PBMAN%BOSS.EQ.0)LEX=.TRUE. IF(LEX)THEN WRITE(*,'(A)') ' >>> Exporting Jacobians to a TXT file' !## number of observations and parameters WRITE(IUPESTJACOBIAN,'(2A10)') 'NBR_OBS','NPARAMS' IPARAM=0; DO IP1=1,SIZE(PEST%PARAM); IF(ABS(PEST%PARAM(IP1)%PACT).NE.1)CYCLE; IPARAM=IPARAM+1; ENDDO !## total number of active parameters WRITE(IUPESTJACOBIAN,'(2I10)') MSR%NOBS,IPARAM ALLOCATE(CPARAM(IPARAM)) IPARAM=0; DO IP1=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(IP1)%PACT).NE.1)CYCLE IPARAM=IPARAM+1; CPARAM(IPARAM)=PEST%PARAM(IP1)%ACRONYM ENDDO WRITE(IUPESTJACOBIAN,'(A3,A32,9999A15)') 'IPF','LABEL','DATE','ACCEPT','OBSERVATION','WEIGHT','HEAD',(CPARAM(I),I=1,IPARAM) DO J=1,MSR%NOBS WRITE(IUPESTJACOBIAN,'(I3,A32,I15,9999F15.7)') MSR%IPF(J),ADJUSTR(MSR%CLABEL(J)),MSR%IDATE(J),MSR%D(J),MSR%O(J),MSR%W(J),MSR%HL(0,J),(MSR%HG(I,J),I=1,IPARAM) ENDDO DEALLOCATE(CPARAM) WRITE(*,'(A)') ' Exporting Finished Jacobians <<< ' ENDIF ENDIF !## construct jqj - NORMAL MATRIX/HESSIAN - can take a while, so only in case ljqj=.true. IF(LJQJ)THEN WRITE(*,'(A)') ' >>> Computing JQJ on '//TRIM(VTOS(UTL_OMP_GET_MAX_THREADS()))//' threads <<<' CALL SYSTEM_CLOCK (I1) ALLOCATE(DJ(NP,MSR%NOBS)); DJ=0.0D0 WRITE(*,'(A)') ' >>> Computing Gradients' DJ=0.0D0 !## compute gradients I=0; IPARAM=0; DO IP1=1,SIZE(PEST%PARAM) IF(ABS(PEST%PARAM(IP1)%PACT).NE.1)CYCLE IPARAM=IPARAM+1; IF(PEST%PARAM(IP1)%PACT.NE.1)CYCLE I=I+1; DF1=PEST%PARAM(IP1)%PDELTA !$OMP PARALLEL PRIVATE(J,H1,H2) SHARED (I,DJ,DF1,MSR) !$OMP DO DO J=1,MSR%NOBS H1=MSR%HG(IPARAM,J); H2=MSR%HL(0,J); DJ(I,J)=(H1-H2)/DF1 ENDDO !$OMP END DO !$OMP END PARALLEL F=REAL(I,8)/REAL(NP,8); WRITE(6,'(A,F10.2,A)') '+ >>> Progress ',F*100.0D0,'%' ENDDO WRITE(*,'(A)') ' End Computing Gradients <<<' ! !### TEMP ! WRITE(IUPESTJACOBIAN,'(2A10)') 'DJ' ! DO J=1,MSR%NOBS ! WRITE(IUPESTJACOBIAN,'(9999F15.7)') (DJ(I,J),I=1,SIZE(DJ,1)) ! ENDDO ! !### TEMP WRITE(*,'(A)') ' >>> Computing Product' JQJ=0.0D0; N=0 !$OMP PARALLEL PRIVATE(I,II,J,DJ1,DJ2,W) SHARED (NP,DJ,MSR,JQJ,N) !$OMP DO DO I=1,NP DO II=1,I DO J=1,MSR%NOBS DJ1=DJ(I,J); DJ2=DJ(II,J); W=MSR%W(J) JQJ(II,I)=JQJ(II,I)+(DJ1*W*DJ2) ENDDO ENDDO N=N+1 F=REAL(N,8)/REAL(NP,8); WRITE(6,'(A,F10.2,A)') '+ >>> Progress ',F*100.0D0,'%' ENDDO !$OMP END DO !$OMP END PARALLEL !## mirror results DO I=1,NP DO J=I,NP JQJ(J,I)=JQJ(I,J) ENDDO ENDDO WRITE(*,'(A)') ' End Computing Product <<<' DEALLOCATE(DJ) ! !### TEMP ! WRITE(IUPESTJACOBIAN,'(2A10)') 'Jqj' ! DO J=1,size(jqj,1) ! WRITE(IUPESTJACOBIAN,'(9999F15.7)') (jqj(I,J),I=1,SIZE(Jqj,2)) ! ENDDO ! !### TEMP CALL SYSTEM_CLOCK (I2,CLOCK_RATE) WRITE(*,'(A)') ' Computing Finished JQJ '//TRIM(VTOS(REAL(I2-I1,8)/REAL(CLOCK_RATE,8),'F',2))//' secs. <<<' BJQJ=JQJ ELSE WRITE(*,'(A)') ' >>> Reusing JQJ' JQJ=BJQJ WRITE(*,'(A)') ' Reusing Finished JQJ <<< ' ENDIF IF(.NOT.LCOV)THEN IF(PEST%PE_REGULARISATION.EQ.0)THEN !## levenberg DO I=1,NP JQJ(I,I)=JQJ(I,I)+MARQUARDT ENDDO ELSE ! WF=IPEST_GLM_BALANCEFACTOR() !## levenberg-marquardt ALLOCATE(DIAG(NP,NP)); DIAG=0.0D0 I=0 DO IP1=1,SIZE(PEST%PARAM) IF(PEST%PARAM(IP1)%PACT.NE.1)CYCLE I=I+1; DIAG(I,I)=1.0D0/PEST%PARAM(IP1)%PARSTD**2.0D0 ENDDO !## levenberg-marquardt (use parameter prior covariance) DO I=1,NP DO J=1,NP JQJ(J,I)=JQJ(J,I)+MARQUARDT*DIAG(J,I) ENDDO ENDDO DEALLOCATE(DIAG) ENDIF RETURN ENDIF IF(ALLOCATED(COR ))DEALLOCATE(COR); ALLOCATE(COR(NP,NP)) IF(ALLOCATED(COV ))DEALLOCATE(COV); ALLOCATE(COV(NP,NP)) IF(ALLOCATED(INDX))DEALLOCATE(INDX); ALLOCATE(INDX(NP)) IF(ALLOCATED(B ))DEALLOCATE(B); ALLOCATE(B(NP,NP)) !## compute inverse of (JQJ)-1 -> B = covariance matrix WRITE(*,'(A)') 'Computing inverse Jacobian ...' CALL IPEST_LUDECOMP_DBL(JQJ,INDX,NP,ISING) !## matrix not singular - compute inverse IF(ISING.EQ.0)THEN B=0.0D0; DO I=1,NP; B(I,I)=1.0D0; ENDDO DO I=1,NP; CALL IPEST_LUBACKSUB_DBL(JQJ,INDX,B(1,I),NP); ENDDO ! !### TEMP ! WRITE(IUPESTJACOBIAN,'(2A10)') 'b' ! DO J=1,size(b,1) ! WRITE(IUPESTJACOBIAN,'(9999F15.7)') (b(I,J),I=1,SIZE(b,2)) ! ENDDO ! !### TEMP !## b1 is variance of each of model response N=MAX(1,MSR%NOBS-NP); B1=MSR%TJ/REAL(N,8) ! !### TEMP ! WRITE(IUPESTJACOBIAN,*) 'n',n,'b1',b1,'tj',msr%tj ! !### TEMP !## parameter covariance matrix scaled to this variance DO I=1,NP; DO J=1,NP; B(I,J)=B1*B(I,J); ENDDO; ENDDO ! !### TEMP ! WRITE(IUPESTJACOBIAN,'(A,f15.7)') 'b_scaled_with_',b1 ! DO J=1,size(b,1) ! WRITE(IUPESTJACOBIAN,'(9999F15.7)') (b(I,J),I=1,SIZE(b,2)) ! ENDDO ! !### TEMP IF(NP.LT.PRINTLIMIT)THEN WRITE(IUPESTOUT,'(/A/)') 'Parameter Covariance Matrix (m2):' CALL IPEST_GLM_WRITEHEADER(' ',NP,IUPESTOUT) WRITE(IUPESTOUT,'(9999A1)') ('-',I=1,15+NP*15) ENDIF I=0 DO IP1=1,SIZE(PEST%PARAM) IF(PEST%PARAM(IP1)%PACT.EQ.1)THEN I=I+1 IF(NP.LT.PRINTLIMIT)THEN WRITE(IUPESTOUT,'(A15,99999G15.7)') PEST%PARAM(IP1)%ACRONYM,(B(I,J),J=1,NP) ENDIF DO J=1,NP; COV(I,J)=B(I,J); ENDDO ENDIF ENDDO !## parameter correlation matrix IF(NP.LT.PRINTLIMIT)THEN WRITE(IUPESTOUT,'(/A/)') 'Parameter Correlation Matrix (-)' CALL IPEST_GLM_WRITEHEADER(' ',NP,IUPESTOUT) WRITE(IUPESTOUT,'(9999A1)') ('-',I=1,15+NP*15) ENDIF COR=0.0D0; I=0 DO IP1=1,SIZE(PEST%PARAM) IF(PEST%PARAM(IP1)%PACT.EQ.1)THEN I=I+1 DO J=1,NP CB=B(I,I)*B(J,J) IF(CB.GT.0.0D0)THEN COR(I,J)=B(I,J)/SQRT(CB) ELSE COR(I,J)=0.0D0 ENDIF ENDDO IF(NP.LT.PRINTLIMIT)THEN WRITE(IUPESTOUT,'(A15,99999F15.7)',IOSTAT=IOS) PEST%PARAM(IP1)%ACRONYM,(COR(I,J),J=1,NP) ENDIF ENDIF ENDDO !## write per parameter highly correlated other parameter DO JJ=1,2 IF(JJ.EQ.2)WRITE(IUPESTOUT,'(/A/)') 'List of Parameter Highly Correlated (correlation > 95%)' IP1=0; N=0 DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PACT.EQ.1)THEN WRITE(BLINE,'(A15)') PEST%PARAM(I)%ACRONYM IP1=IP1+1 IP2=0 II=0 DO J=1,SIZE(PEST%PARAM) IF(PEST%PARAM(J)%PACT.EQ.1)THEN IP2=IP2+1 IF(I.NE.J)THEN WRITE(SLINE,'(A15)') PEST%PARAM(J)%ACRONYM IF(ABS(COR(IP1,IP2)).GE.0.95D0)THEN II=II+1 BLINE=TRIM(BLINE)//','//TRIM(SLINE)//'('//TRIM(VTOS(COR(IP1,IP2)*100.0D0,'F',1))//'%)' ENDIF ENDIF ENDIF ENDDO IF(II.GT.0)THEN N=N+1; IF(JJ.EQ.2)WRITE(IUPESTOUT,'(A)') TRIM(BLINE) ENDIF ENDIF ENDDO IF(N.EQ.0)EXIT ENDDO WRITE(IUPESTOUT,'(/A/)') 'Parameter Variance - Standard Parameter Error (Confidence Limits ~96%)' J=0; IERROR=0 DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PACT.EQ.1)THEN J=J+1 IF(COV(J,J).GT.0.0D0)THEN !## stdev PEST%PARAM(I)%ALPHA_ERROR_STDEV=SQRT(COV(J,J)) ELSE !## error value - should not happen PEST%PARAM(I)%ALPHA_ERROR_STDEV=-999.99D0 WRITE(IUPESTOUT,*) 'Active Parameter#,',J,' Variance ',COV(J,J) IERROR=IERROR+1 ENDIF !## check whether current other parameters belong to this group DO IP1=1,SIZE(PEST%PARAM) !## active and follower of group IF(PEST%PARAM(IP1)%PACT.EQ.2)THEN IF(ABS(PEST%PARAM(IP1)%PIGROUP).EQ.PEST%PARAM(I)%PIGROUP)THEN PEST%PARAM(IP1)%ALPHA_ERROR_STDEV=PEST%PARAM(I)%ALPHA_ERROR_STDEV ENDIF ENDIF ENDDO ENDIF ENDDO IF(IERROR.GT.0)THEN CALL IPEST_GLM_ERROR(IBATCH,'Errors (#'//TRIM(VTOS(IERROR))//') found in the Covariance Matrix, check your matrix, might by singular'); RETURN ENDIF WRITE(IUPESTOUT,'(15X,3A15)') 'Lower_Limit','Average','Upper Limit' WRITE(IUPESTOUT,'(60A1)') ('-',I=1,60) DO I=1,SIZE(PEST%PARAM) IF(PEST%PARAM(I)%PACT.EQ.1)THEN ZW=PEST%PARAM(I)%ALPHA_ERROR_STDEV*1.96D0 IF(PEST%PARAM(I)%PLOG.EQ.1)THEN !## maximize uncertainty ZW=MIN(10.0D0,ZW) Z =EXP(PEST%PARAM(I)%ALPHA(2)) Z1=EXP(PEST%PARAM(I)%ALPHA(2)-ZW) Z2=EXP(PEST%PARAM(I)%ALPHA(2)+ZW) ELSEIF(PEST%PARAM(I)%PLOG.EQ.2)THEN !## maximize uncertainty ZW=MIN(10.0D0,ZW) Z =10.0D0**(PEST%PARAM(I)%ALPHA(2)) Z1=10.0D0**(PEST%PARAM(I)%ALPHA(2)-ZW) Z2=10.0D0**(PEST%PARAM(I)%ALPHA(2)+ZW) ELSE Z= PEST%PARAM(I)%ALPHA(2) Z1=PEST%PARAM(I)%ALPHA(2)-ZW Z2=PEST%PARAM(I)%ALPHA(2)+ZW ENDIF WRITE(BLINE,'(3G15.7)') Z1,Z,Z2 WRITE(IUPESTOUT,'(A15,A)') PEST%PARAM(I)%ACRONYM,TRIM(BLINE) IF(PEST%PARAM(I)%PACT.EQ.-1)THEN WRITE(BLINE,'(3A15)') 'Insens.','Insens.','Insens.' WRITE(IUPESTOUT,'(A15,A)') PEST%PARAM(I)%ACRONYM,TRIM(BLINE) ENDIF ENDIF ENDDO ELSE CALL IPEST_GLM_ERROR(IBATCH,'Singular matrix, cannot compute covariance matrix: skipped.') ENDIF IF(ALLOCATED(COV ))DEALLOCATE(COV) IF(ALLOCATED(COR ))DEALLOCATE(COR) IF(ALLOCATED(INDX))DEALLOCATE(INDX) IF(ALLOCATED(B)) DEALLOCATE(B) END SUBROUTINE IPEST_GLM_JQJ !###==================================================================== SUBROUTINE IPEST_GLM_EIGDECOM(IBATCH,JQJ,EIGW,EIGV,NP,LPRINT) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: NP,IBATCH LOGICAL,INTENT(IN) :: LPRINT REAL(KIND=DP_KIND),DIMENSION(NP,NP),INTENT(IN) :: JQJ REAL(KIND=DP_KIND),DIMENSION(NP,NP),INTENT(OUT) :: EIGV REAL(KIND=DP_KIND),DIMENSION(NP),INTENT(OUT) :: EIGW ! REAL(KIND=DP_KIND) :: DET INTEGER :: I,J,K,IU REAL(KIND=DP_KIND) :: TV,TEV,KAPPA REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: B REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: E IF(ALLOCATED(E))DEALLOCATE(E); ALLOCATE(E(NP)) IF(ALLOCATED(B))DEALLOCATE(B); ALLOCATE(B(NP,NP)) ! !## compute determinant of JQJ ! DET=IPEST_GLM_DET(JQJ,NP) ! IF(LPRINT)THEN ! WRITE(IUPESTOUT,'(A,E15.7)') 'Determinant JQJ + '//TRIM(VTOS(MARQUARDT,'F',7))//' * I = ',DET ! ENDIF !## copy jqj to b for eigenvalue decomposition B=JQJ !## eigenvalue of covariance matrix CALL LUDCMP_TRED2(B,NP,NP,EIGW,E) CALL LUDCMP_TQLI(EIGW,E,NP,NP,B) CALL LUDCMP_EIGSRT(EIGW,B,NP,NP) !## check computation of eigenvalues DO I=1,NP IF(EIGW(I).NE.EIGW(I).OR.EIGW(I).GT. HUGE(1.0D0).OR.EIGW(I).LT.-HUGE(1.0D0))THEN WRITE(*,*) 'Error finding eigenvalues for matrix, see eigenvalues.txt' IU=UTL_GETUNIT(); OPEN(IU,FILE='eigenvalues.txt',STATUS='REPLACE',ACTION='WRITE') DO J=1,SIZE(JQJ,1) WRITE(IU,'(999F15.7)') (JQJ(J,K),K=1,SIZE(JQJ,2)) ENDDO CLOSE(IU) ENDIF ENDDO IF(LPRINT)WRITE(IUPESTOUT,'(/10X,4A15)') 'Eigenvalues','Sing.Values','Variance','Explained Var.' DO I=1,NP; IF(EIGW(I).LE.0.0)EIGW(I)=0.0; ENDDO; TEV=SUM(EIGW) TV=0.0D0 DO I=1,NP TV=TV+(EIGW(I)*100.0D0/TEV) IF(EIGW(I).GT.0.0D0)THEN IF(LPRINT)WRITE(IUPESTOUT,'(I10,G15.7,3F15.7)') I,EIGW(I),SQRT(EIGW(I)),EIGW(I)*100.0D0/TEV,TV ELSE IF(LPRINT)WRITE(IUPESTOUT,'(I10,G15.7,3F15.7)') I,EIGW(I), EIGW(I) ,EIGW(I)*100.0D0/TEV,TV ENDIF ENDDO EIGV= B IF(SUM(EIGW).LT.0.0D0)THEN CALL IPEST_GLM_ERROR(IBATCH,'Warning, there is NO information (no eigenvalues) in parameter perturbation'); RETURN ENDIF ! EIGW=(EIGW*100.0D0)/SUM(EIGW) !## condition number !## get lowest non-zero DO I=NP,1,-1; IF(EIGW(I).GT.0.0D0)EXIT; ENDDO KAPPA=SQRT(EIGW(1))/SQRT(EIGW(I)); KAPPA=LOG10(KAPPA) IF(LPRINT)THEN WRITE(IUPESTOUT,'(/A,F15.7/)') 'Condition Number (kappa):',KAPPA IF(KAPPA.GT.15.0D0)THEN WRITE(IUPESTOUT,'(A)') '>>> If Kappa > 15, inversion is a concern due to parameters that are highly correlated <<<' ELSEIF(KAPPA.GT.30.0D0)THEN WRITE(IUPESTOUT,'(A)') '>>> If Kappa > 30, inversion is highly questionable due to parameters that are highly correlated <<<' ENDIF ENDIF IF(ALLOCATED(E))DEALLOCATE(E) IF(ALLOCATED(B))DEALLOCATE(B) END SUBROUTINE IPEST_GLM_EIGDECOM !###======================================================================== SUBROUTINE IPEST_GLM_WRITEHEADER(TXT,NP,IU) !###======================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: TXT INTEGER,INTENT(IN) :: IU,NP INTEGER :: I,J,N CHARACTER(LEN=15),ALLOCATABLE,DIMENSION(:) :: CTMP ALLOCATE(CTMP(NP)) N=0 DO J=1,SIZE(PEST%PARAM) IF(PEST%PARAM(J)%PACT.EQ.1)THEN N=N+1 WRITE(CTMP(N),'(A15)') PEST%PARAM(J)%ACRONYM ENDIF ENDDO WRITE(IU,'(A15,99999A15)') TXT,(CTMP(I),I=1,NP) DEALLOCATE(CTMP) END SUBROUTINE IPEST_GLM_WRITEHEADER !###======================================================================== DOUBLE PRECISION FUNCTION IPEST_GLM_DET(JQJ,N) !###======================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: N DOUBLE PRECISION,DIMENSION(N,N),INTENT(IN) :: JQJ DOUBLE PRECISION,DIMENSION(:,:),ALLOCATABLE :: MATRIX DOUBLE PRECISION :: M, TEMP INTEGER :: I, J, K, L LOGICAL :: DETEXISTS = .TRUE. ALLOCATE(MATRIX(N,N)) MATRIX=JQJ ! open(99,file='d:\tmp.m',status='REPLACE',action='write') ! write(99,*) '[' ! do j=1,n ! write(99,'(999f15.7)') (jqj(i,j),i=1,n) ! enddo ! close(99) L = 1 !## convert to upper triangular form DO K = 1, N-1 IF (MATRIX(K,K).EQ.0.0D0) THEN DETEXISTS = .FALSE. DO I = K+1, N IF (MATRIX(I,K).NE.0.0D0) THEN DO J = 1, N TEMP = MATRIX(I,J) MATRIX(I,J)= MATRIX(K,J) MATRIX(K,J) = TEMP END DO DETEXISTS = .TRUE. L=-L EXIT ENDIF END DO IF (DETEXISTS .EQV. .FALSE.) THEN IPEST_GLM_DET = 0.0D0 DEALLOCATE(MATRIX) RETURN END IF ENDIF DO J = K+1, N M = MATRIX(J,K)/MATRIX(K,K) DO I = K+1, N MATRIX(J,I) = MATRIX(J,I) - M*MATRIX(K,I) END DO END DO END DO !## calculate determinant by finding product of diagonal elements IPEST_GLM_DET = L DO I = 1, N IPEST_GLM_DET = IPEST_GLM_DET * MATRIX(I,I) END DO DEALLOCATE(MATRIX) END FUNCTION IPEST_GLM_DET !###==================================================================== LOGICAL FUNCTION IPEST_GLM_GETJ_AVG(DIR,IGRAD,IPARAM,CTYPE,IBATCH,ILAMBDA,MNAME) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: DIR,CTYPE,MNAME INTEGER,INTENT(IN) :: IGRAD,IPARAM,IBATCH,ILAMBDA INTEGER :: I,J,JJ,II,NC,NP,NPERIOD,III,K,KK,ILAY,NROWIPFTXT,IUIPFTXT,NCOLIPFTXT,IOS,NAJ,IERROR,SEED,IROW,ICOL,JROW,JCOL,IORG REAL(KIND=DP_KIND) :: X,Y,Z,H,WW,W,MC,MM,DHH,XCOR,YCOR,ZCOR,DHW,NSC,F,D,WF,F1,F2,N,WWGHG,WWGLG,DRES CHARACTER(LEN=256) :: DIRNAME,FNAME,DIRO CHARACTER(LEN=52) :: CID,TXT,CDATE CHARACTER(LEN=52),ALLOCATABLE,DIMENSION(:) :: IPFSTRING CHARACTER(LEN=12) :: CEXT REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: TSNODATA,M,C,MCOPY,CCOPY,STD INTEGER(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: IDATE,ID1,ID2 REAL(KIND=DP_KIND),DIMENSION(3) :: PC,PM,DYN !## percentiles computed/measured REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: AVGOBS,NUMOBS,AVGCOM,AVGW,SUMW,VAR,COV,COR INTEGER(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: VARIDATE INTEGER :: IU,NR,IEXT,IPER LOGICAL :: LYN,LDUAL TYPE(PSTMEASURE),ALLOCATABLE,DIMENSION(:) :: MEASURES !## no measurements IPEST_GLM_GETJ_AVG=.TRUE. !## initiate number of observations MSR%NOBS=0 !## run batch files and get number of observations IF(.NOT.PEST_GLM_BATCHFILES(1,CTYPE,IPARAM))THEN IPEST_GLM_GETJ_AVG=.FALSE.; RETURN ENDIF !## initialise variables MSR%TJ=0.0D0; MSR%RJ=0.0D0 !## write header IF(ASSOCIATED(PEST%MEASURES))THEN IF(TRIM(CTYPE).EQ.'L'.AND.IUPESTRESIDUAL.GT.0)THEN DO I=1,SIZE(MEASURES) J=INDEX(MEASURES(I)%IPFNAME,'\',.TRUE.)+1 FNAME=TRIM(DIRNAME)//'\'//TRIM(MEASURES(I)%IPFNAME(J:)) WRITE(IUPESTRESIDUAL,'(I10,A)') I,','//TRIM(FNAME) ENDDO ENDIF ENDIF IF(TRIM(CTYPE).EQ.'L'.OR.TRIM(CTYPE).EQ.'R')THEN !## write head for the first IF(IUPESTRESIDUAL.GT.0)THEN !## steady-state IF(PRJNPER.EQ.1)THEN WRITE(IUPESTRESIDUAL,'(2A16,A11,6A16,A11,A5)') 'X,','Y,','ILAY,','MSR,','MDL,', & 'MDL-MSR,','MDL-MSR*,','WRESIDUAL,','WEIGHT,','IPF,','LABEL' !## transient ELSE WRITE(IUPESTRESIDUAL,'(2A16,A11,8A16,A11,A5,27X,1X,A15)') 'X,','Y,','ILAY,','WEIGHT,','MSR,','MDL,','MDL-MSR*,', & 'DYNMSR,','DYNMDL,','DYNMSR-DYNMDL,','NASH-SUTCLIFF,','IPF,','LABEL','DATE' ENDIF ENDIF IF(IUPESTPRESIDUAL.GT.0)THEN WRITE(IUPESTPRESIDUAL,'(2A16,A11,4A16,2A11,A33,A14)') 'X,','Y,','LAYER,','MEASURE,','COMPUTED,','WEIGHT,','ACCEPTDIFF,','IPF,','IPOS,','LABEL,','DATE' ENDIF ENDIF IF(ASSOCIATED(PEST%MEASURES))THEN IPEST_GLM_GETJ_AVG=.FALSE.; LYN=.TRUE.; LDUAL=.TRUE. SEED=12345; DIRO=DIR; IF(PBMAN%OUTPUT.NE.'')DIRO=PBMAN%OUTPUT IF(PBMAN%IIES.EQ.1)THEN DIRNAME=TRIM(DIRO)//'\IIES_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IPARAM))//'_L#'//TRIM(VTOS(ILAMBDA))//'\TIMESERIES' ELSE DIRNAME=TRIM(DIRO)//'\IPEST_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IPARAM))//'\TIMESERIES' ENDIF !## create output files in case modflow6 is used (aggregate results from different modflow6- models) IF((PBMAN%IFORMAT.EQ.3.OR.PBMAN%IFORMAT.EQ.6).AND.PBMAN%IPESTP.EQ.1)THEN IF(.NOT.IPEST_GLM_CREATE_RESIDUALSFILES(DIRO,IPARAM,CTYPE,MNAME))RETURN ALLOCATE(MEASURES(1)) MEASURES(1)%IPFNAME=TRIM(DIRNAME)//'\OUTPUT_OBS.IPF' MEASURES(1)%IPFTYPE=1 MEASURES(1)%IXCOL=1; MEASURES(1)%IYCOL=2 MEASURES(1)%ILCOL=3; MEASURES(1)%IMCOL=4 !## check whether imcol needs to be negative IF(PEST%MEASURES(1)%IMCOL.LT.0)MEASURES(1)%IMCOL=-4 !## here ivcol is set to standard stdev input MEASURES(1)%IVCOL=5; MEASURES(1)%IDCOL=0 ELSE NP=SIZE(PEST%MEASURES); ALLOCATE(MEASURES(NP)) !## make temporary copy of measurement settings MEASURES=PEST%MEASURES DO I=1,SIZE(MEASURES) !## ensure a global path here IF(INDEX(MEASURES(I)%IPFNAME,'.\').GT.0)THEN CALL UTL_SUBST(MEASURES(I)%IPFNAME,'.\',TRIM(DIRNAME)//'\') CALL UTL_SUBST(MEASURES(I)%IPFNAME,'\\','\') ENDIF ENDDO ENDIF !## initialise variables II=0 NPERIOD=0; IF(ASSOCIATED(PEST%S_PERIOD))NPERIOD=SIZE(PEST%S_PERIOD) IF(NPERIOD.GT.0)THEN ALLOCATE(ID1(NPERIOD),ID2(NPERIOD)) DO I=1,NPERIOD; READ(PEST%S_PERIOD(I),*) ID1; READ(PEST%E_PERIOD(I),*) ID2; ENDDO ENDIF ELSE ALLOCATE(PEST%MEASURES(0)); MEASURES=PEST%MEASURES ENDIF !## if not yet allocated, allocate msr IF(.NOT.ASSOCIATED(MSR%CLABEL))THEN IF(.NOT.IPEST_GLM_ALLOCATEMSR(0,0,0,MEASURES))RETURN ENDIF IF(SIZE(PEST%MEASURES).LE.0)DEALLOCATE(PEST%MEASURES) IF(SIZE(PEST%MEASURES).GT.0)THEN !## process files DO I=1,SIZE(MEASURES) J=INDEX(MEASURES(I)%IPFNAME,'\',.TRUE.)+1 FNAME=TRIM(DIRNAME)//'\'//TRIM(MEASURES(I)%IPFNAME(J:)) IU=UTL_GETUNIT(); OPEN(IU,FILE=FNAME,STATUS='OLD',ACTION='READ',iostat=ios) IF(IOS.NE.0)THEN WRITE(*,'(/A/)') '>>> Cannot find: '//TRIM(FNAME)//' <<<' ENDIF READ(IU,*) NR; READ(IU,*) NC; DO J=1,NC; READ(IU,*); ENDDO; READ(IU,*) IEXT,CEXT IF(IEXT.EQ.0)THEN ALLOCATE(IPFSTRING(NC)) DO J=1,NR F1=REAL(J,8)/REAL(NR,8)*100.0D0 F2=REAL(I,8)/REAL(SIZE(MEASURES),8)*100.0D0 WRITE(6,'(A)') '+Progressing residuals '//TRIM(VTOS(F1,'F',2))//' for observation file '//TRIM(VTOS(F2,'F',2)) II=II+1 IF(.NOT.UTL_READCSVENTRY(IU,IPFSTRING))EXIT READ(IPFSTRING(1 ),*) MSR%X(II) READ(IPFSTRING(2 ),*) MSR%Y(II) READ(IPFSTRING(3 ),*) MSR%L(II) READ(IPFSTRING(4 ),*) MSR%O(II) READ(IPFSTRING(5 ),*) MSR%W(II) READ(IPFSTRING(6 ),*) MSR%C(II) READ(IPFSTRING(7 ),*) F1 READ(IPFSTRING(8 ),*) F2 READ(IPFSTRING(9 ),*) MSR%D(II) READ(IPFSTRING(10),*) IORG MSR%CLABEL(II)='MEASURE'//TRIM(VTOS(IORG))//'_IPF'//TRIM(VTOS(I)) !## steady-state MSR%IDATE(II)=0 MSR%IPF(II)=I MSR%LOC(II)=J ENDDO DEALLOCATE(IPFSTRING) !## transient ELSE IUIPFTXT=0; ALLOCATE(IPFSTRING(NC)) DO J=1,NR !## read x,y,layer,name,weight,acceptable_error,original_ipf_number IF(.NOT.UTL_READCSVENTRY(IU,IPFSTRING))EXIT READ(IPFSTRING(1),*,IOSTAT=IOS) X; IF(IOS.NE.0)THEN; WRITE(*,'(A)') '>>> ERROR READING X <<<'; STOP; ENDIF READ(IPFSTRING(2),*,IOSTAT=IOS) Y; IF(IOS.NE.0)THEN; WRITE(*,'(A)') '>>> ERROR READING Y <<<'; STOP; ENDIF READ(IPFSTRING(3),*,IOSTAT=IOS) ILAY; IF(IOS.NE.0)THEN; WRITE(*,'(A)') '>>> ERROR READING ILAY <<<'; STOP; ENDIF CID=IPFSTRING(4) READ(IPFSTRING(5),*,IOSTAT=IOS) WW; IF(IOS.NE.0)THEN; WRITE(*,'(A)') '>>> ERROR READING WEIGHT <<<'; STOP; ENDIF READ(IPFSTRING(6),*,IOSTAT=IOS) D; IF(IOS.NE.0)THEN; WRITE(*,'(A)') '>>> ERROR READING D <<<'; STOP; ENDIF READ(IPFSTRING(7),*,IOSTAT=IOS) IORG; IF(IOS.NE.0)THEN; WRITE(*,'(A)') '>>> ERROR READING IORG <<<'; STOP; ENDIF !## usage of gxg values IF(MEASURES(I)%IMCOL.LT.0.AND.MEASURES(I)%IVCOL.EQ.0)THEN READ(IPFSTRING(8),*,IOSTAT=IOS) WWGHG; IF(IOS.NE.0)THEN; WRITE(*,'(A)') '>>> ERROR READING WWGHG <<<'; STOP; ENDIF READ(IPFSTRING(9),*,IOSTAT=IOS) WWGLG; IF(IOS.NE.0)THEN; WRITE(*,'(A)') '>>> ERROR READING WWGLG <<<'; STOP; ENDIF ELSE WWGHG=WW; WWGLG=WW ENDIF LINE=TRIM(DIRNAME)//CHAR(92)//TRIM(CID)//'.'//TRIM(CEXT) IF(IUIPFTXT.EQ.0)IUIPFTXT=UTL_GETUNIT(); OPEN(IUIPFTXT,FILE=LINE,STATUS='OLD',ACTION='READ') READ(IUIPFTXT,*) NROWIPFTXT; READ(IUIPFTXT,*) NCOLIPFTXT ALLOCATE(TSNODATA(MAX(3,NCOLIPFTXT))) DO K=1,NCOLIPFTXT; READ(IUIPFTXT,*) TXT,TSNODATA(K); ENDDO ALLOCATE(M(NROWIPFTXT),C(NROWIPFTXT),IDATE(NROWIPFTXT),MCOPY(NROWIPFTXT),CCOPY(NROWIPFTXT)) IDATE=0; C=0.0D0; M=0.0D0; MCOPY=M; CCOPY=C IF(NCOLIPFTXT.LT.3)TSNODATA(3)=TSNODATA(2) !## get mean measure KK=0 DO K=1,NROWIPFTXT KK=KK+1 READ(IUIPFTXT,*,IOSTAT=IOS) IDATE(KK),M(KK),C(KK) !## error reading, skip it (can be caused by steady-state periods in between) IF(IOS.NE.0)THEN; KK=KK-1; CYCLE; ENDIF !## make double precision dates - if needed IDATE(KK)=UTL_COMPLETEDATE(IDATE(KK)) !## check period (if available) IF(NPERIOD.GT.0)THEN DO III=1,NPERIOD; IF(IDATE(KK).GE.ID1(III).AND.IDATE(KK).LE.ID2(III))EXIT; ENDDO IF(III.GT.NPERIOD)C(KK)=TSNODATA(3) ENDIF IF(M(KK).EQ.TSNODATA(2).OR.C(KK).EQ.TSNODATA(3))KK=KK-1 ENDDO !## add this measurement IF(KK.GT.0)THEN !## compute mean measurement in period XCOR=-9999.99D0 NSC =-9999.99D0 DYN =-9999.99D0 !## mean values MM=SUM(M(1:KK))/DBLE(KK) !## measurements MC=SUM(C(1:KK))/DBLE(KK) !## computed IF(PEST%PE_TARGET(2).GT.0.0D0)THEN DO K=1,KK; MCOPY(K)=M(K); CCOPY(K)=C(K); ENDDO !## percentiles CALL UTL_GETMED(MCOPY,KK,-999.99D0,(/10.0D0,50.0D0,90.0D0/),3,NAJ,PM) CALL UTL_GETMED(CCOPY,KK,-999.99D0,(/10.0D0,50.0D0,90.0D0/),3,NAJ,PC) !## measurements DYN(1)=PM(3)-PM(1) !## computed DYN(2)=PC(3)-PC(1) ENDIF !## compute cross-correlation IF(KK.GT.1)THEN XCOR=0.0D0; YCOR=0.0D0; ZCOR=0.0D0 DO K=1,KK XCOR=XCOR+(MM-M(K))*(MC-C(K)); YCOR=YCOR+(MM-M(K))**2.0D0; ZCOR=ZCOR+(MC-C(K))**2.0D0 ENDDO IF(YCOR.NE.0.0.AND.ZCOR.NE.0.0)XCOR=XCOR/(SQRT(YCOR)*SQRT(ZCOR)) NSC =UTL_NASH_SUTCLIFFE(M,C,KK) NSC=MIN(MAX(-999.9D0,NSC),999.9D0) ENDIF !## compute ghg/glg IF(MEASURES(I)%IMCOL.LT.0)THEN !## compute ghg(1)/glg(2) if ldual=.true. IF(.NOT.UTL_COMPUTE_GXG(M,IDATE,KK,.TRUE.,LDUAL,-999.99D0))THEN WRITE(*,'(/A)') '>>> Cannot compute GXG for measurements for:' WRITE(*,'(A/)') ' '//TRIM(CID)//'.'//TRIM(CEXT)//' <<<' KK=0 !## compute ghg(1)/glg(2) if ldual=.true. ELSE IF(.NOT.UTL_COMPUTE_GXG(C,IDATE,KK,.TRUE.,LDUAL,-999.99D0))THEN WRITE(*,'(/A)') '>>> Cannot compute GXG for computed heads for:' WRITE(*,'(A/)') ' '//TRIM(CID)//'.'//TRIM(CEXT)//' <<<' KK=0 ENDIF ENDIF !## only two values for gxg measurements IF(LDUAL)KK=2 ENDIF !## add observation DO K=1,KK II =II+1 !## random error not yet set IF(MEASURES(I)%IMCOL.LT.0)THEN IF(K.EQ.1)MSR%W(II)=WWGHG IF(K.EQ.2)MSR%W(II)=WWGLG ELSE MSR%W(II)=WW ENDIF !## save information for measurement MSR%CLABEL(II)=TRIM(CID) MSR%X(II)=X MSR%Y(II)=Y MSR%L(II)=ILAY !## dynamics IF(PEST%PE_TARGET(2).GT.0.0D0)THEN MSR%O(II)=DYN(1) MSR%C(II)=DYN(2) !## absolute values ELSE MSR%O(II)=M(K) MSR%C(II)=C(K) ENDIF MSR%D(II)=D MSR%IDATE(II)=IDATE(K) MSR%IPF(II)=I MSR%LOC(II)=J MSR%NS(II)=NSC IF(PEST%PE_TARGET(1).EQ.0.0D0.AND.PEST%PE_TARGET(2).GT.0.0D0)EXIT ENDDO ENDIF DEALLOCATE(TSNODATA,C,M,MCOPY,CCOPY,IDATE); CLOSE(IUIPFTXT) ENDDO DEALLOCATE(IPFSTRING) ENDIF CLOSE(IU) ENDDO MSR%NOBS=II; IF(MSR%NOBS.LE.0)THEN CALL IPEST_GLM_ERROR(IBATCH,'No measurements available within current spatial/temporal space.') ENDIF IERROR=1; IF(PBMAN%IIES.EQ.1)IERROR=IGRAD !## aggregate to mean values per gridcel IF(PBMAN%AVGOBS.EQ.1)THEN !## aggregate per layer ALLOCATE(NUMOBS(PRJIDF%NCOL,PRJIDF%NROW)) !## transient measurements IF(IEXT.NE.0)THEN NP=1; DO I=2,MSR%NOBS; IF(MSR%LOC(I).NE.MSR%LOC(I-1))NP=NP+1; ENDDO; ALLOCATE(AVGOBS(NP,PRJNPER)) DO ILAY=1,PRJNLAY !## aggregate per location NUMOBS=0.0D0; DO II=1,MSR%NOBS !## skip observations that have been processed already IF(MSR%W(II).LE.0.0D0)CYCLE !## skip observations in the wrong modellayer IF(MSR%L(II).NE.ILAY)CYCLE CALL IDFIROWICOL(PRJIDF,IROW,ICOL,MSR%X(II),MSR%Y(II)) !## count number of observations per location NUMOBS(ICOL,IROW)=NUMOBS(ICOL,IROW)+1.0D0 ENDDO !## process values per gridcell DO IROW=1,PRJIDF%NROW; DO ICOL=1,PRJIDF%NCOL !## nothing here IF(NUMOBS(ICOL,IROW).LE.0.0D0)CYCLE; AVGOBS=HNOFLOW; NP=0 DO II=1,MSR%NOBS !## already processed IF(MSR%LOC(II).LT.0)CYCLE !## skip observations in the wrong modellayer IF(MSR%L(II).NE.ILAY)CYCLE !## not in current gridcell CALL IDFIROWICOL(PRJIDF,JROW,JCOL,MSR%X(II),MSR%Y(II)); IF(IROW.NE.JROW.OR.JCOL.NE.ICOL)CYCLE !## get correct time-slots for current measurement NP=NP+1 DO JJ=II,MSR%NOBS !## different location IF(MSR%LOC(JJ).NE.ABS(MSR%LOC(II)).OR.MSR%IPF(JJ).NE.MSR%IPF(II))CYCLE !## make sure not to use this again MSR%LOC(JJ)=-1*ABS(MSR%LOC(JJ)) DO IPER=1,PRJNPER IF(SIM(IPER)%CDATE.NE.'STEADY-STATE')THEN WRITE(CDATE,'(I4.4,6I2.2)') SIM(IPER)%IYR,SIM(IPER)%IMH,SIM(IPER)%IDY,SIM(IPER)%IHR,SIM(IPER)%IMT,SIM(IPER)%ISC IF(TRIM(CDATE).NE.TRIM(VTOS(MSR%IDATE(JJ))))CYCLE ELSE IF(MSR%IDATE(JJ).NE.0D0)CYCLE ENDIF !## get measurement and quit looking AVGOBS(NP,IPER)=MSR%O(JJ); EXIT ENDDO ENDDO ENDDO !## get number of complete timesteps DO K=1,2 JJ=1; DO I=1,PRJNPER II=0; DO J=1,NP IF(AVGOBS(J,I).NE.HNOFLOW)THEN II=II+1 IF(K.EQ.2)THEN VAR(J,JJ)=AVGOBS(J,I) VARIDATE(JJ)=YMDHMSTOITIME(SIM(I)%IYR,SIM(I)%IMH,SIM(I)%IDY,SIM(I)%IHR,SIM(I)%IMT,SIM(I)%ISC) ENDIF ENDIF ENDDO IF(II.EQ.NP)JJ=JJ+1 IF(K.EQ.2.AND.JJ.GT.SIZE(VAR,2))EXIT ENDDO IF(K.EQ.1)THEN JJ=JJ-1; ALLOCATE(VAR(NP,JJ),VARIDATE(JJ)) ENDIF ENDDO !## compute variance/correlation DO I=1,NP H=0.0D0; DO J=1,SIZE(VAR,2) H=H+VAR(I,J) ENDDO H=H/REAL(SIZE(VAR,2),8) DO J=1,SIZE(VAR,2) VAR(I,J)=VAR(I,J)-H ENDDO ENDDO ALLOCATE(STD(SIZE(VAR,1)),COV(SIZE(VAR,1),SIZE(VAR,1)),COR(SIZE(VAR,1),SIZE(VAR,1))) COV=MATMUL(VAR,TRANSPOSE(VAR)) !## compute standard deviation DO I=1,SIZE(VAR,1) STD(I)=SQRT(COV(I,I)/(REAL(SIZE(VAR,2),8)-1.0D0)) ENDDO !## compute correlation COR=0.0D0; I=0 DO I=1,SIZE(VAR,1) DO J=1,SIZE(VAR,1) H=COV(I,I)*COV(J,J) IF(H.GT.0.0D0)THEN COR(I,J)=COV(I,J)/SQRT(H) ELSE COR(I,J)=0.0D0 ENDIF ENDDO ENDDO !## add minimal correlation on measurements H=MINVAL(ABS(COR)) DO I=1,MSR%NOBS IF(MSR%LOC(I).LT.0)THEN MSR%COR(I)=H MSR%LOC(I)=ABS(MSR%LOC(I)) !## if not in current selected dataset (overlap) set weighting to zero DO II=1,SIZE(VARIDATE) IF(MSR%IDATE(I).EQ.VARIDATE(II))EXIT ENDDO IF(II.GT.SIZE(VARIDATE))MSR%W(I)=0.0D0 ENDIF ENDDO DEALLOCATE(VAR,COV,STD,COR,VARIDATE) ENDDO ENDDO; ENDDO DEALLOCATE(AVGOBS) ENDIF ALLOCATE(AVGOBS(PRJIDF%NCOL,PRJIDF%NROW),AVGCOM(PRJIDF%NCOL,PRJIDF%NROW), & AVGW(PRJIDF%NCOL,PRJIDF%NROW), SUMW(PRJIDF%NCOL,PRJIDF%NROW)) !## apply average measurements/computed values DO IPER=1,PRJNPER DO ILAY=1,PRJNLAY AVGOBS=0.0D0; AVGCOM=0.0D0; NUMOBS=0.0D0; AVGW=0.0D0; SUMW=0.0D0 DO II=1,MSR%NOBS IF(MSR%L(II).NE.ILAY)CYCLE IF(SIM(IPER)%CDATE.NE.'STEADY-STATE')THEN WRITE(CDATE,'(I4.4,6I2.2)') SIM(IPER)%IYR,SIM(IPER)%IMH,SIM(IPER)%IDY,SIM(IPER)%IHR,SIM(IPER)%IMT,SIM(IPER)%ISC IF(TRIM(CDATE).NE.TRIM(VTOS(MSR%IDATE(II))))CYCLE ELSE IF(MSR%IDATE(II).NE.0D0)CYCLE ENDIF CALL IDFIROWICOL(PRJIDF,IROW,ICOL,MSR%X(II),MSR%Y(II)) IF(IROW.EQ.0.OR.ICOL.EQ.0)THEN WRITE(*,'(/A/)') 'Assignment of an observation yields a location outside the current modelling domain'; PAUSE; STOP ENDIF H=MSR%C(II) Z=MSR%O(II) W=MSR%W(II) AVGOBS(ICOL,IROW)=AVGOBS(ICOL,IROW)+Z*W AVGCOM(ICOL,IROW)=AVGCOM(ICOL,IROW)+H*W AVGW (ICOL,IROW)=AVGW (ICOL,IROW)+W*W SUMW (ICOL,IROW)=SUMW (ICOL,IROW)+W NUMOBS(ICOL,IROW)=NUMOBS(ICOL,IROW)+1.0D0 ENDDO !## do something with stdev van de metingen binnen 1 cel DO II=1,MSR%NOBS IF(MSR%L(II).NE.ILAY)CYCLE IF(SIM(IPER)%CDATE.NE.'STEADY-STATE')THEN WRITE(CDATE,'(I4.4,6I2.2)') SIM(IPER)%IYR,SIM(IPER)%IMH,SIM(IPER)%IDY,SIM(IPER)%IHR,SIM(IPER)%IMT,SIM(IPER)%ISC IF(TRIM(CDATE).NE.TRIM(VTOS(MSR%IDATE(II))))CYCLE ELSE IF(MSR%IDATE(II).NE.0D0)CYCLE ENDIF CALL IDFIROWICOL(PRJIDF,IROW,ICOL,MSR%X(II),MSR%Y(II)) Z=AVGOBS(ICOL,IROW) H=AVGCOM(ICOL,IROW) W=AVGW (ICOL,IROW) N=SUMW (ICOL,IROW) IF(N.NE.0.0D0)THEN Z=Z/N H=H/N W=W/N N=NUMOBS(ICOL,IROW) W=W/N ELSE W=0.0D0 ENDIF MSR%C(II)=H MSR%O(II)=Z MSR%W(II)=W ENDDO ENDDO ENDDO DEALLOCATE(AVGOBS,AVGCOM,AVGW,SUMW,NUMOBS) ENDIF ENDIF !## postprocess batch files and fill in observations IF(.NOT.PEST_GLM_BATCHFILES(2,CTYPE,IPARAM))THEN IPEST_GLM_GETJ_AVG=.FALSE.; RETURN ENDIF !## write summary of residuals (data) DO II=1,MSR%NOBS H=MSR%C(II) Z=MSR%O(II) H=REAL(INT(H*10000.0D0),8)/10000.0D0 Z=REAL(INT(Z*10000.0D0),8)/10000.0D0 !## random error not yet set CALL IPEST_GLM_ADDNOISE(II,SEED,LYN) !## add random error Z=Z+MSR%E(II) !## calculated - measured DHH=H-Z DRES=PEST%PE_DRES IF(DRES.EQ.0.0D0)DRES=MSR%D(II) !## exclude big residuals IF(DRES.LT.0.0D0)THEN IF(ABS(DHH).GT.ABS(DRES))DHH=0.0D0 ELSE IF(ABS(DHH).LT.DRES)THEN DHH=0.0D0 ELSE IF(DHH.GT. DRES)DHH=DHH-DRES IF(DHH.LT.-DRES)DHH=DHH+DRES ENDIF ENDIF !## save information for measurement SELECT CASE (TRIM(CTYPE)) CASE ('P') MSR%HG(IGRAD,II)=H CASE ('R') MSR%HG(IGRAD,II)=H MSR%IESC(IGRAD,II)=H CASE ('L') MSR%HL(IGRAD,II) =H MSR%DHL(IGRAD,II)=DHH END SELECT GF_H(II)=H GF_O(II)=Z !## compute part of objective function weighted DHW=MSR%W(II)*DHH**2.0D0 !## add to total objective function MSR%TJ=MSR%TJ+DHW IF(MSR%IDATE(II).EQ.0)THEN IF(IUPESTRESIDUAL.GT.0.AND.(TRIM(CTYPE).EQ.'L'.OR.TRIM(CTYPE).EQ.'R'))THEN WRITE(IUPESTRESIDUAL,'(2(F15.2,A1),I10,A1,6(F15.3,A1),I10,A1,A32)') & MSR%X(II),',',MSR%Y(II),',',MSR%L(II),',',MSR%O(II),',',MSR%C(II),',',MSR%C(II)-MSR%O(II),',',DHH,',', & MSR%W(II)*DHH,',',MSR%W(II),',',MSR%IPF(II),',',MSR%CLABEL(II) ENDIF ELSE IF(IUPESTRESIDUAL.GT.0.AND.(TRIM(CTYPE).EQ.'L'.OR.TRIM(CTYPE).EQ.'R'))THEN IF(PEST%PE_TARGET(1).GT.0.0D0)THEN WRITE(IUPESTRESIDUAL,'(2(F15.2,A1),I10,A1,8(F15.3,A1),I10,A1,A32,A1,I15)') & MSR%X(II),',',MSR%Y(II),',',MSR%L(II),',',MSR%W(II),',',MSR%O(II),',',MSR%C(II),',',DHH,',',MSR%COR(II),',',-999.0,',',-999.0,',', & MSR%NS(II),',',MSR%IPF(II),',',MSR%CLABEL(II),',',MSR%IDATE(II) ELSE WRITE(IUPESTRESIDUAL,'(2(F15.2,A1),I10,A1,8(F15.3,A1),I10,A1,A32,A1,I15)') & MSR%X(II),',',MSR%Y(II),',',MSR%L(II),',',MSR%W(II),',',-999.0,',',-999.0,',',-999.0,',',MSR%O(II),',',MSR%C(II),',',DHH,',', & MSR%NS(II),',',MSR%IPF(II),',',MSR%CLABEL(II),',',MSR%IDATE(II) ENDIF ENDIF ENDIF IF(IUPESTPRESIDUAL.GT.0)THEN WRITE(IUPESTPRESIDUAL,'(2(F15.2,A1),I10,A1,4(F15.3,A1),2(I10,A1),A32,A1,I14)') MSR%X(II),',',MSR%Y(II),',',MSR%L(II),',',MSR%O(II),',',MSR%C(II),',', & MSR%W(II),',',MSR%D(II),',',MSR%IPF(II),',',MSR%LOC(II),',',ADJUSTR(MSR%CLABEL(II)),',',MSR%IDATE(II) ENDIF ENDDO IF(ALLOCATED(ID1))DEALLOCATE(ID1); IF(ALLOCATED(ID2))DEALLOCATE(ID2) !## insert regularisation to objective function NP=0; DO I=1,SIZE(PEST%PARAM); IF(PEST%PARAM(I)%PACT.EQ.0.OR.PEST%PARAM(I)%PIGROUP.LE.0)CYCLE; NP=NP+1; ENDDO !## add parameter regularisation (not in gradient-simulations) IF(PEST%PE_REGULARISATION.EQ.1.AND.CTYPE.EQ.'L')THEN !## total weight observations WF=IPEST_GLM_BALANCEFACTOR() DO I=1,SIZE(PEST%PARAM) !## row !## skip inactive IF(PEST%PARAM(I)%PACT.EQ.0)CYCLE !## skip others parts of parameter IF(PEST%PARAM(I)%PIGROUP.LE.0)CYCLE IF(CTYPE.EQ.'L')THEN F1=PEST%PARAM(I)%LALPHA(IGRAD) ELSEIF(CTYPE.EQ.'P')THEN F1=PEST%PARAM(I)%GALPHA(IGRAD) ENDIF F2=PEST%PARAM(I)%PPRIOR IF(PEST%PARAM(I)%PLOG.EQ.1)THEN F1=EXP(F1) F2=EXP(F2) ELSEIF(PEST%PARAM(I)%PLOG.EQ.2)THEN F1=10.0D0**F1 F2=10.0D0**F2 ENDIF IF(F2.LT.F1)THEN F=F1/F2 ELSE F=F2/F1 ENDIF F=(F**2.0D0)-1.0D0 !## compute weighting from stdev WW=1.0D0/PEST%PARAM(I)%PARSTD**2.0D0 !## add balance-weighting WW=WW*WF MSR%RJ=MSR%RJ+F*WW ENDDO ENDIF MSR%TJ=MSR%TJ+MSR%RJ IF(IUPESTRESIDUAL.GT.0) CLOSE(IUPESTRESIDUAL); IUPESTRESIDUAL=0 IF(IUPESTPRESIDUAL.GT.0)CLOSE(IUPESTPRESIDUAL); IUPESTPRESIDUAL=0 IF(ALLOCATED(MEASURES))DEALLOCATE(MEASURES) IPEST_GLM_GETJ_AVG=.TRUE. END FUNCTION IPEST_GLM_GETJ_AVG !###==================================================================== SUBROUTINE PEST_GLM_ADD_BATCHFILES(IU,CTYPE,IPARAM) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: IU,IPARAM CHARACTER(LEN=*),INTENT(IN) :: CTYPE INTEGER :: I,J CHARACTER(LEN=256) :: COMMANDLINE,DIRNAME !## skip if not ipestp activated IF(PBMAN%IPESTP.EQ.0)RETURN !## run batchfile locally IF(.NOT.ASSOCIATED(PEST%B_BATCHFILE))RETURN CALL IOSDIRNAME(DIRNAME) !## goto folder WRITE(IU,'(/A)') 'CD /D '//TRIM(PBMAN%OUTPUT) !## compute what needs to be computed first DO I=1,SIZE(PEST%B_BATCHFILE) WRITE(IU,'(/50A1)') ('=',J=1,50) WRITE(IU,'(A )') 'REM Executing: '//TRIM(PEST%B_BATCHFILE(I)) WRITE(IU,'(50A1/)') ('=',J=1,50) !## add argument of current model; in case the batchfile needs it COMMANDLINE=TRIM(PEST%B_BATCHFILE(I))//' '//TRIM(CTYPE)//'#'//TRIM(VTOS(IPARAM)) WRITE(IU,'(A )') 'START /B /WAIT '//TRIM(COMMANDLINE) WRITE(IU,'(/A)') 'IF %ERRORLEVEL% NEQ 0 (ECHO AN ERROR WAS FOUND %EXIT /B %ERRORLEVEL%)' ENDDO !## got back WRITE(IU,'(/A)') 'CD /D '//TRIM(DIRNAME) END SUBROUTINE PEST_GLM_ADD_BATCHFILES !###==================================================================== LOGICAL FUNCTION PEST_GLM_BATCHFILES(IOPTION,CTYPE,IPARAM) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: IOPTION,IPARAM CHARACTER(LEN=*),INTENT(IN) :: CTYPE INTEGER :: I,II,J,N,NR,IUBAT,IIPF,IFLAGS LOGICAL :: LEX CHARACTER(LEN=256) :: DIRNAME,FNAME !## run batchfile locally PEST_GLM_BATCHFILES=.TRUE.; IF(.NOT.ASSOCIATED(PEST%B_BATCHFILE))RETURN PEST_GLM_BATCHFILES=.FALSE. IFLAGS=PROCCMDPROC+PROCBLOCKED N=SIZE(PEST%B_BATCHFILE); IF(.NOT.ASSOCIATED(PEST%B_NOBS))ALLOCATE(PEST%B_NOBS(N)); PEST%B_NOBS=0 CALL IOSDIRNAME(DIRNAME); CALL IOSDIRCHANGE(PBMAN%OUTPUT) WRITE(*,'(/A)') 'Currently iMOD operates BATCHFILES here: '//TRIM(PBMAN%OUTPUT) !## get maximum number of observations DO I=1,SIZE(PEST%B_BATCHFILE) !## skip this batchfile for processing IF(PEST%B_FRACTION(I).LE.0.0D0)CYCLE !## try to replace $run$ by current model FNAME=UTL_CAP(PEST%B_OUTFILE(I),'U') CALL UTL_SUBST(FNAME,'$RUN$','IPEST_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IPARAM))) INQUIRE(FILE=FNAME,EXIST=LEX) IF(.NOT.LEX)THEN WRITE(*,'(/A)') '>>> Cannot find: '//TRIM(FNAME) WRITE(*,'(A )') ' ioption= '//TRIM(VTOS(IOPTION)) WRITE(*,'(A )') ' ctype = '//TRIM(CTYPE) WRITE(*,'(A )') ' iparam = '//TRIM(VTOS(IPARAM))//' <<<' WRITE(*,'(A/)') ' Model is restarted <<<' RETURN ENDIF IUBAT=UTL_GETUNIT(); OPEN(IUBAT,FILE=FNAME,STATUS='OLD') READ(IUBAT,*) PEST%B_NOBS(I) WRITE(*,'(1X,A)') ' Found '//TRIM(VTOS(PEST%B_NOBS(I)))//' measurements from '//TRIM(FNAME) IF(IOPTION.EQ.2)THEN READ(IUBAT,*) NR WRITE(*,'(1X,A)') ' Reading them ...' !## check whether it is an IPF file IIPF=0; IF(INDEX(UTL_CAP(PEST%B_OUTFILE(I),'U'),'IPF').GT.0)THEN DO J=1,NR+1; READ(IUBAT,*); ENDDO IIPF=1 ENDIF J=MSR%NOBS DO II=1,PEST%B_NOBS(I) J=J+1 IF(IIPF.EQ.1)THEN READ(IUBAT,*) MSR%X(J),MSR%Y(J),MSR%CLABEL(J),MSR%O(J),MSR%W(J),MSR%C(J) ELSE READ(IUBAT,*) MSR%CLABEL(J),MSR%O(J),MSR%W(J),MSR%C(J) MSR%X(J)=-999.99D0; MSR%Y(J)=-999.99D0 ENDIF MSR%D(J) =0.0D0 MSR%L(J) =0 MSR%NS(J) =0.0D0 MSR%IDATE(J)=0 MSR%IPF(J) =I MSR%LOC(J) =II ENDDO MSR%NOBS=J ENDIF CLOSE(IUBAT) ENDDO CALL IOSDIRCHANGE(DIRNAME); WRITE(*,'(A/)') 'Currently iMOD operates MODELS here: '//TRIM(DIRNAME) PEST_GLM_BATCHFILES=.TRUE. END FUNCTION PEST_GLM_BATCHFILES !###==================================================================== SUBROUTINE IPEST_GLM_ADDNOISE(II,SEED,LYN) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(INOUT) :: SEED INTEGER,INTENT(IN) :: II REAL(KIND=DP_KIND) :: SIGMA,E LOGICAL,INTENT(INOUT) :: LYN CHARACTER(LEN=1) :: YN IF(PBMAN%IIES.EQ.1)THEN !## initiate measurement error IF(MSR%E(II).EQ.-999.99D0)THEN IF(MSR%W(II).NE.0.0D0)THEN !## compute stdev from weights again SIGMA=SQRT(1.0D0/MSR%W(II)) IF(SIGMA.GT.1.0D0.AND.LYN)THEN WRITE(*,'(/A$/)') 'Found STDEV for measurement of > 1.0 ('//TRIM(VTOS(SIGMA,'F',3))//'), Are you sure this correct ?' READ(*,'(A1)') YN; IF(YN.NE.'Y'.AND.YN.NE.'y')STOP; LYN=.FALSE. ENDIF E=MSR%E(II) CALL IPEST_NORMAL_MS_SAMPLE(0.0D0,SIGMA,SEED,E) MSR%E(II)=E ELSE MSR%E(II)=0.0D0 ENDIF ENDIF ELSE MSR%E(II)=0.0D0 ENDIF END SUBROUTINE IPEST_GLM_ADDNOISE !###==================================================================== REAL(KIND=DP_KIND) FUNCTION IPEST_GLM_BALANCEFACTOR() !###==================================================================== IMPLICIT NONE REAL(KIND=DP_KIND) :: WO,WP INTEGER :: I,NP,NO WO=0.0D0; NO=0; DO I=1,MSR%NOBS WO=WO+MSR%W(I) NO=NO+1 ENDDO !## total weight parameters WP=0.0D0; NP=0; DO I=1,SIZE(PEST%PARAM) !## row !## skip inactive IF(PEST%PARAM(I)%PACT.EQ.0)CYCLE !## skip others parts of parameter IF(PEST%PARAM(I)%PIGROUP.LE.0)CYCLE WP=WP+1.0D0/PEST%PARAM(I)%PARSTD**2.0D0 NP=NP+1 ENDDO !## determine balance factor IPEST_GLM_BALANCEFACTOR=WO/WP !## determine balance factor as ratio between number of measurements/parameters IPEST_GLM_BALANCEFACTOR=NO/NP ! IPEST_GLM_BALANCEFACTOR=1.0D0 END FUNCTION IPEST_GLM_BALANCEFACTOR !###==================================================================== LOGICAL FUNCTION IPEST_GLM_CREATE_RESIDUALSFILES(DIR,IGRAD,CTYPE,MNAME) !###==================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: DIR,CTYPE,MNAME INTEGER,INTENT(IN) :: IGRAD INTEGER :: IU,JU,KU,ISUB,NPER,NOBS,IOBS,IPER,N,MPER,IOS INTEGER,ALLOCATABLE,DIMENSION(:) :: NMES REAL(KIND=DP_KIND) :: CUMTT,DH,DHW,NODATA CHARACTER(LEN=52) :: TXT,CFOLDER CHARACTER(LEN=256) :: FNAME TYPE HTYPE INTEGER(KIND=DP_KIND) :: ITIME INTEGER :: ILAY REAL(KIND=DP_KIND) :: H,W,M,D END TYPE HTYPE TYPE(HTYPE),ALLOCATABLE,DIMENSION(:,:) :: H REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: X,Y CHARACTER(LEN=40),ALLOCATABLE,DIMENSION(:) :: ONAME IPEST_GLM_CREATE_RESIDUALSFILES=.FALSE. IF(TRIM(CTYPE).EQ.'')THEN CALL UTL_CREATEDIR(TRIM(DIR)//'\TIMESERIES') FNAME=TRIM(DIR)//'\TIMESERIES\OUTPUT_OBS.IPF_' ELSE CALL UTL_CREATEDIR(TRIM(DIR)//'\IPEST_'//CTYPE//'#'//TRIM(VTOS(IGRAD))//'\TIMESERIES') FNAME=TRIM(DIR)//'\IPEST_'//CTYPE//'#'//TRIM(VTOS(IGRAD))//'\TIMESERIES\OUTPUT_OBS.IPF_' ENDIF IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='REPLACE',ACTION='WRITE',FORM='FORMATTED'); IF(IU.EQ.0)RETURN WRITE(IU,'(A)') 'NaN1#' IF(PRJNPER.EQ.1)THEN WRITE(IU,'(A)') '10' ELSE WRITE(IU,'(A)') '7' ENDIF WRITE(IU,'(A)') 'X-COORDINATE' WRITE(IU,'(A)') 'Y-COORDINATE' WRITE(IU,'(A)') 'MODELLAYER' WRITE(IU,'(A)') 'OBSERVATION' WRITE(IU,'(A)') 'WEIGHT' IF(PRJNPER.EQ.1)THEN WRITE(IU,'(A)') 'COMPUTED_HEAD' WRITE(IU,'(A)') 'DIFFERENCE' WRITE(IU,'(A)') 'WEIGHTED_DIFFERENCE' ENDIF WRITE(IU,'(A)') 'ACCEPTABLE_ERROR' WRITE(IU,'(A)') 'ORG_MEASUREMENT_NUMBER' IF(PRJNPER.EQ.1)THEN WRITE(IU,'(A)') '0,TXT' ELSE WRITE(IU,'(A)') '4,TXT' ENDIF N=0; DO ISUB=1,PBMAN%NSUBMODEL IF(PBMAN%IFORMAT.EQ.3)CFOLDER='\GWF_'//TRIM(VTOS(ISUB))//'\MODELOUTPUT\' IF(PBMAN%IFORMAT.EQ.6)CFOLDER='\' !## skip zero IF(TRIM(CTYPE).EQ.'')THEN FNAME=TRIM(DIR)//TRIM(CFOLDER)//'OUTPUT_OBS.TXT' ELSE IF(PBMAN%IFORMAT.EQ.3)THEN FNAME=TRIM(DIR)//TRIM(CFOLDER)//'IPEST_'//TRIM(CTYPE)//'#'// & TRIM(VTOS(IGRAD))//'\OUTPUT_OBS_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IGRAD))//'.TXT' ELSE FNAME=TRIM(DIR)//TRIM(CFOLDER)//'IPEST_'//TRIM(CTYPE)//'#'// & TRIM(VTOS(IGRAD))//'\OUTPUT_OBS_'//TRIM(CTYPE)//'#'//TRIM(VTOS(IGRAD))//'._OS' ENDIF ENDIF JU=UTL_GETUNIT(); CALL OSD_OPEN(JU,FILE=FNAME,STATUS='OLD',ACTION='READ',FORM='FORMATTED'); IF(JU.EQ.0)THEN; CLOSE(IU); RETURN; ENDIF IF(PBMAN%IFORMAT.EQ.3)THEN IF(PBMAN%IFORMAT.EQ.3)CFOLDER='\GWF_'//TRIM(VTOS(ISUB))//'\' FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//TRIM(MNAME)//'.MES6' ELSE FNAME=TRIM(DIR)//TRIM(CFOLDER)//'MODELINPUT\'//TRIM(MNAME)//'.MES7' ENDIF KU=UTL_GETUNIT(); CALL OSD_OPEN(KU,FILE=FNAME,STATUS='OLD',ACTION='READ',FORM='FORMATTED') IF(KU.EQ.0)THEN; CLOSE(JU); CLOSE(IU); RETURN; ENDIF READ(KU,*) READ(KU,*) READ(KU,*) TXT,NPER READ(KU,*) TXT,NOBS READ(KU,*) TXT,NODATA IF(NOBS.GT.0)THEN ALLOCATE(H(NOBS,NPER),ONAME(NOBS),X(NOBS),Y(NOBS),NMES(NOBS)) READ(KU,*) !## read measurements from mes6 file (not used directly by modflow6 and/or seawat) DO IOBS=1,NOBS READ(KU,*) READ(KU,*) ONAME(IOBS),NMES(IOBS),X(IOBS),Y(IOBS) IF(NMES(IOBS).NE.NPER)THEN WRITE(*,'(/A/)') 'Something goes wrong (number of stressperiods in OBS6 and MET6 ne number of stresses'; STOP ENDIF READ(KU,*) DO IPER=1,NMES(IOBS) READ(KU,*,IOSTAT=IOS) H(IOBS,IPER)%ITIME,H(IOBS,IPER)%M,H(IOBS,IPER)%D,H(IOBS,IPER)%W,H(IOBS,IPER)%ILAY !## weight is always specified in stdev H(IOBS,IPER)%W=1.0D0/H(IOBS,IPER)%W**2.0D0 !## error probabily MF6 faces an issue running/closing, try to restart IF(IOS.NE.0)THEN CLOSE(IU); CLOSE(JU); CLOSE(KU); INQUIRE(UNIT=KU,NAME=FNAME) WRITE(*,'(1X,A)') '>>> ERROR READING FILE (IOS='//TRIM(VTOS(IOS))//'): '//TRIM(FNAME)//' <<<' RETURN ENDIF ENDDO ENDDO !## read computed heads from mf6 IF(PBMAN%IFORMAT.EQ.3)THEN READ(JU,*) DO IPER=1,NPER READ(JU,*,IOSTAT=IOS) CUMTT,(H(IOBS,IPER)%H,IOBS=1,NOBS) !## error probabily MF6 faces an issue running/closing, try to restart IF(IOS.NE.0)THEN CLOSE(IU); CLOSE(JU); CLOSE(KU); INQUIRE(UNIT=JU,NAME=FNAME) WRITE(*,'(1X,A)') '>>> ERROR READING FILE (IOS='//TRIM(VTOS(IOS))//'): '//TRIM(FNAME)//' <<<' RETURN ENDIF ENDDO !## read computed heads from seawat ELSE DO IOBS=1,NOBS DO IPER=1,NPER READ(JU,*,IOSTAT=IOS) H(IOBS,IPER)%H !## error probably MF6 faces an issue running/closing, try to restart IF(IOS.NE.0)THEN CLOSE(IU); CLOSE(JU); CLOSE(KU); INQUIRE(UNIT=JU,NAME=FNAME) WRITE(*,'(1X,A)') '>>> ERROR READING FILE (IOS='//TRIM(VTOS(IOS))//'): '//TRIM(FNAME)//' <<<' RETURN ENDIF ENDDO ENDDO ENDIF !## close files CLOSE(JU); CLOSE(KU) !## store them in aggregated IPF-file DO IOBS=1,NOBS N=N+1 IF(PRJNPER.EQ.1)THEN IPER=1 DH =H(IOBS,IPER)%H-H(IOBS,IPER)%M DHW=DH*H(IOBS,IPER)%W WRITE(IU,'(A)') TRIM(VTOS(X(IOBS),'F',3)) //','//TRIM(VTOS(Y(IOBS),'F',3)) //','//TRIM(VTOS(H(IOBS,IPER)%ILAY))//','// & TRIM(VTOS(H(IOBS,IPER)%M,'F',6))//','//TRIM(VTOS(H(IOBS,IPER)%W,'F',3))//','//TRIM(VTOS(H(IOBS,IPER)%H,'F',6))//','// & TRIM(VTOS(DH,'F',6)) //','//TRIM(VTOS(DHW,'F',6)) //','//TRIM(VTOS(H(IOBS,IPER)%D,'F',6))//','// & TRIM(VTOS(IOBS)) ELSE FNAME=ONAME(IOBS) !## find first correct value for w and d DO IPER=1,NPER; IF(H(IOBS,IPER)%M.NE.NODATA)EXIT; ENDDO IF(IPER.LE.NPER)THEN WRITE(IU,'(A)') TRIM(VTOS(X(IOBS),'F',3))//','//TRIM(VTOS(Y(IOBS),'F',3))//','//TRIM(VTOS(H(IOBS,1)%ILAY))//','// & TRIM(FNAME)//','//TRIM(VTOS(H(IOBS,IPER)%W,'F',3))//','//TRIM(VTOS(H(IOBS,IPER)%D,'F',3))//','//TRIM(VTOS(IOBS)) ELSE WRITE(IU,'(A)') TRIM(VTOS(X(IOBS),'F',3))//','//TRIM(VTOS(Y(IOBS),'F',3))//','//TRIM(VTOS(H(IOBS,1)%ILAY))//','// & TRIM(FNAME)//','//TRIM(VTOS(0.0D0,'F',3))//','//TRIM(VTOS(0.0D0,'F',3))//','//TRIM(VTOS(IOBS)) ENDIF IF(TRIM(CTYPE).EQ.'')THEN FNAME=TRIM(DIR)//'\TIMESERIES\'//TRIM(FNAME)//'.TXT' ELSE FNAME=TRIM(DIR)//'\IPEST_'//CTYPE//'#'//TRIM(VTOS(IGRAD))//'\TIMESERIES\'//TRIM(FNAME)//'.TXT' ENDIF JU=UTL_GETUNIT(); CALL OSD_OPEN(JU,FILE=FNAME,STATUS='REPLACE',ACTION='WRITE',FORM='FORMATTED'); IF(JU.EQ.0)RETURN MPER=0; DO IPER=1,NPER IF(H(IOBS,IPER)%ITIME.EQ.0)CYCLE; MPER=MPER+1 ENDDO WRITE(JU,'(A)') TRIM(VTOS(MPER)) WRITE(JU,'(A)') '3,1' WRITE(JU,'(A)') 'IDATE,-9999' WRITE(JU,'(A)') 'MEASURED,'//TRIM(VTOS(HNOFLOW,'G',6)) WRITE(JU,'(A)') 'COMPUTED,'//TRIM(VTOS(HNOFLOW,'G',6)) DO IPER=1,NPER IF(H(IOBS,IPER)%ITIME.EQ.0)CYCLE WRITE(JU,'(A)') TRIM(VTOS(H(IOBS,IPER)%ITIME))//','//TRIM(VTOS(H(IOBS,IPER)%M,'G',6))//','// & TRIM(VTOS(H(IOBS,IPER)%H,'G',6)) ENDDO CLOSE(JU) ENDIF ENDDO DEALLOCATE(H,ONAME,X,Y) ELSE !## close files CLOSE(JU); CLOSE(KU) ENDIF ENDDO CLOSE(IU) !## update number of records in ipf file IF(TRIM(CTYPE).EQ.'')THEN CALL UTL_MF2005_MAXNO(TRIM(DIR)//'\TIMESERIES\OUTPUT_OBS.IPF_',(/N/)) ELSE CALL UTL_MF2005_MAXNO(TRIM(DIR)//'\IPEST_'//CTYPE//'#'//TRIM(VTOS(IGRAD))//'\TIMESERIES\OUTPUT_OBS.IPF_',(/N/)) ENDIF IPEST_GLM_CREATE_RESIDUALSFILES=.TRUE. END FUNCTION IPEST_GLM_CREATE_RESIDUALSFILES !#####================================================================= SUBROUTINE IPEST_GLM_PROGRESS(ITER,JGRAD,IGRAD,CTYPE,LAMBDA) !#####================================================================= IMPLICIT NONE CHARACTER(LEN=1),INTENT(IN) :: CTYPE INTEGER,INTENT(IN) :: JGRAD,IGRAD,ITER REAL(KIND=DP_KIND),INTENT(IN) :: LAMBDA REAL(KIND=DP_KIND) :: X1,X2,ETM,ETP INTEGER :: CLOCK_RATE IF(CTYPE.EQ.'P')THEN !## present parameter value X1=PEST%PARAM(IGRAD)%GALPHA(JGRAD) IF(PEST%PARAM(IGRAD)%PLOG.EQ.1)X1=EXP(X1) IF(PEST%PARAM(IGRAD)%PLOG.EQ.2)X1=10.0D0**X1 X2=PEST%PARAM(IGRAD)%ALPHA(2) IF(PEST%PARAM(IGRAD)%PLOG.EQ.1)X2=EXP(X2) IF(PEST%PARAM(IGRAD)%PLOG.EQ.2)X2=10.0D0**X2 ELSEIF(CTYPE.EQ.'R')THEN !## present lambda X1=LAMBDA ELSEIF(CTYPE.EQ.'L')THEN !## present lambda X1=LAMBDA*PBMAN%LAMBDA_TEST(JGRAD) ENDIF CALL SYSTEM_CLOCK(ETIME(JGRAD),CLOCK_RATE) ETM=REAL(ETIME(JGRAD )-STIME(JGRAD ),8)/REAL(CLOCK_RATE,8) ETP=REAL(ETIMEJ(JGRAD)-STIMEJ(JGRAD),8)/REAL(CLOCK_RATE,8) !## sensitivity IF(MSR%NOBS.GT.0)THEN IF(CTYPE.EQ.'P')THEN IF(INSENS(JGRAD).EQ.0)THEN WRITE(IUPESTPROGRESS,'(I5,A15,4F15.3,30X,2F15.3)') IGRAD,PEST%PARAM(IGRAD)%ACRONYM,X1,MSR%TJ,MSR%TJ-MSR%PJ,X2,ETM,ETP ELSE WRITE(IUPESTPROGRESS,'(I5,A15,F15.3,2A15,F15.5,A15)') IGRAD,PEST%PARAM(IGRAD)%ACRONYM,X2,'NaN','NaN',X2,'> insensitive <' ENDIF !## lambda testing ELSEIF(CTYPE.EQ.'L')THEN IF(ITER.EQ.0)THEN WRITE(IUPESTPROGRESS,'(I5,30X,F15.3,30X,4F15.3)') 0,MSR%TJ,MSR%GOF(JGRAD),MSR%NSC(JGRAD),ETM,ETP ELSE WRITE(IUPESTPROGRESS,'(I5,A15,F15.3,2F15.3,15X,4F15.3)') IGRAD,'LAMBDA'//TRIM(VTOS(JGRAD)),X1,MSR%TJ,MSR%TJ-MSR%PJ, & MSR%GOF(JGRAD),MSR%NSC(JGRAD),ETM,ETP ENDIF !## realization ELSEIF(CTYPE.EQ.'R')THEN IF(ITER.EQ.0)THEN WRITE(IUPESTPROGRESS,'(I5,A15,2F15.3,30X,4F15.3)') IGRAD,'REALS'//TRIM(VTOS(JGRAD)),X1,MSR%TJ,MSR%GOF(JGRAD),MSR%NSC(JGRAD),ETM,ETP ELSE WRITE(IUPESTPROGRESS,'(I5,A15,2F15.3,F15.3,15X,4F15.3)') IGRAD,'REALS'//TRIM(VTOS(JGRAD)),X1,MSR%TJ,MSR%TJ-JE(IGRAD,ITER-1),& MSR%GOF(JGRAD),MSR%NSC(JGRAD),ETM,ETP ENDIF JE(IGRAD,ITER)=MSR%TJ WRITE(IUPESTOUT,'(I5,7X,F15.3)') ITER,MSR%TJ ENDIF ENDIF END SUBROUTINE IPEST_GLM_PROGRESS !###==================================================================== LOGICAL FUNCTION IPEST_GLM_ALLOCATEMSR(NDIM1,NDIM2,NDIM3,MEASURES) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: NDIM1,NDIM2,NDIM3 TYPE(PSTMEASURE),DIMENSION(:),INTENT(IN) :: MEASURES INTEGER :: I,N1,N2,N,M,O,IOS,IU IPEST_GLM_ALLOCATEMSR=.FALSE. CALL IPEST_GLM_DEALLOCATEMSR() IF(NDIM1.EQ.0)THEN !## open files and get total maximum number of observations N=0; DO I=1,SIZE(MEASURES) IU=UTL_GETUNIT(); OPEN(IU,FILE=MEASURES(I)%IPFNAME,STATUS='OLD',ACTION='READ') READ(IU,*) M; CLOSE(IU) !## gxg takes less memory, only 2 per observation IF(MEASURES(I)%IMCOL.LT.0)THEN N=N+2*M ELSE N=N+M*PRJNPER ENDIF ENDDO !## get dimension of possible maximal number of observations M=N !## number of line-searches (-) and number of gradients-simulations (+) N1=0; IF(ALLOCATED(RNL))N1=SIZE(RNL) N2=0; IF(ALLOCATED(RNG))N2=SIZE(RNG) ELSE !## number of observations M =NDIM1 !## number of linesearches N1=NDIM2 !## number of parameters N2=NDIM3 ENDIF IF(M.GT.0)THEN WRITE(*,'(/I10,A)') M ,' Maximal number of observations (from IPF)' ELSE WRITE(*,'(/I10,A)') ABS(M) ,' Maximal number of observations (from LOG_PESTP_RESIDUAL_*.TXT)' ENDIF !## add measurements from batchfiles IF(ASSOCIATED(PEST%B_NOBS).AND.M.GE.0)THEN M=M+SUM(PEST%B_NOBS) WRITE(*,'(I10,A)') SUM(PEST%B_NOBS) ,' Maximal number of observations (from Batchfiles)' WRITE(*,'(60A1)') ('-',I=1,60) WRITE(*,'(I10,A/)') M ,' Maximal number of observations (totally)' ENDIF M=ABS(M) O=ABS(PEST%PE_MXITER) WRITE(*,'(I10,A)') N1,' Maximal number of linesearches' WRITE(*,'(I10,A)') N2,' Maximal number of parameters' WRITE(*,'(I10,A)') O ,' Maximal number of iterations' WRITE(*,*) ALLOCATE(LAMBDAS(O),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR LAMBDAS'; RETURN; ENDIF ALLOCATE(MSR%GOF_H(0:O),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%GOF_H'; RETURN; ENDIF ALLOCATE(MSR%NSC_H(0:O),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%NSC_H'; RETURN; ENDIF ALLOCATE(MSR%TJ_H(0:O),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%TJ_H'; RETURN; ENDIF ALLOCATE(MSR%RJ_H(0:O),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%RJ_H'; RETURN; ENDIF IF(N1.GT.0)THEN ALLOCATE(MSR%GOF(N1),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%GOF'; RETURN; ENDIF ALLOCATE(MSR%NSC(N1),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%NSC'; RETURN; ENDIF ALLOCATE(MSR%DHL_J(N1),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%DHL_J'; RETURN; ENDIF IF(M.GT.0)THEN ALLOCATE(MSR%DHL(0:N1,M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%DHL'; RETURN; ENDIF ALLOCATE(MSR%HL(0:N1,M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%HL'; RETURN; ENDIF ! ALLOCATE(MSR%DHL(0:N1,M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%DHL'; RETURN; ENDIF ENDIF ENDIF IF(PBMAN%IIES.EQ.1)THEN ALLOCATE(MSR%GOF(N2),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%GOF'; RETURN; ENDIF ALLOCATE(MSR%NSC(N2),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%NSC'; RETURN; ENDIF ALLOCATE(MSR%DHL_J(N2),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%DHL_J'; RETURN; ENDIF ALLOCATE(MSR%IESC(N2,M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%IESC'; RETURN; ENDIF ALLOCATE(MSR%DHG(N2,M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%DHG'; RETURN; ENDIF ENDIF IF(N2.GT.0)THEN ALLOCATE(MSR%DHG_J(N2),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%DHG_J'; RETURN; ENDIF IF(M.GT.0)THEN ALLOCATE(MSR%HG(N2,M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%HG'; RETURN; ENDIF ENDIF ENDIF IF(M.GT.0)THEN ALLOCATE(MSR%W(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%W'; RETURN; ENDIF ALLOCATE(MSR%X(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%X'; RETURN; ENDIF ALLOCATE(MSR%Y(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%Y'; RETURN; ENDIF ALLOCATE(MSR%O(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%O'; RETURN; ENDIF ALLOCATE(MSR%C(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%C'; RETURN; ENDIF ALLOCATE(MSR%D(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%D'; RETURN; ENDIF ALLOCATE(MSR%L(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%L'; RETURN; ENDIF ALLOCATE(MSR%NS(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%NS'; RETURN; ENDIF ALLOCATE(MSR%COR(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%COR'; RETURN; ENDIF; MSR%COR=-999.0D0 ALLOCATE(MSR%IDATE(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%IDATE'; RETURN; ENDIF ALLOCATE(MSR%IPF(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%IPF'; RETURN; ENDIF ALLOCATE(MSR%LOC(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%LOC'; RETURN; ENDIF ALLOCATE(MSR%CLABEL(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR MSR%CLABEL'; RETURN; ENDIF ALLOCATE(GF_H(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR GF_H'; RETURN; ENDIF ALLOCATE(GF_O(M),STAT=IOS); IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'CANNOT ALLOCATE MEMORY FOR GF_O'; RETURN; ENDIF ENDIF !## allocate stochastic error per measurement ALLOCATE(MSR%E(M)); MSR%E=-999.99D0 IPEST_GLM_ALLOCATEMSR=.TRUE. END FUNCTION IPEST_GLM_ALLOCATEMSR !###==================================================================== SUBROUTINE IPEST_GLM_DEALLOCATEMSR() !###==================================================================== IMPLICIT NONE IF(ASSOCIATED(MSR%TJ_H)) DEALLOCATE(MSR%TJ_H) IF(ASSOCIATED(MSR%RJ_H)) DEALLOCATE(MSR%RJ_H) IF(ASSOCIATED(MSR%GOF_H)) DEALLOCATE(MSR%GOF_H) IF(ASSOCIATED(MSR%NSC_H)) DEALLOCATE(MSR%NSC_H) IF(ASSOCIATED(MSR%GOF)) DEALLOCATE(MSR%GOF) IF(ASSOCIATED(MSR%NSC)) DEALLOCATE(MSR%NSC) IF(ASSOCIATED(MSR%DHL_J)) DEALLOCATE(MSR%DHL_J) IF(ASSOCIATED(MSR%DHG_J)) DEALLOCATE(MSR%DHG_J) IF(ASSOCIATED(MSR%HL)) DEALLOCATE(MSR%HL) IF(ASSOCIATED(MSR%DHL)) DEALLOCATE(MSR%DHL) IF(ASSOCIATED(MSR%HG)) DEALLOCATE(MSR%HG) IF(ASSOCIATED(MSR%W )) DEALLOCATE(MSR%W) IF(ASSOCIATED(MSR%X )) DEALLOCATE(MSR%X) IF(ASSOCIATED(MSR%Y )) DEALLOCATE(MSR%Y) IF(ASSOCIATED(MSR%O )) DEALLOCATE(MSR%O) IF(ASSOCIATED(MSR%C )) DEALLOCATE(MSR%C) IF(ASSOCIATED(MSR%IESC)) DEALLOCATE(MSR%IESC) IF(ASSOCIATED(MSR%D )) DEALLOCATE(MSR%D) IF(ASSOCIATED(MSR%E )) DEALLOCATE(MSR%E) IF(ASSOCIATED(MSR%L )) DEALLOCATE(MSR%L) IF(ASSOCIATED(MSR%NS)) DEALLOCATE(MSR%NS) IF(ASSOCIATED(MSR%COR)) DEALLOCATE(MSR%COR) IF(ASSOCIATED(MSR%IDATE)) DEALLOCATE(MSR%IDATE) IF(ASSOCIATED(MSR%IPF)) DEALLOCATE(MSR%IPF) IF(ASSOCIATED(MSR%LOC)) DEALLOCATE(MSR%LOC) IF(ASSOCIATED(MSR%CLABEL))DEALLOCATE(MSR%CLABEL) IF(ALLOCATED(GF_H)) DEALLOCATE(GF_H) IF(ALLOCATED(GF_O)) DEALLOCATE(GF_O) IF(ALLOCATED(LAMBDAS)) DEALLOCATE(LAMBDAS) END SUBROUTINE IPEST_GLM_DEALLOCATEMSR !###==================================================================== SUBROUTINE IPEST_LUDECOMP_DBL(AA,IDX,N,ISING) !###==================================================================== IMPLICIT NONE ! INTEGER,PARAMETER :: NMAX=2000 REAL(KIND=DP_KIND),PARAMETER :: TINY=1.0D-20 INTEGER,INTENT(IN) :: N INTEGER,INTENT(OUT) :: ISING INTEGER,DIMENSION(N),INTENT(OUT) :: IDX REAL(KIND=DP_KIND),DIMENSION(N,N),INTENT(INOUT) :: AA REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: VV REAL(KIND=DP_KIND) :: AAMAX,DUM,SUM INTEGER :: I,IMAX,J,K ALLOCATE(VV(N)); VV=0.0D0 DO I=1,N IDX(I)=0 END DO ISING=0 DO I=1,N AAMAX=0.0D0 DO J=1,N IF(ABS(AA(I,J)).GT.AAMAX)AAMAX=ABS(AA(I,J)) ENDDO IF(AAMAX.EQ.0.0D0)THEN WRITE(*,*) 'Matrix is singular' ISING=1 RETURN ENDIF VV(I)=1.0D0/AAMAX ENDDO DO J=1,N DO I=1,J-1 SUM=AA(I,J) DO K=1,I-1 SUM=SUM-AA(I,K)*AA(K,J) ENDDO AA(I,J)=SUM ENDDO AAMAX=0.0D0 DO I=J,N SUM=AA(I,J) DO K=1,J-1 SUM=SUM-AA(I,K)*AA(K,J) ENDDO AA(I,J)=SUM DUM=VV(I)*ABS(SUM) IF(DUM.GE.AAMAX)THEN IMAX=I AAMAX=DUM ENDIF ENDDO IF(J.NE.IMAX)THEN DO K=1,N DUM=AA(IMAX,K) AA(IMAX,K)=AA(J,K) AA(J,K)=DUM ENDDO VV(IMAX)=VV(J) ENDIF IDX(J)=IMAX IF(AA(J,J).EQ.0.0D0)AA(J,J)=TINY IF(J.NE.N)THEN DUM=1.0D0/AA(J,J) DO I=J+1,N AA(I,J)=AA(I,J)*DUM ENDDO ENDIF ENDDO DEALLOCATE(VV) END SUBROUTINE IPEST_LUDECOMP_DBL !###==================================================================== SUBROUTINE IPEST_LUBACKSUB_DBL(AA,IDX,BB,N) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: N REAL(KIND=DP_KIND),DIMENSION(N,N),INTENT(IN) :: AA REAL(KIND=DP_KIND),DIMENSION(N),INTENT(INOUT) :: BB INTEGER,DIMENSION(N),INTENT(IN) :: IDX INTEGER :: I,II,J,LL REAL(KIND=DP_KIND) :: SUM II=0 DO I=1,N LL=IDX(I) SUM=BB(LL) BB(LL)=BB(I) IF(II.NE.0)THEN DO J=II,I-1 SUM=SUM-AA(I,J)*BB(J) ENDDO ELSE IF(SUM.NE.0.0D0)THEN II=I ENDIF BB(I)=SUM ENDDO DO I=N,1,-1 SUM=BB(I) DO J=I+1,N SUM=SUM-AA(I,J)*BB(J) ENDDO BB(I)=SUM/AA(I,I) ENDDO END SUBROUTINE IPEST_LUBACKSUB_DBL !###==================================================================== SUBROUTINE IPEST_ECHELON_DBL(AO,BO,NROW,NCOL) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: NROW,NCOL REAL(KIND=DP_KIND),DIMENSION(NCOL,NROW),INTENT(IN) :: AO INTEGER,DIMENSION(NROW),INTENT(OUT) :: BO REAL(KIND=DP_KIND),PARAMETER :: TINY=1.0D-10 REAL(KIND=DP_KIND),DIMENSION(:,:),ALLOCATABLE :: A INTEGER,DIMENSION(:),ALLOCATABLE :: B INTEGER :: IROW,JROW,KROW,ICOL,PCOL REAL(KIND=DP_KIND) :: AMAX,ATMP,F !## copy matrix ALLOCATE(A(NCOL,NROW),B(NROW)); A=AO PCOL=0 MLOOP: DO IROW=1,NROW !## assume next column to be pivot column PCOL=PCOL+1; IF(PCOL.GT.NCOL)EXIT !## find pivot icol PLOOP: DO DO JROW=IROW+1,NROW IF(ABS(A(PCOL,JROW)).GT.TINY)EXIT PLOOP ENDDO PCOL=PCOL+1 !## finished IF(PCOL.GT.NCOL)EXIT MLOOP ENDDO PLOOP !## find row with largest pivot value KROW=IROW; AMAX=0.0D0; DO JROW=IROW,NROW !+1,NROW IF(ABS(A(PCOL,JROW)).GT.ABS(AMAX))THEN; AMAX=A(PCOL,JROW); KROW=JROW; ENDIF ENDDO !## interchange rows IF(KROW.NE.IROW)THEN DO ICOL=PCOL,NCOL ATMP =A(ICOL,IROW) A(ICOL,IROW)=A(ICOL,KROW) A(ICOL,KROW)=ATMP ENDDO ENDIF !## reduce all rows using the pivot row AMAX=A(PCOL,IROW) DO JROW=IROW+1,NROW F=A(PCOL,JROW)/AMAX DO ICOL=PCOL,NCOL A(ICOL,JROW)=A(ICOL,JROW)-A(ICOL,IROW)*F ENDDO ENDDO ENDDO MLOOP ! DO IROW=1,NROW ! WRITE(*,'(99F10.2)') (A(ICOL,IROW),ICOL=1,NCOL) ! ENDDO !## determine over/bad dimension of matrix B=0; DO IROW=1,NROW DO ICOL=1,NCOL; IF(ABS(A(ICOL,IROW)).GT.TINY)THEN; B(IROW)=ICOL; EXIT; ENDIF; ENDDO ENDDO BO=0; DO IROW=1,NROW; ICOL=B(IROW); IF(ICOL.GT.0)BO(ICOL)=1; ENDDO DEALLOCATE(A,B) END SUBROUTINE IPEST_ECHELON_DBL END MODULE MOD_IPEST_GLM