text:
      SUBROUTINE CALCPR(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT,
     1                  PROB, IUSER, RUSER, MNPCDF)
C
C *** THIS SUBROUTINE CALCULATES A PROBABILITY FOR THE MODEL AND   ***
C *** DATA GIVEN.  FOR MULTINOMIAL PROBIT SOME ADDITIONAL STORAGE  ***
C *** CUSTOMIZATION IS REQUIRED.  THIS APPROACH CAN BE             ***
C *** USED FOR OTHER CHOICE MODELS, WITH APPROPRIATE MODIFICATIONS ***
C *** TO THE ARRAYS USED BELOW.                                    ***
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
      INTEGER NPAR, IERR, ICH, IALT, II, ICDAT(*), IR, IUSER(*)
      REAL X(NPAR), RCDAT(*), PROB, RUSER(*)
      EXTERNAL MNPCDF
C
C *** CALCPR PARAMETER USAGE ***
C
C IALT.... NUMBER OF CHOICES AVAILABLE IN THE CHOICE SET.
C ICDAT... VECTOR OF INTEGER DATA VALUES.
C ICH..... INTEGER INDICATING THE CHOICE.  1 <= ICH <= IALT.
C IERR.... INTEGER FOR PASSING ERROR INFORMATION.
C             IN THIS ROUTINE, IF IERR = 1 ON RETURN THEN THERE WERE
C               NO PROBLEMS.
C             IF IERR = 0 ON RETURN, THEN THE PROBABILITY COULD NOT
C               BE COMPUTED USING THE CURRENT PARAMETERS IN X.
C II...... NUMBER OF INTEGER VALUES STORED IN VECTIR ICDAT.
C IUSER... MODEL-RELATED INTEGER VALUES USED BY CALCPR.  CONTAINS
C             ARRAY POINTERS TO MANAGE DATA STORAGE, AND OTHER
C             PARAMETERS.
C MNPCDF.. SUBROUTINE WHICH CALCULATES THE CDF OF A MULTIVARIATE
C          NORMAL DISTRIBUTION.
C NPAR.... NUMBER OF PARAMETERS IN VECTOR X.
C PROB.... ON RETURN, CHOICE PROBABILITY COMPUTED USING PARAMETERS IN
C             X AND DATA IN ICDAT AND RCDAT.
C RCDAT... VECTOR OF REAL DATA VALUES.
C RUSER... MODEL-RELATED REAL VALUES USED BY CALCPR.  CAN CONTAIN
C             USEFUL PARAMETERS, AND ALSO EXTRA WORK STORAGE.
C
      EXTERNAL CALCP1
      INTEGER ISIGU, IW, NALT, NW
C
      ISIGU = IUSER(16)
      IW = IUSER(15)
      NALT = MAX(1,IUSER(6))
      NW = MAX(1, IUSER(7))
      CALL CALCP1(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT,
     1            PROB, IUSER, RUSER, NALT, RUSER(ISIGU),
     2            NW, RUSER(IW), MNPCDF)
C *** LAST LINE OF CALCPR FOLLOWS ***
      END
      SUBROUTINE CALCP1(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT,
     1                PROB, IUSER, RUSER, NALT, SIGU, NW, W, MNPCDF)
C
C *** THIS SUBROUTINE CALCULATES A PROBABILITY FOR THE MNP MODEL   ***
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
      INTEGER NPAR, IERR, ICH, IALT, II, ICDAT(*), IR, IUSER(*), NALT,
     1        NATTR, NW
      REAL X(NPAR), RCDAT(*), PROB, RUSER(*),
     1                 SIGU(NALT,NALT), W(NW,NALT)
      EXTERNAL MNPCDF
