Skip to content
Snippets Groups Projects
Select Git revision
  • fitacs-lib
  • master default protected
  • multi-isotope-testing
  • program_flow_updates
  • SAMMY-8.2-beta2
  • SAMMY-8.2-beta1
6 results
Blame Permalink
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
mdat5.f 13.59 KiB 
C
C
C --------------------------------------------------------------
C
      SUBROUTINE Adjust_Auxil (Energb, Kdatb)
C
C *** PURPOSE -- add enough points to auxiliary grid Energb to ensure
c ***               that transition from one region to next will be smooth
C
      use ifwrit_m
      use AllocateFunctions_m
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
C
C
C      DIMENSION Energb(Ndatb)
C
      DIMENSION E(1000)
      real(kind=8),allocatable,dimension(:)::X21, X3, X2, X1
      real(kind=8),allocatable,dimension(:)::Energb
      DATA Jmax /1000/, Ratio /2.5d0/, Zero /0.0d0/
C
C
      call allocate_real_data(X21, 30)
      call allocate_real_data(X2 , 60)
      call allocate_real_data(X3 , Kdatb)
      call allocate_real_data(X1 , Kdatb)
      Itx = 0
   10 CONTINUE
      call reallocate_real_data(Energb, max(Kdatb,3))
      A = Energb(2) - Energb(1)
      X1(1) = Energb(1)
      C = X1(1)
      Ax = A
      Cx = C
      M = 1
      I = 1
      DO I=2,Kdatb-1
         B = Energb(I+1) - Energb(I)
C
         IF (A/B.GT.Ratio) THEN
C ***       Here distance between previous two points is much greater
C ***          than distance between current two points.  Ergo add
C ***          points between previous two points.
            CALL Dgrad (C, A, B, X2, Kk)   ! can add at most 30 points to x2
            A = X2(Kk) - C
            IF (Ax/A.GT.Ratio) THEN
C ***          In this case have to rethink M-th (previous) point
C ***             because it's too big relative to this new one
               call reallocate_real_data(X2, kk+30)
               CALL Dgrad (Cx, Ax, A, X2(Kk+1), Kkk) ! can add at most 30 points to x2
               Ii = Kkk
               call reallocate_real_data(X1, M+Kkk+1, 100)
               DO K=1,Kkk
                  X1(M) = X2(Kk+Ii)
                  M = M + 1
                  Ii = Ii - 1
               END DO
               X1(M) = C
            END IF
            Ii = Kk
            call reallocate_real_data(X1, M+Kk+1, 100)
            DO K=1,kk
               M = M + 1
               X1(M) = X2(Ii)
               Ii = Ii - 1
            END DO
            M = M + 1
            X1(M) = Energb(I)
            A = B
            C = X1(M)
            Cx = X1(M-1)
            Ax = C - Cx
         ELSE
C
            IF (B/A.GT.Ratio) THEN
C ***          Here distance between previous two points is much smaller
C ***              than distance between current two points.  Ergo add
C ***              points between current two points.
               M = M + 1
               call reallocate_real_data(X1, M+30)
               X1(m) = Energb(I)
               CALL Ugrad (Energb(I+1), B, A, X2, Kk) ! can add at most 30 points to x2
               call reallocate_real_data(X1, M+kk)
               DO K=1,Kk
                  M = M + 1
                  X1(M) = X2(K)
               END DO
               A = Energb(I+1) - X1(M)
               C = X1(M)
               Cx = X1(M-1)
               Ax = C - Cx
C
            ELSE
C
C ***          Here distance tween previous two points is comparable to
C ***               distance tween current two points.  Ergo accept the
C ***               center point.
               M = M + 1
               call reallocate_real_data(X1, M+100)
               X1(M) = Energb(I)
               Cx = C
               Ax = A
               A = B
               C = X1(M)
            END IF
         END IF
      END DO
C
C
      M = M + 1
      call reallocate_real_data(X1, M)
      X1(M) = Energb(Kdatb)
      Kdatb = M
      call reallocate_real_data(Energb, Kdatb)
      DO I=1,Kdatb
         Energb(I) = X1(I)
      END DO
      Itx = Itx + 1
      IF (Itx.LE.0) GO TO 10
C
C
C
C
      Ntries = 0
      Size = Zero
   70 CONTINUE
      Nn = Kdatb - 3
      Nnp1 = Kdatb - 2
      Nnp2 = Kdatb - 1
