Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • foresee/4C
  • gutsch/4C
2 results
Show changes
Hide whitespace changes
Inline Side-by-side
source_code/version_2.3_windows/timsort.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* SUBROUTINE *!
!* timsort - for sorting of harvested timber to *!
!* different timber qualities *!
!* definition: *!
!* ste - stems *!
!* sg1/sg2 - stem segments *!
!* in1/in2 - industrial wood *!
!* fue - fuelwood *!
!* Subroutine: *!
!* out_tim - generating field sort *!
!* out_timlist *!
!* fuction rabf *!
!* *!
!* Copyright (C) 1996-2018 *!
!* Potsdam Institute for Climate Impact Reserach (PIK) *!
!* Authors and contributors see AUTHOR file *!
!* This file is part of 4C and is licensed under BSD-2-Clause *!
!* See LICENSE file or under: *!
!* http://www.https://opensource.org/licenses/BSD-2-Clause *!
!* Contact: *!
!* https://gitlab.pik-potsdam.de/foresee/4C *!
!* *!
!*****************************************************************!
SUBROUTINE timsort
use data_stand
use data_species
use data_tsort
use data_simul
use data_par
use data_manag
!***************** wpm ************************
use data_wpm
use data_stand
use data_species
!***************** wpm ************************
implicit none
integer i
real h,dbh,db, dcrb, hbo, llazmin, lx,lxz,lldmin ,llasdmin,liszmin,help,lisdmin,llzmin
real h1,h2,dcrb_org,db_org,suml, h_org, sumbio_help,sumvol, h3
real (KIND = dg) :: calcvol
character(4) K
character(2) standt
integer count,count_old,taxid
real, external :: rabf
real :: diam_base=0. ! diameter at basis
llazmin=0.;lx=0.;lldmin=0.;llasdmin=0.; liszmin=0.;lisdmin=0.;llzmin=0.
count = 1
sumbio_help = 0
DO i=1,nspec_tree
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
! Douglasie --- > Kiefer
taxid = zeig%coh%species
IF(taxid .EQ. 10) taxid = 3
IF(taxid .EQ. i) THEN
calcvol = 0.
sumvol = 0.
h = zeig%coh%height -stoh(i) ! stump correction
hbo = zeig%coh%x_hbole
dbh = zeig%coh%diam
dcrb = zeig%coh%dcrb
dcrb_org = dcrb
suml = stoh(i) ! hight of stump; stockhhe
h_org = zeig%coh%height
! calculation of stump biomass for harvesting
! selection of small trees with out dbh
! IF (dbh.eq.0.) THEN
!litter
diam_base= sqrt((zeig%coh%x_ahb+zeig%coh%asapw)*4/pi)
if(hbo.ne.0) then
db = dcrb + (hbo-stoh(i))*(diam_base-dcrb)/hbo
else if (hbo.eq.0)then
db = diam_base*h/(h + stoh(i))
end if
db_org = db
if( hbo.eq.0) then
llzmin=0
lx=0
end if
! stems
help = rabf(i,lzmin(i))
if(db.ge.(lzmin(i) + rabf(i,lzmin(i)))) then
! calculation of lenght at diameter = lzmin
! llzmin > h can occur
if ( dcrb .gt. lzmin(i)+rabf(i,lzmin(i))) then
llzmin = h-(h-hbo)*(lzmin(i)+rabf(i,lzmin(i)))/dcrb
else if (dcrb.le.lzmin(i)+rabf(i,lzmin(i))) then
llzmin = hbo -hbo*(lzmin(i)+rabf(i,lzmin(i))-dcrb)/(db-dcrb)
end if
! calculation of diameter at llzmin/2
if (llzmin/2.lt. hbo) then
lx = dcrb + (db-dcrb)*(hbo-llzmin/2)/hbo
else
lx = dcrb*(h- llzmin/2)/(h-hbo)
end if
! begin of sorting stem lumbre
if( flag_sort.eq.0) then
if( llzmin.ge. lmin(i).and. lx.ge. ldmin(i)+rabf(i,ldmin(i))) then
k = 'ste'
if(zeig%coh%ntreem.ne.0.) then
standt='ab'
call out_tim(count,i,k,llzmin + zug ,lx,lzmin(i), zeig%coh%ntreem,&
standt,h,hbo,db,dcrb,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt='vb'
call out_tim(count,i,k,llzmin + zug ,lx,lzmin(i), zeig%coh%ntreea,&
standt,h,hbo,db,dcrb,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
standt='tb'
flag_deadsort = 1
call out_tim(count,i,k,llzmin + zug ,lx,lzmin(i), zeig%coh%ntreed,&
standt,h,hbo,db,dcrb,calcvol)
end if
end if
suml = suml +llzmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h-llzmin-zug
if(h.lt.0.) h = 0.
db = lzmin(i)
else
! Hhe Hx berechnen, wo lx = ldmin(i)+rabf(i,ldmin(i)) , dann testen, ob 2*Hx >= lmin ist,
!wenn ja, abspeichern in sto
if (dcrb .gt.(ldmin(i)+rabf(i,ldmin(i)))) then
lldmin= h-(h-hbo)*(ldmin(i)+rabf(i,ldmin(i)))/dcrb
else if (dcrb.le.(ldmin(i)+rabf(i,ldmin(i)))) then
lldmin = hbo - hbo*(ldmin(i)+rabf(i,ldmin(i))-dcrb)/(db_org-dcrb)
end if
! calculation of diameter at 2*lldmin and test against lzmin
!(Durchmesser an der Spitze des Stammstckes
if (2*lldmin.lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-lldmin*2)/hbo
else
lx = dcrb*(h- lldmin*2)/(h-hbo)
end if
if (2*lldmin .ge. lmin(i) .and. lx.ge.lzmin(i)+rabf(i,lx)) then
!second test Stammholz!
k = 'ste'
if(zeig%coh%ntreem.ne.0.) then
standt='ab'
call out_tim(count,i,k,2*lldmin + zug,ldmin(i),lx, zeig%coh%ntreem,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt='vb'
call out_tim(count,i,k,2*lldmin + zug,ldmin(i),lx, zeig%coh%ntreea,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
standt='tb'
flag_deadsort = 1
call out_tim(count,i,k,2*lldmin + zug,ldmin(i),lx, zeig%coh%ntreed,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + 2*lldmin + zug
sumvol=sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - 2*lldmin-zug
if(h.lt.0.) h = 0.
db = lx
end if
end if
end if ! db.gt. lzmin(i)
end if ! flag_sortst
! ende test Stammholz
! begin test auf Stammstcke
! calculation of length at diamter lazmin(i)
Do while((h.ge.lasfixl1(i).or.h.ge.lasfixl2(i)).and. hbo.ne.0.)
count_old = count
IF(db .ge.laszmin(i)+rabf(i,laszmin(i)).and. db.gt.lasdmin(i)+ rabf(i,lasdmin(i))) THEN
if (dcrb .eq.0.) then
llazmin = (h_org-suml) -(h_org-stoh(i))*(laszmin(i)+rabf(i,laszmin(i)))/db
else if ( dcrb .gt. laszmin(i)+rabf(i,laszmin(i))) then
llazmin =(h_org-suml)- (h_org-hbo)*(laszmin(i)+rabf(i,laszmin(i)))/dcrb
else if (dcrb.le.laszmin(i)+rabf(i,laszmin(i))) then
llazmin = (hbo-suml) -hbo*(laszmin(i)+rabf(i,laszmin(i))-dcrb)/(db_org-dcrb)
end if
if(llazmin.ge.lasfixl1(i)) then
if(flag_sort.eq.2) then
llazmin = lasfixl2(i)
else
llazmin = lasfixl1(i)
end if
else if(llazmin.ge.lasfixl2(i)) then
llazmin = lasfixl2(i)
end if
! calculation of diameter lx at llazmin/2
if (dcrb .eq.0.) then
lx = db*(h_org-(suml+llazmin/2))/(h_org-stoh(i))
else if ((suml+llazmin/2).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+llazmin/2))/hbo
else
lx = dcrb*(h_org-(suml+ llazmin/2))/(h_org-hbo)
end if
! calculation of diameter at llazmin
if (dcrb .eq.0.) then
lxz = db*(h_org-(suml+llazmin))/(h_org-stoh(i))
else if ((suml+llazmin).lt. hbo) then
lxz = dcrb + (db_org-dcrb)*(hbo-(suml+llazmin))/hbo
else
lxz = dcrb*(h_org-(suml+ llazmin))/(h_org-hbo)
end if
! test
help = lasdmin(i)+rabf(i,lasdmin(i))
! if flag_sort = 2 only lasfixl2 is used
h3 = lasdmin(i)+rabf(i,lasdmin(i))
if (llazmin.ge. lasfixl1(i).and. lx.ge. lasdmin(i)+rabf(i,lasdmin(i)).and. flag_sort.ne.2) then
k = 'sg1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreem,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreea,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreed,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + llazmin+zug
sumvol = sumvol +calcvol
calcvol =0.
count = count + 1
h = h - llazmin-zug
if(h.lt.0.) h = 0.
db = lxz
! test
! if flag_sort = 3 only lasfixl1 is unsed
else if (llazmin.ge.lasfixl2(i).and.llazmin.lt.lasfixl1(i).and. lx.ge. lasdmin(i)+rabf(i,lasdmin(i)).and. flag_sort.ne.3) then
k = 'sg2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreem,standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreea,standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreed,standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + llazmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - llazmin-zug
if(h.lt.0.) h = 0.
db = lxz
else
if (dcrb.eq.0) then
llasdmin = (h_org-suml)- (h_org-stoh(i))*(lasdmin(i)+rabf(i,lasdmin(i)))/db
else if (dcrb .gt. lasdmin(i)+rabf(i,lasdmin(i))) then
llasdmin = (h_org-suml)-(h_org-hbo)*(lasdmin(i)+rabf(i,lasdmin(i)))/dcrb
else if (dcrb.le.lasdmin(i)+rabf(i,lasdmin(i))) then
llasdmin = (hbo-suml)-hbo*(lasdmin(i)+rabf(i,lasdmin(i))-dcrb)/(db_org-dcrb)
end if
if(2*llasdmin.ge.lasfixl1(i)) then
llasdmin = lasfixl1(i)/2.
else if(2*llasdmin.ge.lasfixl2(i)) then
llasdmin = lasfixl2(i)/2.
end if
!calculation lx diameter at 2*llasdmin
if (dcrb .eq.0.) then
lx = db*(h_org-suml-llasdmin*2)/(h_org-stoh(i))
else if ((suml+2*llasdmin).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+2*llasdmin))/hbo
else
lx = dcrb*(h_org- suml-2*llasdmin)/(h_org-hbo)
end if
! if flag_sort = 2 only lasfixl2 is used
if(2*llasdmin.ge.lasfixl1(i).and.lx.ge.laszmin(i)+rabf(i,laszmin(i)).and. flag_sort.ne.2) then
k = 'sg1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
suml = suml + 2*llasdmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
h = h - 2*llasdmin-zug
if(h.lt.0.) h = 0.
db = lx
! if flag_sort = 3 only lasfixl1 is unsed
else if (2*llasdmin.ge.lasfixl2(i).and.2*llasdmin.lt.lasfixl2(i) .and.lx.ge.laszmin(i)+rabf(i,laszmin(i)).and.flag_sort.ne.3) then
k = 'sg2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
suml = suml + 2*llasdmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
h = h - 2*llasdmin-zug
if(h.lt.0.) h = 0.
db = lx
end if
end if
END IF ! db.gt. laszmin(i)
if(count.eq.count_old) exit
END DO
! end test
! assortment LAS1a for pine
Do while((h.ge.lasfixl1(i).or.h.ge.lasfixl2(i)).and.i.eq.3)
count_old = count
IF(db .ge.las1zmin(i)+rabf(i,las1zmin(i)).and. db.gt.las1dmin(i)+ rabf(i,las1dmin(i))) THEN
if (dcrb .eq.0.) then
llazmin = (h_org-suml) -(h_org-stoh(i))*(las1zmin(i)+rabf(i,las1zmin(i)))/db
else if ( dcrb .gt. las1zmin(i)+rabf(i,las1zmin(i))) then
llazmin =(h_org-suml)- (h_org-hbo)*(las1zmin(i)+rabf(i,las1zmin(i)))/dcrb
else if (dcrb.le.las1zmin(i)+rabf(i,las1zmin(i))) then
llazmin = (hbo-suml) -hbo*(las1zmin(i)+rabf(i,las1zmin(i))-dcrb)/(db_org-dcrb)
end if
if(llazmin.ge.lasfixl1(i)) then
if(flag_sort.eq.2) then
llazmin = lasfixl2(i)
else
llazmin = lasfixl1(i)
end if
else if(llazmin.ge.lasfixl2(i)) then
llazmin = lasfixl2(i)
end if
! calculation of diameter lx at llazmin/2
if (dcrb .eq.0.) then
lx = db*(h_org-(suml+llazmin/2))/(h_org-stoh(i))
else if ((suml+llazmin/2).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+llazmin/2))/hbo
else
lx = dcrb*(h_org-(suml+ llazmin/2))/(h_org-hbo)
end if
! calculation of diameter at llazmin
if (dcrb .eq.0.) then
lxz = db*(h_org-(suml+llazmin))/(h_org-stoh(i))
else if ((suml+llazmin).lt. hbo) then
lxz = dcrb + (db_org-dcrb)*(hbo-(suml+llazmin))/hbo
else
lxz = dcrb*(h_org-(suml+ llazmin))/(h_org-hbo)
end if
! if flag_sort = 2 only lasfixl2 is used
help = las1dmin(i)+rabf(i,las1dmin(i))
if (llazmin.ge. lasfixl1(i).and. lx.ge. las1dmin(i)+rabf(i,las1dmin(i)).and.flag_sort.ne.2) then
k = 'sg1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreem,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreea,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,llazmin+zug,lx,lxz,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + llazmin+zug
sumvol = sumvol +calcvol
calcvol =0.
count = count + 1
h = h - llazmin-zug
if(h.lt.0.) h = 0.
db = lxz
! if flag_sort = 3 only lasfixl1 is used
else if (llazmin.ge.lasfixl2(i).and.llazmin.lt.lasfixl1(i).and. lx.ge. las1dmin(i)+rabf(i,las1dmin(i)).and.flag_sort.ne.3) then
k = 'sg2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreem,standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreea,standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,llazmin+zug,lx,lxz,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + llazmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - llazmin-zug
if(h.lt.0.) h = 0.
db = lxz
else
if (dcrb.eq.0) then
llasdmin = (h_org-suml)- (h_org-stoh(i))*(las1dmin(i)+rabf(i,las1dmin(i)))/db
else if (dcrb .gt. las1dmin(i)+rabf(i,las1dmin(i))) then
llasdmin = (h_org-suml)-(h_org-hbo)*(las1dmin(i)+rabf(i,las1dmin(i)))/dcrb
else if (dcrb.le.las1dmin(i)+rabf(i,las1dmin(i))) then
llasdmin = (hbo-suml)-hbo*(las1dmin(i)+rabf(i,las1dmin(i))-dcrb)/(db_org-dcrb)
end if
if(2*llasdmin.ge.lasfixl1(i)) then
llasdmin = lasfixl1(i)/2.
else if(2*llasdmin.ge.lasfixl2(i)) then
llasdmin = lasfixl2(i)/2.
end if
!calculation lx diameter at 2*llasdmin
if (dcrb .eq.0.) then
lx = db*(h_org-suml-llasdmin*2)/(h_org-stoh(i))
else if ((suml+2*llasdmin).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+2*llasdmin))/hbo
else
lx = dcrb*(h_org- suml-2*llasdmin)/(h_org-hbo)
end if
! if flag_sort = 2 only lasfixl2 is used
if(2*llasdmin.ge.lasfixl1(i).and.lx.ge.las1zmin(i)+rabf(i,las1zmin(i)).and.flag_sort.ne.2) then
k = 'sg1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx, zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx, zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
suml = suml + 2*llasdmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
h = h - 2*llasdmin-zug
if(h.lt.0.) h = 0.
db = lx
! if flag_sort = 3 only lasfixl1 is used
else if (2*llasdmin.ge.lasfixl2(i).and.2*llasdmin.lt.lasfixl2(i) .and.lx.ge.las1zmin(i)+rabf(i,las1zmin(i)).and.flag_sort.ne.3) then
k = 'sg2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx, zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx, zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
suml = suml + 2*llasdmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
h = h - 2*llasdmin-zug
if(h.lt.0.) h = 0.
db = lx
end if
end if
END IF ! db.gt. laszmin(i)
if(count.eq.count_old) exit
END DO
! end test LAS1a for pine
! begin test industrial wood
Do while((h.ge.isfixl1(i).or.h.ge.isfixl2(i)).and.hbo.ne.0)
count_old = count
IF(db.gt.iszmin(i)+rabf(i,iszmin(i)).and.db.gt. isdmin(i)+rabf(i,isdmin(i))) THEN
help = iszmin(i)+rabf(i,iszmin(i))
! calculation of length at diameter iszmin(i)
if (dcrb .eq.0.) then
liszmin = h_org -suml -(h_org-stoh(i))*(iszmin(i)+rabf(i,iszmin(i)))/db
else if ( dcrb .gt. iszmin(i)+rabf(i,iszmin(i))) then
liszmin = (h_org-suml)- (h_org-hbo)*(iszmin(i)+rabf(i,iszmin(i)))/dcrb
else if (dcrb.le.(i)+rabf(i,iszmin(i))) then
liszmin = (hbo-suml) -hbo*(iszmin(i)+rabf(i,iszmin(i))-dcrb)/(db_org-dcrb)
end if
if(liszmin.ge.isfixl1(i)) then
liszmin = isfixl1(i)
else if (liszmin.ge.isfixl2(i)) then
liszmin = isfixl2(i)
end if
! calculation of diameter lx at liszmin/2
if (dcrb .eq.0.) then
lx = db*(h_org-suml-liszmin/2)/(h_org-stoh(i))
else if ((suml+liszmin/2).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+liszmin/2))/hbo
else
lx = dcrb*(h_org-suml- liszmin/2)/(h_org-hbo)
end if
! calculation of diameter at liszmin
if (dcrb .eq.0.) then
lxz = db*(h_org-suml-liszmin)/(h_org-stoh(i))
else if ((suml+liszmin).lt. hbo) then
lxz = dcrb + (db_org-dcrb)*(hbo-(suml+liszmin))/hbo
else
lxz = dcrb*(h_org-(suml+ liszmin))/(h_org-hbo)
end if
! test industrial wood Fix length 1
if (liszmin.ge. isfixl1(i).and. lx.ge. isdmin(i)+rabf(i,isdmin(i))) then
k = 'in1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + liszmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - liszmin-zug
if(h.lt.0.) h = 0.
db = lxz
! test industrial wood fix length 2
else if (liszmin.ge.isfixl2(i).and.liszmin.lt.isfixl1(i).and. lx.ge. isdmin(i)+rabf(i,isdmin(i))) then
k = 'in2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + liszmin + zug
sumvol = sumvol + calcvol
calcvol =0.
count = count + 1
h = h - liszmin-zug
if(h.lt.0.) h = 0.
db = lxz
else
if (dcrb.eq.0) then
h1 = isdmin(i)
h2 = rabf(i,isdmin(i))
llasdmin = h_org - suml-(h_org-stoh(i))*(isdmin(i)+rabf(i,isdmin(i)))/db
else if (dcrb .gt. isdmin(i)+rabf(i,isdmin(i))) then
llasdmin = (h_org-suml)-(h_org-hbo)*(isdmin(i)+rabf(i,isdmin(i)))/dcrb
else if (dcrb.le.isdmin(i)+rabf(i,isdmin(i))) then
llasdmin = hbo-suml -hbo*(isdmin(i)+rabf(i,isdmin(i))-dcrb)/(db_org-dcrb)
end if
if(2*llasdmin.ge.isfixl1(i)) then
llasdmin = isfixl1(i)/2.
else if (2*llasdmin.ge.isfixl2(i)) then
llasdmin = isfixl2(i)/2.
end if
!calculation lx diameter at 2*lisdmin
if (dcrb .eq.0.) then
lx = db*(h_org -suml -llasdmin*2)/(h_org -stoh(i))
else if (2*llasdmin.lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-suml-llasdmin)/hbo
else
lx = dcrb*(h_org- llasdmin)/(h_org-hbo)
end if
! test isfixl1
if(2*lisdmin.ge.isfixl1(i).and.lx.ge.iszmin(i)+rabf(i,iszmin(i))) then
k = 'in1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i), zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i), zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i), zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + 2*lisdmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - 2*lisdmin-zug
if(h.lt.0.) h = 0.
db = lx
! test isfixl2
else if (2*lisdmin.ge.isfixl2(i).and.2*lisdmin.lt.isfixl2(i) .and.lx.ge.iszmin(i)+rabf(i,iszmin(i))) then
k = 'in2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i),zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i),zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i),zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + 2*lisdmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - 2*lisdmin-zug
if(h.lt.0.) h = 0.
db = lx
end if
end if
END IF ! db .ge. iszmin
if(count.eq.count_old) exit
END DO
! ende test industrial wood
! begin fuelwood
if (h.ne.0.and. db .ne.0) then
k = 'fue'
lx=0.
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
if (suml.eq.stoh(i)) then
! calculation of fuel wood in the case of total use of stem for fuel wood
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. ! kg DW/tree ---> m/tree
else
! ! calculation of fuelwood volume from all stem segments and total volume of stem, error because stump is not considered
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. - sumvol ! m/tree
end if
call out_tim(count,i,k,h,db, lx, zeig%coh%ntreem, standt,h_org,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
if(suml.eq.stoh(i)) then
! calculation of fuel wood in the case of total use of stem for fuel wood
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. ! kg DW/tree ---> m/tree
help = zeig%coh%x_sap + zeig%coh%x_hrt
else
! calculation of fuelwood volume from all stem segments and total volume of stem, error because stump is not considered
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. - sumvol ! m/tree
end if
standt = 'vb'
call out_tim(count,i,k,h,db, lx, zeig%coh%ntreea, standt,h_org,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
if(suml.eq.stoh(i)) then
! calculation of fuel wood in the case of total use of stem for fuel wood
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. ! kg DW/tree ---> m/tree
help = zeig%coh%x_sap + zeig%coh%x_hrt
else
! calculation of fuelwood volume from all stem segments and total volume of stem, error because stump is not considered
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. - sumvol ! m/tree
end if
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,h,db, lx, zeig%coh%ntreed, standt,h_org,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
end if
end if
zeig => zeig%next
end do
end do
end subroutine timsort
subroutine out_tim(cou,nr,k, len, d,zapf, anz,standt,h,hbo,db,dcrb,calcvol)
use data_tsort
use data_simul
use data_par
!***************** wpm ************************
use data_wpm
use data_stand
use data_species
use data_manag
!***************** wpm ************************
type(mansort_type) :: mansort_ini
integer nr,cou
real len, d, anz,zapf, volume, v1,v2,r,r1,rc,vhelp
real (KIND = dg) :: calcvol
character(4) k
character(2) standt
type(timber) ::tim_ini
tim_ini%year = time
tim_ini%count= cou
tim_ini%ttype = k
tim_ini%specnr = nr
tim_ini%length = len
tim_ini%dia = d
tim_ini%diaor = d -rabf(nr,d)
if(tim_ini%diaor.lt.0) tim_ini%diaor=0
!calculaiton of volume for stem segment, depending on the charcteristics (cone, 2 cones, or frustum of a cone)
! cone: vol=1./3.(pi*h*r)
! frustum: vol = pi*h(r1+r1*r2+r2)/3
r = db*0.5
r1 = zapf*0.5
rc = dcrb*0.5
if(k.eq. 'ste') then
if((len + 10.).lt.hbo) then
volume =( pi*len*(r*r + r*r1 + r1*r1)/3.)/1000000. ! frustum
else
v1 = pi*(hbo-stoh(nr))*(r*r +r*rc + rc*rc)/3.
v2 = pi*(len-stoh(nr)-hbo)*(rc*rc+ rc*r1 + r1*r1)/3.
volume = (v1+v2)/1000000.
end if
else if (k.eq.'fue'.and.hbo.ne.0)then
if( db.gt.dcrb) then
if(len.lt.(h-hbo)) then
volume = ( pi * len* r*r/3 )/1000000.
else
v1 = pi* (len-h+hbo)*(r*r + r*rc + rc*rc)/3. ! frustum
vhelp = pi*hbo*(r*r + r*rc + rc*rc)/3.
v2 = pi* (h-hbo)* rc*rc/3. ! cone
volume = (v1+v2)/1000000.
end if
else
volume = (pi*len*r*r/3.)/1000000.
end if
else if (k.eq.'fue'.and.hbo.eq.0)then
volume = (pi*len*r*r/3.)/1000000.
else
! stem timber or industrial timber
if(hbo .eq.0.) then
volume = (pi*len*r*r/3.)/1000000.
else
volume = ( pi*len*(r*r +r*r1 + r1*r1)/3.)/1000000.
end if
end if
if( volume.lt.0) then
volume = volume
end if
if(calcvol.eq.0.) then
tim_ini%vol = volume
calcvol = volume
else
tim_ini%vol = volume
end if
tim_ini%zapfd = zapf
tim_ini%zapfdor = zapf - rabf(nr,zapf)
if ( tim_ini%zapfdor.lt.0) tim_ini%zapfdor=0.
tim_ini%tnum = anz
tim_ini%stype = standt
tim_ini%hei_tree = h
tim_ini%hbo_tree = hbo
tim_ini%dcrb = dcrb
IF (.not. associated(st%first)) THEN
ALLOCATE (st%first)
st%first%tim = tim_ini
NULLIFY(st%first%next)
anz_list = 1
ELSE
ALLOCATE(ztim)
ztim%tim = tim_ini
ztim%next => st%first
st%first => ztim
anz_list = anz_list +1
END IF
!***************** wpm ************************
! information needed for wpm
if ( flag_wpm > 0 .and. (tim_ini%stype .eq. 'ab'.or.tim_ini%stype .eq. 'tb')) then
if (flag_manreal.eq.1.and.maninf.ne.'tending'.and.maninf.ne.'brushing') then
mansort_ini%year = tim_ini%year
mansort_ini%count = tim_ini%count
mansort_ini%spec = tim_ini%specnr
mansort_ini%typus = tim_ini%ttype
mansort_ini%diam = tim_ini%dia
mansort_ini%diam_wob = tim_ini%diaor
mansort_ini%volume = (tim_ini%vol/kpatchsize)*10000. ! m/patchsize ---> m3/ha
mansort_ini%dw = (tim_ini%vol/kpatchsize)*10000*spar(tim_ini%specnr)%prhos*1000000.*cpart ! m/patchsize ---> kg C/ha
mansort_ini%number = tim_ini%tnum
if (.not. associated(first_mansort)) then
allocate (first_mansort)
first_mansort%mansort = mansort_ini
nullify(first_mansort%next)
else
! build new mansort object
allocate(act_mansort)
act_mansort%mansort = mansort_ini
! chain into the list
act_mansort%next => first_mansort
! set the first pointer to the new object
first_mansort => act_mansort
end if
end if
end if
! information needed for sea or wpm+sea
if ( (flag_wpm == 2 .or. flag_wpm == 3) .and. tim_ini%stype .eq. 'vb') then
mansort_ini%year = tim_ini%year
mansort_ini%count = tim_ini%count
mansort_ini%spec = tim_ini%specnr
mansort_ini%typus = tim_ini%ttype
mansort_ini%diam = tim_ini%dia
mansort_ini%diam_wob = tim_ini%diaor
mansort_ini%volume = (tim_ini%vol/kpatchsize)*10000. ! m/patchsize ---> m3/ha
mansort_ini%dw = (tim_ini%vol/kpatchsize)*10000*spar(tim_ini%specnr)%prhos*1000000.*cpart ! m/patchsize ---> kg C/ha
mansort_ini%number = tim_ini%tnum
if (.not. associated(first_standsort)) then
allocate (first_standsort)
first_standsort%mansort = mansort_ini
nullify(first_standsort%next)
else
! build new mansort object
allocate(act_standsort)
act_standsort%mansort = mansort_ini
! chain into the list
act_standsort%next => first_standsort
! set the first pointer to the new object
first_standsort => act_standsort
end if
end if
!***************** wpm ************************
end subroutine out_tim
subroutine out_timlist
use data_tsort
use data_simul
integer timunit
timunit = getunit()
open (timunit,file = 'timlist.dat', status='unknown')
write( timunit,*) ' year ','count ',' spec','type ',' length',' diameter', 'diam wo bark', 'top diam. ',' top d. wo bark','Volume(m)',' number '
ztim=>st%first
do
IF (.not.ASSOCIATED(ztim)) exit
write(timunit,'(3I6,1x,A5,1x,F8.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,f10.2)') ztim%tim%year, ztim%tim%count, &
ztim%tim%specnr,ztim%tim%ttype,ztim%tim%length,ztim%tim%dia,ztim%tim%diaor,ztim%tim%zapfd,ztim%tim%zapfdor, ztim%tim%vol,ztim%tim%tnum
ztim=>ztim%next
end do
close(timunit)
end subroutine out_timlist
real function rabf(spec, db)
! calculation of rabz i.A.
use data_tsort
use data_species
integer iz, spec
do iz = 1,nspec_tree
if(iz.eq.spec) then
if(db.lt.rabth(spec,1)) then
rabf = rabz(spec,1)
else if (db.ge.rabth(spec,1).and. db.lt.rabth(spec,2)) then
rabf = rabz(spec,2)
else
rabf = rabz(spec,3)
end if
end if
end do
end function rabf
\ No newline at end of file
source_code/version_2.3_windows/tool.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* Subroutines for standard tasks *!
!* *!
!* contains: *!
!* SOLV_QUADR solving quadratic equation, real*4 *!
!* DSOLV_QUADR solving quadratic equation, real*8 *!
!* NEWT Newton method *!
!* TRICOF Harmonic Analysis *!
!* SORT_INDEX Sorts two arrays *!
!* SORT sort an array by quicksort method *!
!* MOMENT Descriptive statistics of a data set *!
!* *!
!* Copyright (C) 1996-2018 *!
!* Potsdam Institute for Climate Impact Reserach (PIK) *!
!* Authors and contributors see AUTHOR file *!
!* This file is part of 4C and is licensed under BSD-2-Clause *!
!* See LICENSE file or under: *!
!* http://www.https://opensource.org/licenses/BSD-2-Clause *!
!* Contact: *!
!* https://gitlab.pik-potsdam.de/foresee/4C *!
!* *!
!*****************************************************************!
SUBROUTINE solv_quadr (meth, p, q, x1, x2, res1, res2, rnum)
! Solution of quadratic equation in normal form
! x*x + p*x + q = 0
IMPLICIT NONE
! Input
integer meth ! solver method
real p, q ! parameter of the quadratic equation
! Output
real x1, x2 ! solutions; initial value of Newton method
real res1, res2 ! residua
integer rnum ! return code
real discr
! Variables of Newton method
real df ! quotient of quadratic function and its first derivative
real :: precision = 1E-5
integer:: maxloop, iloop
! Variables of solver program ZPORC of ISML Library
!external ZPLRC
discr = (p*p/4.)-q
if (discr .lt. 0.) then
rnum = -1 ! no real solution
return
else
select case (meth)
case (1)
! standard solution
discr = SQRT(discr)
x1 = -p/2. + discr
x2 = -p/2. - discr
rnum = 0
case (2)
! Vieta's formulae (root theorem)
discr = SQRT(discr)
x2 = -p/2. - discr
x1 = q / x2
rnum = 0
case (3)
! Newton method
! initial value x2
maxloop = 100
iloop = 1
df = (x2*x2 + p*x2 + q) / (2.* x2 + p)
do while (abs(df) .gt. precision .and. iloop .le. maxloop)
x2 = x2 - df
df = (x2*x2 + p*x2 + q) / (2.* x2 + p)
iloop = iloop + 1
enddo
if (iloop .lt. maxloop) then
rnum = 0
else
rnum = 1
endif
end select
endif ! discr
res1 = x1*x1 + p*x1 + q
res2 = x2*x2 + p*x2 + q
END SUBROUTINE solv_quadr
!**************************************************************
SUBROUTINE dsolv_quadr (meth, p, q, x1, x2, res1, res2, rnum)
! Solution of quadratic equation in normal form
! x*x + p*x + q = 0
! with double precision
IMPLICIT NONE
! Input
integer meth ! solver method
real (kind (0.0D0)) :: p, q ! parameter of the quadratic equation
! Output
real (kind (0.0D0)) :: x1, x2 ! solutions
real (kind (0.0D0)) :: res1, res2 ! residua
integer rnum ! return code
real (kind (0.0D0)) :: discr
! Variables of Newton method
real (kind (0.0D0)) :: df ! quotient of quadratic function and its first derivative
real (kind (0.0D0)) :: precision = 1E-5
integer:: maxloop, iloop
! Variables for solver program ZPORC of ISML Library
real (kind (0.0D0)) :: coeff(3)
real (kind (0.0D0)) :: zero(2)
discr = (p*p/4.)-q
if (discr .lt. 0.) then
rnum = -1 ! no real solution
return
else
select case (meth)
case (1)
! standard solution
discr = DSQRT(discr)
x1 = -p/2. + discr
x2 = -p/2. - discr
rnum = 0
case (2)
! Vieta's formulae (root theorem)
discr = DSQRT(discr)
x2 = -p/2. - discr
x1 = q / x2
rnum = 0
case (3)
! Newton method
! initial value x2
maxloop = 100
iloop = 1
df = (x2*x2 + p*x2 + q) / (2.* x2 + p)
do while (abs(df) .gt. precision .and. iloop .le. maxloop)
x2 = x2 - df
df = (x2*x2 + p*x2 + q) / (2.* x2 + p)
iloop = iloop + 1
enddo
if (iloop .lt. maxloop) then
rnum = 0
else
rnum = 1
endif
end select
endif ! discr
res1 = x1*x1 + p*x1 + q
res2 = x2*x2 + p*x2 + q
END SUBROUTINE dsolv_quadr
!**************************************************************
SUBROUTINE newt (x, f, df, ddf, prec, maxit, rnum)
! Newton method
implicit none
integer :: maxit ! maximum number of iteration
real :: x ! initial value and result
real :: prec ! precision
real :: dx ! quotient of function and its first derivative
real :: hf, hdf, hddf
integer :: rnum ! options: 0 - change of sign allowed,
! 1 - no change of sign
! return code
integer :: i
real, external :: f, df, ddf ! function and its first derivative
hf = f(x)
hdf = df(x)
hddf = ddf(x)
dx = hddf * hf
if (abs(dx) .lt. abs(hdf*hdf)) then ! Test of convergence
! Iteration
i = 1
if (abs(dx) .gt. 0.) then
dx = hf / hdf
endif
do while (abs(hf) .gt. prec .and. i .le. maxit)
if (dx .gt. x .and. rnum .gt. 0) dx = x/2.
x = x - dx
hdf = df(x)
if (abs(hdf) .gt. 0.) then
hf = f(x)
dx = hf/hdf
endif
i = i + 1
enddo
if (i .lt. maxit) then
rnum = 0
else
rnum = 1 ! not enough iteration steps
endif
else
rnum = -1 ! no convergence
endif
END SUBROUTINE newt
!**************************************************************
SUBROUTINE TRICOF(F,NF,A,NE,B,NO,IOP)
! PURPOSE = TO COMPUTE THE COEFFICIENTS IN A TRIGONOMETRIC EXPANSION
! FOR A FUNCTION GIVEN IN EQUIDISTANT POINTS
! PARAMETERS
! F = AN ARRAY USED FOR STORING THE FUNCTION VALUES F(X)
! NF = THE NUMBER OF FUNCTION VALUES IN THE ARRAY F.NF MUST
! HAVE THE STRUCTURE , NF = 2*N+1 , THE GENERAL CASE
! NF = N+1 , THE EVEN CASE
! NF = N-1 , THE ODD CASE
! A = AN ARRAY USED FOR RETURNING THE COEFFICIENTS OF THE COS-
! INE TERMS
! NE = THE NUMBER OF COEFFICIENTS IN THE ARRAY A , NE = N+1
! B = AN ARRAY USED FOR RETURNING THE COEFFICIENTS OF THE SINE
! TERMS
! NO = THE NUMBER OF COEFFICIENTS IN THE ARRAY B , NO = N-1
! IOP = OPTION NUMBER , IOP = 1 , THE GENERAL CASE
! IOP = 2 , THE EVEN CASE
! IOP = 3 , THE ODD CASE
DIMENSION F(NF) , A(NE) , B(NO)
REAL KSI0 , KSI1 , KSIK
DATA ZERO , FOURTH , HALF , ONE , TWO , PI / 0. , .25 , .5 , 1. , 2. , 3.14159265358979 /
! COMPUTE THE NUMBER N (SEE EXPLANATION OF PARAMETERS)
1000 N=0
IF (IOP.EQ.1) N=(NF-1)/2
IF (IOP.EQ.2) N=NF-1
IF (IOP.EQ.3) N=NF+1
IF (N.EQ.0) STOP
! STOP IF IOP DOES NOT HAVE A CORRECT VALUE
IF (IOP.GT.1) GO TO 1030
IF ((2*N-NF+1).NE.0) STOP
! STOP IF NF DOES NOT HAVE THE CORRECT STRUCTURE IN THE GENERAL CASE
! SPLIT THE FUNCTION F(X) IN AN EVEN AND ODD PART
M=N+1
DO 1020 J=1,N
COF1=HALF*(F(M+J)+F(M-J))
COF2=HALF*(F(M+J)-F(M-J))
F(M+J)=COF2
F(M-J)=COF1
1020 CONTINUE
! REWRITE N IN POWERS OF 2 I.E. N=NBASE*2**NEXP
1030 NBASE=N
NEXP =0
1040 NINT =NBASE/2
IF ((NBASE-2*NINT).NE.0) GO TO 1050
NBASE=NINT
NEXP =NEXP+1
GO TO 1040
! DO SOME INITIAL CALCULATIONS
1050 REALN=NBASE
ARG =HALF*PI/REALN
KSI0 =COS(ARG)
ETA0 =SIN(ARG)
! START CALCULATION OF COEFFICIENTS
IF (IOP.EQ.3) GO TO 1160
! ********** EVEN COEFFICIENT CALCULATION **********
! COMPUTE THE BASIC COEFFICIENTS A(K) , K=1(1)(NBASE+1)
! START CALCULATION OF A(1)
NN =NBASE-1
NPOINT=1
NINCRE=2**NEXP
NLOCAL=NINCRE+1
BASEIN=ONE/REALN
A(1) =HALF*(F(1)+F(N+1))
IF (NN.EQ.0) GO TO 1065
DO 1060 J=1,NN
A(1) =A(1)+F(NLOCAL)
NLOCAL=NLOCAL+NINCRE
1060 CONTINUE
1065 A(1) =TWO*BASEIN*A(1)
! START CALCULATION OF A(K) , K=2(1)(NBASE+1)
KSI1=KSI0
KSIK=KSI1
ETA1=ETA0
ETAK=ETA1
CONST=HALF*F(N+1)
DO 1090 K=1,NBASE
COF1=TWO*(TWO*KSIK**2-ONE)
A2 =ZERO
A1 =A2
A0 =CONST
NLOCAL=N+1-NINCRE
DO 1070 J=1,NBASE
A2=A1
A1=A0
A0=F(NLOCAL)+COF1*A1-A2
NLOCAL=NLOCAL-NINCRE
1070 CONTINUE
1080 A(K+1)=BASEIN*(A0-A2)
COF1 =KSIK
COF2 =ETAK
KSIK =KSI1*COF1-ETA1*COF2
ETAK =ETA1*COF1+KSI1*COF2
1090 CONTINUE
! CALCULATION OF THE BASIC EVEN COEFFICIENTS FINISHED
IF (NEXP.EQ.0) GO TO 1145
! CONTINUE CALCULATION OF EVEN COEFFICIENTS
NUMCOF=NBASE
DO 1140 NSTEP=1,NEXP
NINCRE=2**(NEXP-NSTEP)
NPOINT=NINCRE+1
NINCRE=2*NINCRE
NLOCAL=NPOINT
NUMBER=2*NUMCOF+1
! COMPUTE CONSTANT TERM IN MID-POINT APPROXIMATION I.E. K=1
SUM=ZERO
DO 1100 J=1,NUMCOF
SUM=SUM+F(NLOCAL)
NLOCAL=NLOCAL+NINCRE
1100 CONTINUE
SUM =TWO*BASEIN*SUM
COF1=A(1)
A(1)=HALF*(COF1+SUM)
A(NUMBER)=HALF*(COF1-SUM)
IF (NUMCOF.EQ.1) GO TO 1135
! COMPUTE MID-POINT APPROXIMATION FOR K=2(1)NUMCOF
1105 NN =NUMCOF-1
KSIK=KSI1
ETAK=ETA1
DO 1130 K=1,NN
COF1=TWO*(TWO*KSIK**2-ONE)
A2=ZERO
A1=A2
NLOCAL=N+2-NPOINT
A0=F(NLOCAL)
DO 1110 J=1,NN
A2=A1
A1=A0
NLOCAL=NLOCAL-NINCRE
A0=F(NLOCAL)+COF1*A1-A2
1110 CONTINUE
1120 SUM=TWO*BASEIN*(A0-A1)*KSIK
COF1=A(K+1)
A(K+1)=HALF*(COF1+SUM)
A(NUMBER-K)=HALF*(COF1-SUM)
COF1=KSIK
COF2=ETAK
KSIK=KSI1*COF1-ETA1*COF2
ETAK=ETA1*COF1+KSI1*COF2
1130 CONTINUE
1135 A(NUMCOF+1)=HALF*A(NUMCOF+1)
! CALCULATIONS OF MID-POINT APPROXIMATIONS FINISHED
! DO CHANGES RELATED TO HALVING OF THE INTERVAL
ARG =HALF*ARG
COF1=ETA1
ETA1=SIN(ARG)
KSI1=HALF*COF1/ETA1
BASEIN=HALF*BASEIN
NUMCOF=2*NUMCOF
1140 CONTINUE
1145 IF (NEXP.EQ.0) NUMBER=NBASE+1
A(NUMBER)=HALF*A(NUMBER)
! CALULATION OF EVEN COEFFICIENTS FINISHED
1150 IF (IOP.EQ.2) RETURN
! RETURN TO CALLING PROGRAM IF F(X) WAS AN EVEN FUNCTION
! IF IOP=1 CHANGE SIGN OF EACH SECOND COEFFICIENTS
NINT=(N+1)/2
IF (NINT.EQ.0) GO TO 1166
DO 1164 K=1,NINT
A(2*K)=-A(2*K)
1164 CONTINUE
! ********** ODD COEFFICIENT CALCULATION **********
! COMPUTE THE BASIC COEFFICIENTS B(K) , K=1(1)NBASE
1166 ARG=HALF*PI/REALN
1160 IF (IOP.EQ.1) NMAX=2*N+1
IF (IOP.EQ.3) NMAX=N
NINCRE=2**NEXP
NPOINT=NMAX-NINCRE
NLOCAL=NPOINT
BASEIN=ONE/REALN
B(1)=ZERO
IF (NBASE.EQ.1) GO TO 1200
KSI1=TWO*KSI0**2-ONE
KSIK=KSI1
ETA1=TWO*KSI0*ETA0
ETAK=ETA1
NN =NBASE-1
NNN=NN-1
DO 1190 K=1,NN
COF1=TWO*KSIK
A2 =ZERO
A1 =A2
A0 =F(NPOINT)
NLOCAL=NPOINT-NINCRE
IF (NNN.EQ.0) GO TO 1180
DO 1170 J=1,NNN
A2=A1
A1=A0
A0=F(NLOCAL)+COF1*A1-A2
NLOCAL=NLOCAL-NINCRE
1170 CONTINUE
1180 B(K)=TWO*BASEIN*A0*ETAK
COF1=KSIK
COF2=ETAK
KSIK=KSI1*COF1-ETA1*COF2
ETAK=ETA1*COF1+KSI1*COF2
1190 CONTINUE
! CALCULATION OF THE BASIC ODD COEFFICIENTS FINISHED
1200 IF (NEXP.EQ.0) GO TO 1260
! CONTINUE CALCULATION OF ODD COEFFICIENTS
KSI1=KSI0
ETA1=ETA0
NUMCOF=NBASE
DO 1250 NSTEP=1,NEXP
KSIK=KSI1
ETAK=ETA1
NINCRE=2**(NEXP-NSTEP)
NPOINT=NMAX-NINCRE
NINCRE=2*NINCRE
NUMBER=2*NUMCOF
B(NUMCOF)=ZERO
! COMPUTE MID-POINT APPROXIMATIONS FOR K=1(1)NUMCOF
NN =NUMCOF-1
DO 1240 K=1,NUMCOF
COF1=TWO*(TWO*KSIK**2-ONE)
A2 =ZERO
A1 =A2
NLOCAL=NPOINT
A0 =F(NLOCAL)
IF (NN.EQ.0) GO TO 1220
DO 1210 J=1,NN
A2=A1
A1=A0
NLOCAL=NLOCAL-NINCRE
A0=F(NLOCAL)+COF1*A1-A2
1210 CONTINUE
1220 SUM=TWO*BASEIN*(A0+A1)*ETAK
COF1=B(K)
B(K)=HALF*(COF1+SUM)
IF (K.EQ.NUMCOF) GO TO 1230
B(NUMBER-K)=-HALF*(COF1-SUM)
1230 COF1=KSIK
COF2=ETAK
KSIK=KSI1*COF1-ETA1*COF2
ETAK=ETA1*COF1+KSI1*COF2
1240 CONTINUE
! CALCULATION OF MID-POINT APPROXIMATION FINISHED
! DO CHANGES RELATED TO HALVING OF INTERVAL
ARG =HALF*ARG
COF1=ETA1
ETA1=SIN(ARG)
KSI1=HALF*COF1/ETA1
BASEIN=HALF*BASEIN
NUMCOF=2*NUMCOF
1250 CONTINUE
! CALCULATION OF ODD COEFFICIENTS FINISHED
1260 IF (IOP.EQ.3) RETURN
! IF IOP=1 RECOMPUTE FUNCTION VALUES
DO 1270 J=1,N
COF2=F(M+J)
COF1=F(M-J)
F(M+J)=COF1+COF2
F(M-J)=COF1-COF2
1270 CONTINUE
RETURN
END SUBROUTINE TRICOF
!**************************************************************
SUBROUTINE sort_index(n,arr,brr)
! variation of sort2 for integer array
! sorts array arr(1:n) into an ascending order and
! makes the corresponding rearrangement of the array brr(1:n)
INTEGER n,M,NSTACK
Integer arr(n)
INTEGER brr(n)
PARAMETER (M=7,NSTACK=50)
INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
REAL a,b,temp
jstack=0
l=1
ir=n
1 if(ir-l.lt.M)then
do 12 j=l+1,ir
a=arr(j)
b=brr(j)
do 11 i=j-1,1,-1
if(arr(i).le.a)goto 2
arr(i+1)=arr(i)
brr(i+1)=brr(i)
11 continue
i=0
2 arr(i+1)=a
brr(i+1)=b
12 continue
if(jstack.eq.0)return
ir=istack(jstack)
l=istack(jstack-1)
jstack=jstack-2
else
k=(l+ir)/2
temp=arr(k)
arr(k)=arr(l+1)
arr(l+1)=temp
temp=brr(k)
brr(k)=brr(l+1)
brr(l+1)=temp
if(arr(l+1).gt.arr(ir))then
temp=arr(l+1)
arr(l+1)=arr(ir)
arr(ir)=temp
temp=brr(l+1)
brr(l+1)=brr(ir)
brr(ir)=temp
endif
if(arr(l).gt.arr(ir))then
temp=arr(l)
arr(l)=arr(ir)
arr(ir)=temp
temp=brr(l)
brr(l)=brr(ir)
brr(ir)=temp
endif
if(arr(l+1).gt.arr(l))then
temp=arr(l+1)
arr(l+1)=arr(l)
arr(l)=temp
temp=brr(l+1)
brr(l+1)=brr(l)
brr(l)=temp
endif
i=l+1
j=ir
a=arr(l)
b=brr(l)
3 continue
i=i+1
if(arr(i).lt.a)goto 3
4 continue
j=j-1
if(arr(j).gt.a)goto 4
if(j.lt.i)goto 5
temp=arr(i)
arr(i)=arr(j)
arr(j)=temp
temp=brr(i)
brr(i)=brr(j)
brr(j)=temp
goto 3
5 arr(l)=arr(j)
arr(j)=a
brr(l)=brr(j)
brr(j)=b
jstack=jstack+2
if(jstack.gt.NSTACK)pause 'NSTACK too small in sort2'
if(ir-i+1.ge.j-l)then
istack(jstack)=ir
istack(jstack-1)=i
ir=j-1
else
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif
endif
goto 1
END
! (C) Copr. 1986-92 Numerical Recipes Software "!D#+.
!**************************************************************
SUBROUTINE sort(n,arr)
! sort a n-dimensional array arr(1:n) by quicksort method
INTEGER n,M,NSTACK
REAL arr(n)
PARAMETER (M=7,NSTACK=50)
INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
REAL a,temp
jstack=0
l=1
ir=n
1 if(ir-l.lt.M)then
do 12 j=l+1,ir
a=arr(j)
do 11 i=j-1,1,-1
if(arr(i).le.a)goto 2
arr(i+1)=arr(i)
11 continue
i=0
2 arr(i+1)=a
12 continue
if(jstack.eq.0)return
ir=istack(jstack)
l=istack(jstack-1)
jstack=jstack-2
else
k=(l+ir)/2
temp=arr(k)
arr(k)=arr(l+1)
arr(l+1)=temp
if(arr(l+1).gt.arr(ir))then
temp=arr(l+1)
arr(l+1)=arr(ir)
arr(ir)=temp
endif
if(arr(l).gt.arr(ir))then
temp=arr(l)
arr(l)=arr(ir)
arr(ir)=temp
endif
if(arr(l+1).gt.arr(l))then
temp=arr(l+1)
arr(l+1)=arr(l)
arr(l)=temp
endif
i=l+1
j=ir
a=arr(l)
3 continue
i=i+1
if(arr(i).lt.a)goto 3
4 continue
j=j-1
if(arr(j).gt.a)goto 4
if(j.lt.i)goto 5
temp=arr(i)
arr(i)=arr(j)
arr(j)=temp
goto 3
5 arr(l)=arr(j)
arr(j)=a
jstack=jstack+2
if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
if(ir-i+1.ge.j-l)then
istack(jstack)=ir
istack(jstack-1)=i
ir=j-1
else
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif
endif
goto 1
END
! C (C) Copr. 1986-92 Numerical Recipes Software )$!.
!**************************************************************
SUBROUTINE moment(array,n,ave,adev,sdev,var,skew,curt)
! Calculates statistics of array (n-dimensional array of data)
! n - number of observations
! adev - average deviation
! ave - average
! curt - curtosis
! sdev - standard deviation
! skew - skewness
! var - variance
INTEGER n
REAL adev,ave,curt,sdev,skew,var,array(n)
INTEGER j
REAL p,s,ep
if(n.le.1)pause 'n must be at least 2 in moment'
s=0.
do 11 j=1,n
s=s+array(j)
11 continue
ave=s/n
adev=0.
var=0.
skew=0.
curt=0.
ep=0.
do 12 j=1,n
s=array(j)-ave
ep=ep+s
adev=adev+abs(s)
p=s*s
var=var+p
p=p*s
skew=skew+p
p=p*s
curt=curt+p
12 continue
adev=adev/n
var=(var-ep**2/n)/(n-1)
sdev=sqrt(var)
if(var.ne.0.)then
skew=skew/(n*sdev**3)
curt=curt/(n*var**2)-3.
else
skew = -99.
curt = -99.
endif
return
END
! (C) Copr. 1986-92 Numerical Recipes Software )$!.
FUNCTION rtbis(func,x1,x2,xacc)
INTEGER JMAX
REAL rtbis,x1,x2,xacc,func
EXTERNAL func
PARAMETER (JMAX=40)
INTEGER j
REAL dx,f,fmid,xmid
fmid=func(x2)
f=func(x1)
if(f.lt.0.)then
rtbis=x1
dx=x2-x1
else
rtbis=x2
dx=x1-x2
endif
do j=1,JMAX
dx=dx*.5
xmid=rtbis+dx
fmid=func(xmid)
if(fmid.le.0.)rtbis=xmid
if(abs(dx).lt.xacc .or. fmid.eq.0.) return
end do
END function
\ No newline at end of file
source_code/version_2.3_windows/tool1.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for standard tasks *!
!* *!
!* contains: *!
!* DAINTZ Date to day of the year *!
!* TZINDA Day of the year to date *!
!* TAB_INT Table function *!
!* CHARACTER_IN_INTEGER Conversion of character in integer *!
!* INTEGER_IN_CHARACTER Conversion of integer in character *!
!* QUANTILE calculates the 0.05 and 0.95 quantile *!
!* QUANT_CALC calculates a quantile of a sorted array *!
!* *!
!* Copyright (C) 1996-2018 *!
!* Potsdam Institute for Climate Impact Reserach (PIK) *!
!* Authors and contributors see AUTHOR file *!
!* This file is part of 4C and is licensed under BSD-2-Clause *!
!* See LICENSE file or under: *!
!* http://www.https://opensource.org/licenses/BSD-2-Clause *!
!* Contact: *!
!* https://gitlab.pik-potsdam.de/foresee/4C *!
!* *!
!*****************************************************************!
SUBROUTINE DAINTZ(IT,IM,IJ,TZ)
! Umrechnen von Datum in Tageszaehler
implicit none
INTEGER IT, IM, IJ, TZ
INTEGER I, ME
REAL, DIMENSION(12):: MNL
! COMMON /MONTH/ MMM(12),JAHR,JS,IC
DATA MNL /31,28,31,30,31,30,31,31,30,31,30,31/
TZ=IT
IF (IM.EQ.1) RETURN
ME=IM-1
if (mod(IJ,4).EQ.0) MNL(2)=29
if ((IJ .eq. 1900) .or. (IJ .eq. 1800) .or. (IJ .eq. 1700)) MNL(2)=28
DO I=1,ME
TZ=TZ+MNL(I)
enddo
MNL(2)=28
END SUBROUTINE DAINTZ
!***********************************************************************
SUBROUTINE TZINDA(T,M,J,TZ)
! Umrechnen von Tageszaehler in Datum
implicit none
INTEGER MNL(12)
INTEGER T, M, J, TZ
DATA MNL /31,28,31,30,31,30,31,31,30,31,30,31/
if (mod(J,4).EQ.0) MNL(2)=29
if ((J .eq. 1900) .or. (J .eq. 1800) .or. (J .eq. 1700)) MNL(2)=28
T=TZ
M=1
do while (T .gt. MNL(M))
T=T-MNL(M)
M=M+1
if (M .gt. 12) return
enddo
MNL(2)=28
END SUBROUTINE TZINDA
!***********************************************************************
SUBROUTINE tab_int(x,y,idim,arg,val)
! Read a table function with ordered pairs x,y (sortet)
! linear interpolation between
implicit none
! input
integer idim ! dimension of array x, y
real, dimension(idim) :: x, y ! table values
real arg ! argument of function
! output
real val ! result
integer i
if (arg .le. x(1)) then
val = y(1)
else if (arg .ge. x(idim)) then
val = y(idim)
else
i = 2
do while ((i .lt. idim) .and. (arg .gt. x(i)))
i = i+1
enddo
if (arg .eq. x(i)) then
val = y(i)
else
val = y(i) + (y(i)-y(i-1)) * (arg-x(i)) / (x(i)-x(i-1))
endif
endif
END subroutine tab_int
!***********************************************************************
SUBROUTINE character_in_integer(string, vint)
! Conversion of character variable in integer variable
implicit none
integer vint
character (100) string
character (10) help
write(help,'(A)') string
read(help,*) vint
END subroutine character_in_integer
!**************************************************************
SUBROUTINE integer_in_character(vint, string)
! Conversion of integer variable in character variable
implicit none
integer vint
character (10) string
character (10) help
write(help,'(I10)') vint
read(help,*) string
END subroutine integer_in_character
!**************************************************************
SUBROUTINE quantile(idim, arr, quant05, quant95, median)
! sorts and calculates the 0.05 and 0.95 quantile of an array with dimension idim
implicit none
! input
integer idim ! dimension of array arr
real, dimension(idim) :: arr ! array
! output
real quant05, quant95, median ! 0.05 and 0.95 quantile
call sort(idim,arr)
call quant_calc(idim, arr, 0.05, quant05) ! 0.05 quantile
call quant_calc(idim, arr, 0.95, quant95) ! 0.95 quantile
call quant_calc(idim, arr, 0.5, median) ! 0.95 quantile
END SUBROUTINE quantile
!**************************************************************
SUBROUTINE quant_calc(idim, arr, pord, quant)
! calculates a quantile of a sorted array with dimension idim
implicit none
integer idim ! dimension of array arr
real, dimension(idim) :: arr ! array
real quant ! quantile
real pord, help ! order
integer ihelp
help = idim * pord
ihelp = int(help)
if (ihelp*1.0 .lt. help) then
quant = arr(ihelp+1)
else
quant = (arr(ihelp+1) + arr(ihelp)) / 2.
endif
END SUBROUTINE quant_calc
!**************************************************************