C
C  ***  CALCP1 PARAMETER USAGE ***
C
C IALT.... NUMBER OF CHOICES AVAILABLE IN THE CHOICE SET.
C ICDAT... VECTOR OF INTEGER DATA VALUES.
C             IN THIS SUBROUTINE, ICDAT STORES INTEGER INDEXES WHICH
C             DEFINE WHICH OF THE NOMINAL ALTERNATIVES ARE AVAILABLE
C             IN  THE CHOICE SET.  (THIS IS FOR THE CASE WHEN THERE
C             ARE NALT NOMINAL CHOICE ALTERNATIVES, BUT NOT ALL OF
C             THEM NECESSARILY APPEAR IN EVERY SUBSET.  IF ALL NALT
C             ALTERNATIVES APPEAR IN ALL SUBSETS, THEN ICSET = NALT >0
C             SHOULD BE USED WITH IDUM = 1.
C ICH..... INTEGER INDICATING THE CHOICE.  1 <= ICH <= IALT.
C ICOV.... INDICATOR FOR TYPE OF ALTERNATIVE-SPECIFIC ERRORS,
C             = 0 FOR IID ERRORS, = 1 FOR CORRELATED ERRORS.
C             IF ICSET .NE. 0, THEN THE SAME CORRELATION MATRIX IS
C             USED FOR EVERY SUBSET.  OTHERWISE, INTEGER DATA SHOULD
C             BE USED TO IDENTIFY THE ALTERNATIVES IN EACH CHOICE SET.
C             (STORED IN IUSER.)
C IDUM... INDICATOR FOR ALTERNATIVE-SPECIFIC DUMMIES,
C             = 0 FOR NO, = 1 FOR YES.  IF ICSET .NE. 0, THEN
C             THE SAME SET OF DUMMIES IS USED FOR EACH CHOICE SET.
C             OTHERWISE, INTEGER DATA SHOULD BE USED TO IDENTIFY THE
C             ALTERNATIVES IN EACH CHOICE SET (SEE NALT BELOW).
C             (STORED IN IUSER).
C IERR.... INTEGER FOR PASSING ERROR INFORMATION.
C             IN THIS ROUTINE, IF IERR = 1 ON RETURN THEN THERE WERE
C               NO PROBLEMS.
C             IF IERR = 0 ON RETURN, THEN THE PROBABILITY COULD NOT
C               BE COMPUTED USING THE CURRENT PARAMETERS IN X.
C II...... NUMBER OF INTEGER VALUES STORED IN VECTIR ICDAT.
C IUSER... MODEL-RELATED INTEGER VALUES USED BY CALCPR.  THE FIRST
C             PORTION OF IUSER CONTAINS SUCH THINGS AS ARRAY POINTERS.
C             IUSER ALSO CONTAINS STORED VALUES OF NATTR, IDUM, ETC.
C IR...... NUMBER OF REAL VALUES STORED IN VECTOR IRDAT.
C ITASTE.. INDICATOR FOR TASTE VARIATION,
C             = 0 FOR NO TASTE VARIATION, = 1 FOR UNCORRELATED TASTE
C             VARIATION, = 2 FOR CORRELATED TASTE VARIATION.
C             (STORED IN IUSER.)
C NPAR.... NUMBER OF PARAMETERS IN VECTOR X.
C PROB.... ON RETURN, CHOICE PROBABILITY COMPUTED USING PARAMETERS IN
C             X AND DATA IN ICDAT AND RCDAT.
C RCDAT... VECTOR OF REAL DATA VALUES.
C             IN THIS SUBROUTINE, THE NUMBER OF DATA VALUES SHOULD
C             BE = IALT * NATTR SO THAT THE "GENERIC" PART OF THE
C             SCALE VALUE V MAY BE COMPUTED.
C NALT.... TOTAL NUMBER OF NOMINAL CHOICE ALTERNATIVES (IF APPLICABLE).
C             IF ICSET .NE. 0 AND IDUM = 1 OR ICOV = 1 (OR BOTH), THEN
C             NALT SHOULD BE EQUAL TO ICSET.
C             OTHERWISE, NALT SHOULD BE > 0 IF EITHER IDUM OR ICOV
C             (OR BOTH) ARE > 0, AND ICDAT SHOULD BE USED TO PASS
C             INDEX INFORMATION (SEE ICDAT ABOVE).
C NATTR... NUMBER OF ATTRIBUTES (I.E., REAL DATA VARS.) PER ALTERNATIVE.
C NW...... NUMBER OF ROWS IN THE WORK-ARRAY W.
C RUSER... MODEL-RELATED REAL VALUES USED BY CALCPR.  FOR THIS MODEL,
C             IT CONTAINS A CONSTANT FOR THE COVARIANCE MATRIX SCALE,
C             AND INFORMATION USED FOR COMPUTING STEP SIZES IN FINITE-
C             DIFFERENCE CALCULATIONS.
C SIGU.... MATRIX CONTAINING THE "UNPACKED" THE FULL COVARIANCE MATRIX
C             FOR ALL NALT ALTERNATIVE-SPECIFIC ERROR TERMS.  THE
C             MATRIX IS OF DIMENSION 2 TO FACILITATE CODING.  THE
C             NORMALIZATION USED LEAVES A ROW OF ZEROS IN THE LAST
C             (NALT) ROW.  IT IS COMPUTED BEFORE THE CALL TO MINIMIZE
C             WORK WHEN CALLS ARE TO BE REPEATED.
C W....... ARRAY CONTAINING WORKSPACE FOR COVARIANCE COMPUTATIONS.
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
C
      EXTERNAL  L7VML,  V7SCP
C
      INTEGER I, IALTM1, ICOL, ICOV, ICSET, ID, IDUM, IFAULT, IIR,
     1        IOUNIT, IPCOEF, IPDUM, IPT, IPTAST, IROW, ISCALE, ISZ,
     2        ITASTE, IX, J, JP, K, KP
C
      INTEGER MAXALT, MAXAM1, LR
      PARAMETER (MAXALT=20, MAXAM1=MAXALT-1, LR=MAXAM1*(MAXAM1-1)/2)
C
      REAL SCALE, SII
      REAL V(MAXALT), SIGMA(MAXALT,MAXALT)
      REAL Z(MAXAM1), SIGZ(MAXAM1,MAXAM1), R(LR)
C
      REAL ZERO
      PARAMETER (ZERO=0.E0)
C
C  SET UP V AND SIGMA MATRIX FOR MNP SPECIFICATION.
C
C  ALTERNATIVE-SPECIFIC DUMMIES:
C
      IALTM1 = IALT - 1
      IDUM = IUSER(8)
      IF (IDUM.NE.0) THEN
         IPDUM = IUSER(12)
C        CASE 1:  ICSET = 0.
         ICSET = IUSER(5)
         IF (ICSET.EQ.0) THEN
            DO 10 I = 1, IALT
               IX = ICDAT(I)
               IF (IX.NE.NALT) THEN
                  V(I) = X(IX+IPDUM)
               ELSE
                  V(I) = ZERO
               ENDIF
 10         CONTINUE
          ELSE
C         CASE 2:  ICSET.NE.0
            V(IALT) = ZERO
            DO 20 I = 1, IALTM1
                V(I) = X(I+IPDUM)
 20         CONTINUE
          ENDIF
      ELSE
          CALL  V7SCP(IALT, V, ZERO)
      ENDIF
C
C  BETA COEFFICIENTS:
C
      NATTR = IUSER(7)
      IF (NATTR.NE.0) THEN
         IPCOEF = IUSER(11)
         ID = 0
         DO 30 I = 1, IALT
            DO 30 K = 1, NATTR
               ID = ID + 1
               V(I) = V(I) + X(IPCOEF+K)*RCDAT(ID)
 30      CONTINUE
      ENDIF

C
C  ALTERNATIVE-SPECIFIC ERRORS:
C
      ICOV = IUSER(9)
      IF (ICOV.NE.0) THEN
         ICSET = IUSER(5)
         IF (ICSET.EQ.0) THEN
            DO 40 I = 1, IALT
               IROW = ICDAT(I)
               DO 40 J = 1, I
                  ICOL = ICDAT(J)
                  IF (ICOL.LE.IROW) THEN
                     SIGMA(I,J) = SIGU(IROW,ICOL)
                  ELSE
                     SIGMA(I,J) = SIGU(ICOL,IROW)
                  ENDIF
 40         CONTINUE
         ELSE
            DO 50 I = 1, IALT
               DO 50 J = 1, I
                  SIGMA(I,J) = SIGU(I,J)
 50         CONTINUE
         ENDIF
      ELSE
         ISCALE = IUSER(18)
         SCALE = RUSER(ISCALE)
         DO 60 I = 1, IALT
            DO 60 J = 1, I
               IF (I.EQ.J) THEN
                  SIGMA(I,J) = SCALE
               ELSE
                  SIGMA(I,J) = ZERO
               ENDIF
 60      CONTINUE
       ENDIF
C
C  TASTE VARIATION:
C
      ITASTE = IUSER(10)
      IF (ITASTE.EQ.1) THEN
C        UNCORRELATED TASTE VARIATION
C        SET UP W MATRIX:
         ID = 0
         IPTAST = IUSER(14)
         DO 70 J = 1, IALT
            IPT = IPTAST
            DO 70 K = 1, NATTR
               IPT = IPT + 1
               ID = ID + 1
               W(K,J) = X(IPT) * RCDAT(ID)
 70      CONTINUE
      ENDIF
C
      IF (ITASTE.EQ.2) THEN
C        CORRELATED TASTE VARIATION
C        SET UP W MATRIX:
         ID = 1
         IPTAST = IUSER(14) + 1
         DO 80 J = 1, IALT
            CALL  L7VML(NATTR, W(1,J), X(IPTAST), RCDAT(ID))
            ID = ID + NATTR
 80      CONTINUE
      ENDIF

      IF (ITASTE.NE.0) THEN
C        TASTE VARIATION
C        ADD W(**T)W TO SIGMA:
         DO 100 I = 1, IALT
            DO 100 J = 1, I
               DO 90 K = 1, NATTR
                  SIGMA(I,J) = SIGMA(I,J) + W(K,I)*W(K,J)
 90            CONTINUE
 100     CONTINUE
      ENDIF
C
C  SYMMETRIZE SIGMA (MAY NOT BE NECESSARY???)
C
C      IF ((ICOV.NE.0).OR.(ITASTE.NE.0)) THEN
         DO 110 I = 1, IALT
            DO 110 J = 1, I
               SIGMA(J,I) = SIGMA(I,J)
 110     CONTINUE
C      ENDIF
C
C  LOWER DIMENSION VIA STANDARD TRANSFORMATION
C  (REF. PAGE 43 OF DAGANZO OR BUNCH(1991))
      ISZ = 0
      SII = SIGMA(ICH,ICH)
      DO 130 JP = 1, IALT
         IF (JP.LT.ICH) THEN
            J = JP
         ELSE
            J = JP - 1
         ENDIF
         IF (JP.NE.ICH) THEN
            Z(J) = V(JP)-V(ICH)
            DO 120 KP = 1, JP
               IF (KP.LT.ICH) THEN
                  K = KP
               ELSE
                  K = KP - 1
               ENDIF
               IF(KP.NE.ICH) THEN
                  ISZ = ISZ + 1
                  SIGZ(J,K)=SIGMA(JP,KP)-SIGMA(ICH,KP)-SIGMA(ICH,JP)+SII
               ENDIF
 120        CONTINUE
         ENDIF
 130   CONTINUE
C
      IIR = 0
      DO 150 J = 1, IALTM1
         IF (SIGZ(J,J).LE.ZERO) THEN
            IERR = 0
            RETURN
         ENDIF
         SIGZ(J,J) = SQRT(SIGZ(J,J))
         Z(J) = Z(J)/SIGZ(J,J)
         DO 140 K = 1, J-1
               IIR = IIR + 1
               R(IIR) = SIGZ(J,K)/SIGZ(J,J)/SIGZ(K,K)
 140     CONTINUE
 150  CONTINUE
C
      IERR = 1
      CALL MNPCDF(IALTM1, Z, R, PROB, IFAULT)
      IF (IFAULT.NE.0) then
         IERR = 0
         IOUNIT = IUSER(3)
         WRITE(IOUNIT,*) ' Problem evaluating mnpcdf'
      ENDIF
C *** LAST LINE OF CALCP1 FOLLOWS ***
      END
      SUBROUTINE CALCDP(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT,
     1                  PROB0, DP, IUSER, RUSER, MNPCDF)
C
C *** THIS SUBROUTINE CALCULATES FINITE-DIFFERENCE DERIVATIVES FOR ***
C *** CHOICE PROBABILITIES.  THIS VERSION ASSUMES THAT THE CALCPR  ***
C *** BEING CALLED IS THE ONE FOR MULTINOMIAL PROBIT.  HOWEVER,    ***
C *** THE CHANGES REQUIRED FOR OTHER MODELS SHOULD BE MINOR.       ***
C *** NOTE:  THIS SUBROUTINE REQUIRES  S7GRD, AND THE ARRAYS ALPHA ***
C *** AND D SHOULD HAVE THE SAME DIMENSION AS X.                   ***
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
      INTEGER NPAR, IERR, ICH, IALT, II, ICDAT(*), IR, IUSER(*)
      REAL X(NPAR), RCDAT(*), PROB0, DP(NPAR), RUSER(*)
      EXTERNAL MNPCDF
C
      EXTERNAL CALCPR,  S7GRD,  V7SCP
      INTEGER I, ICOV, IETA0, IPCOV, IPP, IPU, IPUP, IRC, ISCALE, ISIGP,
     1        ISIGU, J, NALT, NALTM1, NFC
      REAL ETA, ETA0, PROB, SCALE, XTEMP
C
      INTEGER LX
      REAL ONE, ZERO
      PARAMETER (ZERO=0.E0, ONE=1.E0, LX=60)
C
      REAL ALPHA(LX), D(LX), WRK(6)
C
C ***  PARAMETER USAGE ***
C
C SEE CALCPR AND CALCP1
C
C *** BODY ***
C
      IERR = 1
      ICOV = IUSER(9)
      NALT = IUSER(6)
      NALTM1 = NALT - 1
      ISCALE = IUSER(18)
      SCALE = RUSER(ISCALE)
      IETA0 = IUSER(17)
      ETA0 = RUSER(IETA0)
C
      DO 10 I = 1, NPAR
          ALPHA(I) = ONE
          D(I) = ONE
 10   CONTINUE
C
      ETA = ETA0
C     ETA = ETA0/PROB
      IRC = 0
C
      PROB = PROB0
 20   CONTINUE
      CALL  S7GRD(ALPHA, D, ETA, PROB, DP, IRC,
     1             NPAR, WRK, X)
      IF (IRC.EQ.0) GO TO 40
C        IF ICOV.NE.0, SET UP AN UNPACKED SIGMA MATRIX
         IF (ICOV.NE.0) THEN
C          SQUARE THE CHOLESKY FACTOR TO GET (PACKED) SIGMA:
            IPCOV = IUSER(13)
            XTEMP = X(IPCOV)
            X(IPCOV) = SCALE
            ISIGP = IUSER(15)
            CALL  L7SQR(NALTM1, RUSER(ISIGP), X(IPCOV))
            X(IPCOV) = XTEMP
C          "UNPACK" FOR EASIER ACCESS IN CALCPR:
            IPP = ISIGP - 1
            ISIGU = IUSER(16)
            CALL  V7SCP(NALT*NALT, RUSER(ISIGU), ZERO)
            IPUP = ISIGU - 1
            DO 30 I = 1, NALTM1
               IPU = I + IPUP
               DO 30 J = 1, I
                  IPP = IPP + 1
                  RUSER(IPU) = RUSER(IPP)
                  IPU = IPU + NALT
 30         CONTINUE
         ENDIF
      CALL CALCPR(NPAR, X, NFC, ICH, IALT, II, ICDAT,
     1      IR, RCDAT, PROB, IUSER, RUSER, MNPCDF)
      IF (NFC.EQ.0) THEN
         IERR = 0
         RETURN
      ENDIF
      GO TO 20
 40   CONTINUE
C
C *** LAST LINE OF CALCDP FOLLOWS ***
      END
      SUBROUTINE  S7GRD (ALPHA, D, ETA0, FX, G, IRC, N, W, X)
C
C  ***  COMPUTE FINITE DIFFERENCE GRADIENT BY STWEART*S SCHEME  ***
C  ***  THIS IS SGRAD2 FROM TOMS ALGORITHM 611.
C
C     ***  PARAMETERS  ***
C
      INTEGER IRC, N
      REAL ALPHA(N), D(N), ETA0, FX, G(N), W(6), X(N)
C
C.......................................................................
C
C     ***  PURPOSE  ***
C
C        THIS SUBROUTINE USES AN EMBELLISHED FORM OF THE FINITE-DIFFER-
C     ENCE SCHEME PROPOSED BY STEWART (REF. 1) TO APPROXIMATE THE
C     GRADIENT OF THE FUNCTION F(X), WHOSE VALUES ARE SUPPLIED BY
C     REVERSE COMMUNICATION.
C
C     ***  PARAMETER DESCRIPTION  ***
C
C  ALPHA IN  (APPROXIMATE) DIAGONAL ELEMENTS OF THE HESSIAN OF F(X).
C      D IN  SCALE VECTOR SUCH THAT D(I)*X(I), I = 1,...,N, ARE IN
C             COMPARABLE UNITS.
C   ETA0 IN  ESTIMATED BOUND ON RELATIVE ERROR IN THE FUNCTION VALUE...
C             (TRUE VALUE) = (COMPUTED VALUE)*(1+E),   WHERE
C             ABS(E) .LE. ETA0.
C     FX I/O ON INPUT,  FX  MUST BE THE COMPUTED VALUE OF F(X).  ON
C             OUTPUT WITH IRC = 0, FX HAS BEEN RESTORED TO ITS ORIGINAL
C             VALUE, THE ONE IT HAD WHEN  S7GRD WAS LAST CALLED WITH
C             IRC = 0.
C      G I/O ON INPUT WITH IRC = 0, G SHOULD CONTAIN AN APPROXIMATION
C             TO THE GRADIENT OF F NEAR X, E.G., THE GRADIENT AT THE
C             PREVIOUS ITERATE.  WHEN  S7GRD RETURNS WITH IRC = 0, G IS
C             THE DESIRED FINITE-DIFFERENCE APPROXIMATION TO THE
C             GRADIENT AT X.
C    IRC I/O INPUT/RETURN CODE... BEFORE THE VERY FIRST CALL ON  S7GRD,
C             THE CALLER MUST SET IRC TO 0.  WHENEVER  S7GRD RETURNS A
C             NONZERO VALUE FOR IRC, IT HAS PERTURBED SOME COMPONENT OF
C             X... THE CALLER SHOULD EVALUATE F(X) AND CALL  S7GRD
C             AGAIN WITH FX = F(X).
C      N IN  THE NUMBER OF VARIABLES (COMPONENTS OF X) ON WHICH F
C             DEPENDS.
C      X I/O ON INPUT WITH IRC = 0, X IS THE POINT AT WHICH THE
C             GRADIENT OF F IS DESIRED.  ON OUTPUT WITH IRC NONZERO, X
C             IS THE POINT AT WHICH F SHOULD BE EVALUATED.  ON OUTPUT
C             WITH IRC = 0, X HAS BEEN RESTORED TO ITS ORIGINAL VALUE
C             (THE ONE IT HAD WHEN  S7GRD WAS LAST CALLED WITH IRC = 0)
C             AND G CONTAINS THE DESIRED GRADIENT APPROXIMATION.
C      W I/O WORK VECTOR OF LENGTH 6 IN WHICH  S7GRD SAVES CERTAIN
C             QUANTITIES WHILE THE CALLER IS EVALUATING F(X) AT A
C             PERTURBED X.
C
C     ***  APPLICATION AND USAGE RESTRICTIONS  ***
C
C        THIS ROUTINE IS INTENDED FOR USE WITH QUASI-NEWTON ROUTINES
C     FOR UNCONSTRAINED MINIMIZATION (IN WHICH CASE  ALPHA  COMES FROM
C     THE DIAGONAL OF THE QUASI-NEWTON HESSIAN APPROXIMATION).
C
C     ***  ALGORITHM NOTES  ***
C
C        THIS CODE DEPARTS FROM THE SCHEME PROPOSED BY STEWART (REF. 1)
C     IN ITS GUARDING AGAINST OVERLY LARGE OR SMALL STEP SIZES AND ITS
C     HANDLING OF SPECIAL CASES (SUCH AS ZERO COMPONENTS OF ALPHA OR G).
C
C     ***  REFERENCES  ***
C
C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION
C        METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES,
C        J. ASSOC. COMPUT. MACH. 14, PP. 72-83.
C
C     ***  HISTORY  ***
C
C     DESIGNED AND CODED BY DAVID M. GAY (SUMMER 1977/SUMMER 1980).
C
C     ***  GENERAL  ***
C
C        THIS ROUTINE WAS PREPARED IN CONNECTION WITH WORK SUPPORTED BY
C     THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS76-00324 AND
C     MCS-7906671.
C
C.......................................................................
C
C     *****  EXTERNAL FUNCTION  *****
C
      EXTERNAL  R7MDC
      REAL  R7MDC
C  R7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS.
C
C     ***** LOCAL VARIABLES *****
C
      INTEGER FH, FX0, HSAVE, I, XISAVE
      REAL AAI, AFX, AFXETA, AGI, ALPHAI, AXI, AXIBAR,
     1                 DISCON, ETA, GI, H, HMIN
      REAL C2000, FOUR, HMAX0, HMIN0, H0, MACHEP, ONE, P002,
     1                 THREE, TWO, ZERO
C
      PARAMETER (C2000=2.0E+3, FOUR=4.0E+0, HMAX0=0.02E+0, HMIN0=5.0E+1,
     1     ONE=1.0E+0, P002=0.002E+0, THREE=3.0E+0,
     2     TWO=2.0E+0, ZERO=0.0E+0)
      PARAMETER (FH=3, FX0=4, HSAVE=5, XISAVE=6)
C
C---------------------------------  BODY  ------------------------------
C
      IF (IRC) 50, 10, 110
C
C     ***  FRESH START -- GET MACHINE-DEPENDENT CONSTANTS  ***
C
C     STORE MACHEP IN W(1) AND H0 IN W(2), WHERE MACHEP IS THE UNIT
C     ROUNDOFF (THE SMALLEST POSITIVE NUMBER SUCH THAT
C     1 + MACHEP .GT. 1  AND  1 - MACHEP .LT. 1),  AND  H0 IS THE
C     SQUARE-ROOT OF MACHEP.
C
 10   W(1) =  R7MDC(3)
      W(2) =  SQRT(W(1))
C
      W(FX0) = FX
C
C     ***  INCREMENT  I  AND START COMPUTING  G(I)  ***
C
 20   I =  ABS(IRC) + 1
      IF (I .GT. N) GO TO 120
         IRC = I
         AFX =  ABS(W(FX0))
         MACHEP = W(1)
         H0 = W(2)
         HMIN = HMIN0 * MACHEP
         W(XISAVE) = X(I)
         AXI =  ABS(X(I))
         AXIBAR =   MAX(AXI, ONE/D(I))
         GI = G(I)
         AGI =  ABS(GI)
         ETA =  ABS(ETA0)
         IF (AFX .GT. ZERO) ETA =   MAX(ETA, AGI*AXI*MACHEP/AFX)
         ALPHAI = ALPHA(I)
         IF (ALPHAI .EQ. ZERO) GO TO 80
         IF (GI .EQ. ZERO .OR. FX .EQ. ZERO) GO TO 90
         AFXETA = AFX*ETA
         AAI =  ABS(ALPHAI)
C
C        *** COMPUTE H = STEWART*S FORWARD-DIFFERENCE STEP SIZE.
C
         IF (GI**2 .LE. AFXETA*AAI) GO TO 30
              H = TWO* SQRT(AFXETA/AAI)
              H = H*(ONE - AAI*H/(THREE*AAI*H + FOUR*AGI))
              GO TO 40
 30      H = TWO*(AFXETA*AGI/(AAI**2))**(ONE/THREE)
         H = H*(ONE - TWO*AGI/(THREE*AAI*H + FOUR*AGI))
C
C        ***  ENSURE THAT  H  IS NOT INSIGNIFICANTLY SMALL  ***
C
 40      H =   MAX(H, HMIN*AXIBAR)
C
C        *** USE FORWARD DIFFERENCE IF BOUND ON TRUNCATION ERROR IS AT
C        *** MOST 10**-3.
C
         IF (AAI*H .LE. P002*AGI) GO TO 70
C
C        *** COMPUTE H = STEWART*S STEP FOR CENTRAL DIFFERENCE.
C
         DISCON = C2000*AFXETA
         H = DISCON/(AGI +  SQRT(GI**2 + AAI*DISCON))
C
C        ***  ENSURE THAT  H  IS NEITHER TOO SMALL NOR TOO BIG  ***
C
         H =   MAX(H, HMIN*AXIBAR)
         IF (H .GE. HMAX0*AXIBAR) H = AXIBAR * H0**(TWO/THREE)
C
C        ***  COMPUTE CENTRAL DIFFERENCE  ***
C
         IRC = -I
         GO TO 100
C
 50      H = -W(HSAVE)
         I =  ABS(IRC)
         IF (H .GT. ZERO) GO TO 60
         W(FH) = FX
         GO TO 100
C
 60      G(I) = (W(FH) - FX) / (TWO * H)
         X(I) = W(XISAVE)
         GO TO 20
C
C     ***  COMPUTE FORWARD DIFFERENCES IN VARIOUS CASES  ***
C
 70      IF (H .GE. HMAX0*AXIBAR) H = H0 * AXIBAR
         IF (ALPHAI*GI .LT. ZERO) H = -H
         GO TO 100
 80      H = AXIBAR
         GO TO 100
 90      H = H0 * AXIBAR
C
 100     X(I) = W(XISAVE) + H
         W(HSAVE) = H
         GO TO 999
C
C     ***  COMPUTE ACTUAL FORWARD DIFFERENCE  ***
C
 110     G(IRC) = (FX - W(FX0)) / W(HSAVE)
         X(IRC) = W(XISAVE)
         GO TO 20
C
C  ***  RESTORE FX AND INDICATE THAT G HAS BEEN COMPUTED  ***
C
 120  FX = W(FX0)
      IRC = 0
C
 999  RETURN
C  ***  LAST CARD OF  S7GRD FOLLOWS  ***
      END
      SUBROUTINE PCMRJ(NOBS, NPAR, X, NF, NEED, R, RP, UI, UR, UF)
      INTEGER NOBS, NPAR, NF, NEED(2), UI(*)
      REAL X(NPAR), R(NOBS), RP(NPAR,NOBS), UR(*)
      EXTERNAL UF
C
      EXTERNAL PCMRJ1
C
      INTEGER ICP, IICDAT, IICH, IIIV, IIRV, IIU, INALT, IRCDAT, IRU
C
C *** BODY ***
C
      IIU    = UI(1)
      IICH   = UI(2)
      INALT  = UI(3)
      IIIV   = UI(4)
      IIRV   = UI(5)
      IICDAT = UI(6)
C
      IRU = UI(7)
      ICP = UI(8)
C     IRW = UI(9)
      IRCDAT = UI(10)
C
      CALL PCMRJ1(NOBS, NPAR, X, NF, NEED, R, RP,
     1     UI(IIU), UI(IICH), UI(INALT), UI(IIIV), UI(IIRV), UI(IICDAT),
     2     UR(IRU), UR(ICP), UR(IRCDAT), UF)
 999  RETURN
C *** LAST LINE OF PCMRJ FOLLOWS ***
      END
      SUBROUTINE PCMRJ1(NOBS, NPAR, X, NF, NEED, R, RP,
     1     IUSER, ICHV, NALTV, IIV, IRV, ICDAT,
     2     RUSER, CPROB, RCDAT, UF)
C
C *** THIS SUBROUTINE EXPANDS THE STORAGE IN UI AND UR TO MAKE THEM ***
C *** COMPATIBLE WITH ESTIMATION OF CHOICE MODELS.                  ***
C
      INTEGER NOBS, NPAR, NF, NEED(2), IUSER(*), ICHV(NOBS),
     1        NALTV(NOBS), IIV(NOBS), IRV(NOBS), ICDAT(*)
      REAL X(NPAR), R(NOBS), RP(NPAR,NOBS), RUSER(*),
     1                 CPROB(NOBS,2), RCDAT(*)
      EXTERNAL UF
C
      EXTERNAL CALCDP, CALCPR,  L7SQR,  V7SCP
C
      INTEGER I, IALT, ICH, ICOV, IERR, II, III, IIR, IOBS, IOUNIT,
     1        IPCOV, IPP, IPU, IPUP, IR, ISCALE, ISIGP, ISIGU, J, KS,
     2        NALT, NALTM1, NFC
      REAL PROB, SCALE, XTEMP
C
      INTEGER LX
      REAL ONE, ZERO
      PARAMETER (ZERO=0.E0, ONE=1.E0, LX=60)
C
C ARRAYS:
C
C CPROB... VECTOR FOR STORING CHOICE PROBABILITIES.  CPROB(IOBS,J)
C             FOR J=1,2 STORES CHOICE PROBABILITIES FOR OBSERVATION
C             IOBS.  ONE IS THE CURRENT PROBABILITY, WHILE THE OTHER
C             ONE IS THE PROBABILITY COMPUTED AT THE PREVIOUS TRIAL
C             X.  THE CODE KEEPS TRACK OF WHICH IS WHICH USING THE
C             POINTERS STORED IN IUSER(1) AND IUSER(2).  THIS IS USED
C             IN VARIOUS WAYS TO MAKE COMPUTATION MORE EFFICIENT.
C ICHV.... VECTOR OF LENGTH NOBS.  ICHV(IOBS) CONTAINS THE INDEX OF
C             THE CHOSEN ALTERNATIVE FOR OBSERVATION IOBS.
C IIV..... VECTOR OF LENGHT NOBS.  IIV(IOBS) INDICATES THE NUMBER OF
C             INTEGER DATA VALUES STORED IN ICDAT FOR OBSERVATION IOBS.
C IRV..... VECTOR OF LENGHT NOBS.  IRV(IOBS) INDICATES THE NUMBER OF
C             REAL DATA VALUES STORED IN RCDAT FOR OBSERVATION IOBS.
C NALTV... VECTOR OF LENGHT NOBS.  NALTV(IOBS) INDICATES THE NUMBER OF
C             CHOICES AVAILABLE FOR OBSERVATION IOBS.
C
C *** BODY ***
C
      ICOV = IUSER(9)
      NALT = IUSER(6)
      NALTM1 = NALT - 1
      ISCALE = IUSER(18)
      SCALE = RUSER(ISCALE)
C
      IF (NEED(1).EQ.1) THEN
C
C *** CALCULATE RESIDUAL VECTOR ***
         KS = 1
         IF (NEED(2).EQ.IUSER(1)) KS = 2
         IUSER(KS) = NF
C
C      IF ICOV.NE.0, SET UP AN UNPACKED SIGMA MATRIX
         IF (ICOV.NE.0) THEN
C            SQUARE THE CHOLESKY FACTOR TO GET (PACKED) SIGMA:
             IPCOV = IUSER(13)
             XTEMP = X(IPCOV)
             X(IPCOV) = SCALE
             ISIGP = IUSER(15)
             CALL  L7SQR(NALTM1, RUSER(ISIGP), X(IPCOV))
             X(IPCOV) = XTEMP
C            "UNPACK" FOR EASIER ACCESS IN CALCPR:
             IPP = ISIGP - 1
             ISIGU = IUSER(16)
             CALL  V7SCP(NALT*NALT, RUSER(ISIGU), ZERO)
             IPUP = ISIGU - 1
             DO 10 I = 1, NALTM1
                IPU = I + IPUP
                DO 10 J = 1, I
                   IPP = IPP + 1
                   RUSER(IPU) = RUSER(IPP)
                   IPU = IPU + NALT
 10          CONTINUE
         ENDIF
         III = 1
         IIR = 1
         DO 20 IOBS = 1, NOBS
             ICH = ICHV(IOBS)
             IALT = NALTV(IOBS)
             II = IIV(IOBS)
             IR = IRV(IOBS)
             CALL CALCPR(NPAR, X, NFC, ICH, IALT, II, ICDAT(III),
     1            IR, RCDAT(IIR), PROB, IUSER, RUSER, UF)
             IF ((PROB.LE.ZERO).OR.(PROB.GT.ONE).OR.(NFC.EQ.0)) THEN
                NF = 0
                RETURN
             ENDIF
             R(IOBS) = PROB
             CPROB(IOBS,KS) = PROB
             III = III + II
             IIR = IIR + IR
 20      CONTINUE
      ELSE
C
C *** CALCULATE JACOBIAN OF RESIDUAL VECTOR ***
C
         KS = 1
         IF (IUSER(1).NE.NF) KS = 2
         IF (IUSER(KS).NE.NF) THEN
            IOUNIT = IUSER(3)
            WRITE(IOUNIT,*) ' PROBLEM WITH INITIAL ESTIMATE...'
         ENDIF
C
         III = 1
         IIR = 1
         DO 30 IOBS = 1, NOBS
             ICH = ICHV(IOBS)
             IALT = NALTV(IOBS)
             II = IIV(IOBS)
             IR = IRV(IOBS)
             PROB = CPROB(IOBS,KS)
             CALL CALCDP(NPAR, X, IERR, ICH, IALT, II, ICDAT(III),
     1             IR, RCDAT(IIR), PROB, RP(1,IOBS), IUSER, RUSER, UF)
             IF (IERR.EQ.0) THEN
                NF = 0
                RETURN
             ENDIF
             III = III + II
             IIR = IIR + IR
 30      CONTINUE
      ENDIF
 999  RETURN
C *** LAST LINE OF PCMRJ1 FOLLOWS ***
      END
      SUBROUTINE PCMRHO(NEED, F, NOBS, NF, XN, R, RD, UI, UR, W)
      INTEGER NEED(2), NOBS, NF, UI(*)
      REAL F, XN(*), R(*), RD(NOBS,*), UR(*), W(NOBS)
C
      INTEGER ICP, IOBS, IOUNIT, IRW, WEIGHT, KS
      REAL OOR, VT
C
      REAL NEGONE, ZERO
      PARAMETER (NEGONE=-1.E0, ZERO=0.E0)
C
C *** BODY ***
C
      WEIGHT = UI(14)
      IF (NEED(1).EQ.1) THEN
         VT = ZERO
         IF (WEIGHT.EQ.0) THEN
            DO 10 IOBS = 1, NOBS
                VT = VT - LOG(R(IOBS))
 10         CONTINUE
         ELSE
            IRW = UI(9)
            DO 20 IOBS = 1, NOBS
                VT = VT - UR(IRW) * LOG(R(IOBS))
                IRW = IRW + 1
 20         CONTINUE
         ENDIF
         F = VT
      ELSE
         KS = 1
         IF (UI(11).NE.NF) KS = 2
         IF (UI(10+KS).NE.NF) THEN
            IOUNIT = UI(13)
            WRITE(IOUNIT,*) ' PROBLEM WITH INITIAL POINT...'
            NF = 0
            RETURN
         ENDIF
         ICP = UI(8)
         IF (KS.EQ.2) ICP = ICP + NOBS
         IF (WEIGHT.EQ.0) THEN
            DO 30 IOBS = 1, NOBS
                OOR = NEGONE/UR(ICP)
                R(IOBS) = OOR
                W(IOBS) = R(IOBS) * OOR
                RD(IOBS,1) = W(IOBS)
                ICP = ICP + 1
 30         CONTINUE
         ELSE
            IRW = UI(9)
            DO 40 IOBS = 1, NOBS
                OOR = NEGONE/UR(ICP)
                R(IOBS) = UR(IRW) * OOR
                W(IOBS) = R(IOBS) * OOR
                RD(IOBS,1) = W(IOBS)
                ICP = ICP + 1
                IRW = IRW + 1
 40         CONTINUE
         ENDIF
      ENDIF
 999  RETURN
C *** LAST LINE OF PCMRHO FOLLOWS ***
      END
      SUBROUTINE FPRINT(NOBS, NPAR, X, NF, UI, UR, UF)
      INTEGER NOBS, NPAR, NF, UI(*)
      REAL X(NPAR), UR(*)
      EXTERNAL UF
C
      EXTERNAL FPRNT1
C
      INTEGER ICP, IICDAT, IICH, IIIV, IIRV, IIU, INALT, IRCDAT, IRU,
     1        IRW
C
C *** BODY ***
C
      IIU    = UI(1)
      IICH   = UI(2)
      INALT  = UI(3)
      IIIV   = UI(4)
      IIRV   = UI(5)
      IICDAT = UI(6)
C
      IRU = UI(7)
      ICP = UI(8)
      IRW = UI(9)
      IRCDAT = UI(10)
C
      CALL FPRNT1(NOBS, NPAR, X, NF,
     1     UI(IIU), UI(IICH), UI(INALT), UI(IIIV), UI(IIRV), UI(IICDAT),
     2     UR(IRU), UR(IRCDAT), UR(IRW), UF)
 999  RETURN
C *** LAST LINE OF FPRINT FOLLOWS ***
      END
      SUBROUTINE FPRNT1(NOBS, NPAR, X, NF,
     1     IUSER, ICHV, NALTV, IIV, IRV, ICDAT,
     2     RUSER, RCDAT, WT, UF)
C
C *** THIS SUBROUTINE EXPANDS THE STORAGE IN UI AND UR TO MAKE THEM ***
C *** COMPATIBLE WITH ESTIMATION OF CHOICE MODELS.                  ***
C *** SEE PCMRJ1 DOCUMENTATION ON ARRAYS.                           ***
C
      INTEGER NOBS, NPAR, NF, IUSER(*), ICHV(NOBS), NALTV(NOBS),
     1        IIV(NOBS), IRV(NOBS), ICDAT(*)
      REAL X(NPAR), RUSER(*), RCDAT(*), WT(NOBS)
      EXTERNAL UF
C
      EXTERNAL CALCPR,  L7SQR,  V7SCP
C
      INTEGER I, IALT, ICH, ICOV, ICSET, II, III, IIR, IOBS, IOUNIT,
     1        IPCOV, IPP, IPU, IPUP, IR, ISCALE, ISIGP, ISIGU, J, NALT,
     2        NALTM1, NFC
      REAL FPROB(20), PROB, SCALE, XTEMP
C
      INTEGER LX
      REAL ONE, ZERO
      PARAMETER (ZERO=0.E0, ONE=1.E0, LX=60)
C
C *** BODY ***
C
      ICOV = IUSER(9)
      ICSET = IUSER(5)
      IOUNIT = IUSER(3)
      NALT = IUSER(6)
      NALTM1 = NALT - 1
      ISCALE = IUSER(18)
      SCALE = RUSER(ISCALE)
C
      WRITE(IOUNIT, 10)
 10   FORMAT(//,' FINAL CHOICE SET PROBABILITIES: ',/)
C
C      IF ICOV.NE.0, SET UP AN UNPACKED SIGMA MATRIX
         IF (ICOV.NE.0) THEN
C            SQUARE THE CHOLESKY FACTOR TO GET (PACKED) SIGMA:
             IPCOV = IUSER(13)
             XTEMP = X(IPCOV)
             X(IPCOV) = SCALE
             ISIGP = IUSER(15)
             CALL  L7SQR(NALTM1, RUSER(ISIGP), X(IPCOV))
             X(IPCOV) = XTEMP
C            "UNPACK" FOR EASIER ACCESS IN CALCPR:
             IPP = ISIGP - 1
             ISIGU = IUSER(16)
             CALL  V7SCP(NALT*NALT, RUSER(ISIGU), ZERO)
             IPUP = ISIGU - 1
             DO 20 I = 1, NALTM1
                IPU = I + IPUP
                DO 20 J = 1, I
                   IPP = IPP + 1
                   RUSER(IPU) = RUSER(IPP)
                   IPU = IPU + NALT
 20          CONTINUE
         ENDIF
         III = 1
         IIR = 1
         DO 90 IOBS = 1, NOBS
             ICH = ICHV(IOBS)
             IALT = NALTV(IOBS)
             II = IIV(IOBS)
             IR = IRV(IOBS)
             DO 30 I = 1, IALT
                CALL CALCPR(NPAR, X, NFC, I, IALT, II, ICDAT(III),
     1               IR, RCDAT(IIR), PROB, IUSER, RUSER, UF)
                FPROB(I) = PROB
 30          CONTINUE
             WRITE(IOUNIT, 40) IOBS
 40          FORMAT(/,' IOBS: ',I4)
             IF (ICSET.EQ.0) WRITE(IOUNIT,50) (ICDAT(I),I=1,IALT)
 50          FORMAT('    CHOICE SET: ',20I3)
             WRITE(IOUNIT, 60) IALT, ICH, WT(IOBS)
 60          FORMAT('    NO. OF ALTS: ',I2,'    ICH: ',I2,
     1              '    WT: ',F7.3)
             WRITE(IOUNIT, 70) (FPROB(I),I=1,IALT)
 70          FORMAT('    PROBS: ',8F7.4,/,18X,8F7.4,/,18X,4F7.4)
             WRITE(IOUNIT, 80) FPROB(ICH)
 80          FORMAT('    PROB(ICH): ',F7.4)
             III = III + II
             IIR = IIR + IR
 90      CONTINUE
C
 999  RETURN
C *** LAST LINE OF FPRNT1 FOLLOWS ***
      END

ext: f

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
The ASUCD Senate Office is located in room 354 of the Memorial Union and is usually locked when nobody is around. It also goes by the name of ASUCD Legislative Office.

The purpose of this office is to allow senators, and their interns as well as commission chairs and commissioners, a place to hold office hours, and work on ASUCD business, as well as any of their own personal projects. For the senators use it has three computers, a circular table, a X HP LaserJet 4 printer (as seen in the film Office Space)X a new printer which never has any fucking toner, and three maybe working computers (one is found frequently having loose components or being unplugged.

This office is busiest Mondays through Thursdays, during the afternoon, and is practically empty on Fridays and during Senate Meetings. If you want to meet a senator, or leave a message for him or her, going to the senate office would be the place to go. Many students have admitted feeling intimidated by the third floor of the MU, especially ASUCD Senate Office due to its cliquish vibe, thus making the office the hub of the proverbial ASUCD Ivory Tower.

Though this office is mainly for ASUCD senatorial and commission use, CalPIRG agents and involved students (such as Users/BrentLaabs) use this office as well, for generally good purposes (they do yield computers and chairs to the current senators of course. In addition there has been much controversy on the wiki involving many edits made from the senate offices computers, such as those made under the Users/JoeLevy username, however this is usually limited to the wiki:drama:Daviswiki Deepthroat / Fingercuffs Incident drama surrounding the student government and has little impact to the whole of the Davis community.

Presumably only Senators and Commission Chairs have access to the door code but this is CalPIRG not always the case.

As a random note, the place is often a mess as the senators are quite busy folk, and will not have much time to clean up, then again, it is their office anyway. ~DavePoole

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
! :::::::::::::::::::::::: TRIMBD :::::::::::::::::::::::::::;
!  if used array is completely set (=1.) then return set=true,
!  otherwise return false, along with the dimensions of the smallest
!  rectangle containing all unset points in il,ir,jb,jt.
! ::::::::::::::::::::::::::::::::::::;::::::::::::::::::::::;
!
subroutine trimbd(used,nrow,ncol,set,unset_rect)

    implicit none

    ! Input
    integer, intent(in) :: nrow, ncol
    integer(kind=1), intent(in) :: used(nrow,ncol)

    ! Output
    logical, intent(out) :: set
!     integer, intent(out) :: il, ir, jb, jt
    integer, intent(out) :: unset_rect(4)

    ! Locals
    integer :: i, j, utot
    integer(kind=1) :: check

       utot = 0
        do 100 j = 1,ncol
        do 100 i = 1,nrow
100        utot = utot + used(i,j)

    if (utot .eq. nrow * ncol ) then
        set = .true.
    else
        set = .false.

        check = 1
        do i = 1,nrow
             do j = 1,ncol
                check = min(check,used(i,j))
            enddo
            unset_rect(1) = i
            if (check == 0) exit
        enddo

        check = 1
        do i = 1,nrow
            do j = 1,ncol
                check = min(check,used(nrow - i + 1,j))
            enddo
            unset_rect(2) = nrow - i + 1
            if (check == 0) exit
        enddo

        check = 1
        do j = 1,ncol
            do i = 1,nrow
                check = min(check,used(i,j))
            enddo
           unset_rect(3) = j
           if (check == 0) exit
        enddo

       check = 1
        do j = 1,ncol
            do i = 1,nrow
                check = min(check,used(i,ncol - j + 1))
            enddo
            unset_rect(4) = ncol - j + 1
            if (check == 0) exit
        enddo

    endif

end subroutine trimbd

max_stars_repo_name: ian-r-rose/visualization

max_issues_repo_name: ian-r-rose/visualization

max_forks_repo_name: ian-r-rose/visualization


text:
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
!     HOANG Danh Tai - Laboratoire de Physique Théorique et Modélisation
!     UMR 8089 CNRS-Université de Cergy-Pontoise
!     2, Avenue Adolphe Chauvin, 95032 Cergy-Pontoise Cedex, France
!----------------------------------------------------------------------------------------------------!
!     PROGRAMME: MONTE CARLO TRANSPORT OF SPIN ON MODERN POTTS
!     4/5/2011: sua lai, tinh r 1 lan de chay nhanh hon
!     CHUONG TRINH NAY DA CHUAN, DA TEST, THOI GIAN CHAY LA NHANH NHAT (chieu toi ngay 4/5/2011)
!     13/5/2011: sua lai GS cho he N=24, sua lai cach tinh Mz
!!    17/05/2011: sua lai config_GS cho toan bo gia tri r0 va d. (d=1,2,3,4,... )
!!    19/05/2011: sua thanh phan thu 2 cua dipolar, chi su dung cac lien ket song song
!!    01/06/2011: histogram
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!

      PROGRAM main_transport_modern_potts
      IMPLICIT NONE

      CHARACTER (LEN=150):: CONFIG_INI_DIREC_X
      CHARACTER (LEN=50) :: name
      CHARACTER (LEN=15) :: tmp
      CHARACTER (LEN=3)  :: spinAt
      CHARACTER (256)    :: Ligne20,Ligne126

      REAL    (KIND=8),PARAMETER :: aaa=5.
      REAL    (KIND=8),PARAMETER :: nul=0.

      INTEGER (KIND=4) :: iT,nT,i_case,d,i_layer,n_layer,m
      INTEGER (KIND=4) :: i,j,k,ip,jp,kp,im,jm,km,natx,naty,natz,natx_p,naty_p,natz_p,number

      INTEGER (KIND=4) :: i_loop,i_loop2,n_equi_reseau1,n_equi_reseau2,n_average_thermal,i_times,n_times

      INTEGER (KIND=4) :: deli0,delj0,delk0,r01,iv,n_vois,i1,j1,k1,k2,na

      REAL    (KIND=8) :: T,delT,Tmax,Tmin,rdn_etat,rdn_mtp,tab_tmp,rdn_config,n_atom
      REAL    (KIND=8) :: x,y,z,EN_tmp1,EN_tmp2
      REAL    (KIND=8) :: energy,energy_2,Mz,Mz_2
      REAL    (KIND=8) :: EN_moy,EN_2_moy,Mz_moy,Mz_2_moy,Cv,Ksi
      REAL    (KIND=8) :: r0,r1,D_tmp,EN_D_tmpx,EN_D_tmpy,EN_D_tmpz
      REAL    (KIND =8):: M12,M13,M23,M12x,M12y,M12z,M13x,M13y,M13z,M23x,M23y,M23z

      INTEGER (KIND=8),DIMENSION(3)                  :: clock

      REAL    (KIND=8),DIMENSION(:,:,:,:),ALLOCATABLE  :: spin

      REAL    (KIND=8),DIMENSION(:),      ALLOCATABLE :: deli,delj,delk
      INTEGER (KIND=4),DIMENSION(:),      ALLOCATABLE :: int_deli,int_delj,int_delk
      REAL    (KIND=8),DIMENSION(:),      ALLOCATABLE :: r

       !!! Khai bao phan Histogram
      INTEGER (KIND=4) :: i_histo,n_EN_histo
      REAL    (KIND=8) :: EN_min_histo,EN_max_histo,del_EN_histo

      REAL    (KIND=8),DIMENSION(:),ALLOCATABLE :: EN_histo,H_histo,P_histo

!!!==================================================================================================
!!!==================================================================================================
      CALL system('rm -r config_ini_3D')
      CALL system('mkdir config_ini_3D')
      CALL system('rm -r config_3D')
      CALL system('mkdir config_3D')
      CALL system('rm *.dat*')

      CALL ini_rdm_number()

      CALL read_input_parameter_file()

      natx_p=natx+1
      naty_p=naty+1
      natz_p=natz+1

      n_atom=real(natx*naty*natz)
      na=natx*naty*natz
      r01=int(r0)

      n_layer=natz/(2*d)

      ALLOCATE(spin(natx_p,naty_p,natz_p,3))
      ALLOCATE(EN_histo(n_EN_histo))
      ALLOCATE(H_histo(n_EN_histo))
      ALLOCATE(P_histo(n_EN_histo))

      IF (nT==1) THEN
            delT=0.
      ELSE
            delT=(Tmax-Tmin)/real(nT-1)
      END IF

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!! ====== MAIN PROGRAM === MAIN PROGRAM === MAIN PROGRAM === MAIN PROGRAM === MAIN PROGRAM =====
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      OPEN(unit=20,file='average_thermal.dat')

      CALL load_config_ini()

      CALL number_voisin()

      CALL energy_gs()

      CALL write_config_ini_3D()

      CALL calcul_time_run()
!!!======================================
!!!======================================
      DO iT=1,nT
            WRITE(*,*)'iT = ', iT

            T=Tmin+delT*real(iT-1)

            DO i_times=1,n_times
                  CALL equi_reseau1()
                  CALL average_thermal()
            END DO

            !CALL value_thermal()

            CALL write_config_3D()

      END DO

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      CLOSE(20)

      DEALLOCATE(spin)

      CONTAINS

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!! SUBROUTINE init_rdm_number()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE ini_rdm_number()
      IMPLICIT NONE

      INTEGER (KIND=4) :: i_time,i

      CALL ITIME(clock)
      i_time=(clock(1)+1)*(clock(2)+1)*(clock(3)+1)

      DO i=1,i_time
            CALL random_number(tab_tmp)
      ENDDO

      END SUBROUTINE ini_rdm_number

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!! SUBROUTINE read_input_parameter_file()
!!! OPEN the parameter from file "parameter.in"
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE read_input_parameter_file()
      IMPLICIT NONE

      CHARACTER (LEN=150) :: tamp
      OPEN(11,file='1parameter.in')

      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A30,(A10))')   tamp, CONFIG_INI_DIREC_X
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A30,(I5))')    tamp, natx
      READ(11, '(A30,(I5))')    tamp, naty
      READ(11, '(A30,(I5))')    tamp, natz
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A30,(I5))')    tamp, d
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A30,(F7.4))')  tamp, D_tmp
      READ(11, '(A50)')         tamp
      READ(11, '(A30,(F7.4))')  tamp, r0
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A30,(I8))')    tamp, nT
      READ(11, '(A30,(F7.4))')  tamp, Tmin
      READ(11, '(A30,(F7.4))')  tamp, Tmax
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A30,(I8))')    tamp,n_equi_reseau1
      READ(11, '(A30,(I8))')    tamp,n_equi_reseau2
      READ(11, '(A30,(I8))')    tamp,n_average_thermal
      READ(11, '(A30,(I8))')    tamp,n_times
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp
      READ(11, '(A30,(I8))')    tamp,n_EN_histo
      READ(11, '(A30,(F7.4))')  tamp,EN_min_histo
      READ(11, '(A30,(F7.4))')  tamp,EN_max_histo
      READ(11, '(A50)')         tamp
      READ(11, '(A50)')         tamp

      CLOSE(11)

      END SUBROUTINE read_input_parameter_file

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!! SUBROUTINE calcul_time_run
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE calcul_time_run()

      IMPLICIT NONE

      CHARACTER (8)  :: date
      CHARACTER (10) :: time
      CHARACTER (5)  :: zone
      INTEGER,DIMENSION(8) :: value_time1,value_time2,value_time3,value_time4

      INTEGER (KIND=4) :: n_date_run,n_hour_run,n_minute_run,i
      REAL    (KIND=8) :: n_time_run_total

      CALL system('rm time_run.dat')

      OPEN (90,file='time_run.dat')

