text:
C
C
*+ stack_enqmode

       subroutine stack_enqmode(imap,mode,ipm,ipr,iunit,status)
C      --------------------------------------------------------
C
C Return access mode for the specified map
C
C Input:
C    Map catalogue entry
       integer        imap
C Returned
C    Access mode
       integer        mode
C    Pointer to data
       integer        ipm
C    Pointer to redtape
       integer        ipr
C    Unit number if map is open for sequential access
       integer        iunit
C    Status
       integer        status
C
C The access mode and associated information is returned for map IMAP.
C If the map is not "open" for any access then an error is returned via
C the status word (ill_access).
*-
       include 'mapcat_stack.inc'
       include 'mapcat_errors.inc'

       integer   ib, ib0

C check status on entry
       if (status.ne.0) return

C find entry in active list of maps
       ib0 = -1
       do ib=1,max_active_maps
         if (active_list(1,ib).eq.imap) then
           ib0 = ib
         end if
       end do
       if (ib0.eq.-1) then
          status = ill_active
          goto 999
       end if

       mode = active_list(2,ib0)
       if (mode.eq.mode_sequential) then
         iunit = active_list(3,ib0)
         ipm = active_list(4,ib0)
         ipr = active_list(5,ib0)
       else
         ipm = stack_status(ip_pointer,active_list(6,ib0))
         ipr = active_list(5,ib0)
       end if

999    call mapcat_err(status,'stack_enqmode',' ')

       end

ext: f

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
C ****************************************************************************
C                   VPSC 1-site multiphase  - 29/NOV/2011                    *
C ****************************************************************************
C                        BY R.LEBENSOHN & C.TOME                             *
C ****************************************************************************
C                   MULTIPHASE OR MULTIELEMENT VERSION                       *
C ****************************************************************************
C  RICARDO LEBENSOHN                 |  CARLOS TOME                          *
C  MST-8 - LANL - MS G755            |  MST-8 - LANL - MS G755               *
C  LOS ALAMOS - NM 87545 - USA       |  LOS ALAMOS - NM 87545 - USA          *
C  lebenso@lanl.gov                  |  tome@lanl.gov                        *
C ****************************************************************************
C  GENERAL REFERENCE:                                                        *
C  R.A.Lebensohn & C.N.Tome, Acta metall mater 41, 2611 (1993)               *
C ****************************************************************************
C  COPYRIGHT NOTICE                                                          *
C ****************************************************************************
C Portions of this program were prepared by the Regents of the University of *
C California at Los Alamos National Laboratory (the University) under        *
C Contract No.W-7405-ENG-36 with the U.S. Department of Energy (DOE).        *
C The University has certain rights in the program pursuant to the contract  *
C and the program should not be copied or distributed outside your           *
C organization. All rights in the program are reserved by the DOE and the    *
C University. Neither the U.S. Government nor the University makes any       *
C warranty, express or implied, or assumes any liability or responsibility   *
C for the use of this software.                                              *
C ****************************************************************************
C                              # FEATURES #                                  *
C ****************************************************************************
C  - SINGLE PHASE & MULTIPHASE AGGREGATES.                                   *
C  - SC RELATION FOR MIXED RATE-SENSITIVITY EXPONENTS.                       *
C  - ESHELBY AND PRESSURE TENSORS FOR FULLY INCOMPRESSIBLE INCLUSION & HEM.  *
C  - LOCAL CAUCHY STRESSES                                                   *
C  - MIXED BOUNDARY CONDITIONS ON DISPLACEMENT GRADIENT AND STRESS           *
C  - SYMMETRIC BASIS INTERNAL REPRESENTATION OF TENSORS                      *
C  - OPTION FOR ANALYZING ROLLING FCC COMPONENTS (icubcom=1)                 *
C  - SELF & LATENT HARDENING COEFFICIENTS                                    *
C  - 'VOCE' TYPE HARDENING MODULATION FOR EACH SYSTEM (ihardlaw=0) OR        *
C    'MTS'  HARDENING LAW FOR EACH SYSTEM (ihardlaw=1)                       *
C  - PREDOMINANT TWIN REORIENTATION SCHEME (if ihardlaw=0)                   *
C  - COMPOSITE GRAIN SCHEME (if ihardlaw=2)                                   *
C  - DIFFERENT INDIVIDUAL GRAIN SHAPES.                                      *
C  - GENERAL TREATMENT OF GRAIN SHAPE EVOLUTION.                             *
C  - INCREMENTAL CONTROL: VON MISES STRAIN, SINGLE STRAIN COMPONENT OR TIME  *
C  - POSTMORTEM CAPABILITY: RECOVERS ARRAYS FROM PREVIOUS RUN                *
C  - POLYCRYSTAL YIELD SURFACE AND LANKFORD COEFFICIENTS CALCULATION         *
C  - VARIABLE DEFORMATION HISTORY ROUTINE.                                   *
C  - OPTION FOR ROTATION COUPLING BETWEEN ORIENTATIONS                       *
C  - ALL 'CANNED' SUBROUTINES ARE FROM 'NUMERICAL RECIPIES' LIBRARY          *
C  - INTERNAL Bunge CONVENTION (phi1,Phi,phi2 -> phi,theta,omega) FOR        *
C    EULER ANGLES. BUT READS ALSO Kocks and Roe ANGLES AS INPUT.             *
C                                                                            *
C ****************************************************************************
C                         # INPUT FILES #                                    *
C ****************************************************************************
C ALWAYS NEEDED:                                                             *
C                                                                            *
C       VPSC7.IN      - general input and parameters                         *
C       "FILECRYS"    - deformation modes and CRSS (one per phase)           *
C       "FILETEXT"    - initial texture file (one per phase)                 *
C       "FILEPROC"    - deformation process to be simulated                  *
C                                                                            *
C SOMETIMES NEEDED:                                                          *
C                                                                            *
C  if ishape>1     "FILEAXES"  - initial orientation of ellipsoids           *
C                               (morphologic texture, one per phase)         *
C  if irecover=1   POSTMORT.IN - initial state from previous run             *
C  if icubcom=1    CUBCOMP.IN  - ideal orientations - fcc rolling            *
C  if ivgvar=1     "FILEHIST"  - strain history                              *
C  if ihardlaw=1   MTS parameters for the material at the end of "FILECRYS"  *
C                                                                            *
C ****************************************************************************
C                         # OUTPUT FILES #                                   *
C ****************************************************************************
C DEPENDING ON CASE AND DRIVER:                                              *
C        RERR.OUT     - convergence history                                  *
C        STATS.OUT    - statistics                                           *
C        ACT_PHn.OUT  - activity of deformation modes                        *
C        TEX_PHn.OUT  - cryst. texture for each phase/elem at dif steps      *
C        MOR_PHn.OUT  - morph. texture for each phase/elem at dif steps      *
C        STR_STR.OUT  - macroscopic stress-strain components                 *
C        RUN_LOG.OUT  - a log containing the conditions and input to the run *
C        STAT_AX.OUT  - statistics on grain axes                             *
C                                                                            *
C     if isave=1    POSTMORT.OUT - final state in grains and PX              *
C     if icauchy=1  CAUCHY.OUT - final local cauchy stresses                 *
C     if icubcom=1  CUBCOMPn.OUT - rolling components for fcc                *
C                                                                            *
C   PARAMETER & COMMON FILE (incl. during compilation): VPSC6.DIM            *
C                                                                            *
C*****************************************************************************

      PROGRAM VPSC_DELUXE_SS

      INCLUDE 'vpsc7.dim'

      CHARACTER FLABEL*2,ITERLBL*40
      DIMENSION ROTMAT(3,3)      ! used for rigid sample rotations