C
C *** define X21 = Energb(2) - Energb(1)
      call reallocate_real_data(X21, Nnp2)
      CALL Vsub (Energb(1), Energb(2), X21(1), Nnp2)
C
C *** define X2 = X21 ** 2
      call reallocate_real_data(X2,  Nnp2)
      CALL Vsq (X21(1), X2(1), Nnp2)
C
C *** define X1 = 1./X21
      call reallocate_real_data(X1, Nnp2)
      CALL Vrecip (X21(1), X1(1), Nnp2)
C
C *** define X3(I) = X2(I+2) - X2(I)
      call reallocate_real_data(X3, Nnp2)
      CALL Vsub (X2(1), X2(3), X3(1), NN)
C
C *** redefine X2(I) = X3(I) / X21(I+1) = X3(I) X1(I+1)
      call reallocate_real_data(X21, NN)
      call reallocate_real_data(X1, NN+3)
      CALL Vmul (X3(1), X1(2), X2(1), NN)
C
C *** redefine X1 = 5 * X21
      call reallocate_real_data(X1, Nnp2)
      CALL Vsmul (X21(1), 5.0D0, X1(1), Nnp2)
C
C *** now W(I) = X21(I-2) + X1(I-1) + X1(I) + X21(I+1)
C ***           + X2(I-2) - X2(I-1)
C
C *** here X3 is reused to be weighting factor
      X3(1)    =                        X1(1)    + X21(2)
      X3(2)    =             X1(1)    + X1(2)    + X21(3)
     *                                - X2(1)
      call reallocate_real_data(X3, max(Nn+2, Kdatb))
      DO K=3,Nn+1
         X3(K) = X21(K-2)  + X1(K-1)  + X1(K)    + X21(K+1)
     *                     + X2(K-2)  - X2(K-1)
      END DO
      X3(Nn+2) = X21(Nn)   + X1(Nn+1) + X1(Nn+2)
     *                     + X2(Nn)
      X3(Kdatb)= X21(Nn+1) + X1(Nn+2)
C
C
C
      IF (Ntries.EQ.0) THEN
C ***    First time in, figure "size" which is average size of positive
C ***       weights
         A = Zero
         K = 0
         DO I=2,Nnp2
            IF (X3(I).GT.Zero) K = K + 1
            IF (X3(I).GT.Zero) A = A + X3(I)
         END DO
         Size = -A/dFLOAT(K)
      END IF
C
      A = Zero
      K = 0
      Negati = 0
      Ix = 0
C *** Count how many negative weights there are
      DO I=2,Nnp2
         IF (X3(I).LE.Zero) Ix = Ix + 1
         IF (X3(I).LE.Size) Negati = 1
         IF (X3(I).GT.Zero) A = A + X3(I)
         IF (X3(I).GT.Zero) K = K + 1
      END DO
      IF (Kdebug.NE.0 .AND. Ix.GT.0) WRITE (6,20200) Size, Ix
20200 FORMAT (' Size =', 1PG14.6, ' number of negative weights=', I10)
      IF (Kdebug.NE.0 .AND. Negati.EQ.0 .AND. Ntries.GT.0)
     *   WRITE (6,20000) Ntries
20000 FORMAT (' Number of iterations to get rid of negatives in grid =',
     *   I5)
      IF (Negati.EQ.0) THEN
         deallocate(X21)
         deallocate(X3)
         deallocate(X2)
         deallocate(X1)
         RETURN
      END IF
      Ntries = Ntries + 1
      IF (Ntries.GT.10) THEN
         deallocate(X21)
         deallocate(X3)
         deallocate(X2)
         deallocate(X1)
         RETURN
      END IF
C
      B = A/dFLOAT(K)
      B = -B
C
C
C *** Here there are negative weights so we have to decide how to
C ***    adjust the grid to eliminate negatives
      call reallocate_real_data(X21, Nnp1+1)
      call reallocate_real_data(X2,  max(nnp1, nnp2))
      call reallocate_real_data(X2,  Kdatb)
      call reallocate_real_data(X1, Nnp1)
      call reallocate_real_data(Energb, max(nnp1, nnp2))
      call reallocate_real_data(Energb, Kdatb)
      DO I=1,Nnp1
         X1(I) = X21(I)/X21(I+1)
         IF (X1(I).LT.Zero) THEN
            WRITE (6,30000) I