!!! Tinh thoi gian moi lan thermalisation ==============================

      CALL date_and_time(date,time,zone,value_time1)
      WRITE(90,'(5I6)')value_time1(5),value_time1(6),value_time1(7),value_time1(8)

      DO i=1,1

      CALL equi_reseau()

      END DO

      CALL date_and_time(date,time,zone,value_time2)
      WRITE(90,'(5I6)')value_time2(5),value_time2(6),value_time2(7),value_time2(8)
!!! ===================================================================

!!! Tinh thoi gian moi lan tinh gia tri ==============================

      CALL date_and_time(date,time,zone,value_time3)
      WRITE(90,'(5I6)')value_time3(5),value_time3(6),value_time3(7),value_time3(8)

      DO i=1,1

      CALL value_thermal()

      END DO

      CALL date_and_time(date,time,zone,value_time4)
      WRITE(90,'(5I6)')value_time4(5),value_time4(6),value_time4(7),value_time4(8)
!!! =================================================================

      n_time_run_total = real(nT*(n_equi_reseau1+n_equi_reseau2*n_average_thermal))&
           *(real(value_time2(5)-value_time1(5))*60.+real(value_time2(6)-value_time1(6))&
            +real(value_time2(7)-value_time1(7))/60.+real(value_time2(8)-value_time1(8))/60000.)&

            +real(nT*n_average_thermal)&
           *(real(value_time4(5)-value_time3(5))*60.+real(value_time4(6)-value_time3(6))&
            +real(value_time4(7)-value_time3(7))/60.+real(value_time4(8)-value_time3(8))/60000.)

      n_date_run=int(n_time_run_total/1440.)
      n_hour_run=int(n_time_run_total/60.-n_date_run*24.)
      n_minute_run=int(n_time_run_total-n_date_run*1440.-n_hour_run*60.)

      WRITE(90,*)'n_time_run_total:',n_time_run_total,'mm'
      WRITE(90,*)'n_date_run      :',n_date_run
      WRITE(90,*)'n_hour_run      :',n_hour_run
      WRITE(90,*)'n_minute_run    :',n_minute_run

      CLOSE(90)

      END SUBROUTINE calcul_time_run

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!! Initial position configuration
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE load_config_ini()
      IMPLICIT NONE

