!gfortran -O3 -fopenmp decaimento.f90 -o decaimento
PROGRAM DECAY
  IMPLICIT NONE
  REAL*8::x0,y0,xn,yn,a,b
  INTEGER::AllocateStatus,NPARAM,i,NF,j,k,NT,n
  REAL*8, PARAMETER:: pi=4.0*ATAN(1.0)
  INTEGER, PARAMETER:: NMAX=10000
  REAL*8, DIMENSION(:,:), ALLOCATABLE:: AB0
  !--------------------------------------------------------------------
  OPEN(1, FILE='NoTransport_P1_1E2.dat', STATUS='OLD')
  OPEN(2, FILE='Decay_P1.dat',STATUS='UNKNOWN')
  !--------------------------------------------------------------------
  NPARAM=0
  DO
     READ(1,*,END=10)
     NPARAM=NPARAM+1
  END DO
  10 CLOSE(1)
  !
  ALLOCATE(AB0(NPARAM,2),STAT = AllocateStatus)
  IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
  !
  OPEN(1, FILE='NoTransport_P1_1E2.dat', STATUS='OLD')
  DO i=1,NPARAM
     READ(1,*) AB0(i,1),AB0(i,2)
  END DO
  !--------------------------------------------------------------------
  NF=100
  FINAL_TIME: DO i=1,9
     NF=NF+100
     !
     OPEN(4, FILE='Param.dat')
     !$OMP PARALLEL DO DEFAULT(SHARED) &
     !$OMP PRIVATE(j,k,n,x0,y0,xn,yn,a,b,NT) &
     !$OMP SCHEDULE(dynamic)
     !
     PARAMETER: DO j=1,NPARAM
        a=AB0(j,1)
        b=AB0(j,2)        
        NT=0
        !
        IC: DO k=1,NMAX
           x0=k/(1.0*nmax)
           y0=-10.0
           !
           TIME: DO n=1,NF
              yn=y0-b*SIN(2.0*pi*x0)
              xn=MOD(x0+a*(1.0-yn**2.0),1.0)
              !
              IF(xn.LT.0.0) xn=xn+1.0
              !
              IF(yn.GT.10.0) THEN
                 NT=NT+1
                 EXIT IC
              END IF
              x0=xn;y0=yn
           END DO TIME
        END DO IC
        !
        !$OMP CRITICAL
77      IF(NT.EQ.0) WRITE(4,*) a,b
        !$OMP END CRITICAL
     END DO PARAMETER
     !$OMP END PARALLEL DO
     CLOSE(4)
     !---------------------------------------------------------
     NPARAM=0
     OPEN(4, FILE='Param.dat')
     DO
        READ(4,*,END=100) 
        NPARAM=NPARAM+1
     END DO
100   CLOSE(4)
     WRITE(2,*) NF,NPARAM
     !
     OPEN(4, FILE='Param.dat')
     DO j=1,NPARAM
        READ(4,*)  AB0(j,1),AB0(j,2)
     END DO
     CLOSE(4)
  END DO FINAL_TIME
END PROGRAM DECAY