c      CALL CPU_TIME (START_TIME)

C ***********************************************************************
C     ASSIGNS # TO UNITS AND OPENS I/O FILES.

      UR0= 9     ! VPSC7.IN      (OPEN/CLOSE IN MAIN)
      UR1= 1     ! FILECRYS      (OPEN/READ/CLOSE INSIDE VPSC_INPUT)
      UR2= 2     ! FILETEXT      (OPEN/READ/CLOSE INSIDE VPSC_INPUT)
      UR3= 3     ! FILEAXES      (OPEN/READ/CLOSE INSIDE VPSC_INPUT)
      UR4= 4     ! POSTMORT.IN   (OPEN/CLOSE IN MAIN)
cx      UR5= 5     ! CUBCOMP.IN    (OPEN/CLOSE IN MAIN)
cx      UR6= 6     ! FILEHIST      (OPEN/CLOSE IN MAIN)
      UR5= 98     ! CUBCOMP.IN    (OPEN/CLOSE IN MAIN)
      UR6= 99    ! FILEHIST      (OPEN/CLOSE IN MAIN)

C         10     ! RUN_LOG.OUT   (OPEN IN MAIN)
C         11     ! STATS.OUT     (OPEN IN MAIN)
      UW1=12     ! RERR.OUT
C         13     ! STR_STR.OUT   (OPEN IN MAIN)
C         14     ! PCYS.OUT      (OPEN IN MAIN)
C         15     ! LANKFORD.OUT  (OPEN IN MAIN)
C         16     ! STAT_AX.OUT   (OPEN IN MAIN)
      UW2=19     ! POSTMORT.OUT  (OPEN/CLOSE IN MAIN)
C         83     ! FLUCT.OUT     (OPEN IN VPSC_INPUT IF IFLU=1)
C         84     ! FLCUB.OUT     (OPEN IN VPSC_INPUT IF INTERACT=5)
C         97     ! SO.OUT        (OPEN IN VPSC_INPUT IF INTERACT=5)

      OPEN(10, FILE='RUN_LOG.OUT',STATUS='UNKNOWN')
      OPEN(11, FILE='STATS.OUT'  ,STATUS='UNKNOWN')
      OPEN(UW1,FILE='RERR.OUT'   ,STATUS='UNKNOWN')

C ***************************************************************************
C     CALL SUBROUTINE VPSC_INPUT FOR READING CRYSTAL, GRAIN AND TEXTURE DATA.
C     INITIALIZES ARRAYS.

      OPEN(UR0,FILE='vpsc7.in',STATUS='OLD')

      CALL VPSC_INPUT

      CALL ELSC (0)     ! SC INITIAL ELASTIC MODULI