!!!====================================================
!!! Hinh dang ban dau theo truc x

      IF (CONFIG_INI_DIREC_X == 'YES') THEN

      DO k=1,natz
            DO j=1,naty
                  DO i=1,natx
                  spin(i,j,k,1)=1.
                  spin(i,j,k,2)=0.
                  spin(i,j,k,3)=0.

                  ENDDO
            ENDDO
      ENDDO

      END IF

!!!====================================================
!!! Hinh dang ban dau bat ki

      IF (CONFIG_INI_DIREC_X == 'NO') THEN

      DO k=1,natz
            DO j=1,naty
                  DO i=1,natx

                  CALL random_number(rdn_config)

                  IF (rdn_config<1./3.) THEN
                        spin(i,j,k,1)=1.
                        spin(i,j,k,2)=0.
                        spin(i,j,k,3)=0.
                  ELSE
                        IF (rdn_config<2./3.) THEN
                        spin(i,j,k,1)=0.
                        spin(i,j,k,2)=1.
                        spin(i,j,k,3)=0.

                        ELSE
                        spin(i,j,k,1)=0.
                        spin(i,j,k,2)=0.
                        spin(i,j,k,3)=1.
                        END IF

                  END IF

                  ENDDO
            ENDDO
      ENDDO

      END IF