30000       FORMAT (' For I=', I6)
            Mmin = I - 100
            IF (Mmin.LE.0) Mmin = 1
            Mmax = I + 10
            IF (Mmax.GT.Kdatb) Mmax = Kdatb
            WRITE (6,30100) (M,Energb(M),M=Mmin,Mmax)
30100       FORMAT (('  E=', 4(I5,1PG12.5)))
            WRITE (6,30200) (M,X21(M),M=Mmin,Mmax)
30200       FORMAT ((' dE=', 4(I5,1PG12.5)))
         END IF
         X2(I) = Energb(I)
      END DO
      X2(Nnp2) = Energb(Nnp2)
      X2(Kdatb) = Energb(Kdatb)
C
      Ii = 1
      DO 220 I=1,Kdatb
         call reallocate_real_data(Energb, Ii)
         Energb(Ii) = X2(I)     
         Ii = Ii + 1
         IF (I.LE.4 .OR. I.GE.NN) GO TO 220
C
         IF (X1(I).GE.0.5) THEN
            IF (X1(I).LE.2.0) GO TO 220
C
C ***       Here (e(I+1)-e(I))/(e(I+2)-e(I+1)) > 2.0
C ***            X1(I) = [ X21(I) / X21(I+1) ] > 2.0
            DO J=I-4,I+4
               IF (X3(J).LE.Size) GO TO 130
            END DO
            GO TO 220
  130       CONTINUE
            IF (X1(I+1).GT.6.0D0 .AND. X1(I).GT.4.0D0)
     *         A = X21(I+1)*2.0D0
            IF (X1(I+1).LE.6.0D0 .OR. X1(I).LE.4.0D0)
     *         A = X21(I+1)
            J = 1
            E(J) = X2(I+1) - A
C
            IF (X1(I+1).GT.6.0D0) THEN
C ***          Here we Make further divisions of region i when the next
C ***             region is very small.  Note that array E is ordered in
C ***             reverse.
  140          CONTINUE 
               J = J + 1
               IF (J.GT.Jmax) THEN
                  WRITE (6,10000) J, Jmax
10000             FORMAT (' array size should be', I7,' but is only',I7)
                  WRITE (6,10100) (E(M),M=1,10)
10100             FORMAT (' E( 1...10)=', 1p5g12.4)
                  WRITE (6,10200) (E(M),M=991,1000)
10200             FORMAT (' E(991...1000)=', 1P5G12.4)
                  STOP  '[#2 in Adjust_Auxil in dat/mdat5]'
               END IF
               A = A * 2.0D0
               E(J) = E(J) - A
               IF (Energb(Ii-1).LT.E(J)) GO TO 140
C
C ***          Now we're done adding those extra points
               J = J - 1
               IF (J.GT.Jmax) THEN
                  WRITE (6,10000) J, Jmax
                  STOP '[#3 in Adjust_Auxil in mdat5]'
               END IF
               IF ((E(J)-Energb(II-1)) .LT. A/2.0D0) J = J - 1
            END IF
C
  150       CONTINUE
C ***       Update Energb array
            Jjj = J
            DO 170 Jj=1,J            
               IF (Jjj.GT.Jmax) WRITE (6,10000) Jjj, Jmax              
               IF (Jjj.GT.Jmax) STOP '[#5 in Adjust_Auxil in mdat5]'
               call reallocate_real_data(Energb, Ii, 100)
               Energb(Ii) = E(Jjj)
               Ii = Ii + 1
               Jjj = Jjj - 1
  170       CONTINUE
            GO TO 220
C
C
         ELSE
C
C ***       Here (e(I+1)-e(I))/(e(I+2)-e(I+1)) < 1./2.
C ***             X1(I) = [ X21(I) / X21(I+1) ] < 0.5
            DO J=I-4,I+4
               IF (X3(J).LE.Size) GO TO 190
            END DO
            GO TO 220
  190       CONTINUE
C
C ***       Determine the size of the adjustment to be Made
            IF (X1(I).LT.0.1667D0) A = X21(I)*2.0D0
            IF (X1(I).GE.0.1667D0) A = X21(I)
            call reallocate_real_data(Energb, Ii, 100)
            Energb(Ii) = X2(I+1)           
            Ii = Ii + 1
C
            IF (X1(I).LT.0.1667D0) THEN
