Re: memory leak help!

On Jun 18, 3:38 pm, nos...@xxxxxxxxxxxxx (Richard Maine) wrote:
julesr <julesrei...@xxxxxxxxxxx> wrote:
Thanks Richard - would you suggest posting my whole program? i guess
that way it can be scanned for potential problems?

That's about all I can suggest at the moment. Or a link to the code as
an alternative. I can't guarantee that will elicit answers, but at least
it has a chance.

--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain

Righto, here's my code! I apologise profusely for slapping my whole
program in here, but I guess it's the only way...
Although i'm not supposed to include my own suppositions ;), I am
quite sure that the problem occurs when 'similarity' is called within
a loop from the main program (or similarity's children).

Anyhow, if anyone has the time to have a look and has any bright
ideas, I'd greatly appreciate it!!

thanks Julian

PROGRAM SimilaritySCE_simhyd


IMPLICIT NONE

! FUNCTION TYPES
REAL :: similarity ! function calculates
weights and itself calls RR Model
REAL :: ran1

! Problem-specific inputs:
CHARACTER*70 :: fl_in_Att,fl_in_ID
REAL :: Att(10,7)

CHARACTER*6 :: ids(10)
INTEGER :: ngauge
INTEGER :: num_mod_pars=9

CHARACTER*70 :: fl_in_prob_vars='D:\SENSITIVITY\code\fortran
\input_files\Similarity_SCE_inputs.txt'
CHARACTER*70 :: fl_in_bounds='D:\SENSITIVITY\code\fortran\input_files
\Similarity_SCE_bounds.txt'
CHARACTER*70 :: fl_out_loops='D:\SENSITIVITY\code\fortran\output_files
\Similarity_SCE_loops.txt'
CHARACTER*70 :: fl_out_term='D:\SENSITIVITY\code\fortran\output_files
\Similarity_SCE_term.txt'
CHARACTER*70 :: fl_in_pars='D:\SENSITIVITY\code\fortran\input_files
\D_simhyd_OptPars.dat'

REAL :: modpar(10,9)

REAL :: obj(10)
INTEGER :: nopt=10

REAL :: x0(10),bl(10),bu(10),bound(10)
INTEGER :: iseed=0,iseed1
INTEGER :: iniflg=0
INTEGER :: ngs=7
INTEGER :: npg
INTEGER :: nps
INTEGER :: nspl
INTEGER :: mings
INTEGER :: maxn=20000
REAL :: peps=.01
INTEGER :: kstop=6
REAL :: pcento=0.1

REAL :: x(147,10),cx(21,10),s(11,10)
REAL :: xf(147),cf(21),sf(11)
REAL :: xx(10)
REAL :: bestx(10),worstx(10)
REAL :: bestf,worstf
REAL :: unit1(10),xnstd(10)
INTEGER :: lcs(11),lcs1(11)
REAL, DIMENSION(20) :: criter
REAL :: lpos,ipt,nopt1,icall=1,a,ipcnvg,gnrng,rand,denomi,timeou
INTEGER :: igs,ngs1
INTEGER :: ngs2,nopt2
REAL :: ipr,xname
INTEGER :: npt,npt1,l
INTEGER :: Qdays_all(10),Pdays_all(10)
INTEGER :: nlines
CHARACTER*70 :: filrain,filrun
! Loop variables
INTEGER :: i,j,k1,k2,k3,k,loop
INTEGER :: ierror ! IOSTAT error
INTEGER :: nloop
INTEGER :: arrsize
INTEGER :: ncont

OPEN (9,FILE=fl_in_prob_vars,STATUS='old',IOSTAT=ierror)
READ (9,109,IOSTAT=ierror) fl_in_Att,fl_in_ID,ngauge
CLOSE (9)

OPEN (20,FILE=fl_in_pars,STATUS='old',IOSTAT=ierror)
if (ierror.ne.0) then
WRITE (*,*) 'IOSTAT, model parameters =,',ierror
end if
do i=1,ngauge
READ (20,99,IOSTAT=ierror) obj(i),modpar(i,1),modpar(i,2),modpar(i,
3), &
modpar(i,4),modpar(i,5),modpar(i,6),modpar(i,7),modpar(i,
8),modpar(i,9)
end do
99 format (2x,f14.12,16x,9(2x,f14.12))
CLOSE (20)

OPEN (10,FILE=fl_in_Att,STATUS='old',IOSTAT=ierror)
if (ierror.ne.0) then
WRITE (*,*) 'IOSTAT, attributes =',ierror
end if
OPEN (13,FILE=fl_in_ID,STATUS='old',IOSTAT=ierror)
if (ierror.ne.0) then
WRITE (*,*) 'IOSTAT, id file =',ierror
end if
do i=1,ngauge
READ (10,101,IOSTAT=ierror) Att(i,1:7)
READ (13,113,IOSTAT=ierror) ids(i)
end do
CLOSE (10)
CLOSE (13)

do i=1,ngauge
WRITE (filrain,'(a)') 'D:\SENSITIVITY\daily_precip\fortran\' //
ids(i) // 'precip.dat'
WRITE (filrun,'(a)') 'D:\SENSITIVITY\daily_flow\fortran\' //
ids(i) // 'flow.dat'

CALL filesize(nlines,filrun)
Qdays_all(i)=nlines
CALL filesize(nlines,filrain)
Pdays_all(i)=nlines
END do
OPEN (14,FILE=fl_out_loops,STATUS='new',IOSTAT=ierror)
OPEN (15,FILE=fl_out_term,STATUS='new',IOSTAT=ierror)

! Begin SCE....
iseed1=-ABS(iseed)
! Initialise SCE parameters
do i=1,20
criter(i)=0
end do
nloop=0
loop=0
igs=0
nopt1=8
if (nopt.lt.8) then
nopt1=nopt
end if
nopt2=12
if (nopt2.gt.nopt) then
nopt2=nopt
end if

npg=2*nopt+1
nps=nopt+1
nspl=npg
mings=ngs
npt=npg*ngs

ngs1=ngs
npt1=npt

! Initial guess
OPEN (11,FILE=fl_in_bounds,STATUS='old',IOSTAT=ierror)
DO i=1,nopt
READ (11,100,IOSTAT=ierror) bl(i),bu(i),x0(i)
END DO
CLOSE (11)


!WRITE (*,*) 'Beginning evolution Loop number ',nloop

! Compute the bounds for the parameters being optimised
do i=1,nopt
bound(i)=bu(i)-bl(i)
unit1(i)=1
end do

! Generate npt random points distributed uniformly in the parameter
space,
! and compute the corresponding function values
WRITE (14,'(a)') 'Loop Call F X1 X2 X3 X4
X5 X6 X7 X8 X9 X10'
do i=1,npt
call getpnt(nopt,1,iseed1,xx,bl,bu,unit1,bl)
do j=1,nopt
x(i,j)=xx(j)
end do
ncont=INT(xx(8))

xf(i)=similarity(nopt,xx,Att,ngauge,ids,num_mod_pars,modpar,obj,Qdays_all,Pdays_all,ncont)
icall=icall+1
if (icall.ge.maxn) then
npt1=i
GOTO 45
end if
end do

! if iniflg==1, set x(1,:)=x0
if (iniflg.eq.1) then
do j=1,nopt
x(1,j)=x0(j)
end do
ncont=INT(x0(8))

xf(1)=similarity(nopt,x0,Att,ngauge,ids,num_mod_pars,modpar,obj,Qdays_all,Pdays_all,ncont)
icall=icall+1
end if
! Arrange points in order of increasing function value

45 call hpsort1(npt1,nopt,xf,x)

! Record best and worst points
bestx(1:nopt)=x(1,1:nopt)
worstx(1:nopt)=x(npt1,1:nopt)

bestf=xf(1)
worstf=xf(npt1)

! Compute the parameter range for the initial population
call parstt(npt1,nopt,x,xnstd,bound,gnrng,ipcnvg)

! Print the results for the initial population
WRITE (14,114) nloop,icall,bestf,bestx(1:nopt)

if (icall.ge.maxn) then
GOTO 9000
end if
if (ipcnvg.eq.1) then
GOTO 9200
end if

! Begin the main loop
1000 continue

nloop=nloop+1

do 3000 igs=1,ngs1

! Assign points into complexes

do k1=1,npg
k2=(k1-1)*ngs1+igs
cx(k1,1:nopt)=x(k2,1:nopt)
cf(k1) = xf(k2)
END do

! BEGIN INNER LOOP - RANDOM SELECTION OF SUB-COMPLEXES
---------------
do loop=1,nspl

! CHOOSE A SUB-COMPLEX (nps points) ACCORDING TO A LINEAR
! PROBABILITY DISTRIBUTION
if (nps.eq.npg) then
do k=1,nps
lcs(k)=k
end do
go to 85
end if

rand=ran1(iseed1)
lcs(1)=1
do k=2,nps
60 rand=ran1(iseed1)
lpos=1+int(npg+0.5-sqrt((npg+0.5)**2-npg*(npg+1)*rand))
do k3=1,k-1
if (lpos.eq.lcs(k3)) go to 60
end do
lcs(k)=lpos
end do

! ARRANGE THE SUB-COMPLEX IN ORDER OF INCEASING FUNCTION VALUE
call hpsort3(nps,lcs,lcs1)

! CREATE THE SUB-COMPLEX ARRAYS
85 do k=1,nps
s(k,1:nopt)=cx(lcs1(k),1:nopt)

sf(k)=cf(lcs1(k))
end do

! USE THE SUB-COMPLEX TO GENERATE NEW POINT(S)

call
cce(nopt,nps,s,sf,bl,bu,xnstd,icall,maxn,iseed1,Att,ngauge,ids,num_mod_pars,modpar,obj,Qdays_all,Pdays_all)
! IF THE SUB-COMPLEX IS ACCEPTED, REPLACE THE NEW SUB-COMPLEX
! INTO THE COMPLEX
do k=1,nps
cx(lcs1(k),1:nopt)=s(k,1:nopt)
cf(lcs1(k))=sf(k)
end do

! SORT THE POINTS
call hpsort1(npg,nopt,cf,cx)

! IF MAXIMUM NUMBER OF RUNS EXCEEDED, BREAK OUT OF THE LOOP
if (icall.ge.maxn) go to 2222

! End of inner loop for competitive evolution of simplexes loop=1,nspl
end do



! END OF INNER LOOP ------------
2000 continue
2222 continue

! REPLACE THE NEW COMPLEX INTO ORIGINAL ARRAY x(.,.)
do k3=1,npg
k2=(k3-1)*ngs1+igs
x(k2,1:nopt)=cx(k3,1:nopt)
xf(k2)=cf(k3)
end do
if (icall.ge.maxn) go to 3333

! END LOOP ON COMPLEXES
3000 continue

! RE-SORT THE POINTS
3333 continue
call hpsort1(npt1,nopt,xf,x)

! RECORD THE BEST AND WORST POINTS
bestx(1:nopt)=x(1,1:nopt)
worstx(1:nopt)=x(npt1,1:nopt)
bestf=xf(1)
worstf=xf(npt1)
WRITE (14,114) nloop,icall,bestf,bestx(1:nopt)
! TEST THE POPULATION FOR PARAMETER CONVERGENCE
call parstt(npt1,nopt,x,xnstd,bound,gnrng,ipcnvg)

! PRINT THE RESULTS FOR CURRENT POPULATION
501 continue
601 continue

! TEST IF MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDED
if (icall.ge.maxn) go to 9000

! COMPUTE THE COUNT ON SUCCESSIVE LOOPS W/O FUNCTION IMPROVEMENT
criter(20)=bestf
if (nloop.lt.(kstop+1)) go to 132
denomi=abs(criter(20-kstop)+criter(20))/2
timeou=abs(criter(20-kstop)-criter(20))/denomi
if (timeou.lt.pcento) go to 9100
132 continue
do l=1,19
criter(l)=criter(l+1)
end do

! IF POPULATION IS CONVERGED INTO A SUFFICIENTLY SMALL SPACE
if (ipcnvg.eq.1) go to 9200

! NONE OF THE STOPPING CRITERIA IS SATISFIED, CONTINUE SEARCH

! CHECK FOR COMPLEX NUMBER REDUCTION
if (ngs1.gt.mings) then
ngs2=ngs1
ngs1=ngs1-1
npt1=ngs1*npg
call comp(nopt,npt1,ngs1,ngs2,npg,x,xf,cx,cf)
end if

! END OF MAIN LOOP -----------
go to 1000

! SEARCH TERMINATED

9000 continue
WRITE (*,*) 'terminated: maximum number of function evaluations
has been exceeded'
WRITE (15,'(a)') 'terminated: maximum number of function
evaluations has been exceeded'

go to 9999
9100 continue
WRITE (*,*) 'terminated: percent change in the last ',kstop,'
loops less than ',pcento*100
WRITE (15,*) 'terminated: percent change in the last ',kstop,'
loops less than ',pcento*100

go to 9999
9200 continue
WRITE (*,*) 'terminated: population has converged to a
sufficiently small space'
WRITE (15,*) 'terminated: population has converged to a
sufficiently small space'

9999 continue

! PRINT THE FINAL PARAMETER ESTIMATE AND ITS FUNCTION VALUE

801 continue

! END OF SUBROUTINE SCEUA

100 format (f2.0,1x,f2.0,1x,f4.1)
101 format (7(1x,f15.12))
109 format (a60,/,a60,/,i5)
113 format (a6)
114 format (2x,i2,2x,f4.0,9(2x,f5.3),2(2x,f7.3))
400 format(//,2x,50(1h=),/,2x,'ENTER THE SHUFFLED COMPLEX EVOLUTION
GLOBAL SEARCH',/,2x,50(1h=))
500 format(//,'*** PRINT THE INITIAL POINT AND ITS CRITERION VALUE
***')
510 format(/,' CRITERION',12(5x,a5),/1x,60(1h-))
520 format(g10.3,12f10.3)
530 format(10x,12(5x,a5))
540 format(10x,12f10.3)
600 format(//,1x,'*** PRINT THE RESULTS OF THE SCE SEARCH ***')
610 format(/,1x,'LOOP',1x,'TRIALS',1x,'COMPLXS',2x,'BEST F',
3x,'WORST F',3x,'PAR RNG',1x,8(5x,a5))
620 format(49x,8(5x,a5))
630 format(i5,1x,i5,3x,i5,3g10.3,8(f10.3))
640 format(49x,8(f10.3))
650 format(/,1x,'POPULATION AT LOOP ',i3,/,1x,22(1h-))
660 format(15x,g10.3,20x,8(f10.3))
800 format(//,1x,'*** OPTIMIZATION SEARCH TERMINATED BECAUSE THE
LIMIT ON THE MAXIMUM',/,5x,'NUMBER OF TRIALS ',i5,' EXCEEDED.' &
' SEARCH WAS STOPPED AT',/,5x,'SUB-COMPLEX ',i3,' OF COMPLEX ',i3,'
IN SHUFFLING LOOP ',i3,' ***')
810 format(//,1x,'*** OPTIMIZATION TERMINATED BECAUSE THE CRITERION
VALUE HAS NOT CHANGED ',/,5x,f5.2,' PERCENT IN',i3, &
' SHUFFLING LOOPS ***')
820 format(//,1x,'*** OPTIMIZATION TERMINATED BECAUSE THE POPULATION
HAS CONVERGED INTO ',/,4x,f5.2, &
' PERCENT OF THE FEASIBLE SPACE ***')
830 format(//,'*** PRINT THE FINAL PARAMETER ESTIMATE AND ITS
CRITERION VALUE ***')
! DEALLOCATE
(Att,ids,modpar,obj,x0,bl,bu,bound,x,cx,s,xf,cf,sf,xx,bestx,worstx,lcs,unit1,xnstd,Qdays_all,Pdays_all)
! End SCE-UA
END PROGRAM
!-------------------------------------------------------------------------

!
-------------------------------------------------------------------------
SUBROUTINE
cce(nopt,nps,s,sf,bl,bu,xnstd,icall,maxn,iseed1,Att,ngauge,ids,num_mod_pars,modpar,obj,Qdays_all,Pdays_all)
! ALGORITHM GENERATES A NEW POINT(S) FROM A SUB-COMPLEX

! SUB-COMPLEX VARIABLES
IMPLICIT NONE
INTEGER :: iseed1
INTEGER, INTENT(IN) :: maxn,nps,nopt
REAL :: icall
REAL, INTENT(IN) :: bl(nopt),bu(nopt)
REAL :: s(nps,nopt),sf(nps)
REAL :: xnstd(nopt),npar
REAL :: sw(nopt),sb(nopt),ce(nopt),snew(nopt),alpha,beta,fw
INTEGER :: i,j,ibound,m,n
REAL :: fnew

! to pass on to FUNCTION similarity:
INTEGER, INTENT(IN) :: num_mod_pars
INTEGER :: Qdays_all(ngauge),Pdays_all(ngauge)

INTEGER, INTENT(IN) :: ngauge
CHARACTER*6 :: ids(ngauge)
REAL :: Att(ngauge,7)
REAL :: modpar(ngauge,num_mod_pars)
REAL :: obj(ngauge)
REAL :: similarity
INTEGER :: ncont

! LIST OF LOCAL VARIABLES
! sb(.) = the best point of the simplex
! sw(.) = the worst point of the simplex
! w2(.) = the second worst point of the simplex
! fw = function value of the worst point
! ce(.) = the centroid of the simplex excluding wo
! snew(.) = new point generated from the simplex
! iviol = flag indicating if constraints are violated
! = 1 , yes
! = 0 , no
!


! EQUIVALENCE OF VARIABLES FOR READABILTY OF CODE
n = nps
m = nopt
alpha = 1.0
beta = 0.5

! IDENTIFY THE WORST POINT wo OF THE SUB-COMPLEX s
! COMPUTE THE CENTROID ce OF THE REMAINING POINTS
! COMPUTE step, THE VECTOR BETWEEN wo AND ce
! IDENTIFY THE WORST FUNCTION VALUE fw
do j=1,m
sb(j)=s(1,j)
sw(j)=s(n,j)
ce(j)=0.0
do i=1,n-1
ce(j)=ce(j)+s(i,j)
end do
ce(j)=ce(j)/dble(n-1)
end do
fw=sf(n)
!
! COMPUTE THE NEW POINT snew
!
! FIRST TRY A REFLECTION STEP
do j=1,m
snew(j)=ce(j)+alpha*(ce(j)-sw(j))
end do

! CHECK IF snew SATISFIES ALL CONSTRAINTS
call chkcst(nopt,snew,bl,bu,ibound)


! snew IS OUTSIDE THE BOUND,
! CHOOSE A POINT AT RANDOM WITHIN FEASIBLE REGION ACCORDING TO
! A NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-COMPLEX
! AS MEAN AND STANDARD DEVIATION OF THE POPULATION AS STD
if (ibound.ge.1) call getpnt(nopt,2,iseed1,snew,bl,bu,xnstd,sb)


! COMPUTE THE FUNCTION VALUE AT snew
ncont=INT(snew(8))

fnew=similarity(nopt,snew,Att,ngauge,ids,num_mod_pars,modpar,obj,Qdays_all,Pdays_all,ncont)

icall=icall + 1

! COMPARE fnew WITH THE WORST FUNCTION VALUE fw

! fnew IS LESS THAN fw, ACCEPT THE NEW POINT snew AND RETURN
if (fnew.le.fw) go to 2000
if (icall.ge.maxn) go to 3000


! fnew IS GREATER THAN fw, SO TRY A CONTRACTION STEP
do j=1,m
snew(j)=ce(j)-beta*(ce(j)-sw(j))
end do

! COMPUTE THE FUNCTION VALUE OF THE CONTRACTED POINT
ncont=INT(snew(8))

fnew=similarity(nopt,snew,Att,ngauge,ids,num_mod_pars,modpar,obj,Qdays_all,Pdays_all,ncont)
icall=icall+1

! COMPARE fnew TO THE WORST VALUE fw
! IF fnew IS LESS THAN OR EQUAL TO fw, THEN ACCEPT THE POINT AND
RETURN
if (fnew.le.fw) go to 2000
if (icall.ge.maxn) go to 3000

! IF BOTH REFLECTION AND CONTRACTION FAIL, CHOOSE ANOTHER POINT
! ACCORDING TO A NORMAL DISTRIBUTION WITH BEST POINT OF THE SUB-
COMPLEX
! AS MEAN AND STANDARD DEVIATION OF THE POPULATION AS STD
1000 call getpnt(nopt,2,iseed1,snew,bl,bu,xnstd,sb)
!
! COMPUTE THE FUNCTION VALUE AT THE RANDOM POINT
ncont=INT(snew(8))

fnew=similarity(nopt,snew,Att,ngauge,ids,num_mod_pars,modpar,obj,Qdays_all,Pdays_all,ncont)
icall=icall+1


! REPLACE THE WORST POINT BY THE NEW POINT
2000 continue
do j=1,m
s(n,j)=snew(j)
end do
sf(n)=fnew
3000 continue

! END OF SUBROUTINE CCE
END SUBROUTINE
!--------------------------------------------------------------------------

!--------------------------------------------------------------------------
REAL FUNCTION
similarity(nopt,x,Att,ngauge,ids,num_mod_pars,modpar,obj,Qdays_all,Pdays_all,ncont)

IMPLICIT NONE
CHARACTER*70 :: filrun,filapet
INTEGER, INTENT(IN) :: ncont
INTEGER, INTENT(IN) :: ngauge
CHARACTER*6 :: ids(ngauge)
INTEGER, INTENT(IN) :: nopt
REAL, INTENT(IN) :: x(nopt)
REAL :: obj(ngauge)
REAL :: S(ngauge),S1(ngauge)
INTEGER :: SP(ngauge),SP1(ngauge)
REAL, INTENT(IN) :: Att(ngauge,7)
INTEGER :: ierror,allocstat ! IOSTAT

REAL :: WSim(ncont),WObj(ncont),Weight(ncont)
REAL :: sum_WSim,max_WSim,sum_WObj,sum_Weight
INTEGER, INTENT(IN) :: num_mod_pars
REAL :: pars(ncont,num_mod_pars)
REAL :: modpar(ngauge,num_mod_pars)
REAL :: Qobs
REAL :: apet(12,2)
INTEGER :: Qdays_all(ngauge),Pdays_all(ngauge)
REAL :: NSE_all(ngauge),NSE,NSE_all1(ngauge)
INTEGER :: i,j,k ! loop variables

INTEGER :: nlines,Qdays,Pdays
INTEGER :: arrsize

do i=1,ngauge
WRITE (filrun,'(a)') 'D:\SENSITIVITY\daily_flow\fortran\' //
ids(i) // 'flow.dat'
WRITE (filapet,'(a)') 'D:\SENSITIVITY\apet\' // ids(i) // 'apet.dat'

Pdays=Pdays_all(i)
Qdays=Qdays_all(i)

OPEN (23,FILE=filapet,STATUS='old',IOSTAT=ierror)
if (ierror.ne.0) then
WRITE (*,*) 'filapet, IOSTAT= ',ierror
end if

do k=1,12
READ (23,103,IOSTAT=ierror) apet(k,1),apet(k,2)
end do

103 format (6x,f2.0,7x,f5.2)

do j=1,ngauge
S(j)=(x(1)*ABS(Att(i,1)-Att(j,1))**x(9)+ &
x(2)*ABS(Att(i,2)-Att(j,2))**x(9)+ &
x(3)*ABS(Att(i,3)-Att(j,3))**x(9)+ &
x(4)*ABS(Att(i,4)-Att(j,4))**x(9)+ &
x(5)*ABS(Att(i,5)-Att(j,5))**x(9)+ &
x(6)*ABS(Att(i,6)-Att(j,6))**x(9)+ &
x(7)*ABS(Att(i,7)-Att(j,7))**x(9))**(1/x(9))

SP(j)=j
end do

! Sort similarity using heapsort
call hpsort(ngauge,S,SP,S1,SP1)

sum_WSim=0
max_WSim=0
sum_WObj=0
sum_Weight=0
do j=1,ncont
WSim(j)=S1(1+j)
if (WSim(j).gt.max_WSim) then
max_WSim=WSim(j)
end if
WObj(j)=obj(SP1(j+1))**x(10)
sum_WObj=sum_WObj+WObj(j)
end do
do j=1,ncont
WSim(j)=1.1-WSim(j)/max_WSim
sum_WSim=sum_WSim+WSim(j)
WObj(j)=WObj(j)/sum_WObj
end do
do j=1,ncont
WSim(j)=WSim(j)/sum_WSim
Weight(j)=WSim(j)*WObj(j)
sum_Weight=sum_Weight+Weight(j)
end do
do j=1,ncont
Weight(j)=Weight(j)/sum_Weight
pars(j,1:num_mod_pars)=modpar(SP1(j+1),1:num_mod_pars)
end do

call
RRMod(NSE,pars,Weight,ncont,Qdays,Pdays,apet,num_mod_pars,ids(i),filrun)

NSE_all(i)=NSE

! end ngauge loop,i
end do


! find median
call hpsort2(ngauge,NSE_all,NSE_all1)
if (MOD(ngauge,2)==0) then
similarity=1-(NSE_all1(ngauge/2)+NSE_all1(ngauge/2+1))/2
else
similarity=1-NSE_all1(ngauge/2)
end if
CLOSE (22)
CLOSE (23)
! End subroutine similarity
END

!--------------------------------------------------------------------

!---------------------------------------------------------------------
SUBROUTINE
RRMod(NSE,pars,Weight,ncont,Qdays,Pdays,apet,num_mod_pars,id,filrun)

IMPLICIT NONE
INTEGER, INTENT(IN) :: ncont,Qdays,Pdays,num_mod_pars
REAL, INTENT(IN) :: pars(ncont,num_mod_pars)
REAL, INTENT(IN) :: Weight(ncont)
REAL, INTENT(IN) :: apet(12,2)
REAL :: P(Pdays,4)
INTEGER :: i,z
REAL :: sms,gw
REAL :: et
REAL :: dt,c0,c1,c2,fip,fop,fdr
REAL :: pa,aet1,rinfcap,runsur,rinfil,runint,rech,gw1,runbas,aet2
REAL :: rtot,fin,fon
INTEGER :: ix
REAL :: Q1
REAL :: Q(Qdays)
CHARACTER*70 :: filrain
CHARACTER*70, INTENT(IN) :: filrun
CHARACTER*6 :: id
INTEGER :: ierror, arrsize,allocstat
REAL, INTENT(OUT) :: NSE

WRITE (filrain,'(a)') 'D:\SENSITIVITY\daily_precip\fortran\' //
id // 'precip.dat'
OPEN (21,FILE=filrain,STATUS='old',IOSTAT=ierror)
if (ierror.ne.0) then
WRITE (*,*) 'filrain, IOSTAT= ',ierror
end if

do i=1,ncont
sms=0.3*pars(i,4)
gw=0

dt=1
c0=(-2.0*pars(i,8)*pars(i,9)+dt)/(2.0*pars(i,8)*(1.0-pars(i,9))+dt);
c1=(2.0*pars(i,8)*pars(i,9)+dt)/(2.0*pars(i,8)*(1.0-pars(i,9))+dt);
c2=(2.0*pars(i,8)*(1.0-pars(i,9))-dt)/(2.0*pars(i,8)*(1.0-pars(i,9))
+dt);
fip=0.1;
fop=0.05;

do z=1,Pdays
if (i.eq.1) then
READ (21,121,IOSTAT=ierror) P(z,1:4)
end if

121 format (4(2x,f14.12))
et=apet(INT(P(z,2)),2)
if (P(z,4)>pars(i,1)) then
pa=P(z,4)
else
pa=pars(i,1)
end if

if (P(z,4).gt.pars(i,1)) then
aet1=pars(i,1)
else
aet1=P(z,4)
end if
if (aet1.gt.et) then
aet1=et
end if

rinfcap=pars(i,2)*EXP(-pars(i,3)*sms/pars(i,4))
if (pa-rinfcap.gt.0) then
runsur=pa-rinfcap
rinfil=rinfcap
else
runsur=0
rinfil=pa
end if
runint=pars(i,5)*sms/pars(i,4)*rinfil
rinfil=rinfil-runint
rech=pars(i,6)*sms/pars(i,4)*rinfil
rinfil=rinfil-rech
sms=sms+rinfil

gw=gw+rech
if (sms-pars(i,4).gt.0) then
gw1=sms-pars(i,4)
else
gw1=0
end if
sms=sms-gw1
gw=gw+gw1
runbas=pars(i,7)*gw
gw=gw-runbas

aet2=sms/pars(i,4)*10
if (aet2.gt.et-aet1) then
aet2=et-aet1
end if
sms=sms-aet2
rtot=runsur+runint+runbas

fin=rtot/24;
fdr=0.0;
do ix=1,24
fon=c0*fin+c1*fip+c2*fop;
fdr=fdr+fon;
fop=fon;
fip=fin;
end DO

if (z.gt.Pdays-Qdays) then

Q1=Weight(i)*fdr
if (i.eq.1) then
Q(z-(Pdays-Qdays))=Q1
else
Q(z-(Pdays-Qdays))=Q(z-(Pdays-Qdays))+Q1
END if
end if
arrsize=SIZE(Q)
! end daily loop,z
end do
! end ncont loop, i
end do

call objective(NSE,Q,Qdays,filrun)

! end function RRMod
CLOSE (21)

END
!--------------------------------------------------------------------

!--------------------------------------------------------------------
SUBROUTINE filesize (ndays,filrun)

IMPLICIT NONE

INTEGER, INTENT(OUT) :: ndays
INTEGER :: ierror ! i/o status
CHARACTER*70, INTENT(IN) :: filrun
REAL :: value
OPEN (3,FILE=filrun,STATUS='old',IOSTAT=ierror)
if (ierror.ne.0) then
WRITE (*,*) 'IOSTATUS= ',ierror
end if
ndays=0
IF (ierror==0) THEN
DO
READ (3,*,IOSTAT=ierror) value
IF (ierror/=0) THEN
EXIT
END IF
ndays=ndays+1
END DO
IF (ierror>0) THEN
WRITE (*,*) 'An error occured when reading line ',ndays+1
ELSE
END IF
ELSE
WRITE (*,*) 'Error opening
file'
END IF

! Close file
CLOSE (3)
END SUBROUTINE
!------------------------------------------------------------------

!------------------------------------------------------------------
! Uses the heapsort algorithm, from Numerical Recipes
SUBROUTINE hpsort(n,ra,rap,outi,outip)
IMPLICIT NONE
INTEGER :: n ! size of array ra
REAL :: ra(n) ! array to be sorted
REAL, INTENT(OUT) :: outi(n) ! sorted rearrangement
INTEGER :: rap(n) ! pointing array
INTEGER, INTENT(OUT) :: outip(n) ! pointing array rearrangement
INTEGER :: i,ir,j,l
REAL :: rra,rra1
INTEGER :: arrsize
if (n.lt.2) then
return
end if
l=n/2+1
ir=n
10 continue
if (l.gt.1) then
l=l-1
rra=ra(l)
rra1=rap(l)
else
rra=ra(ir)
rra1=rap(ir)
ra(ir)=ra(l)
outi(ir)=ra(l)
rap(ir)=rap(l)
outip(ir)=rap(l)
ir=ir-1
if (ir.eq.1) then
ra(1)=rra
outi(1)=rra
rap(1)=rra1
outip(1)=rra1
return
end if
end if
i=l
j=l+l
20 if (j.le.ir) then
if (j.lt.ir) then
if (ra(j).lt.ra(j+1)) then
j=j+1
end if
end if
if (rra.lt.ra(j)) then
ra(i)=ra(j)
outi(i)=ra(j)
rap(i)=rap(j)
outip(i)=rap(j)
i=j
j=j+j
else
j=ir+1
end if
GOTO 20
end if
ra(i)=rra
outi(i)=rra
rap(i)=rra1
outip(i)=rra1
GOTO 10
END
!------------------------------------------------------------------

!------------------------------------------------------------------
! Uses the heapsort algorithm, from Numerical Recipes
SUBROUTINE hpsort1(n,nopt,ra,ra1)
IMPLICIT NONE
INTEGER :: nopt
INTEGER :: n ! size of array ra
REAL :: ra(n),ra1(n,nopt) ! array to be sorted - replaced upon
output by its sorted rearrangement
! ra contains function values, ra1
contains coordinate values. sorted and the replaced on output
INTEGER :: i,ir,j,l,k
REAL :: rra,rra1(nopt)

if (n.lt.2) then
return
end if
l=n/2+1
ir=n
10 continue
if (l.gt.1) then
l=l-1
rra=ra(l)
rra1(1:nopt)=ra1(l,1:nopt)
else
rra=ra(ir)
rra1(1:nopt)=ra1(ir,1:nopt)
ra(ir)=ra(l)
ra1(ir,1:nopt)=ra1(l,1:nopt)
ir=ir-1
if (ir.eq.1) then
ra(1)=rra
ra1(1,1:nopt)=rra1(1:nopt)
return
end if
end if
i=l
j=l+l
20 if (j.le.ir) then
if (j.lt.ir) then
if (ra(j).lt.ra(j+1)) then
j=j+1
end if
end if
if (rra.lt.ra(j)) then
ra(i)=ra(j)
ra1(i,1:nopt)=ra1(j,1:nopt)
i=j
j=j+j
else
j=ir+1
end if
GOTO 20
end if
ra(i)=rra
ra1(i,1:nopt)=rra1(1:nopt)
GOTO 10
END
!----------------------------------------------------------------------

!----------------------------------------------------------------------
! Uses the heapsort algorithm, from Numerical Recipes
SUBROUTINE hpsort2(n,ra,outi)
IMPLICIT NONE
INTEGER :: n ! size of array ra
REAL :: ra(n) ! array to be sorted - replaced upon output by its
sorted rearrangement
REAL, INTENT(OUT) :: outi(n)
INTEGER :: i,ir,j,l
REAL :: rra

if (n.lt.2) then
return
end if
l=n/2+1
ir=n
10 continue
if (l.gt.1) then
l=l-1
rra=ra(l)

else
rra=ra(ir)
ra(ir)=ra(l)
outi(ir)=ra(l)
ir=ir-1
if (ir.eq.1) then
ra(1)=rra
outi(1)=rra
return
end if
end if
i=l
j=l+l
20 if (j.le.ir) then
if (j.lt.ir) then
if (ra(j).lt.ra(j+1)) then
j=j+1
end if
end if
if (rra.lt.ra(j)) then
ra(i)=ra(j)
outi(i)=ra(j)
i=j
j=j+j
else
j=ir+1
end if
GOTO 20
end if
ra(i)=rra
outi(i)=rra
GOTO 10
END
!------------------------------------------------------------------
!------------------------------------------------------------------
! Uses the heapsort algorithm, from Numerical Recipes
SUBROUTINE hpsort3(n,ra,outi)
IMPLICIT NONE
INTEGER :: n ! size of array ra
INTEGER :: ra(n) ! array to be sorted - replaced upon output by its
sorted rearrangement
INTEGER, INTENT(OUT) :: outi(n)
INTEGER :: i,ir,j,l
INTEGER :: rra

if (n.lt.2) then
return
end if
l=n/2+1
ir=n
10 continue
if (l.gt.1) then
l=l-1
rra=ra(l)

else
rra=ra(ir)
ra(ir)=ra(l)
outi(ir)=ra(l)
ir=ir-1
if (ir.eq.1) then
ra(1)=rra
outi(1)=rra
return
end if
end if
i=l
j=l+l
20 if (j.le.ir) then
if (j.lt.ir) then
if (ra(j).lt.ra(j+1)) then
j=j+1
end if
end if
if (rra.lt.ra(j)) then
ra(i)=ra(j)
outi(i)=ra(j)
i=j
j=j+j
else
j=ir+1
end if
GOTO 20
end if
ra(i)=rra
outi(i)=rra
GOTO 10
END
!------------------------------------------------------------------
!------------------------------------------------------------------
SUBROUTINE getpnt(nopt,idist,iseed,x,bl,bu,std,xi)

IMPLICIT NONE

INTEGER :: nopt
REAL :: xi(nopt),x(nopt),bl(nopt),bu(nopt),std(nopt)
INTEGER :: j
REAL :: rand
INTEGER :: iseed
INTEGER :: idist
REAL :: gasdev
INTEGER :: ibound
REAL :: ran1
! This subroutine generates a new point within feasible region

! x(.) = new point
! xi(.) = focal point
! bl(.) = lower bound
! bu(.) = upper bound
! std(.) = standard deviation of probability distribution
! idist = probability flag
! = 1 - uniform distribution
! = 2 - Gaussian distribution


1 do j=1,nopt

2 if (idist.eq.1) then
rand=ran1(iseed)
END if
if (idist.eq.2) then
rand=gasdev(iseed)
END if
x(j)=xi(j)+std(j)*rand*(bu(j)-bl(j))

! Check explicit constraints

call chkcst(1,x(j),bl(j),bu(j),ibound)
if (ibound .ge. 1) go to 2
end do

! Check implicit constraints

call chkcst(nopt,x,bl,bu,ibound)
if (ibound.ge.1) then
goto 1
END if
END SUBROUTINE
!-------------------------------------------------------------

!-------------------------------------------------------------
REAL FUNCTION gasdev(idum)

INTEGER :: idum
! Uses ran1
! Returns a normally distributed deviate with zero mean and unit
variance, using ran1(idum) as the
! source of the uniform deviates

INTEGER :: iset
REAL :: fac,gset,rsq,v1,v2,ran1
SAVE :: iset,gset
DATA iset/0/

if (iset.eq.0) then ! we don't have any extra deviate
handy,
1 v1=2*ran1(idum)-1 ! so pick two uniform numbers in the
square
v2=2*ran1(idum)-1 ! extending from -1 to +1 in each
direction
rsq=v1**2+v2**2 ! See if they are in the unit
circle, if not, try again
if (rsq.ge.1.or.rsq.eq.0) then ! now make the box-muller
transformation to get two
GOTO 1 ! normal deviates - return one
and save the other for next time
end if
fac=SQRT(-2*LOG(rsq)/rsq)
gset=v1*fac
gasdev=v2*fac
iset=1 ! set flag
else ! we have an extra deviate handy
gasdev=gset ! so return it
iset=0 ! and unset the flag
endif

! END OF FUNCTION GASDEV
END
!------------------------------------------------------------

!------------------------------------------------------------
REAL FUNCTION ran1(idum)

INTEGER :: idum,IA,IM,IQ,IR,NTAB,NDIV
INTEGER :: arrsize
REAL :: AM,EPS,RNMX
PARAMETER (IA=16807,IM=2147483647,AM=1./
IM,IQ=127773,IR=2836,NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-
EPS)
! "minimal" random number generator of park and miller with bays-
durham shuffle and added safeguards.
! returns a uniform random deviate between 0.0 and 1.0 (exclusive of
the endpoint values). Call with idum
! a negative integer to initialise; thereafter, do not alter idum
between successive deviates in a sequence. RNMX should approximate
! the largest floating value that is less than 1.
INTEGER :: j,k,iv(NTAB),iy
SAVE iv,iy
DATA iv /NTAB*0/, iy /0/
if (idum.le.0) then ! initialise
idum=MAX(-idum,1) ! prevent idum = 0

do j=NTAB+8,1,-1 ! load the shuffle table after 8 warm-ups
k=idum/IQ

idum=IA*(idum-k*IQ)-IR*k

if (idum.lt.0) then
idum=idum+IM

end if
if (j.le.NTAB) then
iv(j)=idum

end if
end do
iy=iv(1)
END if
k=idum/IQ ! start here when not initialising
idum=IA*(idum-k-IQ)-IR*k ! compute idum=mod(IA*idum,IM) without
overflows by Schrage's method
if (idum.lt.0) then
idum=idum+IM
end if
j=1+iy/NDIV ! Will be in the range 1:NTAB
iy=iv(j) ! Output previously stored value and
refill the shuffle table
iv(j)=idum
ran1=MIN(AM*iy,RNMX) ! Because users don't expect endpoint
values
return
END
!---------------------------------------------------------------------

!---------------------------------------------------------------------
SUBROUTINE parstt(npt,nopt,x,xnstd,bound,gnrng,ipcnvg)

! SUBROUTINE CHECKING FOR PARAMETER CONVERGENCE
IMPLICIT NONE
INTEGER :: nopt,npt
REAL ::
x(npt,nopt),xmax(nopt),xmin(nopt),xmean(nopt),xnstd(nopt),bound(nopt)
REAL :: delta,peps
REAL :: gnrng,ipcnvg,gsum,xsum1,xsum2
INTEGER :: i,k

parameter (delta=1.0d-20,peps=1.0d-3)


! COMPUTE MAXIMUM, MINIMUM AND STANDARD DEVIATION OF PARAMETER VALUES
gsum = 0.d0
do k = 1, nopt
xmax(k) = -1.0d+20
xmin(k) = 1.0d+20
xsum1 = 0.d0
xsum2 = 0.d0
do i = 1, npt
xmax(k) = max(x(i,k), xmax(k))
xmin(k) = min(x(i,k), xmin(k))
xsum1 = xsum1 + x(i,k)
xsum2 = xsum2 + x(i,k)*x(i,k)
end do
xmean(k) = xsum1 / dble(npt)
xnstd(k) = (xsum2 / dble(npt) - xmean(k)*xmean(k))
if (xnstd(k) .le. delta) xnstd(k) = delta
xnstd(k) = sqrt(xnstd(k))
xnstd(k) = xnstd(k) / bound(k)
gsum = gsum + log( delta + (xmax(k)-xmin(k))/bound(k) )
end do
gnrng = dexp(gsum/dble(nopt))

! CHECK IF NORMALIZED STANDARD DEVIATION OF PARAMETER IS <= eps
ipcnvg = 0
if (gnrng .le. peps) then
ipcnvg = 1
end if

! END OF SUBROUTINE PARSTT
return
end
!---------------------------------------------------------------------

!---------------------------------------------------------------------
SUBROUTINE comp(nopt,npt,ngs1,ngs2,npg,a,af,b,bf)

! THIS SUBROUTINE REDUCE INPUT MATRIX a(n,ngs2*npg) TO MATRIX
! b(n,ngs1*npg) AND VECTOR af(ngs2*npg) TO VECTOR bf(ngs1*npg)
IMPLICIT NONE
INTEGER :: npg,npt,nopt
REAL :: a(npt,nopt),af(npt),b(npt,nopt),bf(npt)
INTEGER :: k1,k2,ngs1,ngs2
INTEGER :: i,igs,ipg,j
do igs=1,ngs1
do ipg=1,npg
k1=(ipg-1)*ngs2+igs
k2=(ipg-1)*ngs1+igs
do i=1,nopt
b(k2,i)=a(k1,i)
end do
bf(k2)=af(k1)
end do
end do

do j=1,npt
do i=1,nopt
a(j,i)=b(j,i)
end do
af(j)=bf(j)
end do

! END OF SUBROUTINE COMP
return
end
!---------------------------------------------------------------------

!---------------------------------------------------------------------

SUBROUTINE chkcst(nopt,x,bl,bu,ibound)

! This subroutine check if the trial point satisfies all
! constraints.

! ibound - violation indicator
! = -1 initial value
! = 0 no violation
! = 1 violation
! nopt = number of optimizing variables
! ii = the ii'th variable of the arrays x, bl, and bu
IMPLICIT NONE

! implicit real*8 (a-h,o-z)
! dimension x(nopt),bl(nopt),bu(nopt)
INTEGER :: ii,nopt,ibound
REAL :: x(nopt),bl(nopt),bu(nopt)

ibound = -1

! Check if explicit constraints are violated

do ii=1, nopt
if (x(ii).lt.bl(ii).or.x(ii).gt.bu(ii)) then
go to 10
endif
end do
if (nopt.eq.1) then
go to 9
END if

! Check if implicit constraints are violated
! (no implicit constraints for this function)

! No constraints are violated

9 ibound = 0
go to 20

! At least one of the constraints are violated

10 ibound = 1
20 return
! end subroutine chkcst
end
!------------------------------------------------------------------------

!------------------------------------------------------------------------
SUBROUTINE objective(NSE,Q,Qdays,filrun)

IMPLICIT NONE
CHARACTER*70, INTENT(IN) :: filrun
INTEGER, INTENT(IN) :: Qdays
REAL, INTENT(IN) :: Q(Qdays)
REAL :: NSE1,NSE2,NSE
INTEGER :: ierror,j
REAL :: mean_obs
REAL :: Qobs

OPEN (22,FILE=filrun,STATUS='old',IOSTAT=ierror)
if (ierror.ne.0) then
WRITE (*,*) 'filrun, IOSTAT= ',ierror
end if

NSE1=0
mean_obs=0
do j=1,Qdays
READ (22,102,IOSTAT=ierror) Qobs
102 format (49x,f15.12)
if (Qobs.lt.0) then
Qobs=Q(j)
end if
! WRITE (24,124) Q(j),Qobs
NSE1=NSE1+(Q(j)-Qobs)**2
mean_obs=mean_obs+Qobs
end do
mean_obs=mean_obs/Qdays
NSE2=0
do j=1,Qdays
NSE2=NSE2+(Q(j)-mean_obs)**2
end do
NSE=1-NSE1/NSE2

END SUBROUTINE
.

Relevant Pages

  • Re: Pointer-valued function to access inner components
    ... implicit none ... Return a F90 pointer which points to the diagonal of T. ... a copyin/copyout on array A during the procedure call. ... subroutine dirty_pointer_1d ...
    (comp.lang.fortran)
  • Re: Is this a compiler bug?
    ... implicit none ... integer, allocatable, dimension:: array ... subroutine sub1 ... allocate ) ...
    (comp.lang.fortran)
  • Re: how to use 1D array as a multidimensional array
    ... >> But i guess it could be called array aliasing. ... implicit none ... end program xequiv ... subroutine print_mat ...
    (comp.lang.fortran)
  • Re: rolling dice
    ... implicit none ... subroutine next_in_grid ... find the next vector in an N-dimensional grid (do loop) ... dot_product is an intrinsic function introduced in Fortran 90. ...
    (comp.lang.fortran)
  • Re: Neq Max Limit
    ... you are passing incorrect arguments to a subroutine. ... One may be the array size and another an array. ... dimension y,dy ... PRINT *,i at the top of the outer loop may tell you if the program is ...
    (comp.lang.fortran)