!!!====================================================
!!! Hinh dang ban dau theo GS
      IF (CONFIG_INI_DIREC_X == 'GS') THEN

      spin(:,:,:,:)=0.

      DO i_layer=1,n_layer
            DO m=1,d
                  k1 = (i_layer-1)*2*d+m
                  k2= (i_layer-1)*2*d+d+m
                  spin(:,:,k1,1)=1.
                  spin(:,:,k2,2)=1.

            END DO
      END DO

      END IF

      END SUBROUTINE load_config_ini

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!! WRITE initial position configuration in 3D
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE write_config_ini_3D()
      IMPLICIT NONE

      OPEN(unit=12,file='config_ini_3D/config_ini_3D_POTTS.pdb')

      DO k=1,natz
            DO j=1,naty
                  DO i=1,natx

                        x=real(i-1)*aaa
                        y=real(j-1)*aaa
                        z=real(k-1)*aaa

                        IF (int(spin(i,j,k,1))==1) THEN
                              spinAt='Au'
                              WRITE(12,'("ATOM",9X,A3,14X,3F8.3,1X,"1.00",1X,F8.3)') &
                                   spinAt,x,y,z,nul

                        ELSE

                              IF (int(spin(i,j,k,2))==1) THEN
                              spinAt='Cu'
                              WRITE(12,'("ATOM",9X,A3,14X,3F8.3,1X,"1.00",1X,F8.3)') &
                                  spinAt,x,y,z,nul

                              ELSE

                              spinAt='Mn'
                              WRITE(12,'("ATOM",9X,A3,14X,3F8.3,1X,"1.00",1X,F8.3)') &
                                   spinAt,x,y,z,nul

                              END IF
                        END IF

                  END DO
            END DO
      END DO

      CLOSE(12)

      END SUBROUTINE write_config_ini_3D

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!! WRITE initial position configuration in 3D
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE write_config_3D()
      IMPLICIT NONE

      number=iT+10000000
      WRITE(tmp,'(I8)') number

      name='_POTTS_'//TRIM(tmp)

      OPEN(unit=13,file='config_3D/config_3D'//trim(name)//'.pdb')

      DO k=1,natz
            DO j=1,naty
                  DO i=1,natx

                        x=real(i-1)*aaa
                        y=real(j-1)*aaa
                        z=real(k-1)*aaa

                       IF (int(spin(i,j,k,1))==1) THEN
                              spinAt='Au'
                              WRITE(13,'("ATOM",9X,A3,14X,3F8.3,1X,"1.00",1X,F8.3)') &
                                   spinAt,x,y,z,nul

                        ELSE

                              IF (int(spin(i,j,k,2))==1) THEN
                              spinAt='Cu'
                              WRITE(13,'("ATOM",9X,A3,14X,3F8.3,1X,"1.00",1X,F8.3)') &
                                  spinAt,x,y,z,nul

                              ELSE

                              spinAt='Mn'
                              WRITE(13,'("ATOM",9X,A3,14X,3F8.3,1X,"1.00",1X,F8.3)') &
                                   spinAt,x,y,z,nul

                              END IF
                        END IF

                  END DO
            END DO
      END DO

      CLOSE(13)

      END SUBROUTINE write_config_3D

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE equi_reseau1()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE equi_reseau1()
      IMPLICIT NONE

      DO i_loop=1,n_equi_reseau1
            CALL equi_reseau()
      END DO

      END SUBROUTINE equi_reseau1

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE equi_reseau()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE equi_reseau()
      IMPLICIT NONE

      DO k=1,natz
            kp=k+1-(k/natz)*natz
            km=k-1+(1/k)*natz
            DO j=1,naty
                  jp=j+1-(j/naty)*naty
                  jm=j-1+(1/j)*naty
                  DO i=1,natx
                        ip=i+1-(i/natx)*natx
                        im=i-1+(1/i)*natx

                        i_case=int(spin(i,j,k,1))+2*int(spin(i,j,k,2))+3*int(spin(i,j,k,3))

                        SELECT CASE (i_case)
!!!!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!!!!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

                        !!! Truong hop spin nam doc theo truc x
                        !!=================================================================================
                        CASE(1)

                        CALL random_number(rdn_etat)
                        IF (rdn_etat<0.5) THEN

                        CALL value_EN_D_tmp12()

                        EN_tmp1=-(spin(ip,j,k,1)+spin(im,j,k,1)+spin(i,jp,k,1)&
                                 +spin(i,jm,k,1)+spin(i,j,kp,1)+spin(i,j,km,1))+D_tmp*EN_D_tmpx

                        EN_tmp2=-(spin(ip,j,k,2)+spin(im,j,k,2)+spin(i,jp,k,2)&
                                  +spin(i,jm,k,2)+spin(i,j,kp,2)+spin(i,j,km,2))+D_tmp*EN_D_tmpy

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,1)=0.
                                    spin(i,j,k,2)=1.

                              END IF

                        ELSE

                        CALL value_EN_D_tmp13()

                        EN_tmp1=-(spin(ip,j,k,1)+spin(im,j,k,1)+spin(i,jp,k,1)&
                                 +spin(i,jm,k,1)+spin(i,j,kp,1)+spin(i,j,km,1))+D_tmp*EN_D_tmpx

                        EN_tmp2=-(spin(ip,j,k,3)+spin(im,j,k,3)+spin(i,jp,k,3)&
                                 +spin(i,jm,k,3)+spin(i,j,kp,3)+spin(i,j,km,3))+D_tmp*EN_D_tmpz

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,1)=0.
                                    spin(i,j,k,3)=1.

                              END IF

                        END IF

                        !!! Truong hop spin nam doc theo truc y
                        !!=================================================================================
                        CASE(2)
                        CALL random_number(rdn_etat)

                        IF (rdn_etat<0.5) THEN
                        CALL value_EN_D_tmp12()

                        EN_tmp1=-(spin(ip,j,k,2)+spin(im,j,k,2)+spin(i,jp,k,2)&
                                 +spin(i,jm,k,2)+spin(i,j,kp,2)+spin(i,j,km,2))+D_tmp*EN_D_tmpy

                        EN_tmp2=-(spin(ip,j,k,1)+spin(im,j,k,1)+spin(i,jp,k,1)&
                                 +spin(i,jm,k,1)+spin(i,j,kp,1)+spin(i,j,km,1))+D_tmp*EN_D_tmpx

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,2)=0.
                                    spin(i,j,k,1)=1.

                              END IF

                        ELSE

                        CALL value_EN_D_tmp23()
                        EN_tmp1=-(spin(ip,j,k,2)+spin(im,j,k,2)+spin(i,jp,k,2)&
                                 +spin(i,jm,k,2)+spin(i,j,kp,2)+spin(i,j,km,2))+D_tmp*EN_D_tmpy

                        EN_tmp2=-(spin(ip,j,k,3)+spin(im,j,k,3)+spin(i,jp,k,3)&
                                 +spin(i,jm,k,3)+spin(i,j,kp,3)+spin(i,j,km,3))+D_tmp*EN_D_tmpz

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,2)=0.
                                    spin(i,j,k,3)=1.

                              END IF

                        END IF

                        !!! Truong hop spin nam doc theo truc z
                        !!=================================================================================
                        CASE(3)
                        CALL random_number(rdn_etat)
                        IF (rdn_etat<0.5) THEN
                        CALL value_EN_D_tmp13()

                        EN_tmp1=-(spin(ip,j,k,3)+spin(im,j,k,3)+spin(i,jp,k,3)&
                                 +spin(i,jm,k,3)+spin(i,j,kp,3)+spin(i,j,km,3))+D_tmp*EN_D_tmpz

                        EN_tmp2=-(spin(ip,j,k,1)+spin(im,j,k,1)+spin(i,jp,k,1)&
                                 +spin(i,jm,k,1)+spin(i,j,kp,1)+spin(i,j,km,1))+D_tmp*EN_D_tmpx

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,3)=0.
                                    spin(i,j,k,1)=1.

                              END IF

                        ELSE

                        CALL value_EN_D_tmp23()
                        EN_tmp1=-(spin(ip,j,k,3)+spin(im,j,k,3)+spin(i,jp,k,3)&
                                 +spin(i,jm,k,3)+spin(i,j,kp,3)+spin(i,j,km,3))+D_tmp*EN_D_tmpz
                        EN_tmp2=-(spin(ip,j,k,2)+spin(im,j,k,2)+spin(i,jp,k,2)&
                                 +spin(i,jm,k,2)+spin(i,j,kp,2)+spin(i,j,km,2))+D_tmp*EN_D_tmpy

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,3)=0.
                                    spin(i,j,k,2)=1.

                              END IF

                        END IF

                        !!=================================================================================

                        END SELECT

