c Threshold recognition program version 4, Oct 12, 2004 c Developer: Lianhong Gu, Environmental Sciences Division c Oak Ridge National Laboratory c lianhong-gu@ornl.gov c PROGRAM main implicit double precision (a-h,l,o-z) dimension fmet(36),startday(4),endday(4), & dayfrac(8000),timefrac(8000),sinele_ori(8000), & soil_temp(8000),air_temp(8000),ustmedian0(8000), & co2flux(8000),fricv(8000),fneenorm(8000),par(8000) dimension ustarbinmedian(1000), & ratiobinmean(1000),confidbinradius(1000) dimension iflag(8000),ialltimeflag(8000) character*16 fparam(6),fpred(10),seasons(4),fbin(6) character*7 rawfile(13),procfile(13) character filepath1*35,filepath2*37,paramfile*12, & binfile*12 parameter(numpoints=25) parameter(ivarassump=1) c ivarassump = 1, unequal variance assumption c ivarassump = 2, equal variance assumption c c ------iflag or ialltimeflag------ c = 0, passed all three tests (outlier test, u* star test, calm night test) c = 1, passed outlier test but identified as affected by pressure pumping c = 2, passed outlier test but identified as affected by stable condition c = 3, outlier (>3*sd) c = 4, outlier (<-3sd) c = 5, passed outlier test and u* test but failed calm night test c = 6, ustar can not be used as a criterion c = 9999, daytime data c cutust_c = 9999.0d0 normalized flux alway increases with u* c u* correction method cannot be used. c cutust_c=0. No u* dependence. No cutting of data is necessary c c cutust_p= -9999.0d0 normalized flux alway increases with u* c u* correction method cannot be used. c cutust_p= 9999.0d0 No pressure pumping effect observed c fparam(1)='Season' fparam(2)='custar_c' fparam(3)='custar_p' fparam(4)='rsq_temp_co2' fparam(5)='%remaining' fparam(6)='year' c fpred(1)='day' fpred(2)='time' fpred(3)='orinee' fpred(4)='normnee' fpred(5)='soiltemp' fpred(6)='airtemp' fpred(7)='ustar' fpred(8)='iflag' fpred(9)='PAR' fpred(10)='sinele' c startday(1)=334.5d0 endday(1)=59.5d0 seasons(1)='DJF' startday(2)=59.5d0 endday(2)=151.5d0 seasons(2)='MAM' startday(3)=151.5d0 endday(3)=243.5d0 seasons(3)='JJA' startday(4)=243.5d0 endday(4)=334.5d0 seasons(4)='SON' fbin(1)='Year' fbin(2)='StartDay' fbin(3)='EndDay' fbin(4)='ustar' fbin(5)='neenorm' fbin(6)='confid' latit=42.535720d0 longit=-72.171997d0 filepath1= & '/usr/users/l2g/harvarddata/rawdata/' filepath2= & '/usr/users/l2g/harvarddata/processed/' rawfile(1)='hf91gen' rawfile(2)='hf92gen' rawfile(3)='hf93gen' rawfile(4)='hf94gen' rawfile(5)='hf95gen' rawfile(6)='hf96gen' rawfile(7)='hf97gen' rawfile(8)='hf98gen' rawfile(9)='hf99gen' rawfile(10)='hf00gen' rawfile(11)='hf01gen' rawfile(12)='hf02gen' rawfile(13)='hf03gen' procfile(1)='hf91pro' procfile(2)='hf92pro' procfile(3)='hf93pro' procfile(4)='hf94pro' procfile(5)='hf95pro' procfile(6)='hf96pro' procfile(7)='hf97pro' procfile(8)='hf98pro' procfile(9)='hf99pro' procfile(10)='hf00pro' procfile(11)='hf01pro' procfile(12)='hf02pro' procfile(13)='hf03pro' paramfile='hfustarparam' binfile= 'hfustarbined' OPEN(UNIT=3,FILE=filepath2//paramfile) write(3,255)fparam open(unit=4,file=filepath2//binfile) write(4,265)fbin c do 10 iyear=2,13 write(*,*)'The year is ',iyear open(unit=2,file=filepath2//procfile(iyear)) write(2,245)fpred do 40 k=1,4 Nsamp=0 OPEN(UNIT=5,FILE=filepath1//rawfile(iyear)) read(5,*) read(5,*) read(5,*) read(5,*) 20 read(5,*,end=100)fmet locday=dint(fmet(1)) loctime=(fmet(1)-locday)*240.0d0 loctime=dnint(loctime)/10.0d0 locddy=locday+loctime/24.0d0 gmttime=loctime+5.0d0 gmtday=locday if(gmttime.ge.24.0d0)then gmttime=gmttime-24.0d0 gmtday=gmtday+1.0d0 endif year=dble(iyear+1990) call sun(year,gmtday,gmttime,latit,longit, &sinele,sunrise,sunset) if(k.ne.1)then if(locddy.le. & startday(k).or.locddy.gt.endday(k))goto 20 endif if(k.eq.1)then if(locddy.ge. & endday(k).and.locddy.lt.startday(k))goto 20 endif ! Nsamp=Nsamp+1 sinele_ori(Nsamp)=sinele dayfrac(Nsamp)=locddy timefrac(Nsamp)=loctime co2flux(Nsamp)=fmet(11) if(co2flux(Nsamp).lt.-100.0d0.or. & co2flux(Nsamp).gt.200.0d0)then co2flux(Nsamp)=-9999.0d0 else strco2=fmet(34) if(strco2.gt.-100.0d0.and.strco2.lt.100.0d0)then co2flux(Nsamp)=co2flux(Nsamp)+strco2 endif endif par(Nsamp)=fmet(19) air_temp(Nsamp)=fmet(25) if(air_temp(Nsamp).lt.-100.0d0.or. & air_temp(Nsamp).gt.100.0d0)then air_temp(Nsamp)=fmet(26) endif soil_temp(Nsamp)=fmet(30) if(soil_temp(Nsamp).lt.-100.0d0.or. & soil_temp(Nsamp).gt.100.0d0)then soil_temp(Nsamp)=fmet(32) endif ustar=dabs(fmet(12)) if(dabs(ustar-9999.0d0).lt.0.1d-5)then ustar=-9999.0d0 else ustar=dsqrt(ustar)/100.0d0 endif fricv(Nsamp)=ustar goto 20 100 close(5) if(Nsamp.lt.5)goto 40 call ustarfilter(dayfrac,timefrac,sinele_ori,par, & air_temp,soil_temp,co2flux,fricv,Nsamp,numpoints, & ivarassump,fneenorm, & iflag,ialltimeflag,ustmedian0,cutust_c,cutust_p, & r_ust_normco2,premain,ustarbinmedian, & ratiobinmean,confidbinradius,nbins,ierr) if(ierr.ne.0)goto 40 write(3,275)seasons(k),cutust_c, &cutust_p,r_ust_normco2,premain,iyear+1990 write(2,*)seasons(k) do 110 I=1,Nsamp write(2,55)dayfrac(I),timefrac(I),co2flux(I), & fneenorm(I),soil_temp(I),air_temp(I),fricv(I), & ialltimeflag(I),par(I),sinele_ori(I) 110 continue write(2,*) write(4,*)seasons(k) do 1115 I=1,nbins write(4,365)iyear+1990,startday(k),endday(k), & ustarbinmedian(I), & ratiobinmean(I),confidbinradius(I) 1115 continue write(4,*) 40 continue close(2) 10 continue close(3) close(4) 55 FORMAT(1X,7f12.5,i12,2f12.5) 245 format(1x,10a12) 255 format(1x,6a16) 265 format(1x,6a16) 365 format(1x,i5,5f12.5) 275 format(1x,a12,4f12.5,i5) END c c subroutine ustarfilter(dayfrac_ori,time_ori,sinele_ori, & solarradiation,air_temp_ori,soil_temp_ori, & co2flux_ori,fricv_ori,N_ori,numpoints,ivarassump, & fneenorm, & iflag,ialltimeflag,ustmedian0, & cutust_c,cutust_p,r_ust_normco2,premain,ustarbinmedian, & ratiobinmean,confidbinradius,nbins,ierr) implicit double precision (a-h,l,o-z) ! ! --------- Inputs --------------------------------- ! dayfrac_ori(1:N_ori): day of year with fraction ! time_ori(1:N_ori): time of day ! sinele_ori(1:N_ori): sine of solar elevation angle ! solarradiation(1:N_ori): solar radiation ! air_temp_ori(1:N_ori): air temperature (oC) ! soil_temp_ori(1:N_ori): soil temperature (oC) ! co2flux_ori(1:N_ori): net ecosystem exchange of CO2 ! fricv_ori(1:N_ori): friction velocity (ms-1) ! N_ori: number of data points ! numpoints: the number of data points in the moving sample, also the ! number of data points in a bin ! ivarassump: variance assumption indicator ! 1: unequal variance assumption ! 2: equal variance assumption ! dimension dayfrac_ori(N_ori),time_ori(N_ori),sinele_ori(N_ori), & air_temp_ori(N_ori),soil_temp_ori(N_ori), & co2flux_ori(N_ori),fricv_ori(N_ori), & solarradiation(N_ori) ! ! -------- Outputs ----------------------------------- ! fneenorm(1:N): Normalized CO2 flux ! iflag(1:N): datum category indicator for nighttime data only ! = 0, passed all three tests (outlier test, u* star test, calm night test) ! = 1, passed outlier test but identified as affected by pressure pumping ! = 2, passed outlier test but identified as affected by stable condition ! = 3, outlier (>3*sd) ! = 4, outlier (<-3sd) ! = 5, passed outlier test and u* test but failed calm night test ! = 6, ustar can not be used as a criterion ! ialltimeflag(1:N_ori): datum category indicator for all data ! flag indicator same as iflag except including 9999 fordaytime data ! ! ustmedian0(1:N_ori): the median ustar of the night in which the datum belongs, ! for daytime, the value set to 9999 ! cutust_c: Lower ustar threshold ! cutust_p: Upper ustar threshold ! cutust_c = 9999.0d0 normalized flux alway increases with u* ! u* correction method cannot be used. ! cutust_c=0. No u* dependence. No cutting of data is necessary ! cutust_p= -9999.0d0 normalized flux alway increases with u* ! u* correction method cannot be used. ! cutust_p= 9999.0d0 No pressure pumping effect observed ! r_ust_normco2: correlation coefficent of the normalized biological fluxes ! premain: percentage of the biolological fluxes in the oringinal dataset ! ustarbinmedian(1:nbins): median ustar in a bin ! ratiobinmean(1:nbins): mean normalized flux in a bin ! confidbinradius(1:nbins): the total number of bins ! ierr: error message. =0, ok; =1, too few effective data points dimension fneenorm(N_ori),ustmedian0(N_ori), & ustarbinmedian(N_ori), & ratiobinmean(N_ori),confidbinradius(N_ori) dimension iflag(N_ori),ialltimeflag(N_ori) ! ! -------------- Local variables -------------------------------------- dimension tempsoil(N_ori),tair(N_ori),conratio(N_ori), & conustar(N_ori),ratgood(N_ori),ustgood(N_ori), & conustme(N_ori),scrndnee(N_ori),scrndts(N_ori), & scrndta(N_ori),ustH_mem(N_ori),ustL_mem(N_ori), & day(N_ori),time(N_ori),ratio(N_ori), & co2flux(N_ori),fricv(N_ori) parameter(confidlevel=90.0d0) N=0 nsamenight=0 do 2000 i=1,N_ori if(solarradiation(i).lt.10.0d0)goto 2005 ! if(solarradiation(i).lt.20.0d0.and. ! & sinele_ori(i).lt.0.0d0)goto 2005 if(nsamenight.eq.0)goto 2000 ! night ends do 2045 kk=N-nsamenight+1,N ustmedian0(kk)=ustarmet(fricv,N,nsamenight) 2045 continue nsamenight=0 goto 2000 2005 if(nsamenight.eq.0)then starttime=dayfrac_ori(i) endif deltasamenight=dabs(dayfrac_ori(i)-starttime) if(deltasamenight.gt.1.0d0)then do 2035 kk=N-nsamenight+1,N ustmedian0(kk)=ustarmet(fricv,N,nsamenight) 2035 continue nsamenight=0 endif if(soil_temp_ori(i).lt.-100.0d0.or. & soil_temp_ori(i).gt.100.0d0)goto 2000 if(air_temp_ori(i).lt.-100.0d0.or. & air_temp_ori(i).gt.100.0d0)goto 2000 if(fricv_ori(i).lt.0.0d0)goto 2000 if(co2flux_ori(i).lt.0.00001d0.or. & co2flux_ori(i).gt.100.0d0)goto 2000 N=N+1 nsamenight=nsamenight+1 day(N)=dayfrac_ori(i) time(N)=time_ori(i) tempsoil(N)=soil_temp_ori(i)+273.15d0 co2flux(N)=co2flux_ori(i) fricv(N)=fricv_ori(i) tair(N)=air_temp_ori(i)+273.15d0 2000 continue if(nsamenight.ne.0)then do 2055 kk=N-nsamenight+1,N ustmedian0(kk)=ustarmet(fricv,N,nsamenight) 2055 continue nsamenight=0 endif c c u* threshold search starts c ustL_pre=-99999.0d0 ustH_pre=99999.0d0 nrepeat=0 ustH_mem(1)=99999.0d0 ustL_mem(1)=-99999.0d0 ierr=0 if(N.lt.5)then ierr=1 return endif do 45 I=1,N ratio(I)= & co2flux(I)/(2.0d0**((tempsoil(I)-283.15d0)/10.0d0)) 45 continue c 995 call outliers(N,ratio,iflag) icond=1 do 1000 I=1,N if(iflag(I).eq.0)then conratio(icond)=ratio(I) conustar(icond)=fricv(I) conustme(icond)=ustmedian0(I) icond=icond+1 endif 1000 continue icond=icond-1 call cutoffustar(conustar,conratio,conustme,icond, & numpoints,cutust_p,cutust_c,ivarassump) if(cutust_p.lt.0.0d0.or.cutust_c.gt.100.0d0)goto 885 deltap=dabs(cutust_p-ustH_pre) deltac=dabs(cutust_c-ustL_pre) if(deltap.lt.1.0d-5.and.deltac.lt.1.0d-5)goto 885 c------------------------------------ c Detect if a cycle is occurring if(nrepeat.le.1)goto 40 icyc=nrepeat 10 deltac=dabs(cutust_c-ustL_mem(icyc)) deltap=dabs(cutust_p-ustH_mem(icyc)) if(deltap.lt.1.0d-5.and.deltac.lt.1.0d-5)goto 885 if(icyc.eq.1)goto 40 icyc=icyc-1 goto 10 c-------------------------------------- 40 nrepeat=nrepeat+1 ustH_mem(nrepeat)=ustH_pre ustL_mem(nrepeat)=ustL_pre ustH_pre=cutust_p ustL_pre=cutust_c icond=1 do 865 I=1,N if(iflag(I).eq.0.and.ustmedian0(I).ge.cutust_c.and. &fricv(I).ge.cutust_c.and.fricv(I).le.cutust_p)then scrndnee(icond)=co2flux(I) scrndta(icond)=tair(I) scrndts(icond)=tempsoil(I) icond=icond+1 endif 865 continue icond=icond-1 call regression(scrndnee,scrndta,scrndts,icond,param1, & param2,param3,param4,param5) do 845 I=1,N ratio(I)= & co2flux(I)/ecoresp(tair(I),tempsoil(I),param1,param2, & param3,param4,param5) 845 continue goto 995 c u* thresholds search ends c 885 nfinal=1 do 1110 I=1,icond if((conustar(I).ge.cutust_c) & .and.conustar(I).le.cutust_p.and. & conustme(I).ge.cutust_c)then ratgood(nfinal)=conratio(I) ustgood(nfinal)=conustar(I) nfinal=nfinal+1 endif 1110 continue nfinal=nfinal-1 if(nfinal.le.numpoints)then r_ust_normco2=-9999.0d0 else r_ust_normco2=correlation(ratgood,ustgood,nfinal) endif premain=100.0d0*dble(nfinal)/dble(N) do 1100 I=1, N if(cutust_c.gt.100.0d0.or.cutust_p.lt.0.0d0)then iflag(I)=6 else if(iflag(I).eq.0)then if(fricv(I).lt.cutust_c)then iflag(I)=2 else if(fricv(I).gt.cutust_p)then iflag(I)=1 endif endif endif if(ustmedian0(I).lt.cutust_c)then iflag(I)=5 endif endif 1100 continue N0126=0 do 1200 I=1,N if(iflag(I).eq.0.or.iflag(I).eq.1.or.iflag(I).eq.2. & or.iflag(I).eq.6)then N0126=N0126+1 ratgood(N0126)=ratio(I) ustgood(N0126)=fricv(I) endif 1200 continue ! nbinpoints=10 call average_bin(ustgood,ratgood,N0126,nbinpoints, & confidlevel,ustarbinmedian,ratiobinmean, & confidbinradius,nbins) ! do 1300 I=1,N_ori J=1 1400 delta=dabs(day(J)-dayfrac_ori(I)) if(delta.gt.1.0d-6)goto 1500 ialltimeflag(I)=iflag(J) fneenorm(I)=ratio(J) goto 1300 1500 if(J.eq.N)then ! daytime data ialltimeflag(I)=9999 fneenorm(I)=-9999.0d0 goto 1300 endif J=J+1 goto 1400 1300 continue return END c c subroutine outliers(Nsamp,datasample,iflag) c c Outlier detection. iflag=3, >3*sd; c iflags=4, <-3sd; else iflag=0 c double precision datasample(Nsamp) double precision sum, fmean, sd, dif, crit integer Nsamp,i,iflag(Nsamp) c sum=0.0d0 do 10 i=1, Nsamp sum=sum+datasample(i) 10 continue fmean=sum/dble(Nsamp) sum=0.0d0 do 20 i=1, Nsamp sum=sum+(datasample(i)-fmean)*(datasample(i)-fmean) 20 continue sd=dsqrt(sum/dble(Nsamp-1)) crit=3.0d0*sd do 30 i=1, Nsamp dif=datasample(i)-fmean if(dif.gt.crit)then iflag(i)=3 else if(dif.lt.-crit)then iflag(i)=4 else iflag(i)=0 endif endif 30 continue return end c subroutine average_moving(ustar0,ratio0,nsamp,nummov0) implicit double precision (a-h,l,o-z) dimension ustar(nsamp),ustar0(nsamp),ratio0(nsamp),ratio(nsamp) c do 5, i=1, nsamp ustar(i)=ustar0(i) ratio(i)=ratio0(i) 5 continue do 10 i=1, nsamp do 20 j=i+1, nsamp if(ustar(j)sample2, compmeans=1 c else compmeans=0 c sum=0.0d0 do 10 i=1, nsamp1 sum=sum+sample1(i) 10 continue fmean1=sum/dble(nsamp1) c sum=0.0d0 do 20 i=1, nsamp2 sum=sum+sample2(i) 20 continue fmean2=sum/dble(nsamp2) c sum=0.0d0 do 30 i=1, nsamp1 sum=sum+(sample1(i)-fmean1)*(sample1(i)-fmean1) 30 continue sd1=sum/dble(nsamp1-1) c sum=0.0d0 do 40 i=1, nsamp2 sum=sum+(sample2(i)-fmean2)*(sample2(i)-fmean2) 40 continue sd2=sum/dble(nsamp2-1) c if(ivarassump.eq.1)then ! Unequal variance assumption se1=sd1/dble(nsamp1) se2=sd2/dble(nsamp2) se=se1+se2 ndf=idint(se*se & /(se1*se1/dble(nsamp1-1)+se2*se2/dble(nsamp2-1))) else if(ivarassump.eq.2)then ! Equal variance assumption ndf=nsamp1+nsamp2-2 sp2=(sd1*dble(nsamp1-1)+ & sd2*dble(nsamp2-1))/dble(ndf) se=sp2/dble(nsamp1)+sp2/dble(nsamp2) else write(*,*)'ivarassump must be either:' write(*,*)'1 for unequal variance assumption or ' write(*,*)'2 for equal variance assumption ' stop endif endif tvalue=(fmean1-fmean2)/dsqrt(se) t0=tcritic(ndf,alpha) if(tvalue.lt.-t0)then icompmeans=-1 else if(tvalue.gt.t0)then icompmeans=1 else icompmeans=0 endif endif return end c c integer function meansdif(sample1, nsamp1, & sample2, nsamp2) implicit double precision (a-h,l,o-z) dimension sample1(nsamp1), sample2(nsamp2) parameter(critperc=0.95d0) c c This function compares the means of sample1 and sample2. c sample1 is the target sample, sample2 is the reference sample. c if sample1105%*sample2, mdif_pp=-1 c else compmeans=0 c sum=0.0d0 do 10 i=1, nsamp1 sum=sum+sample1(i) 10 continue fmean1=sum/dble(nsamp1) c sum=0.0d0 do 20 i=1, nsamp2 sum=sum+sample2(i) 20 continue fmean2=sum/dble(nsamp2) c tvalue=critperc*fmean2 if(fmean1.ge.tvalue)then mdif_pp=-1 else mdif_pp=0 endif return end c c double precision function tcritic(nsamp,alpha) c c this function reads the upper tail probability c double precision alpha integer nsamp,nsamp0,ndf,ilevel double precision t0(5) c t90, t925, t95, t975, t99; c c alpha =0.1, 0.075, 0.05, 0.025, 0.01 c if(dabs(alpha-0.1d0).lt.0.0001d0)then ilevel=1 else if(dabs(alpha-0.075d0).lt.0.0001d0)then ilevel=2 else if(dabs(alpha-0.05d0).lt.0.0001d0)then ilevel=3 else if(dabs(alpha-0.025d0).lt.0.0001d0)then ilevel=4 else if(dabs(alpha-0.01d0).lt.0.0001d0)then ilevel=5 else write(*,*)'Reset confidence level' stop endif endif endif endif endif c nsamp0=nsamp if(nsamp.gt.1000)then c write(*,*)"Degree of freedom exceeds maximum 1000" c write(*,*)"t-critical is taken as that for 1000 samples" nsamp0=1000 endif open(unit=6,file='t_tab_uptail.txt') read(6,*) 10 read(6,*)ndf,t0 if(ndf.ne.nsamp0)goto 10 tcritic=t0(ilevel) close(6) return end double precision function ran2() implicit double precision(a-h,o-z) parameter(im1=2147483563,im2=2147483399,am=1./im1, &imm1=im1-1,ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1= &12211,ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2d-7, &rnmx=1.0d0-eps) dimension iv(ntab) save iv,iy,idum2,idum data idum2/123456789/,iv/ntab*0/,iy/0/,idum/-1/ if(idum.le.0) then idum=max0(-idum,1) idum2=idum do 11 j=ntab+8,1,-1 k=idum/iq1 idum=ia1*(idum-k*iq1)-k*ir1 if(idum.lt.0) idum=idum+im1 if(j.le.ntab) iv(j)=idum 11 continue iy=iv(1) end if k=idum/iq1 idum=ia1*(idum-k*iq1)-k*ir1 k=idum2/iq2 idum2=ia2*(idum2-k*iq2)-k*ir2 if(idum2.lt.0) idum2=idum2+im2 j=1+iy/ndiv iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+imm1 ran2=dmin1(am*dble(iy),rnmx) return end c c This subroutine calculates R2 and root mean square error and index of cagreement subroutine rsq_rms(y1,y2,n,rsq,rms,agrind) implicit double precision (a-h,l,o-z) dimension y1(n),y2(n) sum=0.0d0 do 10 i=1,n sum=sum+(y1(i)-y2(i))*(y1(i)-y2(i)) 10 continue rms=dsqrt(sum/dble(n)) ymean1=0.0d0 ymean2=0.0d0 do 20 i=1,n ymean1=ymean1+y1(i) ymean2=ymean2+y2(i) 20 continue ymean1=ymean1/dble(n) ymean2=ymean2/dble(n) sum1=0.0d0 sum2=0.0d0 sum3=0.0d0 sum4=0.0d0 sum5=0.0d0 do 30 i=1,n sum1=(y1(i)-ymean1)*(y2(i)-ymean2)+sum1 sum2=(y1(i)-ymean1)*(y1(i)-ymean1)+sum2 sum3=(y2(i)-ymean2)*(y2(i)-ymean2)+sum3 sum4=(y1(i)-y2(i))*(y1(i)-y2(i))+sum4 sum5=(dabs(y2(i)-ymean1)+dabs(y1(i)-ymean1))* &(dabs(y2(i)-ymean1)+dabs(y1(i)-ymean1))+sum5 30 continue rsq=sum1/dsqrt(sum2*sum3) rsq=rsq*rsq agrind=1.0d0-sum4/sum5 return end c double precision function correlation(sample1,sample2,num1) implicit double precision (a-h,l,o-z) dimension sample1(num1),sample2(num1) c c This function calculates the correlation coefficient between sample1 c and sample2. lxy=0 lxx=0 lyy=0 sumx=0 sumy=0 k=0 do 10 i=1, num1 dif1=dabs(dabs(sample1(i))-9999.0d0) dif2=dabs(dabs(sample2(i))-9999.0d0) if((dif1.gt.0.001d0).and.(dif2.gt.0.001d0))then sumx=sumx+sample1(i) sumy=sumy+sample2(i) k=k+1 endif 10 continue if(k>2)then fmeanx=sumx/dble(k) fmeany=sumy/dble(k) do 20 i=1, num1 dif1=dabs(dabs(sample1(i))-9999) dif2=dabs(dabs(sample2(i))-9999) if((dif1.gt.0.001d0).and.(dif2.gt.0.001d0))then lxx=lxx+(sample1(i)-fmeanx)*(sample1(i)-fmeanx) lyy=lyy+(sample2(i)-fmeany)*(sample2(i)-fmeany) lxy=lxy+(sample1(i)-fmeanx)*(sample2(i)-fmeany) endif 20 continue r=lxy/dsqrt(lxx*lyy) if((k-2).ge.1000)then r0=0.062d0 else open(unit=7,file='/usr/users/l2g/ustar/corrtable95') read(7,*) 30 read(7,*)ndeg2,r2 if((k-2).eq.ndeg2)then r0=r2 ngetout=1 else if((k-2).gt.ndeg2)then r1=r2 ndeg1=ndeg2 ngetout=0 else r0=r2+dble(ndeg2-k+2)*(r1-r2)/dble(ndeg2-ndeg1) ngetout=1 endif endif if(ngetout.eq.0)goto 30 close(7) endif if(dabs(r) FCN NAME OF THE USER SUPPLIED FUNCTION SUBROUTINE C ==> N NUMBER OF OBSERVATIONS C ==> M COLUMNS OF DATA IN THE EXPLANATORY VARIABLE C ==> NP NUMBER OF PARAMETERS C ==> NQ NUMBER OF RESPONSES PER OBSERVATION C <==> BETA FUNCTION PARAMETERS C ==> Y RESPONSE VARIABLE C ==> LDY LEADING DIMENSION OF ARRAY Y C ==> X EXPLANATORY VARIABLE C ==> LDX LEADING DIMENSION OF ARRAY X C ==> WE "EPSILON" WEIGHTS C ==> LDWE LEADING DIMENSION OF ARRAY WE C ==> LD2WE SECOND DIMENSION OF ARRAY WE C ==> WD "DELTA" WEIGHTS C ==> LDWD LEADING DIMENSION OF ARRAY WD C ==> LD2WD SECOND DIMENSION OF ARRAY WD C ==> IFIXB INDICATORS FOR "FIXING" PARAMETERS (BETA) C ==> IFIXX INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X) C ==> LDIFX LEADING DIMENSION OF ARRAY IFIXX C ==> JOB TASK TO BE PERFORMED C ==> NDIGIT GOOD DIGITS IN SUBROUTINE FUNCTION RESULTS C ==> TAUFAC TRUST REGION INITIALIZATION FACTOR C ==> SSTOL SUM OF SQUARES CONVERGENCE CRITERION C ==> PARTOL PARAMETER CONVERGENCE CRITERION C ==> MAXIT MAXIMUM NUMBER OF ITERATIONS C ==> IPRINT PRINT CONTROL C ==> LUNERR LOGICAL UNIT FOR ERROR REPORTS C ==> LUNRPT LOGICAL UNIT FOR COMPUTATION REPORTS C ==> STPB STEP SIZES FOR FINITE DIFFERENCE DERIVATIVES WRT BETA C ==> STPD STEP SIZES FOR FINITE DIFFERENCE DERIVATIVES WRT DELTA C ==> LDSTPD LEADING DIMENSION OF ARRAY STPD C ==> SCLB SCALE VALUES FOR PARAMETERS BETA C ==> SCLD SCALE VALUES FOR ERRORS DELTA IN EXPLANATORY VARIABLE C ==> LDSCLD LEADING DIMENSION OF ARRAY SCLD C <==> WORK DOUBLE PRECISION WORK VECTOR C ==> LWORK DIMENSION OF VECTOR WORK C <== IWORK INTEGER WORK VECTOR C ==> LIWORK DIMENSION OF VECTOR IWORK C <== INFO STOPPING CONDITION C PARAMETERS SPECIFYING MAXIMUM PROBLEM SIZES HANDLED BY THIS DRIVER C MAXN MAXIMUM NUMBER OF OBSERVATIONS C MAXM MAXIMUM NUMBER OF COLUMNS IN EXPLANATORY VARIABLE C MAXNP MAXIMUM NUMBER OF FUNCTION PARAMETERS C MAXNQ MAXIMUM NUMBER OF RESPONSES PER OBSERVATION C PARAMETER DECLARATIONS AND SPECIFICATIONS INTEGER LDIFX,LDSCLD,LDSTPD,LDWD,LDWE,LDX,LDY,LD2WD,LD2WE, + LIWORK,LWORK,MAXM,MAXN,MAXNP,MAXNQ PARAMETER (MAXM=25,MAXN=10000,MAXNP=30,MAXNQ=1, + LDY=MAXN,LDX=MAXN, + LDWE=1,LD2WE=1,LDWD=1,LD2WD=1, + LDIFX=MAXN,LDSTPD=1,LDSCLD=1, + LWORK=18 + 11*MAXNP + MAXNP**2 + MAXM + MAXM**2 + + 4*MAXN*MAXNQ + 6*MAXN*MAXM + 2*MAXN*MAXNQ*MAXNP + + 2*MAXN*MAXNQ*MAXM + MAXNQ**2 + + 5*MAXNQ + MAXNQ*(MAXNP+MAXM) + LDWE*LD2WE*MAXNQ, + LIWORK=20+MAXNP+MAXNQ*(MAXNP+MAXM)) C VARIABLE DECLARATIONS INTEGER I,INFO,IPRINT,J,JOB,L,LUNERR,LUNRPT,M,MAXIT,N, + NDIGIT,NP,NQ,icond INTEGER IFIXB(MAXNP),IFIXX(LDIFX,MAXM),IWORK(LIWORK) DOUBLE PRECISION PARTOL,SSTOL,TAUFAC DOUBLE PRECISION BETA(MAXNP),SCLB(MAXNP),SCLD(LDSCLD,MAXM), + STPB(MAXNP),STPD(LDSTPD,MAXM), + WD(LDWD,LD2WD,MAXM),WE(LDWE,LD2WE,MAXNQ), + WORK(LWORK),X(LDX,MAXM),Y(LDY,MAXNQ) c double precision scrndnee(icond),scrndta(icond),scrndts(icond), +param1,param2,param3,param4,param5 EXTERNAL FCN c C SPECIFY DEFAULT VALUES FOR DODRC ARGUMENTS WE(1,1,1) = -1.0D0 WD(1,1,1) = -1.0D0 IFIXB(1) = -1 IFIXX(1,1) = -1 JOB = 00023 NDIGIT = -1 TAUFAC = -1.0D0 SSTOL = -1.0D0 PARTOL = -1.0D0 MAXIT = -1 c IPRINT = -1 IPRINT=0 LUNERR = -1 LUNRPT = -1 STPB(1) = -1.0D0 STPD(1,1) = -1.0D0 SCLB(1) = -1.0D0 SCLD(1,1) = -1.0D0 MAXIT = 20000 C SET UP ODRPACK REPORT FILES LUNERR = 9 LUNRPT = 9 c M=2 NP=5 NQ=1 BETA(1)=1.0d0 BETA(2)=0.07d0 BETA(3)=0.5d0 BETA(4)=1.0d0 BETA(5)=0.07d0 N=icond do 885 I=1,N X(I,1)=scrndta(I) X(I,2)=scrndts(I) Y(I,1)=scrndnee(I) 885 continue C READ PROBLEM DATA, AND SET NONDEFAULT VALUE FOR ARGUMENT IFIXX DO 10 I=1,N DO 15 J=1, M IFIXX(I,J) = 1 15 CONTINUE 10 CONTINUE 60 CALL DODRC(FCN, + N,M,NP,NQ, + BETA, + Y,LDY,X,LDX, + WE,LDWE,LD2WE,WD,LDWD,LD2WD, + IFIXB,IFIXX,LDIFX, + JOB,NDIGIT,TAUFAC, + SSTOL,PARTOL,MAXIT, + IPRINT,LUNERR,LUNRPT, + STPB,STPD,LDSTPD, + SCLB,SCLD,LDSCLD, + WORK,LWORK,IWORK,LIWORK, + INFO) i1=mod(INFO,10) i2=(mod(INFO,100)-i1)/10 i3=(mod(INFO,1000)-mod(INFO,100))/100 i4=(mod(INFO,10000)-mod(INFO,1000))/1000 i5=(INFO-mod(INFO,10000))/10000 if(i5.gt.0)then write(*,*)'Fatal error in regression',i5 stop endif c write(*,*)INFO if(INFO.lt.4.or.INFO.eq.1001)goto 75 BETA(1)=1.0d0*(ran2()+0.5d0) BETA(2)=0.07d0*(ran2()+0.5d0) BETA(3)=ran2() BETA(4)=1.0d0*(ran2()+0.5d0) BETA(5)=0.07d0*(ran2()+0.5d0) goto 60 c 75 param1=BETA(1) param2=BETA(2) param3=BETA(3) param4=BETA(4) param5=BETA(5) return END c SUBROUTINE FCN(N,M,NP,NQ, + LDN,LDM,LDNP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + IDEVAL,F,FJACB,FJACD, + ISTOP) C SUBROUTINE ARGUMENTS C ==> N NUMBER OF OBSERVATIONS C ==> M NUMBER OF COLUMNS IN EXPLANATORY VARIABLE C ==> NP NUMBER OF PARAMETERS C ==> NQ NUMBER OF RESPONSES PER OBSERVATION C ==> LDN LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING N C ==> LDM LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING M C ==> LDNP LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING NP C ==> BETA CURRENT VALUES OF PARAMETERS C ==> XPLUSD CURRENT VALUE OF EXPLANATORY VARIABLE, I.E., X + DELTA C ==> IFIXB INDICATORS FOR "FIXING" PARAMETERS (BETA) C ==> IFIXX INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X) C ==> LDIFX LEADING DIMENSION OF ARRAY IFIXX C ==> IDEVAL INDICATOR FOR SELECTING COMPUTATION TO BE PERFORMED C <== F PREDICTED FUNCTION VALUES C <== FJACB JACOBIAN WITH RESPECT TO BETA C <== FJACD JACOBIAN WITH RESPECT TO ERRORS DELTA C <== ISTOP STOPPING CONDITION, WHERE C 0 MEANS CURRENT BETA AND X+DELTA WERE C ACCEPTABLE AND VALUES WERE COMPUTED SUCCESSFULLY C 1 MEANS CURRENT BETA AND X+DELTA ARE C NOT ACCEPTABLE; ODRPACK SHOULD SELECT VALUES C CLOSER TO MOST RECENTLY USED VALUES IF POSSIBLE C -1 MEANS CURRENT BETA AND X+DELTA ARE C NOT ACCEPTABLE; ODRPACK SHOULD STOP C INPUT ARGUMENTS, NOT TO BE CHANGED BY THIS ROUTINE: INTEGER I,IDEVAL,ISTOP,L,LDIFX,LDM,LDN,LDNP,M,N,NP,NQ DOUBLE PRECISION BETA(NP),XPLUSD(LDN,M) INTEGER IFIXB(NP),IFIXX(LDIFX,M) C OUTPUT ARGUMENTS: DOUBLE PRECISION F(LDN,NQ),FJACB(LDN,LDNP,NQ),FJACD(LDN,LDM,NQ) c local variables DOUBLE PRECISION TERM1,TERM2,zerotemp C LOCAL VARIABLES INTRINSIC dexp,dlog parameter(zerotemp=273.15d0) c C CHECK FOR UNACCEPTABLE VALUES FOR THIS PROBLEM c IF(BETA(1).LT.0.0D0.or.BETA(2).lt.0.0d0.or.BETA(3) &.lt.0.0d0.or.BETA(4).lt.0.0d0.or.BETA(5).lt.0.0d0)then ISTOP = 1 RETURN ELSE if(BETA(3).gt.1.0d0)then ISTOP=1 RETURN else ISTOP = 0 endif END IF C COMPUTE PREDICTED VALUES c XPLUSD(I,1): air temperature c XPLUSD(I,2): soil temperature c IF (MOD(IDEVAL,10).GE.1) THEN DO 110 L = 1,NQ DO 100 I = 1,N TERM1=BETA(3)*(XPLUSD(I,1)-zerotemp)+ &(1.0d0-BETA(3))*(XPLUSD(I,2)-zerotemp) TERM1=BETA(2)*TERM1 TERM2=BETA(5)*(XPLUSD(I,2)-zerotemp) F(I,L)=BETA(1)*dexp(TERM1)+BETA(4)*dexp(TERM2) 100 CONTINUE 110 CONTINUE END IF C COMPUTE DERIVATIVES WITH RESPECT TO BETA IF (MOD(IDEVAL/10,10).GE.1) THEN DO 210 L = 1,NQ DO 200 I = 1,N TERM1=BETA(3)*(XPLUSD(I,1)-zerotemp)+ &(1.0d0-BETA(3))*(XPLUSD(I,2)-zerotemp) TERM1=BETA(2)*TERM1 TERM2=BETA(5)*(XPLUSD(I,2)-zerotemp) FJACB(I,1,L)=dexp(TERM1) FJACB(I,4,L)=dexp(TERM2) FJACB(I,2,L)=BETA(1)*dexp(TERM1)* &(BETA(3)*(XPLUSD(I,1)-zerotemp)+ &(1.0d0-BETA(3))*(XPLUSD(I,2)-zerotemp)) FJACB(I,5,L)= &BETA(4)*dexp(TERM2)*(XPLUSD(I,2)-zerotemp) FJACB(I,3,L)=BETA(1)*dexp(TERM1)*BETA(2)*( &XPLUSD(I,1)-XPLUSD(I,2)) 200 CONTINUE 210 CONTINUE END IF RETURN END c double precision function ecoresp(airtemp,soiltemp, ¶m1,param2,param3,param4,param5) implicit double precision (a-h,l,o-z) parameter(zerotemp=273.15d0) c c air temperature (airtemp) and soil temperature (soiltemp) in K c TERM1=param3*(airtemp-zerotemp)+ &(1.0d0-param3)*(soiltemp-zerotemp) TERM1=param2*TERM1 TERM2=param5*(soiltemp-zerotemp) ecoresp=param1*dexp(TERM1)+param4*dexp(TERM2) return end