DOUBLE PRECISION FUNCTION SYS_SS_POT(R,V,NV,LSETUP) IMPLICIT NONE DOUBLE PRECISION R,V INTEGER NV LOGICAL LSETUP C C include 'PhysicalConst.h' C C ... Local Variables C double precision c6,c8,c10,CA,beta,gamma,as,bs,at,bt, > re1,re3,shift1,shift3,rdatas,vdatas,rdatat,vdatat save c6,c8,c10,CA,beta,gamma,as,bs,at,bt, > re1,re3,shift1,shift3,rdatas,vdatas,rdatat,vdatat double precision hartree,bohra save hartree,bohra double precision correction_ss double precision rs1,vs1,rs2,vs2,rs3,vs3,rt1,vt1,rt2,vt2, > rt3,vt3,res,des,ret,det,rdat,vdats,vdatt double precision vs,vt,vav,vdiff,vex,vdisp,amp,vshift, > rfit,rjoin,rfinal,r0,rl,ri,r1,r2,v1,v2 integer np,ns,nt,ndatas,ndatat save np,ns,nt,ndatas,ndatat integer i,j,k,ni,nsl,ntl,ndat C integer ns1,ns2,ns3,nt1,nt2,nt3,nmax parameter (ns1=3,ns2=42,ns3=41,nt1=8,nt2=12,nt3=11,nmax=500) CFLA parameter (ns1=3,ns2=42,ns3=41,nt1=8,nt2=13,nt3=12,nmax=500) C dimension rs1(ns1),vs1(ns1),rs2(ns2),vs2(ns2),rs3(ns3),vs3(ns3), > rt1(nt1),vt1(nt1),rt2(nt2),vt2(nt2),rt3(nt3),vt3(nt3) dimension rdatas(nmax),vdatas(nmax),rdatat(nmax),vdatat(nmax), > rdat(nmax),vdats(nmax),vdatt(nmax) logical do_setup data do_setup /.true./ save do_setup logical order,norepeat C C output flags for potential routine: C output_data - output raw potential data C output_fitex - exchange fitting constant A as a function of R C output_v - output final potentials to be used in calculations C output_vmesh - output interpolated mesh of final potentials + C correction included C logical output_data, output_fitex, output_v, output_vmesh data output_data /.false./ data output_fitex /.false./ data output_v /.false./ data output_vmesh /.false./ c logical output_data=.false., output_fitex=.false., c > output_v=.false., output_vmesh=.false. C C=====Open output files================================================= C if (LSETUP .and. do_setup) then !open(unit=55,status='old',file='testdat') if (output_data) then open(unit=50,status='unknown',file='vdata/v_singlet.raw') open(unit=51,status='unknown',file='vdata/v_triplet.raw') open(unit=52,status='unknown',file='vdata/v_st.raw') endif if (output_fitex) then open(unit=60,status='unknown',file='vdata/fit_ex.dat') open(unit=61,status='unknown',file='vdata/fit_disp.dat') open(unit=62,status='unknown',file='vdata/fit_vs.dat') open(unit=63,status='unknown',file='vdata/fit_vt.dat') open(unit=64,status='unknown',file='vdata/fit_A.dat') endif if (output_v) then open(unit=70,status='unknown',file='vdata/v_singlet.dat') open(unit=71,status='unknown',file='vdata/v_triplet.dat') open(unit=72,status='unknown',file='vdata/v_singlet.sdat') open(unit=73,status='unknown',file='vdata/v_triplet.sdat') endif if (output_vmesh) then open(unit=80,status='unknown',file='vdata/vmesh_s.dat') open(unit=81,status='unknown',file='vdata/vmesh_t.dat') open(unit=82,status='unknown',file='vdata/vmesh_s.sdat') open(unit=83,status='unknown',file='vdata/vmesh_t.sdat') endif endif C C=====Data for ^7Li2==================================================== C C Physical constants used in potential routines C DATA HARTREE /219474.630680D0/, BOHRA /0.5291772490D0/ C C np is the degree of polynomial used in the akima spline interpolation C np=3 C C rfit= point at which the exchange form is fit to the data and at C which is calculated the shift in the data required to match the C short and long range parts C rjoin= point at which the short range data is smoothly connected C to the long range form C rfinal= last point for short range data set used in the calculations C DATA rfit /15.5D0/, rjoin /17.D0/,rfinal /20.0d0/ C C --- Long range: C Yan, Babb, Dalgarno and Drake, Phys. Rev. A 54, 2824 (1996) C (in atomic units) C DATA C6 /1393.39D0/, C8 /83425.8D0/, C10 /7372100.0D0/ C C --- Long range exchange: C Cote, Dalgarno, and Jamieson, Phys. Rev. A 50, 399 (1994) C (in atomic units) C DATA beta /1.25902D0/, gamma /4.55988D0/ C C --- Short-range singlet potential data C C De determined by Zemke and Stwalley, J. Phys. Chem. 97, 2053 (1993) C (in units of cm-1) C DATA Des /8516.61D0/ C C Section with E>0 : small R C 1st 2 pts from Konowalow and Olson J. Chem. Phys. 71, 450 (1979) C 3rd pt from Schmidt-Mink et al Chem. Phys. 92, 263 (1985) C (in atomic units) C DATA rs1 / 2.75, 3.00, 3.25 / C DATA vs1 / 0.061866D0, 0.036411D0, 0.013546D0 / C C Section with E<0 : RKR potential C from Barakat et al Chem. Phys. 102, 215 (1986) C note: misprints in data corrected, see Cote et al C Phys. Rev A 50 399 (1994) C (in units of angstroms and cm-1) C DATA res / 2.673247D0 / C DATA rs2 / 1.823201D0, 1.823462D0, 1.824057D0, & 1.825118D0, 1.826727D0, 1.828918D0, & 1.831708D0, 1.834867D0, 1.838851D0, & 1.843462D0, 1.848639D0, 1.854376D0, & 1.860716D0, 1.867431D0, 1.874904D0, & 1.882864D0, 1.891368D0, 1.900497D0, & 1.910173D0, 1.920282D0, 1.931173D0, & 1.942717D0, 1.954984D0, 1.967947D0, & 1.981710D0, 1.996354D0, 2.011906D0, & 2.028431D0, 2.046073D0, 2.064963D0, & 2.085263D0, 2.107177D0, 2.130967D0, & 2.156977D0, 2.185675D0, 2.217724D0, & 2.254114D0, 2.296439D0, 2.347580D0, & 2.413867D0, 2.517207D0, 2.673247D0 / C DATA vs2 / 8514.7777D0, 8507.8421D0, 8492.0440D0, & 8463.9645D0, 8421.6125D0, 8364.3068D0, & 8292.0294D0, 8205.2323D0, 8104.4730D0, & 7990.4162D0, 7863.7083D0, 7724.9765D0, & 7574.8736D0, 7413.8431D0, 7242.5556D0, & 7061.4199D0, 6870.8931D0, 6671.3979D0, & 6463.3140D0, 6246.9482D0, 6022.6578D0, & 5790.7056D0, 5551.3992D0, 5304.9322D0, & 5051.5343D0, 4791.4274D0, 4524.7756D0, & 4251.7309D0, 3972.4624D0, 3687.1094D0, & 3395.7978D0, 3090.6412D0, 2795.7419D0, & 2487.1914D0, 2173.0721D0, 1853.4573D0, & 1528.4128D0, 1197.9974D0, 862.2642D0, & 521.2611D0, 175.0320D0, 0.0000D0 / C DATA rs3 / 12.639382D0, 9.863785D0, 8.495553D0, & 7.657218D0, 7.090991D0, 6.675227D0, & 6.352341D0, 6.089904D0, 5.870023D0, & 5.680896D0, 5.514781D0, 5.366746D0, CX & 5.580896D0, 5.514781D0, 5.366746D0, & 5.232506D0, 5.109851D0, 4.996566D0, & 4.890860D0, 4.791595D0, 4.697788D0, CX & 4.490860D0, 4.791595D0, 4.697788D0, & 4.608426D0, 4.522846D0, 4.440679D0, & 4.361370D0, 4.284400D0, 4.209387D0, & 4.135091D0, 4.064160D0, 3.993263D0, & 3.923163D0, 3.853594D0, 3.784295D0, & 3.714965D0, 3.645290D0, 3.574908D0, & 3.503385D0, 3.430172D0, 3.354534D0, & 3.275416D0, 3.191166D0, 3.098852D0, & 2.992096D0, 2.848951D0 / C DATA vs3 / 8514.7777D0, 8507.8421D0, 8492.0440D0, & 8463.9645D0, 8421.6125D0, 8364.3068D0, & 8292.0294D0, 8205.2323D0, 8104.4730D0, & 7990.4162D0, 7863.7083D0, 7724.9765D0, & 7574.8736D0, 7413.8431D0, 7242.5556D0, & 7061.4199D0, 6870.8931D0, 6671.3979D0, & 6463.3140D0, 6246.9482D0, 6022.6578D0, & 5790.7056D0, 5551.3992D0, 5304.9322D0, & 5051.5343D0, 4791.4274D0, 4524.7756D0, & 4251.7309D0, 3972.4624D0, 3687.1094D0, & 3395.7978D0, 3090.6412D0, 2795.7419D0, & 2487.1914D0, 2173.0721D0, 1853.4573D0, & 1528.4128D0, 1197.9974D0, 862.2642D0, & 521.2611D0, 175.0320D0 / C C --- Short-range triplet potential data C C De determined by Linton et al, J. Mol. Spect. 196, 20 (1999) C (in units of cm-1) C DATA Det /333.69D0/ C C Section with E>0 : small R C 1st pt from Konowalow and Olson J. Chem. Phys. 81, 4534 (1984) C last 7 pts from Schmidt-Mink et al Chem. Phys. 92, 263 (1985) C (in atomic units) C DATA rt1 / 3.00D0, 3.25D0, 3.50D0, & 4.00D0, 4.50D0, 5.00D0, & 5.50D0, 6.00D0 / C DATA vt1 / 0.104471D0, 0.082074D0, 0.063725D0, & 0.036778D0, 0.020006D0, 0.010053D0, & 0.004370D0, 0.001258D0 / C C Section with E<0 : RKR potential C from Linton et al J. Mol. Spect. 196, 20 (1999) C (in units of angstroms and cm-1) C DATA ret / 4.173D0 / C DATA rt2 / 3.356D0, 3.358D0, 3.365D0, & 3.377D0, 3.395D0, 3.422D0, & 3.458D0, 3.505D0, 3.571D0, & 3.668D0, 3.846D0, 4.173D0 / CFLA & 3.668D0, 3.846D0, 4.159D0,4.173D0 / C DATA vt2 / 333.269D0, 330.170D0, 322.155D0, & 308.098D0, 287.665D0, 260.837D0, & 227.679D0, 188.240D0, 142.523D0, CFLA & 90.453D0, 31.857D0, 0.024D0, 0.000D0 / & 90.453D0, 31.857D0, 0.000D0 / C DATA rt3 / 16.052D0, 11.392D0, 9.441D0, & 8.297D0, 7.501D0, 6.885D0, & 6.373D0, 5.922D0, 5.503D0, & 5.092D0, 4.630D0 / CFLA & 5.092D0, 4.630D0 , 4.185D0/ C DATA vt3 / 333.269D0, 330.170D0, 322.155D0, & 308.098D0, 287.665D0, 260.837D0, & 227.679D0, 188.240D0, 142.523D0, & 90.453D0, 31.857D0/ CFLA & 90.453D0, 31.857D0, 0.024D0 / C C=============Initial setup============================================= C IF(LSETUP) THEN C if (do_setup) then !namelist / Adjust_Pot / shift1,shift3,c6 ! read(55,Adjust_Pot) ! write(66,990) ! write(66,Adjust_Pot) shift1=0.0d0 shift3=-7d-7 c6=1393.39d0 write(*,*)'--------------------' write(*,*) shift1,shift3,c6 C C --- Calculate total number of potential data points C ns=ns1+ns2+ns3 nt=nt1+nt2+nt3-2 if (ns.gt.nmax.or.nt.gt.nmax) then write(*,1000) ns,nt,nmax endif C C --- Order each set of data, in ascending order, into an array C and smooth data by multiplying by R^6 for interpolation. C Convert data to atomic units. C re1=res/bohra re3=ret/bohra C do i=1,ns1 rdatas(i)=rs1(i) vdatas(i)=vs1(i) vdatas(i)=vdatas(i)*rdatas(i)**6 enddo do i=1,ns2 j=ns1+i rdatas(j)=rs2(i)/bohra vdatas(j)=(vs2(i)-des)/hartree vdatas(j)=vdatas(j)*rdatas(j)**6 enddo do i=1,ns3 j=ns1+ns2+i rdatas(j)=rs3(i)/bohra vdatas(j)=(vs3(i)-des)/hartree vdatas(j)=vdatas(j)*rdatas(j)**6 enddo C order=.true. do i=2,ns if (rdatas(i).eq.rdatas(i-1)) then write(*,1010) i-1,i stop elseif (rdatas(i).lt.rdatas(i-1)) then order=.false. endif enddo if (.not. order) then call sort (rdatas,vdatas,ns) endif C do i=1,nt1 rdatat(i)=rt1(i) vdatat(i)=vt1(i) vdatat(i)=vdatat(i)*rdatat(i)**6 enddo C C Note: 1st and 2nd points of the RKR data had to be omitted C to prevent oscillations in the spline fit between C the two sets of data C do i=3,nt2 j=nt1-2+i rdatat(j)=rt2(i)/bohra vdatat(j)=(vt2(i)-det)/hartree vdatat(j)=vdatat(j)*rdatat(j)**6 enddo do i=1,nt3 j=nt1+nt2-2+i rdatat(j)=rt3(i)/bohra vdatat(j)=(vt3(i)-det)/hartree vdatat(j)=vdatat(j)*rdatat(j)**6 enddo C order=.true. do i=2,nt if (rdatat(i).eq.rdatat(i-1)) then write(*,1010) i-1,i stop elseif (rdatat(i).lt.rdatat(i-1)) then order=.false. endif enddo if (.not. order) then call sort (rdatat,vdatat,nt) endif C if (output_data) then do i=1,ns write(50,1100) rdatas(i),vdatas(i)/rdatas(i)**6 enddo do i=1,nt write(51,1100) rdatat(i),vdatat(i)/rdatat(i)**6 enddo C do i=1,ns rdat(i)=rdatas(i) vdats(i)=vdatas(i) do j=1,nt if (rdatas(i).eq.rdatat(j)) then vdatt(i)=vdatat(j) else call akimasv(np,nt,rdatat,vdatat,1,rdat(i),vt) vdatt(i)=vt endif enddo enddo C k=ns do i=1,nt norepeat=.true. do j=1,ns if (rdatat(i).eq.rdatas(j)) then norepeat=.false. endif enddo if (norepeat) then k=k+1 rdat(k)=rdatat(i) vdatt(k)=vdatat(i) call akimasv(np,ns,rdatas,vdatas,1,rdat(k),vs) vdats(k)=vs endif enddo ndat=k call sort2 (rdat,vdats,vdatt,ndat) C do i=1,ndat write(52,1110) rdat(i),vdats(i),vdatt(i) enddo C endif C C --- Fit Smirnov and Chibisov exchange form: A exp(-beta R) R^gamma C to the exchange energy determined from the potential data to C obtain the amplitude A C C Calculate A at the chosen fitting point C C Interpolate potential data using akima spline interpolation C call akimasv(np,ns,rdatas,vdatas,1,rfit,vs) call akimasv(np,nt,rdatat,vdatat,1,rfit,vt) C C Determine fitting constant A C vt=vt/rfit**6 vs=vs/rfit**6 vex=log((vt-vs)/2.0d0) CA=vex+beta*rfit-gamma*log(rfit) CA=exp(CA) C C --- Calculate the average of vs and vt, which should correspond C with vdisp C vav=(vt+vs)/2.0d0 vdisp=-c6/rfit**6-c8/rfit**8-c10/rfit**10 vshift=vdisp-vav C C --- Shift short range data by calculated difference to match the C long range form C nsl=0 do i=1,ns vdatas(i)=vdatas(i)+vshift*rdatas(i)**6 if(nsl.eq.0.and.rdatas(i).ge.rfit) nsl=i enddo ntl=0 do i=1,nt vdatat(i)=vdatat(i)+vshift*rdatat(i)**6 if(ntl.eq.0.and.rdatat(i).ge.rfit) ntl=i enddo if(nsl.eq.0.or.ntl.eq.0) then write(*,1020) rfit stop endif C C --- Connect short range data onto long range form using akima spline C to smooth across the connection C call akimasv(np,ns,rdatas,vdatas,1,rfit,vs) call akimasv(np,nt,rdatat,vdatat,1,rfit,vt) rdatas(nsl)=rfit rdatat(ntl)=rfit vdatas(nsl)=vs vdatat(ntl)=vt C ni=int((rfinal-rjoin)/0.1d0)+1 ndatas=nsl+ni ndatat=ntl+ni ri=rjoin do i=1,ni vex=CA*exp(-beta*ri)*ri**gamma vdisp=-c6/ri**6-c8/ri**8-c10/ri**10 C rdatas(nsl+i)=ri rdatat(ntl+i)=ri vdatas(nsl+i)=(vdisp-vex)*ri**6 vdatat(ntl+i)=(vdisp+vex)*ri**6 C ri=ri+0.1d0 enddo rfinal=rdatas(ndatas) C C Extrapolate short range data using the form: aexp(-br) C r1=rdatas(1) r2=rdatas(1)+0.1d0 call akimasv(np,ndatas,rdatas,vdatas,1,r1,v1) call akimasv(np,ndatas,rdatas,vdatas,1,r2,v2) if (v1*v2.lt.0.0d0) then write(*,1030) stop endif v1=v1/r1**6 v2=v2/r2**6 bs=log(v1/v2)/(r2-r1) as=v2*exp(bs*r2) C r1=rdatat(1) r2=rdatat(1)+0.1d0 call akimasv(np,ndatat,rdatat,vdatat,1,r1,v1) call akimasv(np,ndatat,rdatat,vdatat,1,r2,v2) if (v1*v2.lt.0.0d0) then write(*,1030) stop endif v1=v1/r1**6 v2=v2/r2**6 bt=log(v1/v2)/(r2-r1) at=v2*exp(bt*r2) C C Output potential details/parameters C write(*,1040) c6,c8,c10 write(*,1050) rfit,CA,beta,gamma write(*,1060) vshift,vshift*hartree write(*,1070) rjoin,rfinal write(*,1080) as,bs,at,bt C C Output test data for constructing the potentials C if (output_fitex) then C r0=rdat(1) rfinal=rdat(ndat) ri=r0 ni=int((rfinal-r0)/0.01d0)+1 do i=1,ni C call akimasv(np,ndat,rdat,vdats,1,ri,vs) call akimasv(np,ndat,rdat,vdatt,1,ri,vt) vs=vs/ri**6 vt=vt/ri**6 C vav=(vt+vs)/2.0d0 vdiff=(vt-vs)/2.0d0 C vdisp=-c6/ri**6-c8/ri**8-c10/ri**10 vex=CA*exp(-beta*ri)*ri**gamma C amp=log(abs(vdiff))+beta*ri-gamma*log(ri) amp=exp(amp) C write(60,1120) ri,vdiff,vex,vex/vdiff write(61,1120) ri,vav,vdisp,vdisp/vav write(62,1120) ri,vs,vdisp-vex,(vdisp-vex)/vs write(63,1120) ri,vt,vdisp+vex,(vdisp+vex)/vt write(64,1100) ri,amp C ri=ri+0.01d0 C enddo endif C C Output potentials used in the calculations C if (output_v) then do i=1,ndatas write(70,1100) rdatas(i),vdatas(i)/rdatas(i)**6 write(72,1100) rdatas(i),vdatas(i) enddo do i=1,ndatat write(71,1100) rdatat(i),vdatat(i)/rdatat(i)**6 write(73,1100) rdatat(i),vdatat(i) enddo endif C C Output potentials as a function of R C if (output_vmesh) then r0=1.0d0 rl=rfinal+10.0d0 ri=r0 ni=int((rl-r0)/0.01d0)+1 do i=1,ni call mpot(np,ndatas,rdatas,vdatas,c6,c8,c10, > CA,beta,gamma,as,bs,1,ri,vs) call mpot(np,ndatat,rdatat,vdatat,c6,c8,c10, > CA,beta,gamma,at,bt,2,ri,vt) C vs=vs+correction_ss(ri,shift1,c6,re1) vt=vt+correction_ss(ri,shift3,c6,re3) C write(80,1100) ri,vs write(81,1100) ri,vt write(82,1100) ri,vs*ri**6 write(83,1100) ri,vt*ri**6 C ri=ri+0.01d0 C enddo endif C if(output_data) then close(50) close(51) endif if(output_fitex) then close(60) close(61) close(62) close(63) close(64) endif if(output_v) then close(70) close(71) close(72) close(73) endif if(output_vmesh) then close(80) close(81) close(82) close(83) endif C endif C do_setup=.not.do_setup SYS_SS_POT=0.0d0 C RETURN ENDIF C C======================================================================= C IF((NV.LE.0).OR.(NV.GT.2)) THEN WRITE(*,900) STOP END IF C GO TO (10,20),NV C C======================================================================= C C --- SECTION FOR NV=1 -- Calculate singlet potential C 10 call mpot(np,ndatas,rdatas,vdatas,c6,c8,c10, > CA,beta,gamma,as,bs,nv,r,v) C v=v+correction_ss(r,shift1,c6,re1) SYS_SS_POT=V C RETURN C C======================================================================= C C --- SECTION FOR NV=2 -- Calculate triplet potential C 20 call mpot(np,ndatat,rdatat,vdatat,c6,c8,c10, > CA,beta,gamma,at,bt,nv,r,v) C v=v+correction_ss(r,shift3,c6,re3) SYS_SS_POT=V C RETURN C C======================================================================= C C --- FORMAT STATEMENTS ------------------------------------------------ C 900 FORMAT(//,' Program HALTED in Potential Subroutine -- NV OUT ', X 'OF RANGE') C 990 FORMAT(//,'MOLECULAR POTENTIAL DATA') 1000 format(//,'ERROR in Potential Subroutine -- no. of data points', X 'exceeds array size:',3i6) 1010 format(//,'ERROR in Potential Subroutine -- identical data for', X 'pair:',2i6) 1020 format(//,'ERROR in Potential Subroutine -- rfit (au)=',d14.6, X 'exceeds short range data') 1030 format(//,'ERROR in Potential Subroutine -- data points used in', X 'short range extrapolation of data cross a zero') 1040 format(//,'Long range dispersion coefficients (au):',/, X ' C6=',d14.6,' C8=',d14.6,' C10=',d14.6) 1050 format(//,'Long range exchange, Aexp(-beta R)R^gamma (au):',/, X ' fit to short range data at R=',d14.6,/, X ' A=',d14.6,' beta=',d14.6,' gamma=',d14.6) 1060 format(/,'Shift applied to short range data to match to', X ' long range form:',/, X d14.6,' au',' or',d14.6,' cm-1') 1070 format(//,'Long range form assumed at R (au)=',d14.6,/, X 'Short range data extended to R (au)=',d14.6) 1080 format(//,'Short range extrapolation, aexp(-bR) (au)',/, X 'Singlet: a=',d14.6,' b=',d14.6,/, X 'Triplet: a=',d14.6,' b=',d14.6) 1100 format(2d16.8) 1110 format(3d16.8) 1120 format(4d16.8) C END C======================================================================= C======================================================================= C subroutine mpot(np,ndata,rdata,vdata,c6,c8,c10, > CA,beta,gamma,a,b,nv,r,v) C double precision rdata,vdata,c6,c8,c10,CA,beta,gamma,a,b,r,v integer np,ndata,nv dimension rdata(ndata),vdata(ndata) C C ... Local Variables C double precision vex,vdisp C if(r.gt.0.0d0 .and. r.lt.rdata(1)) then C C Extrapolate short range data using the form: aexp(-br) C v=a*exp(-b*r) elseif(r.ge.rdata(1) .and. r.le.rdata(ndata)) then C C Interpolate short range data C call akimasv(np,ndata,rdata,vdata,1,r,v) v=v/r**6 elseif(r.gt.rdata(ndata)) then C C Calculate long range form C vex=CA*exp(-beta*r)*r**gamma vdisp=-c6/r**6-c8/r**8-c10/r**10 if(nv.eq.1) then v=vdisp-vex elseif(nv.eq.2) then v=vdisp+vex endif endif C return end C C======================================================================= C======================================================================= double precision function > correction_ss(r,shift,c6,re) implicit none double precision r,shift,c6,re C correction_ss=0.0d0 if (r.lt.re) correction_ss=shift*(r-re)**2 C return end C======================================================================= C======================================================================= subroutine sort(a,b,n) implicit none integer n,i,j real*8 a(n),atem,b(n),btem do 10 j=1,n do 20 i=1,n-j if (a(i).gt.a(i+1)) then atem=a(i+1) a(i+1)=a(i) a(i)=atem btem=b(i+1) b(i+1)=b(i) b(i)=btem end if 20 continue 10 continue return end C======================================================================= C======================================================================= subroutine sort2(a,b,c,n) implicit none integer n,i,j real*8 a(n),atem,b(n),btem,c(n),ctem do 10 j=1,n do 20 i=1,n-j if (a(i).gt.a(i+1)) then atem=a(i+1) a(i+1)=a(i) a(i)=atem btem=b(i+1) b(i+1)=b(i) b(i)=btem ctem=c(i+1) c(i+1)=c(i) c(i)=ctem end if 20 continue 10 continue return end C======================================================================= C=======================================================================