!!!!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!!!!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

                   ENDDO
             ENDDO
      ENDDO

      END SUBROUTINE equi_reseau

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE value_thermal()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE value_thermal()

      IMPLICIT NONE

      energy=0.

      DO k=1,natz
            kp=k+1-(k/natz)*natz
            km=k-1+(1/k)*natz
            DO j=1,naty
                  jp=j+1-(j/naty)*naty
                  jm=j-1+(1/j)*naty
                  DO i=1,natx
                        ip=i+1-(i/natx)*natx
                        im=i-1+(1/i)*natx

                        i_case=int(spin(i,j,k,1))+2*int(spin(i,j,k,2))+3*int(spin(i,j,k,3))

                        SELECT CASE (i_case)
!!!!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!!!!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

                        !!! Truong hop spin nam doc theo truc x
                        !!=================================================================================
                        CASE(1)

                        CALL random_number(rdn_etat)
                        IF (rdn_etat<0.5) THEN

                        CALL value_EN_D_tmp12()

                        EN_tmp1=-(spin(ip,j,k,1)+spin(im,j,k,1)+spin(i,jp,k,1)&
                                 +spin(i,jm,k,1)+spin(i,j,kp,1)+spin(i,j,km,1))+D_tmp*EN_D_tmpx

                        EN_tmp2=-(spin(ip,j,k,2)+spin(im,j,k,2)+spin(i,jp,k,2)&
                                  +spin(i,jm,k,2)+spin(i,j,kp,2)+spin(i,j,km,2))+D_tmp*EN_D_tmpy

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,1)=0.
                                    spin(i,j,k,2)=1.
                                    EN_tmp1=EN_tmp2
                              END IF

                        ELSE

                        CALL value_EN_D_tmp13()

                        EN_tmp1=-(spin(ip,j,k,1)+spin(im,j,k,1)+spin(i,jp,k,1)&
                                 +spin(i,jm,k,1)+spin(i,j,kp,1)+spin(i,j,km,1))+D_tmp*EN_D_tmpx

                        EN_tmp2=-(spin(ip,j,k,3)+spin(im,j,k,3)+spin(i,jp,k,3)&
                                 +spin(i,jm,k,3)+spin(i,j,kp,3)+spin(i,j,km,3))+D_tmp*EN_D_tmpz

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,1)=0.
                                    spin(i,j,k,3)=1.
                                    EN_tmp1=EN_tmp2
                              END IF

                        END IF

                        !!! Truong hop spin nam doc theo truc y
                        !!=================================================================================
                        CASE(2)
                        CALL random_number(rdn_etat)

                        IF (rdn_etat<0.5) THEN
                        CALL value_EN_D_tmp12()

                        EN_tmp1=-(spin(ip,j,k,2)+spin(im,j,k,2)+spin(i,jp,k,2)&
                                 +spin(i,jm,k,2)+spin(i,j,kp,2)+spin(i,j,km,2))+D_tmp*EN_D_tmpy

                        EN_tmp2=-(spin(ip,j,k,1)+spin(im,j,k,1)+spin(i,jp,k,1)&
                                 +spin(i,jm,k,1)+spin(i,j,kp,1)+spin(i,j,km,1))+D_tmp*EN_D_tmpx

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,2)=0.
                                    spin(i,j,k,1)=1.
                                    EN_tmp1=EN_tmp2
                              END IF

                        ELSE

                        CALL value_EN_D_tmp23()
                        EN_tmp1=-(spin(ip,j,k,2)+spin(im,j,k,2)+spin(i,jp,k,2)&
                                 +spin(i,jm,k,2)+spin(i,j,kp,2)+spin(i,j,km,2))+D_tmp*EN_D_tmpy

                        EN_tmp2=-(spin(ip,j,k,3)+spin(im,j,k,3)+spin(i,jp,k,3)&
                                 +spin(i,jm,k,3)+spin(i,j,kp,3)+spin(i,j,km,3))+D_tmp*EN_D_tmpz

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,2)=0.
                                    spin(i,j,k,3)=1.
                                    EN_tmp1=EN_tmp2
                              END IF

                        END IF

                        !!! Truong hop spin nam doc theo truc z
                        !!=================================================================================
                        CASE(3)
                        CALL random_number(rdn_etat)
                        IF (rdn_etat<0.5) THEN
                        CALL value_EN_D_tmp13()

                        EN_tmp1=-(spin(ip,j,k,3)+spin(im,j,k,3)+spin(i,jp,k,3)&
                                 +spin(i,jm,k,3)+spin(i,j,kp,3)+spin(i,j,km,3))+D_tmp*EN_D_tmpz

                        EN_tmp2=-(spin(ip,j,k,1)+spin(im,j,k,1)+spin(i,jp,k,1)&
                                 +spin(i,jm,k,1)+spin(i,j,kp,1)+spin(i,j,km,1))+D_tmp*EN_D_tmpx

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,3)=0.
                                    spin(i,j,k,1)=1.
                                    EN_tmp1=EN_tmp2
                              END IF

                        ELSE

                        CALL value_EN_D_tmp23()
                        EN_tmp1=-(spin(ip,j,k,3)+spin(im,j,k,3)+spin(i,jp,k,3)&
                                 +spin(i,jm,k,3)+spin(i,j,kp,3)+spin(i,j,km,3))+D_tmp*EN_D_tmpz
                        EN_tmp2=-(spin(ip,j,k,2)+spin(im,j,k,2)+spin(i,jp,k,2)&
                                 +spin(i,jm,k,2)+spin(i,j,kp,2)+spin(i,j,km,2))+D_tmp*EN_D_tmpy

                              CALL random_number(rdn_mtp)

                              IF (exp(-(EN_tmp2-EN_tmp1)/T) > rdn_mtp) THEN
                                    spin(i,j,k,3)=0.
                                    spin(i,j,k,2)=1.
                                    EN_tmp1=EN_tmp2
                              END IF

                        END IF

                        !!=================================================================================

                        END SELECT

                        energy=energy+EN_tmp1