C ***          Here we need to subdivide region i several times
  200          CONTINUE
               call reallocate_real_data(Energb, Ii, 100)
               Energb(Ii) = Energb(Ii-1) + A             
               Ii = Ii + 1
               A = A*2.0D0
               IF (Energb(Ii-1)+A .LT. X2(I+2)) GO TO 200
C
C ***          and we're done, so update Energb
               IF (X2(I+2)-Energb(Ii-1) .LT. A/2.0D0) Ii = Ii - 1
            END IF
C
            call reallocate_real_data(X2,  I+2, 100)
            X2(I+1) = Energb(Ii-1)
            Ii = Ii - 1
            IF (I+1.GE.Nnp1) GO TO 220            
            X21(I+1) = X2(I+2) - X2(I+1)
            call reallocate_real_data(X1, I+1, 100)
            X1(I+1) = X21(I+1)/X21(I+2)
C
         END IF
C
  220 CONTINUE
      Kdatb = Ii - 1      
      Size = B
C
C
C *** Finally, double-check to be sure we did not introduce spurious
C ***      points (as in two at same energy)
      N = Kdatb - 1
      Isame = 0
      Imix = 0
      DO J=1,N
         call reallocate_real_data(Energb, J+1, 100)
         IF (Energb(J+1).EQ.Energb(J))THEN
C ***    Here two points have same energy so drop one
            Kdatb = Kdatb - 1
            call reallocate_real_data(Energb, Kdatb + 1, 100)
            DO K=J+1,Kdatb
               Energb(K) = Energb(K+1)
            END DO
            Isame = 1 + Isame
         ELSE IF (Energb(J+1).LT.Energb(J)) THEN
C ***       Here Energb(J+1) .LT. Energb(J) which can't happen
            Imix = Imix + 1
            WRITE (6,89999) J, Energb(J), Energb(J+1)
89999       FORMAT(' Energb(J=',I6,',J+1)=', 1P2G14.6)
         END IF
      END DO
      IF (Imix.GT.0) WRITE (6,99999) Imix, Isame
99999 FORMAT (' Problems in Energy-adjustment algorithm: Imix, Isame =',
     *     I5, I5)
      IF (Imix.GT.0) STOP '[#8 in Adjust_Auxil in mdat5]'
      IF (Isame.NE.0) WRITE (6,99998) Isame
99998 FORMAT (' Two identical energies were chosen', I3, ' times.')
      GO TO 10
      END
C
C
C ______________________________________________________________________
C
      SUBROUTINE Dgrad (E1, Diff, D, Store, Kk)
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      DIMENSION Store(*), X(30)
      K = 0
      A = Diff
   10 CONTINUE
         A = A/2.0
         IF (A.LT.D) GO TO 20
         K = K + 1
         IF (K.GT.30) STOP '[STOP in Dgrad in dat/mdat5]'
         X(K) = A
      GO TO 10
C
   20 CONTINUE
      Kk = K
      DO J=2,kk
         X(J) = X(J-1) + X(J)
         K = K - 1
      END DO
      IF (x(Kk).GT.Diff-d) Kk = Kk - 1
      K = Kk
      DO I=1,kk
         Store(I) = E1 + X(K)
         K = K - 1
      END DO
      RETURN
      END
C
C
C ______________________________________________________________________
C
      SUBROUTINE Ugrad (E1, Diff, D, Store, Kk)
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      DIMENSION Store(*), X(30)
      K = 0
      A = Diff
   10 CONTINUE
         A = A/2.0d0
         IF (A.LT.D) GO TO 20
         K = K + 1
         IF (K.GT.30) STOP '[STOP in Ugrad in dat/mdat5]'
         X(K) = A
      GO TO 10
C
   20 CONTINUE
      IF (K.GT.30) THEN
         WRITE (6,10100) E1, Diff, D
10100    FORMAT(' E1, Diff, D=', 1P5G20.10)
         WRITE (6,10000) (X(I),I=1,K)
10000    FORMAT (' X=', 1p5G14.6)
         STOP '[STOP #1 in Ugrad in mdat5]'
      END IF
      Kk = K
      DO K=2,Kk
         X(K) = X(K) + X(K-1)
      END DO
      IF (X(Kk).GT.Diff-D) Kk = Kk - 1
      Ii = Kk
      DO I=1,Kk
         Store(Ii) = E1 - X(I)
         Ii = Ii - 1
      END DO
      RETURN
      END