c-------------------------------------------------------------------- c Fit halo profiles, c use simple minimization c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (Nmax = 500000) ! PARAMETER (Nbinmax = 120) ! PARAMETER (Nqual = 13) ! COMMON /HALOS/hals(Nqual,Nmax),Nhbin(Nmax), & Prof(7,Nbinmax,Nmax),Nsmall(Nbinmax,Nmax) COMMON /Binn/ Radout(Nbinmax) COMMON /APPROX/rsh(Nmax),Vmaxh(Nmax) DIMENSION Nin(Nbinmax),Radden(Nbinmax),Vc(Nbinmax),Dens(Nbinmax) DIMENSION VVmax(Nmax),Rs(Nmax),r1(Nmax),rho1(Nmax),x2c(Nmax), & a2c(Nmax),error(Nmax),slope(Nmax) Character Header*100 Open(1,file='Catalog.unbound.DAT') Open(2,file='CatFit.unbound.dat') Open(3,file='CatOut.unbound.dat') c Open(1,file='CatalogB.486.DAT') c Open(2,file='CatFitB.486.dat') c Open(3,file='CatOut.486.dat') c Open(1,file='CatalogA.9999.DAT') c Open(2,file='CatFitA.9999.dat') c Open(3,file='CatOut.9999.dat') Lines =14 ! number of lines in the header preceeding radial bins lbins =100 ! number of bins along radius Lskip =5 ! number of header lines between bins and catalog aMfit = 1.e12 ! limit for mass for fitting NFW or Sersic profiles Do i=1,Lines read(1,'(a)') Header write (*,'(a)') Header write (2,'(a)') Header EndDo read(1,*) (Radout(i),i=1,lbins) write(2,'(10f8.2)') (Radout(i),i=1,lbins) write (*,'(10f8.2)') (Radout(i),i=1,lbins) Do i=1,Lskip read(1,'(a)') Header if(i.lt.4)write (2,'(a)') Header if(i.eq.4)write (2,'(a,a)') Header, & ' C Nbin error Slope nLarge' write (*,'(a)') Header EndDo Nh =0 100 Nh =Nh+1 read(1,*,err=400,end=400) (hals(i,Nh),i=1,Nqual),Nhbin(Nh) c write(*,'(" Halo=",i8," coords=",3f8.3," bins=",i3)') c & Nh,(hals(i,Nh),i=1,3),Nhbin(Nh) do i=1,Nhbin(Nh) read(1,*) (Prof(j,i,Nh),j=1,7),Nsmall(i,Nh) EndDo If(Nh.lt.Nmax)goto100 c If(Nh.lt.1000)goto100 Nh =Nh +1 400 Nh =Nh -1 write (*,*) ' Number of halos read =',Nh factor = 8.23e10 Nchunk = max(Nh/20,100) C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (ih,Umax,Rumax,Rvir,Nbin,Radden) C$OMP+PRIVATE (Nin,Vc,Dens) C$OMP+SCHEDULE(DYNAMIC,Nchunk) Do ih=1,Nh If(mod(ih,Nchunk/2).eq.0)write(*,*) ' halo=',ih Umax = hals(10,ih) Rumax= hals(12,ih) Rvir = hals(8,ih) Nbin = Nhbin(ih) do i=1,Nbin Radden(i) = Prof(1,i,ih) Nin(i) = Prof(2,i,ih) Vc(i) = Prof(6,i,ih) Dens(i) = Prof(7,i,ih)/factor EndDo if(hals(7,ih).lt.aMfit)Then CALL HaloFit(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, & Rs(ih),VVmax(ih),Nbin,error(ih)) r1(ih) = 0. rho1(ih) =0. x2c(ih) =0. a2c(ih) =0. slope(ih) =1. Else CALL HaloFitSER(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, & Rs(ih),slope(ih),VVmax(ih),Nbin,error(ih)) c CALL HaloFitExp(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, c & Rs(ih),VVmax(ih),r1(ih),rho1(ih),x2c(ih), c & a2c(ih),Nbin,error(ih)) EndIf EndDo Do ih=1,Nh Rvir = hals(8,ih) Conc = Rvir/Rs(ih) error(ih) =min(error(ih),1.d0) Do i=1,Nhbin(ih) If(Prof(1,i,ih).ge.Rvir)Then nn =INT(Prof(2,i,ih)) nL =INT(Prof(2,i,ih))-Nsmall(i,ih) goto 18 EndIf EndDo 18 write(2,20) (hals(i,ih),i=1,9), & VVmax(ih),INT(hals(11,ih)),Rs(ih)*2.15,Conc,Nhbin(ih), & sqrt(error(ih)),1./slope(ih),nL c & ,r1(ih),rho1(ih),x2c(ih),a2c(ih) 20 format(3f8.3,3F8.1,g11.3,f8.2,2f8.1,i9,f8.1,f8.2,i4,f8.3 & f7.3,i8) c If((Conc.lt.4..or.Conc.gt.25.).and.hals(7,ih).gt.1.e12)Then c write(3,'(" Halo=",i8," coords=",3f8.3," M=",g11.4, c & " C=",f8.3," R1,Rs,Rvir=",3f8.3)') c & ih,(hals(i,ih),i=1,3),hals(7,ih), c & Conc,r1(ih),Rs(ih),hals(8,ih) c write (3,'(5x,"r Nin Dens Dens_4exp", c & 6x,"Vc Vc_4exp")') c do i=1,Nhbin(ih) c Radd = Prof(1,i,ih) c Ni = Prof(2,i,ih) c Vcc = Prof(6,i,ih) c Den = Prof(7,i,ih)/factor c c write(3,50) Radd,Ni,Den, c & RhoExp(Radd,r1(ih),rho1(ih),x2c(ih),a2c(ih)), c & Vcc,VcExp(Radd,r1(ih),rho1(ih),x2c(ih),a2c(ih)) c EndDo c EndIf enddo 50 format(f8.3,i7,2g11.4,2f9.2) stop end c-------------------------------------------------------------------- c Fit halo profile SUBROUTINE HaloFitSer(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, & Rs,slope,VVmax,Nbin,error) c output: Rs = estimate of NFW core radius, c Vvmax = estimate of Vmax c slope = 1/n slope for Sersic fit c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (Nbinmax = 120) ! PARAMETER (facV = 2.0) ! max variation for Vmax PARAMETER (facR = 3.0) ! max variation for rs c PARAMETER (facRvir = 2.0) ! max radius for fit in units of Rvir PARAMETER (facRvir = 0.5) ! max radius for fit in units of Rvir PARAMETER (nStepS = 50) ! number of bins in slope PARAMETER (nStepR = 50) ! number of bins in radius PARAMETER (nStepV = 50) ! number of bins in Vmax DIMENSION Radden(Nbin),Nin(Nbin),Dens(Nbin),Vc(Nbin), & Radout(Nbinmax) DATA iFlag/0/ SAVE iFlag rs = Rumax/2.15 ! initial guess for rs VVmax = Umax Vcmax = Umax Vmin = Umax/facV dV = Umax/nStepV*(facV -1./facV) ! step in Vmax Rmin = rs/facR dR = rs/nStepR*(facR -1./facR) ! step in rs Smin = 0.08 Smax = 0.350 dS = (Smax-Smin)/nStepS ! step in slope iFlag = iFlag +1 error = 1.e10 istart= 1 Do istart =1,Nbin If(Nin(istart).gt.200)goto20 EndDo write (*,*)'not enough particles(first,last=)',Nin(1),Nin(Nbin) return ! not enough particles; get back 20 Rmax = min(Rvir/facRvir,Radout(Nbin)) ! max radius for the fit do i=Nbin,1,-1 if(Radout(i).lt.Rmax)goto30 EndDo 30 imax = i if(imax.le.istart)Then write (*,*)'not enough particles. Bins=',istart,imax return EndIf Do is =0,nStepS slope = Smin +is*dS Do iv =0,nStepV Vmaxx = Vmin +iv*dV Do ir =0,nStepR rss = Rmin +ir*dR rho_0 = 1.022e3*(Vmaxx/rss)**2 ee =0. ww =0. do i=istart,imax w = (Radout(i)-Radout(i-1))/Radden(i) ee=ee+w*(log10(RhoSER(Radden(i),rho_0,rss,slope)/Dens(i)))**2 ww = ww + w EndDo c ee =ee/(imax-istart) ee =ee/ww If(ee.lt.error)then error = ee Rs = rss Slp = slope VVmax = Vmaxx ivv = iv irr = ir rho = rho_0 c write (*,'(10x,"V=",f8.3," r/slope=",2f8.3," err=",3g11.3)') c & Vmaxx,rss,Slp,ee,rho_0 EndIf EndDo EndDo EndDo write (*,'(10x,"V=",f8.3," r/Slp=",2f8.3," err=",7g11.3)') & VVmax,Rs,1./Slp,error,rho,Rvir,Rvir/Rs c write(*,*) ' after ----' c do i=1,Nbin c x = Radden(i)/rs c F = log(1.+x)-x/(1.+x) c write(*,50) i,Radden(i),Dens(i),rho/x/(1.+x)**2, c & Vc(i),VVmax*sqrt(4.62*F/x) c enddo slope = Slp 50 format(i3,' r=',f8.3,' dens=',2g11.4,' v=',3f8.2) rho_0 = 1.022e3*(VVmax/Rs)**2 Do i=1,Nbin write (3,'(4g11.4,i9)')Radout(i),Dens(i),Vc(i), & RhoSER(Radden(i),rho_0,Rs,slope),Nin(i) EndDo c if(iFlag.ge.10)stop Return End c-------------------------------------------------------------------- c Sersic density profile FUNCTION RhoSER(r,rho_0,rs,slope) c c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) x = r/rs RhoSER = rho_0/exp(2./slope*(x**slope-1.))+1. Return End c-------------------------------------------------------------------- c Fit halo profile SUBROUTINE HaloFit(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, & Rs,VVmax,Nbin,error) c output: Rs = estimate of NFW core radius, c Vvmax = estimate of Vmax c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (Nbinmax = 120) ! PARAMETER (facV = 1.25) ! max variation for Vmax PARAMETER (facR = 3.0) ! max variation for rs c PARAMETER (facRvir = 2.0) ! max radius for fit in units of Rvir PARAMETER (facRvir = 1.0) ! max radius for fit in units of Rvir PARAMETER (nStepR = 50) ! number of bins in radius PARAMETER (nStepV = 20) ! number of bins in Vmax DIMENSION Radden(Nbin),Nin(Nbin),Dens(Nbin),Vc(Nbin), & Radout(Nbinmax) rs = Rumax/2.15 ! initial guess for rs VVmax = Umax Vcmax = Umax Vmin = Umax/facV dV = Umax/nStepV*(facV -1./facV) ! step in Vmax Rmin = rs/facR dR = rs/nStepR*(facR -1./facR) ! step in rs error = 1.e10 istart= 1 Do istart =1,Nbin If(Nin(istart).gt.200)goto20 EndDo return ! not enough particles; get back 20 Rmax = min(Rvir/facRvir,Radout(Nbin)) ! max radius for the fit do i=Nbin,1,-1 if(Radout(i).lt.Rmax)goto30 EndDo 30 imax = i if(imax.le.istart)Then c write (*,*)'not enough particles. Bins=',istart,imax return EndIf Do iv =0,nStepV Vmaxx = Vmin +iv*dV Do ir =0,nStepR rss = Rmin +ir*dR rho_0 = 1.022e3*(Vmaxx/rss)**2 ee =0. ww =0. do i=istart,imax w = (Radout(i)-Radout(i-1))/Radden(i) ee = ee + w*(log10(RhoNFW(Radden(i),rho_0,rss)/Dens(i)))**2 ww = ww + w EndDo c ee =ee/(imax-istart) ee =ee/ww If(ee.lt.error)then error = ee Rs = rss VVmax = Vmaxx ivv = iv irr = ir rho = rho_0 c write (*,'(10x,"V=",f8.3," r=",f8.3," err=",3g11.3)') c & Vmaxx,rss,ee,rho_0 EndIf EndDo EndDo c write (*,'(10x,"V=",f8.3," r=",f8.3," err=",3g11.3)') c & VVmax,Rs,error,rho c write(*,*) ' after ----' c do i=1,Nbin c x = Radden(i)/rs c F = log(1.+x)-x/(1.+x) c write(*,50) i,Radden(i),Dens(i),rho/x/(1.+x)**2, c & Vc(i),VVmax*sqrt(4.62*F/x) c enddo 50 format(i3,' r=',f8.3,' dens=',2g11.4,' v=',3f8.2) Return End c-------------------------------------------------------------------- c NFW density profile FUNCTION RhoNFW(r,rho_0,rs) c output: Rs = estimate of NFW core radius, c Vvmax = estimate of Vmax c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) x = r/rs RhoNFW = rho_0/x/(1.+x)**2 Return End c-------------------------------------------------------------------- c Double Exp-Fit for halo profile SUBROUTINE HaloFitExp(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, & Rs,VVmax,r1,rho1,x2c,a2c,Nbin,error) c output: Rs = estimate of NFW core radius, c Vvmax = estimate of Vmax c Uses 4 parameters: c rho_1, r_1, a_2, x_2 c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (Nbinmax = 120) ! PARAMETER (facRho = 10.) ! max variation for density PARAMETER (facR = 10.0) ! max variation for r1 c PARAMETER (facRvir = 2.0) ! max radius for fit in units of Rvir PARAMETER (facRvir = 0.5) ! max radius for fit in units of Rvir DIMENSION Radden(Nbin),Nin(Nbin),Dens(Nbin),Vc(Nbin), & Radout(Nbinmax) Logical Converge rs = Rumax/2.15 ! initial guess for rs VVmax = Umax r1 = min(max(Rumax/2.15,Radden(2)),Rvir/5.) ! initial guess for r1 dR = r1/20. c rho1 = 1.022e3*(Umax/r1)**2 ! initial guess for rho_1 rho1 = 200./(r1/Rvir)**2 dRho = rho1/10. x2c = 1.5 dx2 = x2c/10. error = 100. c write (*,*) c write(*,'(" Guess: R,rho,a2,x2=",4g11.3)')r1,rho1,a2c,x2c c write(*,'(" steps: R,rho,a2,x2=",4g11.3)')dR,dRho,da2,dx2 istart= 1 Do istart =1,Nbin If(Nin(istart).gt.100)goto20 EndDo return ! not enough particles; get back 20 Rmax = min(Rvir/facRvir,Radout(Nbin)) ! max radius for the fit do i=Nbin,1,-1 if(Radout(i).lt.Rmax)goto30 EndDo 30 imax = i c write (*,*) ' max bin=',imax,' Rout=',Radout(imax),Rmax if(imax.le.istart.or.Nin(imax).lt.1000)Then c write (*,*)'not enough particles. Bins=',istart,imax return EndIf Do iter =1,100 r1_old = r1 rho1_old = rho1 x2c_old = x2c e_old = error c write (*,*) ' ==== iteration=== ',iter,error CALL LocalFit(Radden,Radout,Dens,istart,imax,Nbin, & r1,rho1,x2c,dR,dRho,dx2,error,iFlag) c write (*,'(10x,"r1=",f8.4," rho1=",7g11.4," err=",g11.4,i2)') c & r1,rho1,a2c,x2c,dR,dRho,dx2,da2,error,iFlag quality = (r1/r1_old-1.)**2+(rho1/rho1_old-1.)**2+ & (x2c/x2c_old-1)**2 Converge = (quality.lt.1.e-4).and. & (abs(error/e_old-1.).lt.1.e-3).and. & (iter.gt.10).and.(iFlag.eq.0) dRho =max(min(dRho,rho1/10.),rho1/30.) If(Converge)Goto40 If(error.ge.e_old.and.iter.gt.10)Then r1 = r1_old rho1= rho1_old x2c = x2c_old dR =dR/1.21 dRho =dRho/1.21 dx2 =dx2/1.21 c write (*,*) ' restore center, reduce steps' EndIf If(iFlag.eq.0)Then dR =min(dR/1.41,r1/20.) dRho =min(dRho/1.41,rho1/20.) dx2 =dx2/1.41 c write (*,*) ' changed steps' c Else c dR =dr/1.01 c dRho =drho/1.01 EndIf EndDo 40 Rs =r_1 c write (*,'(10x,"r1=",f8.3," rho1=",7g11.4," err=",g11.3,2i4)') c & r1,rho1,a2c,x2c,dR,dRho,dx2,da2,error,iFlag,iter CALL FindVcMax(rmax,Vcmax,Rvir,r_s,r1,rho1,x2c,0.05) Rout = Radden(imax) c r_s =min(r_s,r_s2) c r_s =r_s2 c If(Rvir/r_s.lt.4.)Then c write (*,'(" rmax,Vmax=",2f8.3," rs=",3g11.3)')rmax,Vcmax, c & r_s,Rout,Dens(imax) c do i=istart,imax c write(*,50) i,Radden(i),Dens(i), c & RhoExp(Radden(i),r1,rho1,x2c,a2c), c & rhs*Rout/Radden(i)*(1.+Rout/Radden(i))**2, c & Vc(i),VcExp(Radden(i),r1,rho1,x2c,a2c) c enddo c EndIf c 50 format(i3,' r= ',f8.3,' dens=',3g11.4,' v=',3f8.2) Rs = r_s VVmax = Vcmax Return End c-------------------------------------------------------------------- c Fit halo profile SUBROUTINE LocalFit(Radden,Radout,Dens,istart,imax,Nbin, & Rc,Rhoc,x2c,dR,dRho,dx2,error,iFlag) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (Nbinmax = 120) ! PARAMETER (nStepR = 6) ! number of bins in radius PARAMETER (nStepRho= 6) ! number of bins in Rho PARAMETER (nStepx2 = 6) ! number of bins in x2 DIMENSION Radden(Nbin),Dens(Nbin),Radout(Nbinmax) Logical on_boundary Rmin = max(Rc -dR*nStepR/2.,Rc/10.) Rhomin = max(Rhoc -dRho*nStepRho/2.,Rhoc/10.) x2min = max(x2c - dx2*nStepx2/2.,1.d-2) c write (*,'(5x," Rmin/max =",4g11.4)') c & Rmin,Rmin +nStepR*dR,Rc c write (*,'(5x," Rhomin/max=",4g11.4)') c & Rhomin,Rhomin +nStepRho*dRho,Rhoc c write (*,'(5x," xmin/max =",4g11.4)') c & x2min,x2min +nStepx2*dx2,x2c error =1.e+10 Do ir =0,nStepR r11 = Rmin +ir*dR e_old = error Do irho =0,nStepRho rho11 = Rhomin +irho*dRho Do ix2 = 0,nStepx2 x2 = x2min +ix2*dx2 c write (*,'(4g11.3,4i4)') r11,rho11,x2,a2,ir,irho,ix2,ia2 ee =0. ww =0. do i=istart,imax w = (Radout(i)-Radout(i-1))/Radden(i) ee = ee + w*(log10(RhoExp( Radden(i),r11,rho11,x2)/ & Dens(i)))**2 ww = ww + w EndDo c ee =ee/(imax-istart) ee =ee/ww If(ee.lt.error)then error = ee r_1 = r11 rho_1 = rho11 x_2 = x2 ir_1 = ir irho_1= irho ix_1 = ix2 EndIf EndDo EndDo c write (*,'(10x,"r1=",f8.3," rho1=",3g11.3," err=",g11.3,4i4)') c & r_1,rho_1,a_2,x_2,error,ir_1,irho_1,ix_1,ia_1 c If(abs(error-e_old).lt.1.e-5)Goto100 EndDo 100 on_boundary = .false. if(ir_1.eq.0.or.ir_1.eq.nStepR)on_boundary = .true. if(irho_1.eq.0.or.irho_1.eq.nStepRho)on_boundary = .true. if(ix_1.eq.0.or.ix_1.eq.nStepx2)on_boundary = .true. iFlag = 0 If(on_boundary)iFlag =1 Rc = r_1 Rhoc= rho_1 x2c = x_2 c write (*,'(10x,"r1=",f8.3," rho1=",3g11.3," err=",g11.3,4i4)') c & r_1,rho_1,a_2,x_2,error,ir_1,irho_1,ix_1,ia_1 Return End c-------------------------------------------------------------------- c LogNorm-fit for density profile FUNCTION RhoExp(r,r1,rho_1,x2) c c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) x = r/r1 RhoExp = rho_1/x**2*exp(-0.5*(x/x2)**2)+1. Return End c-------------------------------------------------------------------- c dimensionless integrant for mass profile FUNCTION AxExp(x,x2) c c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) AxExp = exp(-0.5*(LOG(x)/x2)**2) Return End c-------------------------------------------------------------------- c Exp-fit mass profile: dimensionless c FUNCTION FExp(x,x2) c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) REAL*8 INTG EXTERNAL AxEXP xmin = min(1.d-4,x/10.d0) FExp = INTG(AxEXP,xmin,x,x2) Return End c-------------------------------------------------------------------- c Exp-fit velocity profile FUNCTION VcExp(r,r1,rho_1,x2) c output: Vcirc in km/s c c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) x = r/r1 VcExp = 6.698e-2*r1*sqrt(rho_1/x*FExp(x,x2)) Return End c-------------------------------------------------------------------- c Vmax and effective r_s for Exp-fit velocity profile SUBROUTINE FindVcMax(r,Vcmax,Rvir,r_s,r_1,rho_1,x2,dr) c output: Vcirc in km/s c c-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Vmax =0. iter =0 r = r_1 Vv0 = VcExp(r,r_1,rho_1,x2) Vv1 = VcExp(r+dr,r_1,rho_1,x2) If(Vv1.gt.Vv0)Then r = r-dr Else r = 0 EndIf 10 iter =iter +1 ! Find Vcmax using the approximation r = r +dr Vv= VcExp(r,r_1,rho_1,x2) c write (*,*) ' VV,r=',Vv,r,Vmax If(Vv.gt.Vmax.and.r.lt.1.e4.and.r.lt.Rvir)Then Vmax = Vv GoTo 10 EndIf Vcmax = Vmax rmax = r Return End C-------------------------------------------- Simpson integration REAL*8 FUNCTION INTG(FUNC,A,B,sig) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (EPS=5.0d-4, JMAX=22) EXTERNAL FUNC OST=-1.D30 OS= -1.D30 ST =0. DO 11 J=1,JMAX CALL TRAPZD(FUNC,A,B,ST,J,sig) INTG=(4.0d0*ST-OST)/3.0d0 IF (ABS(INTG-OS).Le.EPS*ABS(OS)) RETURN OS=INTG OST=ST 11 CONTINUE WRITE (*,*)'Integration did not converge' RETURN END C---------------------------------------------- SUBROUTINE TRAPZD(FUNCC,A,B,S,N,sig) C--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) SAVE IT EXTERNAL FUNCC IF (N.EQ.1) THEN S=0.5d0*(B-A)*(FUNCC(A,sig)+FUNCC(B,sig)) IT=1 ELSE TNM=IT DEL=(B-A)/TNM X=A+0.5D0*DEL SUM=0.0D0 DO 11 J=1,IT SUM=SUM+FUNCC(X,sig) X=X+DEL 11 CONTINUE S=0.5D0*(S+(B-A)*SUM/TNM) IT=2*IT ENDIF RETURN END