!!!!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!!!!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

                   ENDDO
             ENDDO
      ENDDO

      CALL value_Mz()

      energy=energy/(2.*n_atom)
      energy_2=energy**2.
      Mz_2=Mz**2.

      !WRITE(*,*)'EN=',energy,'Mz=',Mz,'Mz1=',Mz1,'Mz2=',Mz2,'Mz3=',Mz3,'Mz1+Mz2+Mz3=',Mz1+Mz2+Mz3

      END SUBROUTINE value_thermal

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE value_Mz()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE value_Mz()
      IMPLICIT NONE

      !!! Tinh theo cac mat phang xy
      M12=0.
      DO k= 1,natz

            M12x=0.; M12y=0.; M12z=0.

            DO j=1,naty
                  DO i=1,natx
                  M12x=M12x+spin(i,j,k,1)
                  M12y=M12y+spin(i,j,k,2)
                  M12z=M12z+spin(i,j,k,3)
                  END DO
            END DO

            M12=M12+abs(3.*max(M12x,M12y,M12z)/real(natx*naty)-1.)/2.
      END DO

      M12=M12/real(natz)

      !!! Tinh theo cac mat phang xz
      M13=0.
      DO j= 1,naty

            M13x=0.; M13y=0.; M13z=0.

            DO k=1,natz
                  DO i=1,natx
                  M13x=M13x+spin(i,j,k,1)
                  M13y=M13y+spin(i,j,k,2)
                  M13z=M13z+spin(i,j,k,3)
                  END DO
            END DO

            M13=M13+abs(3.*max(M13x,M13y,M13z)/real(natx*natz)-1.)/2.
      END DO

      M13=M13/real(naty)

      !!! Tinh theo cac mat phang yz
      M23=0.
      DO i= 1,natx

            M23x=0.; M23y=0.; M23z=0.

            DO k=1,natz
                  DO j=1,naty
                  M23x=M23x+spin(i,j,k,1)
                  M23y=M23y+spin(i,j,k,2)
                  M23z=M23z+spin(i,j,k,3)
                  END DO
            END DO

            M23=M23+abs(3.*max(M23x,M23y,M23z)/real(naty*natz)-1.)/2.
      END DO

      M23=M23/real(natx)

      Mz=max(M12,M13,M23)

      END SUBROUTINE value_Mz
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE average_thermal()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE average_thermal()

      IMPLICIT NONE

      EN_moy=0.
      EN_2_moy=0.
      Mz_moy=0.
      Mz_2_moy=0.

      !!!------------------------------------------
      OPEN(unit=126,file='value_P_at_To.dat')

      H_histo(:)=0.
      EN_histo(:)=0.

      del_EN_histo = (EN_max_histo-EN_min_histo)/real(n_EN_histo-1)

      !! END khai bao bien cho histogram

      DO i_loop=1,n_average_thermal

            DO i_loop2=1,n_equi_reseau2
                  CALL equi_reseau()
            END DO

            CALL value_thermal()

            EN_moy=EN_moy+energy
            EN_2_moy=EN_2_moy+energy_2
            Mz_moy=Mz_moy+Mz
            Mz_2_moy=Mz_2_moy+Mz_2

            IF ((EN_min_histo <= energy) .and. (energy <= EN_max_histo)) THEN
                  i_histo=int((energy-EN_min_histo)/del_EN_histo)+1
                  H_histo(i_histo)=H_histo(i_histo)+1.

            END IF

      END DO

      DO i_histo=1,n_EN_histo
            P_histo(i_histo)=H_histo(i_histo)/real(n_average_thermal)
            EN_histo(i_histo)=EN_min_histo+(real(i_histo)-0.5)*del_EN_histo

            WRITE(Ligne126,*)T,EN_histo(i_histo),P_histo(i_histo)
            WRITE(126,'(a)') trim(Ligne126)
      ENDDO

      CLOSE(126)

      EN_moy=EN_moy/real(n_average_thermal)
      EN_2_moy=EN_2_moy/real(n_average_thermal)
      Mz_moy=Mz_moy/real(n_average_thermal)
      Mz_2_moy=Mz_2_moy/real(n_average_thermal)

      Cv = n_atom*(EN_2_moy-EN_moy**2.)/(T**2.)
      Ksi= n_atom*(Mz_2_moy-Mz_moy**2.)/T

      WRITE(Ligne20,*) T,EN_moy,Mz_moy,Cv,Ksi,EN_2_moy,Mz_2_moy

      WRITE(20,'(a)') trim(Ligne20)

      WRITE(*,*)'i_times=',i_times,'EN_moy=',EN_moy,'Mz_moy=',Mz_moy,'Cv=',Cv,'Ksi=',Ksi

      END SUBROUTINE average_thermal

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE value_EN_D_tmp12()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE value_EN_D_tmp12()
      IMPLICIT NONE

      EN_D_tmpx=0.
      EN_D_tmpy=0.

      DO iv=1,n_vois

            i1=i+int_deli(iv)
            j1=j+int_delj(iv)
            k1=k+int_delk(iv)

            IF (i1<1)    i1=natx+i1
            IF (i1>natx) i1=i1-natx
            IF (j1<1)    j1=naty+j1
            IF (j1>naty) j1=j1-naty
            IF (k1<1)    k1=natz+k1
            IF (k1>natz) k1=k1-natz

            EN_D_tmpx=EN_D_tmpx+spin(i1,j1,k1,1)*(1./r(iv)**3.-3.*(deli(iv)**2.)/(r(iv)**5.))
            EN_D_tmpy=EN_D_tmpy+spin(i1,j1,k1,2)*(1./r(iv)**3.-3.*(delj(iv)**2.)/(r(iv)**5.))

      END DO

      END SUBROUTINE value_EN_D_tmp12

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE value_EN_D_tmp13()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE value_EN_D_tmp13()
      IMPLICIT NONE

      EN_D_tmpx=0.
      EN_D_tmpz=0.

      DO iv=1,n_vois

            i1=i+int_deli(iv)
            j1=j+int_delj(iv)
            k1=k+int_delk(iv)

            IF (i1<1)    i1=natx+i1
            IF (i1>natx) i1=i1-natx
            IF (j1<1)    j1=naty+j1
            IF (j1>naty) j1=j1-naty
            IF (k1<1)    k1=natz+k1
            IF (k1>natz) k1=k1-natz

            EN_D_tmpx=EN_D_tmpx+spin(i1,j1,k1,1)*(1./r(iv)**3.-3.*(deli(iv)**2.)/(r(iv)**5.))
            EN_D_tmpz=EN_D_tmpz+spin(i1,j1,k1,3)*(1./r(iv)**3.-3.*(delk(iv)**2.)/(r(iv)**5.))

      END DO

      END SUBROUTINE value_EN_D_tmp13

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE value_EN_D_tmp23()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE value_EN_D_tmp23()
      IMPLICIT NONE

      EN_D_tmpy=0.
      EN_D_tmpz=0.

      DO iv=1,n_vois

            i1=i+int_deli(iv)
            j1=j+int_delj(iv)
            k1=k+int_delk(iv)

            IF (i1<1)    i1=natx+i1
            IF (i1>natx) i1=i1-natx
            IF (j1<1)    j1=naty+j1
            IF (j1>naty) j1=j1-naty
            IF (k1<1)    k1=natz+k1
            IF (k1>natz) k1=k1-natz

            EN_D_tmpy=EN_D_tmpy+spin(i1,j1,k1,2)*(1./r(iv)**3.-3.*(delj(iv)**2.)/(r(iv)**5.))
            EN_D_tmpz=EN_D_tmpz+spin(i1,j1,k1,3)*(1./r(iv)**3.-3.*(delk(iv)**2.)/(r(iv)**5.))

      END DO

      !WRITE(*,*)'ENDx=',EN_D_tmpx,'ENDy=',EN_D_tmpy,'ENDz=',EN_D_tmpz

      END SUBROUTINE value_EN_D_tmp23

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE value_EN_D_tmp()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE value_EN_D_tmp()
      IMPLICIT NONE

      EN_D_tmpx=0.
      EN_D_tmpy=0.
      EN_D_tmpz=0.

      DO iv=1,n_vois

            i1=i+int_deli(iv)
            j1=j+int_delj(iv)
            k1=k+int_delk(iv)

            IF (i1<1)    i1=natx+i1
            IF (i1>natx) i1=i1-natx
            IF (j1<1)    j1=naty+j1
            IF (j1>naty) j1=j1-naty
            IF (k1<1)    k1=natz+k1
            IF (k1>natz) k1=k1-natz

            EN_D_tmpx=EN_D_tmpx+spin(i1,j1,k1,1)*(1./r(iv)**3.-3.*(deli(iv)**2.)/(r(iv)**5.))
            EN_D_tmpy=EN_D_tmpy+spin(i1,j1,k1,2)*(1./r(iv)**3.-3.*(delj(iv)**2.)/(r(iv)**5.))
            EN_D_tmpz=EN_D_tmpz+spin(i1,j1,k1,3)*(1./r(iv)**3.-3.*(delk(iv)**2.)/(r(iv)**5.))

      END DO

      END SUBROUTINE value_EN_D_tmp

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE number_voisin()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE number_voisin()

      IMPLICIT NONE
!!!! Tinh so voisin (n_vois) cua moi atom, de su dung cho khai bao cac tab
      iv=0

      DO deli0=-r01,r01
      DO delj0=-r01,r01
      DO delk0=-r01,r01

            r1=sqrt(real(deli0)**2.+real(delj0)**2.+real(delk0)**2.)

            IF ((0.<r1).and.(r1<=r0)) THEN

                  iv=iv+1
            END IF

      END DO
      END DO
      END DO

      n_vois=iv

!!!===============================================================================

!!! Tinh deli, delj, delk, ban kinh r (co gia tri bang nhau cho cac nguyen tu)

      iv=0
      ALLOCATE(deli(n_vois),delj(n_vois),delk(n_vois))
      ALLOCATE(int_deli(n_vois),int_delj(n_vois),int_delk(n_vois))
      ALLOCATE(r(n_vois))

      OPEN(14,file='voisin.dat')

      DO deli0=-r01,r01
      DO delj0=-r01,r01
      DO delk0=-r01,r01

            r1=sqrt(real(deli0)**2.+real(delj0)**2.+real(delk0)**2.)

            IF ((0.<r1).and.(r1<=r0)) THEN

                  iv=iv+1

                  deli(iv)=real(deli0)
                  delj(iv)=real(delj0)
                  delk(iv)=real(delk0)
                  r(iv)=r1

                  int_deli(iv)=deli0
                  int_delj(iv)=delj0
                  int_delk(iv)=delk0

            WRITE(14,*)iv,deli(iv),delj(iv),delk(iv),r(iv)

            END IF

      END DO
      END DO
      END DO

      WRITE(*,*)'Xong buoc tinh del,r,i1,j1,k1'

      CLOSE(14)

      END SUBROUTINE number_voisin

!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!   SUBROUTINE energy_gs()
!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE energy_gs()
      IMPLICIT NONE

      REAL (KIND=8) :: EN_gs

      DO k=1,natz
            kp=k+1-(k/natz)*natz
            km=k-1+(1/k)*natz
            DO j=1,naty
                  jp=j+1-(j/naty)*naty
                  jm=j-1+(1/j)*naty
                  DO i=1,natx
                        ip=i+1-(i/natx)*natx
                        im=i-1+(1/i)*natx

                        CALL value_EN_D_tmp12()

                        i_case=int(spin(i,j,k,1))+2*int(spin(i,j,k,2))

                        SELECT CASE (i_case)

                        CASE(1)

                        EN_tmp1=EN_tmp1+spin(i,j,k,1)*&
                                    (-(spin(ip,j,k,1)+spin(im,j,k,1)+spin(i,jp,k,1)&
                                 +spin(i,jm,k,1)+spin(i,j,kp,1)+spin(i,j,km,1))+D_tmp*EN_D_tmpx)

                        CASE(2)

                        EN_tmp2=EN_tmp2+spin(i,j,k,2)*&
                                    (-(spin(ip,j,k,2)+spin(im,j,k,2)+spin(i,jp,k,2)&
                                  +spin(i,jm,k,2)+spin(i,j,kp,2)+spin(i,j,km,2))+D_tmp*EN_D_tmpy)

                        END SELECT

                  END DO
            END DO

      END DO

      EN_gs = (EN_tmp1+EN_tmp2)/(2.*n_atom)

      WRITE(*,*)'EN_gs=',EN_gs

      END SUBROUTINE energy_gs

      END PROGRAM main_transport_modern_potts

max_stars_repo_name: danhtaihoang/molecular-crystals

max_issues_repo_name: danhtaihoang/molecular-crystals

max_forks_repo_name: danhtaihoang/molecular-crystals


text:
module initialize_module

  use BoxLib
  use ml_boxarray_module
  use ml_layout_module
  use define_bc_module
  use restart_module
  use init_module
  use box_util_module
  use make_new_grids_module
  use tag_boxes_module
  use multifab_fill_ghost_module

  implicit none

  private

  public :: initialize_from_restart, initialize_with_fixed_grids, initialize_with_adaptive_grids, &
            make_new_state