C ***********************************************************************
C     READS SX & PX STATE FROM PREVIOUS RUN (sav,xmsec,sg,ag,crss,..)
C     WHEN IRECOVER=1.

      IF(IRECOVER.EQ.1) THEN
        OPEN(UR4,file='POSTMORT.IN',form='UNFORMATTED',
     #           access='SEQUENTIAL',status='OLD')
          CALL POSTMORTEM (1)
        CLOSE(UNIT=UR4)
      ENDIF

C *** SET TO ZERO ACCUMULATED STRAIN COMPONENTS
      EPSACU=0.
      EPSVM =0.
      DO I=1,3
        DO J=1,3
          EPSTOT(I,J)=0.
        ENDDO
      ENDDO

      IF(ICUBCOM.EQ.1) THEN
        OPEN(UR5,file='CUBCOMP.IN',status='OLD')
          CALL CUBCOMP(0,0)          ! reads ideal rolling components
        CLOSE(UNIT=UR5)
      ENDIF

C ***********************************************************************
C     OPEN OUTPUT FILES
C ***********************************************************************

      OPEN(13,file='STR_STR.OUT' ,status='UNKNOWN')     ! may be redundant

      IF(ICAUCHY.EQ.1) OPEN(18,file='CAUCHY.OUT',status='UNKNOWN')

      IF(NELEM.EQ.1) THEN
        NFILES=NPH
        FLABEL='PH'
      ELSE IF(NELEM.GT.1) THEN
        NFILES=NELEM
        FLABEL='EL'
      ENDIF

      DO I=1,NFILES
        IF(IVGVAR.LE.1) THEN     ! may be redundant
          IUNIT=30+I
          OPEN(IUNIT,FILE='TEX_'//FLABEL//CHAR(48+I)//'.OUT',
     #        STATUS='UNKNOWN')
        ENDIF

        IF(ISHAPE(I).NE.0) THEN
          OPEN(16,FILE='STAT_AX.OUT',STATUS='UNKNOWN')
          IUNIT=40+I
          OPEN(IUNIT,FILE='MOR_'//FLABEL//CHAR(48+I)//'.OUT',
     #        STATUS='UNKNOWN')
        ENDIF

        IF(ICUBCOM.EQ.1) THEN
          IUNIT=60+I
          OPEN(IUNIT,FILE='CUBCOMP'//CHAR(48+I)//'.OUT',
     #        STATUS='UNKNOWN')
cwx        CALL CUBCOMP (0,1)         ! cubcomp analysis for initial texture
        ENDIF
c
          IUNIT=50+I
          OPEN(IUNIT,FILE='ACT_'//FLABEL//CHAR(48+I)//'.OUT',
     #        STATUS='UNKNOWN')
      ENDDO
c
c     cubcomp analysis for initial texture
c
      IF(ICUBCOM.EQ.1) CALL CUBCOMP (0,1)

C *******************************************************************
C     DO LOOP OVER PROCESSES
C *******************************************************************

C *** READS NUMBER OF PROCESSES TO RUN SEQUENTIALLY. IT CAN BE ANY
C     COMBINATION OF UNIFORM LOAD (IVGVAR=0), VARIABLE LOAD (IVGVAR=1),
C     PCYS (IVGVAR=2), LANKFORD (IVGVAR=3), RIGID ROTATION (IVGVAR=4)

      READ(UR0,'(A)') PROSA
      READ(UR0,*)     NPROC
      READ(UR0,'(A)') PROSA

      DO 3000 IPROC=1,NPROC

C *** READS LOAD CONDITIONS: STRESS, VELOCITY GRADIENT, TEMP, INCR, NSTEPS.
C *** CALCULATES SYMMETRIC STRAIN RATE 'DSIM' AND 'DBAR'

      READ(UR0,*) IVGVAR

      IPCYSOPEN=0
      IPCYSSKIP=0
      ILANKOPEN=0

      IF(IVGVAR.EQ.0) THEN
        READ(UR0,'(A)') FILEHIST
        OPEN(UR6,file=FILEHIST,status='OLD')
        CALL LOAD_CONDITIONS (UR6)
        CLOSE(UNIT=UR6)
        NSTEPSX=NSTEPS+1      ! EXTRA STEP REQUIRED FOR UPDATING MACRO STRESS
      ELSE IF(IVGVAR.EQ.1) THEN
        READ(UR0,'(A)') FILEHIST
        OPEN(UR6,file=FILEHIST,status='UNKNOWN')
        CALL VAR_VEL_GRAD(0)
        NSTEPSX=NSTEPS
      ELSE IF(IVGVAR.EQ.2) THEN
        IF(IPCYSOPEN.EQ.0) THEN
          OPEN(14,file='PCYS.OUT',status='UNKNOWN')
          IPCYSOPEN=1
        ENDIF
        READ(UR0,*) INDX,INDY              ! components of 2D subspace
        CALL PCYS(IDUM,INDX,INDY,0)        ! generates 2D-PCYS probes
        NSTEPSX=NSTEPS
      ELSE IF(IVGVAR.EQ.-2) THEN
        IF(IPCYSOPEN.EQ.0) THEN
          OPEN(14,file='PCYS.OUT',status='UNKNOWN')
          IPCYSOPEN=1
        ENDIF
        CALL PCYS_IT(IDUM,IPCYSSKIP,0)     ! generates 5D-PCYS probes
        NSTEPSX=NSTEPS
      ELSE IF(IVGVAR.EQ.3) THEN
        IF(ILANKOPEN.EQ.0) then
          OPEN(15,file='LANKFORD.OUT',status='UNKNOWN')
          ILANKOPEN=1
        ENDIF
        READ(UR0,*) DELTALANK              ! angular increment for probes
        CALL LANKFORD(ISTEP,DELTALANK,0)   ! initializes arrays
        NSTEPSX=NSTEPS
      ELSE IF(IVGVAR.EQ.4) THEN
        READ(UR0,'(A)') FILEHIST
        OPEN(UR6,file=FILEHIST,status='OLD')
          READ(UR6,'(A)') PROSA
          READ(UR6,*) ((ROTMAT(I,J),J=1,3),I=1,3)
        CLOSE(UNIT=UR6)
        CALL TEXTURE_ROTATION(ROTMAT)
        CALL WRITE_TEXTURE
        GO TO 3000
      ENDIF

C *******************************************************************
C     DO LOOP OVER DEFORMATION STEPS
C *******************************************************************

      IF(IVGVAR.LE.1 .AND. NSTEPS.LT.NWRITE) THEN
        WRITE(*,'(/,'' WARNING *** NWRITE='',I3,''  > NSTEPS='',I3,
     #        '' --> WILL NOT WRITE TEXTURE !!!'',/)') NWRITE,NSTEPS
        write(*,'("1")')
        read(*,*)
      ENDIF

      IF(INTERACTION.EQ.0) ITERLBL=' *** FULL CONSTRAINT CALCULATION'
      IF(INTERACTION.EQ.1) ITERLBL=' *** AFFINE CALCULATION'
      IF(INTERACTION.EQ.2) ITERLBL=' *** SECANT CALCULATION'
      IF(INTERACTION.EQ.3) ITERLBL=' *** NEFF=10 CALCULATION'
      IF(INTERACTION.EQ.4) ITERLBL=' *** TANGENT CALCULATION'
      IF(INTERACTION.EQ.5) ITERLBL=' *** SECOND ORDER CALCULATION'

      DO 2500 ISTEP=1,NSTEPSX
cw
        if(iflu.eq.1) then
          write(83,'(a,i7)') ' STEP =',istep
          if(icubcom.eq.1) write(84,'(a,i7)') ' STEP =',istep
        endif
cw
        if(interaction.eq.5) then
          write(97,'(a,i7)') ' STEP =',istep
          write(97,'(a)') '    IRS  ITSO  ERRASO      ERRESO'
        endif

        WRITE(UW1,'(/,''*******   STEP'',I4,5X,A)') ISTEP,ITERLBL
        WRITE( 10,'(/,''*******   STEP'',I4,5X,A)') ISTEP,ITERLBL
        WRITE(  *,'(/,''*******   STEP'',I4,5X,A)') ISTEP,ITERLBL
        WRITE(  *,'(''   ITSGR    SIGGR      SIGAV        DAV'',
     #             ''    ITTAN     MTAN    ITSEC     MSEC'',/)')

C *** IMPOSE VELOCITY GRADIENT AT EACH STEP:
C       IF IMPOSING MONOTONIC LOADING                 (IVGVAR=0)
C       IF IMPOSING A HISTORY OF DEFORMATION          (IVGVAR=1)
C       IF PROBING THE POLYCRYSTAL YIELD SURFACE      (IVGVAR=2)
C       IF PROBING FOR LANKFORD COEFFICIENTS          (IVGVAR=3)

        IF(IVGVAR.EQ.1) CALL VAR_VEL_GRAD(1)
        IF(IVGVAR.EQ.2) CALL PCYS(ISTEP,INDX,INDY,1)      ! 2D PCYS
        IF(IVGVAR.EQ.-2) THEN                             ! 5D PCYS
          CALL PCYS_IT(ISTEP,IPCYSSKIP,1)
          IF(IPCYSSKIP.EQ.1) GO TO 2500       ! SKIP ALREADY CALCULATED STATE
        ENDIF
        IF(IVGVAR.EQ.3) CALL LANKFORD(ISTEP,DELTALANK,1)

C *******************************************************************

        CALL UPDATE_SCHMID

C *******************************************************************
C     FOR VOCE NORMALIZES ALL ARRAYS WITH UNITS OF STRESS, STRAIN RATE &
C     COMPLIANCE TO WORK WITH n POWERS OF ORDER ONE AND IMPROVE CONVERGENCE.
C     RENORMALIZES AT THE END OF EACH STEP.
C     SCALES GAMD0 IF IRATESENS=0 TO ELIMINATE RATE SENSITIVITY ASSOCIATED .
C     WITH 'n'.

        IF(IHARDLAW.EQ.0) THEN
          IF(IRATESENS.EQ.0) CALL SCALE_3 (ISTEP)   ! rate-insensitive response
C         IF(IRATESENS.EQ.1) CALL NORMALIZE (0,ISTEP)
        ENDIF
C *******************************************************************
C     FOR MTS SCALE GAMD0 TO WORK WITH n POWERS OF ORDER ONE, AND
C     CALCULATE CRRS's ASSOCIATED WITH TEMP AND RATE FOR THIS STEP.

        IF(IHARDLAW.EQ.1) THEN
          CALL SCALE_3 (ISTEP)
          EDOT=SQRT(2./3.)*VNORM(DBAR,5)
          CALL UPDATE_CRSS_MTS (EDOT,TEMPERAT,2)
        ENDIF
C *******************************************************************
C     FOR DISL DENS LAW SCALE GAMD0 TO ELIMINATE n-ASSOCIATED RATE SENSITIVITY.
C     CALCULATE INITIAL CRRS's ASSOCIATED WITH TEMP FOR FIRST STEP.

        IF(IHARDLAW.GE.20 .AND. ISTEP.EQ.1) THEN
          CALL SCALE_3 (ISTEP)
          IF(IHARDLAW.EQ.20) CALL UPDATE_CRSS_DD (1)
          IF(IHARDLAW.EQ.22) CALL UPDATE_CRSS_DD_CE2 (1)
        ENDIF
C *******************************************************************
C     IF IMPOSING STRAIN MAKES A TAYLOR GUESS FOR THE FIRST STEP WHEN
C     IRECOVER=0 OR STARTS FROM RECOVERED STRESS STATE 'SG' WHEN IRECOVER=1.
C     IF IMPOSING STRESS MAKES A SACHS GUESS FOR EVERY STEP.

        IF(ISTEP.EQ.1.AND.IRECOVER.EQ.0 .OR. STRAIN_CONTROL.EQ.0) THEN
          CALL INITIAL_STATE_GUESS
        ENDIF
        IF(ISTEP.EQ.1.AND.IRECOVER.EQ.1) THEN
          KGX=1
          DO IPH=IPHBOT,IPHTOP
            IPHEL=IPH-IPHBOT+1
          DO KKK=NGR(IPH-1)+1,NGR(IPH)
            CALL GRAIN_RATE_AND_MODULI (1,1,KGX,KKK,IPHEL,IPH)
            KGX=KGX+1
          ENDDO
          ENDDO
        ENDIF

C *************************************************************************
C     VPSC CALCULATION (OR TAYLOR CALCULATION WHEN INTERACTION=0) OF
C     STRESS AND STRAIN-RATE FOR EVERY GRAIN.

        CALL VPSC (ISTEP)      !!!     THIS IS THE CORE OF THE CODE

C ***************************************************************************
C     VPSC PROVIDES THE MACROSCOPIC STRAIN RATE 'DSIM'.
C     CALCULATES MACROSCOPIC ROTATION RATE AND MACROSCOPIC VELOCITY GRADIENT
C     DEPENDING ON THE BOUNDARY CONDITIONS.

      DO I=1,3
        ROTBAR(I,I)=0.
        DO J=1,3
          IF(IUDOT(I,J).EQ.1.AND.IUDOT(J,I).EQ.1) THEN
            ROTBAR(I,J)=(UDOT(I,J)-UDOT(J,I))/2.
          ELSE IF(IUDOT(I,J).EQ.1) THEN
            ROTBAR(I,J)=UDOT(I,J)-DSIM(I,J)
            ROTBAR(J,I)=-ROTBAR(I,J)
          ENDIF
          UDOT(I,J)=DSIM(I,J)+ROTBAR(I,J)
        ENDDO
      ENDDO

C     VON MISES STRAIN-RATE & STRESS.
C     IF INTERACT=0: DBAR IS IMPOSED AND WE DEFINE SBAR=SAV.
C     IF INTERACT>0: DBAR & SBAR FOLLOW FROM THE CALL TO SUBR. STATE6x6.

      SVM=0.
      DVM=0.
      DO I=1,5
        SVM=SVM+SBAR(I)*SBAR(I)
        DVM=DVM+SBAR(I)*DBAR(I)
      ENDDO
      SVM=SQRT(SVM*3./2.)
      DVM=DVM/SVM

C *************************************************************************
C     RENORMALIZE ALL ARRAYS WITH UNITS OF STRESS, STRAIN RATE & COMPLIANCE

C       IF(IHARDLAW.NE.1) THEN
C         IF(IRATESENS.EQ.1) CALL NORMALIZE (1,ISTEP)
C       ENDIF

C *************************************************************************
C     WRITES OR STORES STRESS-STRAIN STATE FOR THE STEP.

      IF(IVGVAR.LE.1) CALL WRITE_STRESS_STRAIN (ISTEP)
      IF(IVGVAR.EQ.2) CALL PCYS (ISTEP,INDX,INDY,2)      ! 2D PCYS
      IF(IVGVAR.EQ.-2) THEN                              ! 5D PCYS
        CALL PCYS_IT (ISTEP,IPCYSSKIP,2)
        GO TO 2500
      ENDIF
      IF(IVGVAR.EQ.3) CALL LANKFORD (ISTEP,DELTALANK,2)

C *******************************************************************
C     STATISTICS & OUTPUT
C *******************************************************************

cwx      IF(ICUBCOM.EQ.1) CALL CUBCOMP (ISTEP,1)
      IF(ICUBCOM.EQ.1.AND.ISTEP.GT.1) CALL CUBCOMP (ISTEP,1)

      CALL STAT_SHEAR_ACTIVITY
      IF(IHARDLAW.NE.2) CALL STAT_GRAIN_SHAPE
      IF(IHARDLAW.NE.2) CALL STAT_STRESS_STRAIN

      CALL WRITE_STAT(ISTEP)
      CALL WRITE_ACTIV(ISTEP)

      IF(IVGVAR.GE.2) GO TO 2500

      IF(ISTEP.EQ.NSTEPSX) GO TO 2000    ! LAST STEP IS ONLY TO UPDATE STRESS

C ***********************************************************************
C     GIVEN THE OVERALL STRAIN RATE AND STRESS THEN:
C     *  IF ICTRL=0 CALCULATES TIME INCREMENT NECESSARY TO ACHIEVE
C        THE IMPOSED VON MISES STRAIN INCREMENT.
C     *  IF ICTRL=1,6 CALCULATES TIME INCREMENT NECESSARY TO ACHIEVE
C        THE IMPOSED STRAIN COMPONENT INCREMENT.
C     *  IF ICTRL=7 AND STRAIN_CONTROL=1 A TIME INCREMENT IS IMPOSED.
C        --> CALCULATES THE VON MISES STRAIN INCR ASSOCIATED WITH THE
C            IMPOSED STRAIN-RATE.
C     *  IF STRAIN_CONTROL=0 IT CORRESPONDS TO A CREEP SIMULATION
C        --> CALCULATES THE STRAIN INCREMENT ASSOCIATED WITH THE IMPOSED
C            STRESS COMPONENT AND TIME INCREMENT
C ***********************************************************************

      IF(ICTRL.EQ.0) THEN
        TINCR  =EVMINCR/DVM
        EPSACU =EPSACU+EVMINCR
      ELSE IF(ICTRL.GE.1.AND.ICTRL.LE.6 .AND. STRAIN_CONTROL.EQ.1) THEN
        TINCR  =EIJINCR/ABS(DSIMCTRL)
        EPSACU =EPSACU+EIJINCR
      ELSE IF(ICTRL.EQ.7 .OR. STRAIN_CONTROL.EQ.0) THEN
        EVMINCR=DVM*TINCR
        EPSACU =EPSACU+EVMINCR
      ENDIF

      EPSVM=EPSVM+DVM*TINCR
      DO I=1,3
      DO J=1,3
        EPSTOT(I,J)=EPSTOT(I,J)+DSIM(I,J)*TINCR
      ENDDO
      ENDDO

C ***********************************************************************
C     UPDATES ORIENTATION, HARDENING AND SHAPE OF EVERY GRAIN MAKING
C     A LINEAR EXTRAPOLATION OF THE CALCULATED RATE TIMES 'TINCR'.
C     (SKIPS UPDATE FOR PCYS OR LANKFORD RUN)
C ***********************************************************************

        CALL UPDATE_FIJ(0)        ! UPDATES DEFORMATION TENSOR OF ELEMENT
        CALL UPDATE_SHAPE(0)      ! UPDATES SHAPE OF ELEMENT

        DO IPH=IPHBOT,IPHTOP
          IPHEL=IPH-IPHBOT+1
          IF(NTWMOD(IPHEL).NE.0) CALL UPDATE_TWINNING (IPH)
        ENDDO

        IF(IUPDORI.EQ.1) CALL UPDATE_ORIENTATION

        DO IPH=IPHBOT,IPHTOP
          CALL UPDATE_FIJ(IPH)    ! UPDATES DEF TENSOR OF PHASE & GRAINS
          IF(IUPDSHP.EQ.1) CALL UPDATE_SHAPE (IPH)
        ENDDO

        IF(IUPDHAR.EQ.1) THEN

C *** VOCE HARDENING PLUS PREDOMINANT TWIN REORIENTATION SCHEME
          IF(IHARDLAW.EQ.0)  CALL UPDATE_CRSS_VOCE (2)
C *** HARDENING FOR IRRADIATED MATERIAL
CCC         IF(IHARDLAW.EQ.-1) CALL UPDATE_CRSS_IRRAD(2)
C *** MECHANICAL THRESHOLD STRESS HARDENING (NO TWINNING, ONLY ONE PHASE)
          IF(IHARDLAW.EQ.1) THEN
            TEMPX=TEMPERAT
            EDOT=SQRT(2./3.)*VNORM(DBAR,5)
            CALL UPDATE_CRSS_MTS (EDOT,TEMPX,3)
            TEMPERAT=TEMPX
          ENDIF
C *** DISLOCATION DENSITY STRESS HARDENING (ONLY ONE PHASE)
          IF(IHARDLAW.EQ.20) CALL UPDATE_CRSS_DD (2)
          IF(IHARDLAW.EQ.22) CALL UPDATE_CRSS_DD_CE2 (2)

        ENDIF

 2000 CONTINUE

C *******************************************************************
C *** END OF FORWARD UPDATE OF STRAIN, SHAPE, ORIENTATION & CRSS
C *******************************************************************

C *** CALCULATES evm/svm/taylor/work FOR EACH GRAIN USING sg/dg/dbar

      CALL GRAIN_INFO

C *******************************************************************
C     WRITE TEXTURE FILES FOR EACH PHASE
C *******************************************************************

      IF(IVGVAR.LE.1) THEN
        IWRITE=0
        IF(NWRITE.EQ.0) THEN
          IF(ISTEP.EQ.NSTEPS) IWRITE=1
        ELSE IF(NWRITE.NE.0) THEN
          IF(MOD(ISTEP,NWRITE).EQ.0) IWRITE=1
        ENDIF
        IF(IWRITE.EQ.1) CALL WRITE_TEXTURE
      ENDIF

C     WRITES GRAIN'S & PX STATES INTO BINARY FILE 'POSTMORT.OUT'

      IF(ISAVE.EQ.ISTEP) THEN
        OPEN(UW2,file='POSTMORT.OUT',form='UNFORMATTED',
     #           access='SEQUENTIAL',status='UNKNOWN')
        CALL POSTMORTEM (2)
        CLOSE(UNIT=UW2)
        ISAVE=-1
      ENDIF

CCC   CALL ELSC (1)     ! SC CALCULATION OF ELASTIC MODULI AT EACH STEP

 2500 CONTINUE      ! END OF DO LOOP OVER DEFORMATION INCREMENTS

 3000 CONTINUE      ! END OF DO LOOP OVER PROCESSES

C *******************************************************************

c      CALL CPU_TIME (END_TIME)
c      DELTIME=END_TIME-START_TIME
c      WRITE(*,*)
c      WRITE(*,'('' TOTAL RUN TIME:'',F10.2,'' secs'')') DELTIME

      STOP
      END

C*************************************************************************
C     INCLUDE HERE THE FILE 'VPSC7.SUB' WHICH CONTAINS THE SUBROUTINES
C     FOR PERFORMING THE SELF-CONSISTENT CALCULATION AND THE REORIENTATION
C
C     INCLUDE ALSO THE LIBRARY OF CANNED SUBROUTINES 'LIBRARY.SUB'
C
C*************************************************************************

      INCLUDE 'vpsc7.sub'
      INCLUDE 'library7.sub'

ext: for


PROGRAM p13

    IMPLICIT NONE
    INTEGER, PARAMETER :: st=8
    INTEGER :: r,s

    LOGICAL, DIMENSION(10,10) :: k1, k2

    LOGICAL, DIMENSION(11:20, 11:20) :: m

    COMPLEX, DIMENSION(1:10, 2+ST) :: z,c

    CHARACTER(LEN=10), DIMENSION(1:3) :: napisy

    DO r=1,10
        DO s=1,10
            IF(r==s) THEN
                k1(r,s) = .TRUE.
            ELSE
                k1(r,s) = .FALSE.
            END IF
        END DO
    END DO

    k2 = k1
    m=k2

    PRINT *, m

    PRINT *
    DO s=11,20
        PRINT *, (m(s,r), r=11,20,1)
    END DO

END PROGRAM p13


module PSERemeshModule

use NumberKindsModule
use OutputWriterModule
use LoggerModule
use ParticlesModule
use EdgesModule
use FacesModule
use PolyMesh2dModule
use FieldModule
use MPISetupModule
use PlaneGeomModule
use PlanarIncompressibleModule
use BetaPlaneMeshModule
use RefinementModule
use PSEDirectSumModule
use PlanarSWEModule

implicit none

include 'mpif.h'

private

!
!----------------
! Module types and public declarations
!----------------
!

!
!----------------
! Module interfaces
!----------------
!

!
!----------------
! Logging
!----------------
!
logical(klog), save :: logInit = .FALSE.
type(Logger) :: log
character(len=28), save :: logKey = 'PSERemesh'
integer(kint), parameter :: logLevel = DEBUG_LOGGING_LEVEL
character(len=MAX_STRING_LENGTH) :: logString

contains

!
!----------------
! public methods
!----------------
!

end module


text:
        SUBROUTINE MA_SC1D  ( marrpt, iret )
C************************************************************************
C* MA_SC1D						                *
C*								        *
C* This subroutine calls the routines which decode the groups in        *
C* section 1 of the drifting buoy report.  The section contains         *
C* meteorological and other non-marine data.                            *
C*                                                                      *
C* MA_SC1D  ( MARRPT, IRET )                                            *
C*                                                                      *
C* Input parameters:						        *
C*      MARRPT          CHAR*           Report array                    *
C*	LSEC1           INTEGER         Length of section 1 in report   *
C*      ISEC1           INTEGER         Pointer to start of section 1   *
C*								        *
C* Output parameters:						        *
C*	RIVALS(IRQDS1) 	REAL		Quality control indicator for   *
C*	                                section 1. WMO Code table 3334. *
C*	RIVALS(IRQXS1) 	REAL		Indicator of position of group  *
C*	IRET	 	INTEGER	    	Return code		        *
C*				   	  0 = Normal return 	        *
C*                                        1 = Problems                  *
C**								        *
C* Log:								        *
C* R. Hollern/NCEP      6/96                                            *
C* D. Kidwell/NCEP      4/97	Changed interface, reorganized header   *
C*				and comments                            *
C* D. Kidwell/NCEP     10/97	Changed interface                       *
C* R. Hollern/NCEP      7/99	Decoded quality control data            *
C************************************************************************
        INCLUDE         'GEMPRM.PRM'
        INCLUDE  	'macmn.cmn'
C*
        CHARACTER*(*)   marrpt
C*
        INTEGER   	jflg(7)
C*
        CHARACTER       fld1*1
C------------------------------------------------------------------------
        iret = 0
        ip = isec1
C
C*      A group should only appear once in a section.  This array will
C*      be used to flag those groups that are decoded in the section.
C
        DO i = 1, 7
            jflg ( i ) = 0
        END DO
C
        ipt = ip
C
C*      Get the quality control indicator for section 1.
C
        fld1 = marrpt (ipt:ipt)
        CALL  ST_INTG ( fld1, ival, ier )
        IF ( ier .eq. 0 ) rivals ( irqds1 ) = FLOAT ( ival )
C
        ipt = ipt + 1
C
C*      Get the indicator of position of group
C
        fld1 = marrpt (ipt:ipt)
        CALL  ST_INTG ( fld1, ival, ier )
        IF ( ier .eq. 0 ) rivals ( irqxs1 ) = FLOAT ( ival )
C
        iend = isec1 + lsec1 - 1
C
        DO WHILE ( ipt .lt. iend )
C
            ipt = ipt + 1
C
            IF ( marrpt ( ipt:ipt+1 ) .eq. ' 0' .and.
     +           jflg ( 1 ) .eq. 0 ) THEN
C
C*              Decode the 0ddff group.
C
                ipt = ipt + 2
                CALL  MA_DDFF ( marrpt, ipt, jret)
                jflg (1) = 1
              ELSE IF ( marrpt (ipt:ipt+1) .eq. ' 1' .and.
     +                  jflg (2) .eq. 0 ) THEN
C
C*              Decode air temperature group.
C
                ipt = ipt + 2
                iparam = 1
                CALL  MA_TEMP ( marrpt, iparam, ipt, jret )
                jflg (2) = 1
              ELSE IF ( marrpt (ipt:ipt+2) .eq. ' 29' .and.
     +                  jflg (3) .eq. 0 ) THEN
C
C*              Decode relative humidity group.
C
                ipt = ipt + 3
                CALL  MA_RELH ( marrpt, ipt, jret )
                jflg (3) = 1
              ELSE IF ( marrpt (ipt:ipt+1) .eq. ' 2' .and.
     +                  jflg (4) .eq. 0 ) THEN
C
C*              Decode dew point temperature group.
C
                ipt = ipt + 2
                iparam = 2
                CALL  MA_TEMP ( marrpt, iparam, ipt, jret )
                jflg (4) = 1
              ELSE IF ( marrpt (ipt:ipt+1) .eq. ' 3' .and.
     +                  jflg (5) .eq. 0 ) THEN
C
C*              Decode the station pressure group.
C
                iparam = 0
                ipt = ipt + 2
                CALL  MA_PRES ( marrpt, iparam, ipt, jret )
                jflg (5) = 1
              ELSE IF ( marrpt (ipt:ipt+1) .eq. ' 4' .and.
     +               jflg (6) .eq. 0 ) THEN
C
C*              Decode the mean sea level pressure group.
C
                iparam = 1
                ipt = ipt + 2
                CALL  MA_PRES ( marrpt, iparam, ipt, jret )
                jflg (6) = 1
              ELSE IF ( marrpt (ipt:ipt+1) .eq. ' 5' .and.
     +                  jflg (7) .eq. 0 ) THEN
C
C*              Decode the station 3 hourly pressure change.
C
                ipt = ipt + 2
                CALL  MA_PRS3 ( marrpt, ipt, jret )
                jflg (7) = 1
            END IF
        END DO
C*
	RETURN
        END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
C    @(#)swpanl.f	20.3 2/13/96
      subroutine swpanl   (m,n)

      include 'ipfinc/parametr.inc'
      include 'ipfinc/anlys.inc'
      save
C
      entry swapan (m,n)

      kps = lstown(m)
      lstown(m) = lstown(n)
      lstown(n) = kps
      return

      entry swpan2  (m,n)

      kps = lstvlt(m)
      lstvlt(m) = lstvlt(n)
      lstvlt(n) = kps
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      program t
      character(len=6)::fname='infile'
      integer::ierr
      logical::p
      p=.false.
      inquire(file=fname,exist=p,iostat=ierr)
      write (*,fmt='(l)',iostat=ierr) .true.
      open(unit=77,file=fname,status='old',iostat=ierr)
      if (ierr.ne.0) then
         stop 'unable to open file: infile'
      endif
      p=.false.
      inquire(unit=77,opened=p,iostat=ierr)
      write (*,fmt='(l)',iostat=ierr) .true.
      close(77)
      end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SMLIB Module
! Author: Brian Krauth, Simon Fraser University
!
! Machine-specific code for SMLE program.  By default, the
! program uses the standard random number generator associated
! with the random_number subroutine.  The standard random number
! generator is not always all that good; so this file can be altered
! to substitute a preferred algorithm.
!
! Public procedures defined in this module are:
!
!   runif:	        Draw random numbers from uniform distribution
!   setseed: 		Reseed random number generator based on system
!			clock.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module bklib ! contains any machine-specific code
implicit none
public :: setseed, runif, read_buffer, get_parmfile
private :: runif0s,runif1s,runif2s,runif0d,runif1d,runif2d
integer, parameter, public :: SP=KIND(1.0),DP=selected_real_kind(9,99), strlen=32767
interface runif
   module procedure runif0s,runif1s,runif2s,runif0d,runif1d,runif2d
end interface
! external g05faf,g05ccf ! These functions are in the NAG library, which is unavailable on this machine

contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SETSEED subroutine
!
! Reseeds random number generator based on system clock.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine setseed()
  integer :: i,n
  real :: j
  call system_clock(n)
  n = modulo(n,40000)
  do i=1,n
     call random_number(j)
  end do
!   call g05ccf
end subroutine setseed

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! RUNIF, RUNIFD functions
!
! Draw pseudorandom number(s) from uniform(0,1) distribution.
!
! RUNIF gives single-precision numbers, RUNIFD gives double.
!
! Format: RUNIF() gives a single number, RUNIF(n) gives a
! vector of length n, and RUNIF(r,c) gives an r-by-c matrix.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine runif0s(x)
  real(kind=SP), intent(inout) :: x
  call random_number(x)
end subroutine runif0s
subroutine runif1s(x)
  real(kind=SP), dimension(:), intent(inout) :: x
  call random_number(x)
end subroutine runif1s
subroutine runif2s(x)
  real(kind=SP), dimension(:,:), intent(inout) :: x
  call random_number(x)
end subroutine runif2s
subroutine runif0d(x)
  real(kind=DP), intent(inout) :: x
  call random_number(x)
end subroutine runif0d
subroutine runif1d(x)
  real(kind=DP), dimension(:), intent(inout) :: x
  call random_number(x)
end subroutine runif1d
subroutine runif2d(x)
  real(kind=DP), dimension(:,:), intent(inout) :: x
  call random_number(x)
end subroutine runif2d

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! READ_BUFFER subroutine
!
! Format: read_buffer(buffer)
!
! Reads a line of data into a character buffer. This subroutine
! was written to work around a problem in the IBM XLFortran
! compiler: when reading in a line of text from a DOS/Windows
! text file (i.e., one with lines ending in a CR-LF rather than a LF)
! the program would include the CR as a character rather than
! dropping it.  This leads to all sorts of problems, so we
! just read in each line into a big string buffer, and change
! any CR at the end of the record into a space.
!
! Note that this workaround is only necessary for the IBM
! compiler; other compilers seem to strip the CR automatically.
! It also comes at the cost of imposing a fixed limit on
! the size of a record (right now I have it at 1000 characters).
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine read_buffer(unit,buffer)
    integer, intent(in) :: unit
    character(len=*), intent(inout) :: buffer
    integer :: ios
    ! Read the data into the text buffer
    read (unit=unit,iostat=ios,fmt="(a)") buffer
    ! If the last character is anything but a blank, then we might
    ! have records that are too big.
    if (len(buffer) == len_trim(buffer)) then
       stop "Error: Records in data file are larger than maximum buffer size"
    else
       ! Here's the important part.  If the lat non-blank character is a CR, replace it with a space.
       if (iachar(buffer(len_trim(buffer):len_trim(buffer))) == 13) then
          buffer(len_trim(buffer):len_trim(buffer)) = " "
       end if
    end if
  end subroutine read_buffer

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! GET_PARMFILE subroutine
!
! Format: get_parmfile(parmfile)
!
! Retrieves the name of the input file if provided in the command
! arguments.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine get_parmfile(parmfile)
    character(len=*), intent(inout) :: parmfile
    character(len=strlen) :: buffer
    integer :: status
    status = command_argument_count()
    if (status > 0) then
      call get_command_argument(1,buffer)
      parmfile = trim(buffer)
    end if
  end subroutine get_parmfile

end module bklib

max_stars_repo_name: bvkrauth/smle

max_issues_repo_name: bvkrauth/smle

max_forks_repo_name: bvkrauth/smle


text:
! { dg-do compile }
! { dg-shouldfail "VOLATILE not part of F95" }
! { dg-options "-std=f95" }
! Test whether volatile statements and attributes are rejected
! with -std=f95.
! PR fortran/29601
program volatile_test
  implicit none
  real, volatile :: foo ! { dg-error "VOLATILE attribute" }
  real :: l
  volatile :: l         ! { dg-error "VOLATILE statement" }
  l   = 4.0
  foo = 3.0             ! { dg-error "no IMPLICIT type" }
end program volatile_test

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
module GAMESS_read_module

use M_strings
use hess_data
implicit none

contains
!------------------
subroutine GAMESS_read_number_of_atoms(file_id)

implicit none
integer, intent(in)                                    :: file_id
!-----------------------------------------------------------------
character(len=80)                                      :: line
character(len=3)                                       :: point_g
!-----------------------------------------------------------------

                                                  ! Need to shift position to the line containing $DATA

read(file_id,'(A)') line                          ! Check from first line

do while (index(line, '$DATA') == 0)
 read(file_id,'(A)') line
end do                                            ! Now on the $DATA line

                                                  ! Now need to shift position to the line containing $END
read(file_id,'(A)') line                          ! The line next to $DATA: mol. spec.
read(file_id,'(A)') line                          ! The line next to mol. spec.:
                                                  ! the point group of the molecule
point_g = trim(adjustl(line))                     ! Assign the point group of the molecule
if (point_g /= 'C1') then                         ! Check the point group
  read(file_id,'(A)') line                        ! Skip one line if not C1
end if

number_of_atoms = 0

do while (index(line, '$END') == 0)
 read(file_id,'(A)') line
 number_of_atoms = number_of_atoms + 1
end do

number_of_atoms = number_of_atoms - 1             ! Looks like DO WHILE counts the extra $END
                                                  ! line which is not the atom containing line
rewind file_id                                    ! Positions the file associated with the
                                                  ! specified unit to its initial point
end subroutine GAMESS_read_number_of_atoms
!--------------------------------------------------------------------------
subroutine GAMESS_read_atomic_data(file_id)

implicit none
integer, intent(in)                                             :: file_id
!-------------------------------------------------------------------------
character(len=80)                                               :: line
character(len=3)                                                :: point_g
character(len=80), allocatable                                  :: array(:)
integer                                                         :: i
!-------------------------------------------------------------------------
allocate  (symbol (number_of_atoms) )
allocate  (charge (number_of_atoms) )
allocate  (coord (number_of_atoms,3) )
!-------------------------------------------------------------------------
                                                   ! Need to shift position to the
                                                   ! line containing $DATA

do while (index(line, '$DATA') == 0)
 read(file_id,'(A)') line
end do                                             ! Now on the $DATA line

read(file_id,'(A)') line                           ! The line next to $DATA: mol. spec.

read(file_id,'(A)') line                           ! The line next to mol. spec.:
                                                   ! the point group of the molecule
point_g = trim(adjustl(line))                      ! Assign the point group of the molecule

if (point_g /= 'C1') then                          ! Check the point group
 read(file_id,'(A)') line                          ! Skip one line if not C1
end if

do i =1, number_of_atoms
 read(file_id,'(A)') line
 call split(line, array)
 read(array(1),*) symbol(i)
 read(array(2),*) charge(i)
 read(array(3),*) coord(i,1)
 read(array(4),*) coord(i,2)
 read(array(5),*) coord(i,3)
end do

end subroutine GAMESS_read_atomic_data
!---------------------------------------------------------------------------------
subroutine GAMESS_read_grad_hess(file_id,grad,Hess)

implicit none
integer, intent(in)                                             :: file_id
double precision, dimension(:), allocatable, intent(out)        :: grad
double precision, dimension(:,:), allocatable, intent(out)      :: Hess
!----------------------------------------------------------------------------------
character(len=80)                                               :: line
character(len=80), allocatable                                  :: array(:)
integer                                                         :: i, j, k, l, last_number
!----------------------------------------------------------------------------------
allocate  (grad    (3*number_of_atoms) )
allocate  (Hess (3*number_of_atoms,3*number_of_atoms) )
!----------------------------------------------------------------------------------
! Read gradient
                                                   ! Need to shift position to the
                                                   ! line containing $GRAD

do while (index(line, '$GRAD') == 0)
 read(file_id,'(A)') line
end do                                             ! Now on the $GRAD line

read(file_id,'(A)') line                           ! Skip extra line containing
                                                   ! converged energy

do i = 1, number_of_atoms
 read(file_id,'(A)') line
 call split(line, array)
 read(array(3),*) grad(i*3-2)
 read(array(4),*) grad(i*3-1)
 read(array(5),*) grad(i*3)
end do

rewind file_id
!------------------------------------------------------------------------------------
! Read Hessian
                                                   ! Need to shift position to the
                                                   ! line containing $HESS

do while (index(line, '$HESS') == 0)
 read(file_id,'(A)') line
end do                                             ! Now on the $HESS line

read(file_id,'(A)') line                           ! Now on the line of
                                                   ! converged energy

do i = 1, 3*number_of_atoms
 do j = 1, ceiling(3*real(number_of_atoms)/5)
  read(file_id,'(A)') line
  if (j /= ceiling(3*real(number_of_atoms)/5)) then
   last_number = 5
  else
   last_number = 5-5*ceiling(3*real(number_of_atoms)/5) &
                   + 3*number_of_atoms
  end if

  do l = 1, last_number
   k = j*5 + l - 5
   read(line((15*(l-1)+6):(15*l+5)),*) Hess(i,k)
  end do

 end do
end do

end subroutine GAMESS_read_grad_hess
!-------------------------------------------------------------------------------------
end module GAMESS_read_module

max_stars_repo_name: svarganov/NAST

max_issues_repo_name: svarganov/NAST

max_forks_repo_name: svarganov/NAST


module stppsmod

!$$$ module documentation block
!           .      .    .                                       .
! module:   stppsmod    module for stpps and its tangent linear stpps_tl
!  prgmmr:
!
! abstract: module for stpps and its tangent linear stpps_tl
!
! program history log:
!   2005-05-18  Yanqiu zhu - wrap stpps and its tangent linear stpps_tl into one module
!   2005-11-16  Derber - remove interfaces
!   2008-12-02  Todling - remove stpps_tl
!   2009-08-12  lueken - update documentation
!   2010-05-13  todling - uniform interface across stp routines
!   2013-10-28  todling - reame p3d to prse
!   2014-04-12       su - add non linear qc from Purser's scheme
!   2015-02-26       su - add njqc as an option to chose new non linear qc
!   2016-05-18  guo     - replaced ob_type with polymorphic obsNode through type casting
!
! subroutines included:
!   sub stpps
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

implicit none

PRIVATE
PUBLIC stpps

contains

subroutine stpps(pshead,rval,sval,out,sges,nstep)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    stpps       calculate penalty and contribution to
!                             stepsize for sfcp, using nonlinear qc
!   prgmmr: derber           org: np23                date: 1991-02-26
!
! abstract: calculate penalty and contribution to stepsize for
!           surface pressure with nonlinear qc.
!
! program history log:
!   1991-02-26  derber
!   1997-12-14  weiyu yang
!   1999-08-24  derber, j., treadon, r., yang, w., first frozen mpp version
!   2004-08-02  treadon - add only to module use, add intent in/out
!   2004-10-08  parrish - add nonlinear qc option
!   2005-04-11  treadon - merge stpps and stpps_qc into single routine
!   2005-08-02  derber  - modify for variational qc parameters for each ob
!   2005-09-28  derber  - consolidate location and weight arrays
!   2005-10-21  su      - modify for variational qc
!   2006-07-28  derber  - modify to use new inner loop obs data structure
!   2006-09-18  derber  - modify to output of b1 and b3
!   2007-06-04  derber  - use quad precision to get reproducability over number of processors
!   2008-12-03  todling - changed handling of ptr%time
!   2010-01-04  zhang,b - bug fix: accumulate penalty for multiple obs bins
!   2010-05-13  todling  - update to use gsi_bundlemod
!   2015-12-21  yang    - Parrish's correction to the previous code in new varqc.
!
!   input argument list:
!     pshead
!     rp       - search direction for ps
!     sp       - analysis increment for ps
!     sges     - step size estimates (nstep)
!     nstep    - number of stepsize estimates  (==0 means use outer iteration values)
!
!   output argument list:
!     out(1:nstep)   - contribution to penalty for surface pressure - sges(1:nstep)
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind,r_quad
  use qcmod, only: nlnqc_iter,varqc_iter,njqc,vqc
  use constants, only: half,one,two,tiny_r_kind,cg_term,zero_quad,r3600
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use m_obsNode, only: obsNode
  use m_psNode , only: psNode
  use m_psNode , only: psNode_typecast
  use m_psNode , only: psNode_nextcast
  implicit none

! Declare passed variables
  class(obsNode), pointer             ,intent(in   ) :: pshead
  integer(i_kind)                     ,intent(in   ) :: nstep
  real(r_quad),dimension(max(1,nstep)),intent(inout) :: out
  type(gsi_bundle)                    ,intent(in   ) :: rval,sval
  real(r_kind),dimension(max(1,nstep)),intent(in   ) :: sges

! Declare local variables
  integer(i_kind) j1,j2,j3,j4,kk,ier,istatus
  real(r_kind) val,val2,w1,w2,w3,w4
  real(r_kind) cg_ps,ps,wgross,wnotgross,ps_pg
  real(r_kind),dimension(max(1,nstep))::pen
  real(r_kind),pointer,dimension(:) :: sp
  real(r_kind),pointer,dimension(:) :: rp
  type(psNode), pointer :: psptr

  out=zero_quad

!  If no ps data return
  if(.not. associated(pshead))return
! Retrieve pointers
! Simply return if any pointer not found
  ier=0
  call gsi_bundlegetpointer(sval,'prse',sp,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(rval,'prse',rp,istatus);ier=istatus+ier
  if(ier/=0)return

  psptr => psNode_typecast(pshead)
  do while (associated(psptr))
     if(psptr%luse)then
        if(nstep > 0)then
           j1 = psptr%ij(1)
           j2 = psptr%ij(2)
           j3 = psptr%ij(3)
           j4 = psptr%ij(4)
           w1 = psptr%wij(1)
           w2 = psptr%wij(2)
           w3 = psptr%wij(3)
           w4 = psptr%wij(4)
           val =w1* rp(j1)+w2* rp(j2)+w3* rp(j3)+w4* rp(j4)
           val2=w1* sp(j1)+w2* sp(j2)+w3* sp(j3)+w4* sp(j4)-psptr%res
           do kk=1,nstep
              ps=val2+sges(kk)*val
              pen(kk)=ps*ps*psptr%err2
           end do
        else
           pen(1)=psptr%res*psptr%res*psptr%err2
        end if

!  Modify penalty term if nonlinear QC

        if (vqc .and. nlnqc_iter .and. psptr%pg > tiny_r_kind .and.  &
                             psptr%b  > tiny_r_kind) then
           ps_pg=psptr%pg*varqc_iter
           cg_ps=cg_term/psptr%b
           wnotgross= one-ps_pg
           wgross =ps_pg*cg_ps/wnotgross
           do kk=1,max(1,nstep)
              pen(kk) = -two*log((exp(-half*pen(kk))+wgross)/(one+wgross))
           end do
        endif

!   for Dr. Jim purser' non liear quality control
        if(njqc  .and. psptr%jb > tiny_r_kind .and. psptr%jb <10.0_r_kind) then
           do kk=1,max(1,nstep)
              pen(kk) = two*two*psptr%jb*log(cosh(sqrt(pen(kk)/(two*psptr%jb))))
           enddo
           out(1) = out(1)+pen(1)*psptr%raterr2
           do kk=2,nstep
              out(kk) = out(kk)+(pen(kk)-pen(1))*psptr%raterr2
           end do
        else
           out(1) = out(1)+pen(1)*psptr%raterr2
           do kk=2,nstep
              out(kk) = out(kk)+(pen(kk)-pen(1))*psptr%raterr2
           end do
        endif

     end if

     psptr => psNode_nextcast(psptr)
  end do

  return
end subroutine stpps

end module stppsmod


text:
/*
vim:syntax=sml
*/

/* (Heap m) <: (Heap n) iff n <= m */
datatype Heap (min:Int) =
      Empty
    | Node of (x:Int. min <= x) (c:Int) (l:(Heap x)) (r:(Heap x));;

let n1 : (Heap 1) = (Empty 0);; /* BUG: bad arg */
let n2 : (Heap 1) = Node 1 1 6
                       (Node 1 2 2 (Empty 2) (Empty 2))
                       (Node 1 4 2 (Empty 4) (Empty 4));;

let foo (h:(Heap 0)) : Bool = true;;

foo n2;;

let min (a:Int) (b:Int) : (c:Int.(a>=c) && (b>=c)) =
  if[(c:Int.(a>=c) && (b>=c))] (a >= b) then
    b
  else
    a;;

let max = fn (a:Int) (b:Int) =>
  if[c:Int.(c>=a) && (c>=b)] (a >= b) then
    a
  else
    b;;

let heap_children (min:Int) (h:(Heap min)) : Int =
  caseHeap min h Int
    (fn u => 0)
    (fn x c l r => c);;

let rec insert (hmin:Int)
               (h:(Heap hmin))
               (v:Int)
               (rmin:Int.(rmin <= hmin) && (rmin <= v)) :
               (Heap rmin) =
  caseHeap hmin h (Heap rmin)
    (fn u =>
      Node v v 2 (Empty v) (Empty v))
    (fn x c l r =>
      if[Heap rmin] v < x
      then
        let newchildren : Int = (heap_children hmin h) + 1 in
        (Node rmin v newchildren (Node v x c l r) (Empty v))
      else /* v >= x */
        let lchildren : Int = heap_children x l in
        let rchildren : Int = heap_children x r in
        let newmin : (m:Int.(m <= x) && (m <= rmin)) = min x rmin in
        if[Heap rmin] lchildren < rchildren
        then
          Node rmin x (c + 1) (insert x l v x) r
        else
          Node rmin x (c + 1) l (insert x r v x))
;;

let extract_min (min:Int) (h:(Heap min)) : (r:Int.(min <= r)) =
  caseHeap min h (r:Int.(min <= r))
    (fn u => min) /* This seems weird */
    (fn x c l r => x)
;;

let h1 = insert 0 (Empty 0) 1 0;;
let h2 = insert 0 h1 2 0;;
extract_min 0 h2;;

max_stars_repo_name: ucsc-proglang/sage

max_issues_repo_name: ucsc-proglang/sage

max_forks_repo_name: ucsc-proglang/sage


text:
C    @(#)swprt2.f	20.3 2/13/96
      subroutine swprt2(i,j)

      include 'ipfinc/cgrate.inc'

      do 10 k=1,3
         itemp=lrat(k,i)
         lrat(k,i)=lrat(k,j)
         lrat(k,j)=itemp

         temp=rat(k,i)
         rat(k,i)=rat(k,j)
         rat(k,j)=temp
10    continue

      return
      end


text:
C    @(#)kmpkol.f	20.3 2/13/96
      function kmpkol (i,j)
C
      include 'ipfinc/parametr.inc'
      include 'ipfinc/renum.inc'
C
      kmpkol = kolum1(i) - kolum1(j)
      return
      end