Files
piscal/dataassim/math/othersupmath/targetgridsampling.f
2016-02-03 18:52:05 +00:00

72 lines
2.2 KiB
FortranFixed

subroutine targetgridsampling(sampfunc,nparams0,
& msect0,bestguess,yatguess,bmax0,bmin0,iflag)
implicit none
!
!iflag(i)=0, fix bestguess(i) at the input value
!iflag(i)=1, bestguess(i) is not provided an input value. use grid search
integer nparams0,msect0,iflag(nparams0)
double precision bestguess(nparams0),
& bmax0(nparams0),bmin0(nparams0),
& params(nparams0,msect0+1),yatguess,
& beta(nparams0),fatbeta,bmax(nparams0),bmin(nparams0)
integer i,j,k,n,msect,m,nparams,itag(nparams0)
double precision tiny,x1,delta,eps
parameter(eps=1.0d-8)
external sampfunc
!
nparams=0
do i=1,nparams0
if(iflag(i).eq.1)then
nparams=nparams+1
bmax(nparams)=bmax0(i)
bmin(nparams)=bmin0(i)
itag(nparams)=i
endif
enddo
msect=msect0
do i=1,nparams
tiny=(bmax(i)-bmin(i))*eps
x1=(bmax(i)-bmin(i)-2.0d0*tiny)/dble(msect)
params(i,1)=bmin(i)+tiny
do j=2,msect+1
params(i,j)=params(i,j-1)+x1
params(i,j)=dmax1(params(i,j),bmin(i))
params(i,j)=dmin1(params(i,j),bmax(i))
enddo
enddo
msect=msect+1
yatguess=1.0d+100
do i=1,msect**nparams
do j=1,nparams
!the size of the larger repeated unit is msect**(nparams-j+1)
k=i/(msect**(nparams-j+1))
n=mod(i,(msect**(nparams-j+1)))
if(n.eq.0)k=k-1
!k is the number of repeated units before i (not include the unit i is in)
k=i-k*(msect**(nparams-j+1))
!now k is the position in the larger repeated unit
!the size of the smaller repeated unit is (msect**(nparams-j+1))/msect
m=(msect**(nparams-j+1))/msect
n=k/m
if(mod(k,m).ne.0)n=n+1
beta(j)=params(j,n)
enddo
do j=nparams,1,-1
beta(itag(j))=beta(j)
enddo
do j=1,nparams0
if(iflag(j).eq.0)then
beta(j)=bestguess(j)
endif
enddo
call sampfunc(nparams0,beta,fatbeta)
if(fatbeta.lt.yatguess)then
yatguess=fatbeta
do j=1,nparams0
bestguess(j)=beta(j)
enddo
endif
enddo
return
end