contains

  subroutine initialize_from_restart(mla,restart,time,dt,dx,pmask,uold,sold,gp,p,the_bc_tower)

     use probin_module, only : dim_in, nlevs, nscal, ng_cell, ng_grow, nodal

     type(ml_layout),intent(out)   :: mla
     integer       , intent(in   ) :: restart
     real(dp_t)    , intent(  out) :: time,dt
     real(dp_t)    , pointer       :: dx(:,:)
     logical       , intent(in   ) :: pmask(:)
     type(multifab), pointer       :: sold(:),uold(:),gp(:),p(:)
     type(bc_tower), intent(  out) :: the_bc_tower

     type(ml_boxarray)         :: mba
     type(multifab), pointer   :: chkdata(:)
     type(multifab), pointer   :: chk_p(:)
     type(layout)              :: la

     integer :: n,dm

     type(bl_prof_timer), save :: bpt

     call build(bpt,"initialize_from_restart")

     dm = dim_in

     call fill_restart_data(restart,mba,chkdata,chk_p,time,dt)

     call ml_layout_build(mla,mba,pmask)

     nlevs = mba%nlevel

     allocate(uold(nlevs),sold(nlevs),gp(nlevs),p(nlevs))

     do n = 1,nlevs
        call multifab_build(   uold(n), mla%la(n),    dm, ng_cell)
        call multifab_build(   sold(n), mla%la(n), nscal, ng_cell)
        call multifab_build(     gp(n), mla%la(n),    dm, ng_grow)
        call multifab_build(      p(n), mla%la(n),     1, ng_grow, nodal)
     end do
     do n = 1,nlevs
        call multifab_copy_c(uold(n),1,chkdata(n),1         ,dm)
        call multifab_copy_c(sold(n),1,chkdata(n),1+dm      ,nscal)
        call multifab_copy_c(  gp(n),1,chkdata(n),1+dm+nscal,dm)
        call multifab_copy_c(   p(n),1,  chk_p(n),1         ,1)
        !
        ! The layouts for chkdata and chk_p are built standalone, level
        ! by level, and need to be destroy()d as such as well.
        !
        la = get_layout(chkdata(n))
        call multifab_destroy(chkdata(n))
        call destroy(la)
        la = get_layout(chk_p(n))
        call multifab_destroy(chk_p(n))
        call destroy(la)
     end do
     deallocate(chkdata,chk_p)

     call initialize_dx(dx,mba,nlevs)

     call initialize_bc(the_bc_tower,nlevs,pmask)
     do n = 1,nlevs
        call bc_tower_level_build( the_bc_tower,n,mla%la(n))
     end do

     call destroy(mba)

     call destroy(bpt)

  end subroutine initialize_from_restart

  subroutine initialize_with_fixed_grids(mla,pmask,dx,uold,sold,gp,p,the_bc_tower)

     use probin_module, only : dim_in, nlevs, nscal, ng_cell, ng_grow, nodal, fixed_grids

     type(ml_layout),intent(out)   :: mla
     logical       , intent(in   ) :: pmask(:)
     real(dp_t)    , pointer       :: dx(:,:)
     type(multifab), pointer       :: uold(:),sold(:),gp(:),p(:)
     type(bc_tower), intent(  out) :: the_bc_tower

     type(ml_boxarray)         :: mba

     integer :: n,dm

     type(bl_prof_timer), save :: bpt

     call build(bpt,"initialize_with_fixed_grids")

     dm = dim_in

     call read_a_hgproj_grid(mba, fixed_grids)
     call ml_layout_build(mla,mba,pmask)

     ! check for proper nesting
     if (.not. ml_boxarray_properly_nested(mla%mba, ng_cell, pmask)) &
         call bl_error('fixed_grids not properly nested')

     nlevs = mla%nlevel
     allocate(uold(nlevs),sold(nlevs),p(nlevs),gp(nlevs))

     do n = 1,nlevs
        call multifab_build(   uold(n), mla%la(n),    dm, ng_cell)
        call multifab_build(   sold(n), mla%la(n), nscal, ng_cell)
        call multifab_build(     gp(n), mla%la(n),    dm, ng_grow)
        call multifab_build(      p(n), mla%la(n),     1, ng_grow, nodal)
     end do

     ! Set initial pressure to zero just in case ... to make sure ghost cell
     !     values are initialized
     do n = 1,nlevs
        call setval(p(n), 0.d0, all=.true.)
     end do

     call initialize_dx(dx,mba,nlevs)

     call initialize_bc(the_bc_tower,nlevs,pmask)

     do n = 1,nlevs
        call bc_tower_level_build( the_bc_tower,n,mla%la(n))
     end do

     call initdata(nlevs,uold,sold,dx,the_bc_tower%bc_tower_array,mla)

     call destroy(mba)

     call destroy(bpt)

  end subroutine initialize_with_fixed_grids

  subroutine initialize_with_adaptive_grids(mla,pmask,dx,uold,sold,gp,p,the_bc_tower)

     use probin_module, only : dim_in, nlevs, nscal, ng_cell, ng_grow, nodal, &
                               n_cellx, n_celly, n_cellz, &
                               regrid_int, amr_buf_width, max_grid_size, &
                               ref_ratio, max_levs, verbose

     type(ml_layout),intent(inout)  :: mla
     logical       , intent(in   )  :: pmask(:)
     real(dp_t)    , pointer        :: dx(:,:)
     type(multifab), pointer        :: uold(:),sold(:),gp(:),p(:)
     type(bc_tower), intent(  out)  :: the_bc_tower

     type(layout)                   :: la_array(max_levs)
     type(box)                      :: bxs
     type(ml_boxarray)              :: mba

     logical :: new_grid
     integer :: lo(dim_in), hi(dim_in)
     integer :: n, nl, dm, ng_buffer

     type(bl_prof_timer), save :: bpt

     call build(bpt,"initialize_with_adaptive_grids")

     dm = dim_in

     ! set up hi & lo to carry indexing info
     lo = 0
     hi(1) = n_cellx-1
     if (dm > 1) then
        hi(2) = n_celly - 1
        if (dm > 2)  then
           hi(3) = n_cellz -1
        endif
     endif

     ! mba is big enough to hold max_levs levels
     call ml_boxarray_build_n(mba,max_levs,dm)
     do n = 1, max_levs-1
        mba%rr(n,:) = ref_ratio
     enddo

     allocate(uold(max_levs),sold(max_levs),p(max_levs),gp(max_levs))

       ! Build the level 1 boxarray
     call box_build_2(bxs,lo,hi)
     call boxarray_build_bx(mba%bas(1),bxs)
     call boxarray_maxsize(mba%bas(1),max_grid_size)

     ! build pd(:)
     mba%pd(1) = bxs
     do n = 2, max_levs
        mba%pd(n) = refine(mba%pd(n-1),mba%rr((n-1),:))
     enddo

     ! Need to build pd before making dx
     call initialize_dx(dx,mba,max_levs)

     ! Initialize bc's.
     call initialize_bc(the_bc_tower,max_levs,pmask)

     ! Build the level 1 layout.
     call layout_build_ba(la_array(1),mba%bas(1),mba%pd(1),pmask)

     ! Define bc_tower at level 1.
     call bc_tower_level_build(the_bc_tower,1,la_array(1))

     nlevs = 1

     if (max_levs > 1) then

        ! Build and initialize the level 1 data only.
        ! We're going to use it to tag and refine boxes.
        ! Note: fill_boundary() and physbc() are called in initdata_on_level
        call make_new_state(la_array(1),uold(1),sold(1),gp(1),p(1))
        call initdata_on_level(uold(1),sold(1),dx(1,:),the_bc_tower%bc_tower_array(1))

        new_grid = .true.
        nl = 1

        ! Choose 4 because it can accommodate the
        ! largest stencil currently used in VARDEN
        ng_buffer = 4

        do while ( (nl .lt. max_levs) .and. (new_grid) )

           if (tagging_needs_ghost_cells .and. nl.ge.2) then
              call multifab_fill_ghost_cells(sold(nl),sold(nl-1),sold(nl)%ng,mba%rr((nl-1),:), &
                   the_bc_tower%bc_tower_array(nl-1), &
                   the_bc_tower%bc_tower_array(nl), &
                   1,dm+1,nscal)
           end if

           ! Do we need finer grids?
           call make_new_grids(new_grid,la_array(nl),la_array(nl+1),sold(nl),dx(nl,1), &
                               amr_buf_width,ref_ratio,nl,max_grid_size)

           if (new_grid) then

              call copy(mba%bas(nl+1),get_boxarray(la_array(nl+1)))

              ! Need to enforce proper nesting within the grid creation procedure
              ! so that we can fillpatch the new levels.
              if (nl .ge. 2) then

                 ! Test on whether grids are already properly nested
                 if (.not. ml_boxarray_properly_nested(mba, ng_buffer, pmask, 2, nl+1)) then

                    do n = 2,nl
                       ! Delete old multifabs so that we can rebuild them.
                       call destroy(  sold(n))
                       call destroy(  uold(n))
                       call destroy(    gp(n))
                       call destroy(     p(n))
                    end do

                    call enforce_proper_nesting(mba,la_array,max_grid_size)

                    ! Loop over all the lower levels which we might have changed when we enforced proper nesting.
                    do n = 2,nl

                       ! This makes sure the boundary conditions are properly defined everywhere
                       call bc_tower_level_build(the_bc_tower,n,la_array(n))

                       ! Rebuild the lower level data again if it changed.
                       call make_new_state(la_array(n),uold(n),sold(n),gp(n),p(n))
                       call initdata_on_level(uold(n),sold(n),dx(n,:),the_bc_tower%bc_tower_array(n))

                    end do
                 end if
              end if

              ! Define bc_tower at level nl+1.
              call bc_tower_level_build(the_bc_tower,nl+1,la_array(nl+1))

              ! Build the level nl+1 data only.
              call make_new_state(la_array(nl+1),uold(nl+1),sold(nl+1),gp(nl+1),p(nl+1))

              ! fills the physical region of nl+1 with problem data
              call initdata_on_level(uold(nl+1),sold(nl+1),dx(nl+1,:),the_bc_tower%bc_tower_array(nl+1))

              nlevs = nl+1
              nl = nl + 1

           endif ! if (new_grid)

        enddo

      do n = 1,nlevs
         call destroy(sold(n))
         call destroy(uold(n))
         call destroy(gp(n))
         call destroy(p(n))
      end do

      nlevs = nl

      if (nlevs .ge. 3) then
          ! check for proper nesting
          call enforce_proper_nesting(mba,la_array,max_grid_size)
      end if

   end if ! end if (maxlev > 1)

   call ml_layout_restricted_build(mla,mba,nlevs,pmask)

  ! this makes sure the boundary conditions are properly defined everywhere
  do n = 1,nlevs
     call bc_tower_level_build(the_bc_tower,n,mla%la(n))
  end do

   nlevs = mla%nlevel

   do n = 1, nlevs
      call destroy(la_array(n))
   end do

   do n = 1,nlevs
      call make_new_state(mla%la(n),uold(n),sold(n),gp(n),p(n))
      call setval(gp(n),0.d0,all=.true.)
      call setval( p(n),0.d0,all=.true.)
   end do

   call initdata(nlevs,uold,sold,dx,the_bc_tower%bc_tower_array,mla)

   call destroy(mba)

   call destroy(bpt)

  end subroutine initialize_with_adaptive_grids

  subroutine make_new_state(la_loc,uold_loc,sold_loc,gp_loc,p_loc)

    use probin_module, only : dim_in, nscal, ng_cell, ng_grow, nodal
    use bl_constants_module

    type(layout),intent(in   ) :: la_loc
    type(multifab ),intent(inout) :: uold_loc,sold_loc,gp_loc,p_loc

    integer :: dm

    dm = dim_in

    call multifab_build(   uold_loc, la_loc,    dm, ng_cell)
    call multifab_build(   sold_loc, la_loc, nscal, ng_cell)
    call multifab_build(     gp_loc, la_loc,    dm, ng_grow)
    call multifab_build(      p_loc, la_loc,     1, ng_grow, nodal)

    call setval(  uold_loc,ZERO, all=.true.)
    call setval(  sold_loc,ZERO, all=.true.)
    call setval(    gp_loc,ZERO, all=.true.)
    call setval(     p_loc,ZERO, all=.true.)

  end subroutine make_new_state

  subroutine initialize_bc(the_bc_tower,num_levs,pmask)

     use bc_module

     use probin_module, only : dim_in, &
                               bcx_lo, bcx_hi, bcy_lo, bcy_hi, bcz_lo, bcz_hi

     type(bc_tower), intent(  out) :: the_bc_tower
     integer       , intent(in   ) :: num_levs
     logical       , intent(in   ) :: pmask(:)

     integer :: domain_phys_bc(dim_in,2)

     integer :: dm

     dm = dim_in

     ! Define the physical boundary conditions on the domain
     ! Put the bc values from the inputs file into domain_phys_bc
     domain_phys_bc(1,1) = bcx_lo
     domain_phys_bc(1,2) = bcx_hi
     if (pmask(1)) then
        domain_phys_bc(1,:) = BC_PER
        if (bcx_lo .ne. -1 .or. bcx_hi .ne. -1) &
             call bl_error('MUST HAVE BCX = -1 if PMASK = T')
     end if
     if (dm > 1) then
        domain_phys_bc(2,1) = bcy_lo
        domain_phys_bc(2,2) = bcy_hi
        if (pmask(2)) then
           domain_phys_bc(2,:) = BC_PER
           if (bcy_lo .ne. -1 .or. bcy_hi .ne. -1) &
                call bl_error('MUST HAVE BCY = -1 if PMASK = T')
        end if
     end if
     if (dm > 2) then
        domain_phys_bc(3,1) = bcz_lo
        domain_phys_bc(3,2) = bcz_hi
        if (pmask(3)) then
           domain_phys_bc(3,:) = BC_PER
           if (bcz_lo .ne. -1 .or. bcz_hi .ne. -1) &
                call bl_error('MUST HAVE BCZ = -1 if PMASK = T')
        end if
     end if

     ! Initialize the_bc_tower object.

     call bc_tower_init(the_bc_tower,num_levs,dm,domain_phys_bc)

  end subroutine initialize_bc

  subroutine initialize_dx(dx,mba,num_levs)

     use probin_module, only : dim_in, prob_lo, prob_hi

     real(dp_t)       , pointer     :: dx(:,:)
     type(ml_boxarray), intent(in ) :: mba
     integer          , intent(in ) :: num_levs

     integer :: i,n,dm

     dm = dim_in

     allocate(dx(num_levs,dm))

     do i = 1,dm
        dx(1,i) = (prob_hi(i)-prob_lo(i)) / float(mba%pd(1)%hi(i)-mba%pd(1)%lo(i)+1)
     end do
     do n = 2,num_levs
        dx(n,:) = dx(n-1,:) / mba%rr(n-1,:)
     end do

  end subroutine initialize_dx

end module initialize_module

max_stars_repo_name: BoxLib-Codes/VARDEN

max_issues_repo_name: BoxLib-Codes/VARDEN

max_forks_repo_name: BoxLib-Codes/VARDEN