check.f

Go to the documentation of this file.

      SUBROUTINE check(DE)
*
*
*       Error check and restart.
*       ------------------------
*
      include 'common6.h'
*
*
*       See whether output intervals should be increased (at main output).
      IF (kz(32).GT.0.AND.time.GE.tnext) THEN
*       Check current energy level (factor of 2) for possible increase.
          k = kz(32)
          ecrit = 0.25/2.0**k
          IF (abs(e(3)).LT.ecrit) THEN
*       Define dynamical crossing time in case energy is near zero.
              tdyn = 2.0*rscale/sqrt(2.0*zkin/zmass)
              IF (2.0*dtadj.GT.tdyn.OR.time.LE.0.0d0) go to 5
              dtadj = 2.0*dtadj
              deltat = 2.0*deltat
*       Increase SMAX similarly and re-initialize hierarchical step array.
              smax = 2.0*smax
              CALL iblock
              qe = sqrt(2.0)*qe
              kz(32) = kz(32) + 1
              WRITE (6,1)  dtadj, deltat, qe, smax
    1         FORMAT (/,5x,'NEW INTERVALS:   DTADJ =',f6.2,
     &                     '  DELTAT =',f6.2,'  QE =',1p,e8.1,
     &                     '  SMAX =',0p,f7.3)
          END IF
      END IF
*
*       Perform automatic error check & restart (option 2).
    5 de = abs(de)
      etacor = 1.0
*
*       Skip error check after recent mass loss (otherwise reduce KZ(19)).
      IF (kz(19).EQ.2) THEN
          kz(19) = kz(19) - 1
          de = 0.0
          go to 70
      END IF
*
*       Check restart for large errors (two attempts permitted).
      IF (de.LT.5.0*qe) go to 30
*
*       Terminate run if no further restart is allowed.
      IF (kz(2).LE.1.OR.ndump.GE.2) THEN
          WRITE (6,10)
   10     FORMAT (/,9x,'CALCULATIONS HALTED * * *')
*       Increase NDUMP to prevent 3rd restart (safety check in routine MAIN).
          ndump = 2
          IF (kz(1).NE.0.AND.kz(2).GE.1) CALL mydump(1,1)
          stop
      END IF
*
*       Repeat the previous interval with reduced time-step parameters.
      tcomp = cpu
      ntemp = ndump
      CALL mydump(0,2)
      cpu = tcomp
      ndump = ntemp + 1
*       Control variable NDUMP used to prevent a third restart.
      etacor = 0.5
      etai = etacor*etai
      etar = etacor*etar
      IF (kz(17).GT.1) etau = etacor*etau
      dtmin = sqrt(etacor)*dtmin
      smin = sqrt(etacor)*smin
      WRITE (6,20)  time+toff, etai, etar, etau
   20 FORMAT (/,9x,'RESTART * * *   TIME =',f8.2,'  ETAI =',f7.3,
     &                           '  ETAR =',f7.3,'  ETAU =',f7.3)
      CALL mydump(1,2)
      go to 50
*
*       Reset counter and check optional modification of accuracy parameters.
   30 ndump = 0
      IF (kz(17).EQ.0) go to 50
*
      IF (de.GT.qe) THEN
*       Continue calculation but reduce the time-step parameters.
          etacor = sqrt(qe/de)
          etai = etacor*etai
          etar = etacor*etar
          IF (kz(17).GT.1) etau = etacor*etau
          dtmin = sqrt(etacor)*dtmin
          smin = sqrt(etacor)*smin
          IF (etacor.LT.0.99) WRITE (6,40)  etai, etar, etau
   40     FORMAT (8x,'ETAI =',f7.3,'  ETAR =',f7.3,'  ETAU =',f7.3)
      ELSE IF (de.LT.0.2*qe) THEN
*       Increase the time-step parameters (up to initial value only).
          IF (time.GT.0.0d0) THEN
              etacor = min(1.2d0,eta0/etai)
              etai = etacor*etai
              etar = etacor*etar
              IF (kz(17).GT.1) etau = etacor*etau
              dtmin = sqrt(etacor)*dtmin
              smin = sqrt(etacor)*smin
              IF (etacor.GT.1.01) WRITE (6,40)  etai, etar, etau
          END IF
      END IF
*
*       See whether the time-steps should be reduced (Note: KZ(2) > 2).
   50 IF (etacor.LT.1.0.AND.kz(2).GT.2) THEN
          etacor = sqrt(etacor)
          DO 60 i = ifirst,ntot
          IF (dmod(t0(i),0.5d0*step(i)).EQ.0.0d0) THEN
              step(i) = 0.5*step(i)
              tnew(i) = t0(i) + step(i)
          END IF
          IF (dmod(t0r(i),0.5d0*stepr(i)).EQ.0.0d0) THEN
              stepr(i) = 0.5*stepr(i)
          END IF
   60     CONTINUE
*
*       Set IPHASE = -1 to ensure new time-step list in INTGRT.
          iphase = -1
      END IF
*
   70 RETURN
*
      END