text:
      SUBROUTINE second0(stime)
      USE stel_kinds
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec), INTENT(out) :: stime
      INTEGER :: cnt, cnt_rate
C-----------------------------------------------
!     USE CPU_TIME IF ACTUAL CPU USAGE DESIRED
!      CALL CPU_TIME(stime)
!      RETURN

      CALL SYSTEM_CLOCK(count=cnt, count_rate=cnt_rate)
      IF (cnt_rate .ne. 0) THEN
         stime = REAL(cnt, rprec)/cnt_rate
      ELSE
         stime = 0
      END IF

      END SUBROUTINE second0


text:
integer x(g),z,i,j
string s,s_str_(10)
data x/(5.:20.),(4.:-2.),(0.3:0.)/,z/.true./
common /com1/i,j
data s_str_/"string1", *"string2", "string3"/
$ no comment
y(z) = .not. re(x) .and. .not. x(0x10F,i,g(j)).gt.y .or. z .eq. k**3
100 call a
s = "sh"
im(x(1))=-2.1
if (b .and. (x .gt. a(i**k))) stop
return
do x=1,x(10,g(y+im(x(2)))),3
real list a
real i
data i/.314159e-31/,b/5E2/
1100 read (x,i=1,x(i)),z
goto 1000
cdr(a) = [0.1, 6., -3.2, 0.]
if (car(a).gt.0B.01.and.b.lt.x(i+2.)) then
integer z
read z,a(z)
goto z,(100,101,102,103,104,105)
else
logical l(x)
data l/.false.,.true.,.true./
do i=1,N,2
l(i) = x(i+g(i)).le.a
enddo
endif
call try_me_(x(u),y,l)
if (y(g(y(z-2)))+a(1,i)) 100,1000,10
y = .not. x(i)
1000 continue
enddo
end
subroutine try_me_(integer n,a(n))
do i=1,n
if (a(i) .gt. 0) a(i) = a(i) - i
enddo
return
end

ext: f

max_stars_repo_name: konmaz/Compiler-Design-Project

max_issues_repo_name: konmaz/Compiler-Design-Project

max_forks_repo_name: konmaz/Compiler-Design-Project


text:
program slash_descriptor
implicit none
character :: text*100, word*10
print "('cat'/'dog')"
print "(a)","cat" // new_line("") // "dog"
text = "3/20/2022"
read (text,*) word ! stops at '3' because / follows
print "(a)","word = '" // trim(word) // "'"
end program slash_descriptor
! output
! cat
! dog
! cat
! dog
! word = '3'

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


text:
      SUBROUTINE RSVFVM(NEM1,NEM2)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    RSVFVM
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 1994-01-06
C
C ABSTRACT: THIS SUBROUTINE STEPS THROUGH THE "FOLLOWING VALUE"
C   MNEMONIC NEM1 AND, FOR EACH "." CHARACTER ENCOUNTERED (EXCEPT FOR
C   THE INITIAL ONE), OVERWRITES IT WITH THE NEXT CORRESPONDING
C   CHARACTER FROM NEM2 (SEE REMARKS).
C
C PROGRAM HISTORY LOG:
C 1994-01-06  J. WOOLLEN -- ORIGINAL AUTHOR
C 2003-11-04  J. ATOR    -- ADDED DOCUMENTATION
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- UNIFIED/PORTABLE FOR WRF; ADDED HISTORY
C                           DOCUMENTATION
C
C USAGE:    CALL RSVFVM (NEM1, NEM2)
C   INPUT ARGUMENT LIST:
C     NEM1     - CHARACTER*8: "FOLLOWING VALUE" MNEMONIC
C     NEM2     - CHARACTER*8: MNEMONIC IMMEDIATELY FOLLOWING NEM1
C                WITHIN USER DICTIONARY TABLE
C
C   OUTPUT ARGUMENT LIST:
C     NEM1     - CHARACTER*8: COPY OF INPUT NEM1 WITH ALL "."
C                CHARACTERS (EXCEPT INITIAL ONE) OVERWRITTEN WITH
C                CORRESPONDING CHARACTERS FROM NEM2
C
C REMARKS:
C    FOR EXAMPLE:
C     if, on input:    NEM1 = ".DTH...."
C                      NEM2 = "MXTM    "
C     then, on output: NEM1 = ".DTHMXTM"
C
C
C    THIS ROUTINE CALLS:        None
C    THIS ROUTINE IS CALLED BY: NEMTBD   SEQSDX
C                               Normally not called by any application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      CHARACTER*8 NEM1,NEM2

      DO I=1,LEN(NEM1)
      IF(I.EQ.1) THEN

C        Skip initial "." and initialize J.

         J = 1
      ELSE
         IF(NEM1(I:I).EQ.'.') THEN
            NEM1(I:I) = NEM2(J:J)
            J = J+1
         ENDIF
      ENDIF
      ENDDO

      RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak


text:
!***********************************************************************
!  Integrated Water Flow Model (IWFM)
!  Copyright (C) 2005-2021
!  State of California, Department of Water Resources
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU General Public License
!  as published by the Free Software Foundation; either version 2
!  of the License, or (at your option) any later version.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  (http://www.gnu.org/copyleft/gpl.html)
!
!  You should have received a copy of the GNU General Public License
!  along with this program; if not, write to the Free Software
!  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!
!  For tecnical support, e-mail: IWFMtechsupport@water.ca.gov
!***********************************************************************
MODULE Class_Bypass
  USE MessageLogger                , ONLY: SetLastMessage                          , &
                                           EchoProgress                            , &
                                           MessageArray                            , &
                                           iFatal
  USE GeneralUtilities             , ONLY: StripTextUntilCharacter                 , &
                                           CleanSpecialCharacters                  , &
                                           IntToText                               , &
                                           GetArrayData                            , &
                                           AllocArray                              , &
                                           ConvertID_To_Index                      , &
                                           LocateInList
  USE IOInterface                  , ONLY: GenericFileType
  USE Package_Misc                 , ONLY: PairedDataType                          , &
                                           FlowDestinationType                     , &
                                           f_iFlowDest_Outside                     , &
                                           f_iFlowDest_StrmNode                    , &
                                           f_iFlowDest_Lake
  USE Class_RechargeZone           , ONLY: RechargeZoneType                        , &
                                           RechargeZone_New
  USE Class_StrmReach              , ONLY: StrmReachType                           , &
                                           StrmReach_GetReachNumber                , &
                                           StrmReach_GetReaches_InUpstrmNetwork    , &
                                           StrmReach_CompileReachNetwork
  USE Package_ComponentConnectors  , ONLY: StrmLakeConnectorType                   , &
                                           f_iBypassToLakeFlow
  IMPLICIT NONE

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** VARIABLE DEFINITIONS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- PUBLIC ENTITIES
  ! -------------------------------------------------------------
  PRIVATE
  PUBLIC :: BypassType                 ,  &
            Bypass_New                 ,  &
            f_iDestTypes

  ! -------------------------------------------------------------
  ! --- BYPASS DATA TYPE
  ! -------------------------------------------------------------
  TYPE,EXTENDS(FlowDestinationType) :: BypassType
      INTEGER                :: ID              = 0
      CHARACTER(LEN=20)      :: cName           = ''
      INTEGER                :: iNode_Exp       = 0
      INTEGER                :: iColBypass      = 0
      TYPE(PairedDataType)   :: RatingTable
      REAL(8)                :: FracRecvLoss    = 0.0
      REAL(8)                :: FracNonRecvLoss = 0.0
      REAL(8)                :: Bypass_Out      = 0.0
      REAL(8)                :: Bypass_Received = 0.0
      REAL(8)                :: RecvLoss        = 0.0
      REAL(8)                :: NonRecvLoss     = 0.0
      TYPE(RechargeZoneType) :: Recharge
  END TYPE BypassType

  ! -------------------------------------------------------------
  ! --- MISC. ENTITIES
  ! -------------------------------------------------------------
  INTEGER,PARAMETER                   :: f_iDestTypes(3) = [f_iFlowDest_Outside , f_iFlowDest_StrmNode , f_iFlowDest_Lake]
  INTEGER,PARAMETER                   :: ModNameLen = 14
  CHARACTER(LEN=ModNameLen),PARAMETER :: ModName    = 'Class_Bypass::'

CONTAINS

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** CONSTRUCTORS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- INSTANTIATE A SET OF BYPASSES FROM A FILE
  ! -------------------------------------------------------------
  SUBROUTINE Bypass_New(cFileName,NStrmNodes,iStrmNodeIDs,iElemIDs,iLakeIDs,Reaches,StrmLakeConnector,TUnitStrmFlow,TUnitBypass,Bypasses,iStat)
    CHARACTER(LEN=*),INTENT(IN)         :: cFileName
    INTEGER,INTENT(IN)                  :: NStrmNodes,iStrmNodeIDs(NStrmNodes),iElemIDs(:),iLakeIDs(:)
    TYPE(StrmReachType)                 :: Reaches(:)
    TYPE(StrmLakeConnectorType)         :: StrmLakeConnector
    CHARACTER(LEN=6),INTENT(OUT)        :: TUnitStrmFlow,TUnitBypass
    TYPE(BypassType),TARGET,ALLOCATABLE :: Bypasses(:)
    INTEGER,INTENT(OUT)                 :: iStat

    !Local variables
    CHARACTER(LEN=ModNameLen+10) :: ThisProcedure = ModName // 'Bypass_New'
    INTEGER                      :: NBypass,ErrorCode,indxBypass,ID,iNum,NPoints,iReach_Exp,iReach_Imp, &
                                    indxBypass1,iNode_Exp_ID,iDest_ID
    TYPE(GenericFileType)        :: InFile
    REAL(8)                      :: FactFlow,FactBypass,DummyArray(7)
    CHARACTER                    :: ALine*2000
    TYPE(BypassType),POINTER     :: pBypass
    INTEGER,ALLOCATABLE          :: iBypassIDs(:),iBypassOutReachIDs(:),iBypassInReachIDs(:),iReachesUpNetwork(:)
    REAL(8),ALLOCATABLE          :: Dummy2DRealArray(:,:)

    !Initialize
    iStat = 0

    !Return if no filename is supplied
    IF (cFileName .EQ. '') THEN
        ALLOCATE (Bypasses(0))
        RETURN
    END IF

    !Echo progress
    CALL EchoProgress('Instantiating bypasses')

    !Open file
    CALL InFile%New(FileName=TRIM(ADJUSTL(cFileName)),InputFile=.TRUE.,IsTSFile=.FALSE.,Descriptor='Bypass specifications data file',iStat=iStat)
    IF (iStat .EQ. -1) RETURN

    !Number of bypasses
    CALL InFile%ReadData(NBypass,iStat)  ;  IF (iStat .EQ. -1) RETURN
    ALLOCATE (Bypasses(NBypass) , iBypassIDs(NBypass) , iBypassOutReachIDs(NBypass) , iBypassInReachIDs(NBypass) , STAT=ErrorCode)
    IF (ErrorCode .NE. 0) THEN
        CALL SetLastMessage('Error in allocating memory for bypasses!',iFatal,ThisProcedure)
        iStat = -1
        RETURN
    END IF

    !Return if no bypasses are specified
    IF (NBypass .EQ. 0) THEN
        CALL InFile%Kill()
        RETURN
    END IF

    !Read basic data
    CALL InFile%ReadData(FactFlow,iStat)  ;  IF (iStat .EQ. -1) RETURN
    CALL InFile%ReadData(ALine,iStat)  ;  IF (iStat .EQ. -1) RETURN ; CALL CleanSpecialCharacters(ALine)
    TUnitStrmFlow = ADJUSTL(StripTextUntilCharacter(ALine,'/'))
    CALL InFile%ReadData(FactBypass,iStat)  ;  IF (iStat .EQ. -1) RETURN
    CALL InFile%ReadData(ALine,iStat)  ;  IF (iStat .EQ. -1) RETURN ; CALL CleanSpecialCharacters(ALine)
    TUnitBypass = ADJUSTL(StripTextUntilCharacter(ALine,'/'))

    !Read bypass data
    DO indxBypass=1,NBypass
        CALL InFile%ReadData(ALine,iStat)  ;  IF (iStat .EQ. -1) RETURN
        CALL GetArrayData(ALine,DummyArray,'by-pass number '//TRIM(IntToText(indxBypass)),iStat)
        IF (iStat .EQ. -1) RETURN

        pBypass             => Bypasses(indxBypass)
        pBypass%cName       =  ALine(1:20)
        ID                  =  INT(DummyArray(1))
        pBypass%ID          =  ID
        iNode_Exp_ID        =  INT(DummyArray(2))
        pBypass%iDestType   =  INT(DummyArray(3))
        iDest_ID            =  INT(DummyArray(4))
        iNum                =  INT(DummyArray(5))
        IF (iNum .GT. 0) THEN
            Bypasses(indxBypass)%iColBypass = iNum
        ELSEIF (iNum .LT. 0) THEN
            NPoints = -iNum
        END IF
        pBypass%FracRecvLoss    = DummyArray(6)
        pBypass%FracNonRecvLoss = DummyArray(7)

        !Make sure same ID is not used more than once
        DO indxBypass1=1,indxBypass-1
            IF (ID .EQ. Bypasses(indxBypass1)%ID) THEN
                CALL SetLastMessage('Bypass ID '//TRIM(IntToText(ID))//' is used more than once!',iFatal,ThisProcedure)
                iStat = -1
                RETURN
            END IF
        END DO

        !Make sure that iNode_Exp_ID is a legit stream node
        IF (iNode_Exp_ID .GT. 0) THEN
            CALL ConvertID_To_Index(iNode_Exp_ID,iStrmNodeIDs,pBypass%iNode_Exp)
            IF (pBypass%iNode_Exp .EQ. 0) THEN
                MessageArray(1) = 'For bypass ID '//TRIM(IntToText(ID))//', stream node number ('//TRIM(IntToText(iNode_Exp_ID))//') where the '
                MessageArray(2) = 'bypass originates from is not in the model!'
                CALL SetLastMessage(MessageArray(1:2),iFatal,ThisProcedure)
                iStat = -1
                RETURN
            END IF
        ELSE
            pBypass%iNode_Exp = 0
        END IF

        !Make sure that a rating table is not specified if the bypass is coming from outside the model area
        IF (iNode_Exp_ID .EQ. 0) THEN
            IF (pBypass%iColBypass .EQ. 0) THEN
                MessageArray(1) = 'A rating table for bypass number '//TRIM(IntToText(ID))//' ,which originates '
                MessageArray(2) = 'outside the model area, cannot be specifed!'
                CALL SetLastMessage(MessageArray(1:2),iFatal,ThisProcedure)
                iStat = -1
                RETURN
            END IF
        END IF

        !Make sure destination type is recognized
        IF (.NOT. ANY(pBypass%iDestType.EQ.f_iDestTypes)) THEN
            CALL SetLastMessage('Destination type for bypass number '//TRIM(IntToText(ID))//' is not recognized!',iFatal,ThisProcedure)
            iStat = -1
            RETURN
        END IF

        !Make sure that destination location is simulated
        SELECT CASE (pBypass%iDestType)
            CASE (f_iFlowDest_StrmNode)
                CALL ConvertID_To_Index(iDest_ID,iStrmNodeIDs,pBypass%iDest)
                IF (pBypass%iDest .EQ. 0) THEN
                    CALL SetLastMessage('Stream node '//TRIM(IntToText(iDest_ID))//' that receives water from bypass '//TRIM(IntToText(ID))//' is not in the model!',iFatal,ThisProcedure)
                    iStat = -1
                    RETURN
                END IF
            CASE (f_iFlowDest_Lake)
                CALL ConvertID_To_Index(iDest_ID,iLakeIDs,pBypass%iDest)
                IF (pBypass%iDest .EQ. 0) THEN
                    CALL SetLastMessage('Lake '//TRIM(IntToText(iDest_ID))//' that receives water from bypass '//TRIM(IntToText(ID))//' is not in the model!',iFatal,ThisProcedure)
                    iStat = -1
                    RETURN
                END IF
        END SELECT

        !Check if the exporting stream node is upstream from the importing stream node
        IF (pBypass%iNode_Exp .GT. 0) THEN
            IF (pBypass%iDestType .EQ. f_iFlowDest_StrmNode) THEN
                iReach_Exp = StrmReach_GetReachNumber(pBypass%iNode_Exp,Reaches)  !Stream reach that the exporting stream node belongs to
                iReach_Imp = StrmReach_GetReachNumber(pBypass%iDest,Reaches)      !Stream reach that the importing stream node belongs to

                !If exporting and importing nodes are in the same reach
                IF (iReach_Exp .EQ. iReach_Imp) THEN
                    IF (pBypass%iNode_Exp .GE. pBypass%iDest) THEN
                        CALL SetLastMessage('Exporting stream node for by-pass '//TRIM(IntToText(ID))//' must be upstream from the receiving node!',iFatal,ThisProcedure)
                        iStat = -1
                        RETURN
                    END IF

                !If exporting and importing nodes are in different reaches (checks only if the nodes are in the same reach network)
                ELSE IF (iReach_Exp .GT. iReach_Imp) THEN
                    DEALLOCATE (iReachesUpNetwork,STAT=ErrorCode)
                    CALL StrmReach_GetReaches_InUpstrmNetwork(Reaches,iReach_Exp,iReachesUpNetwork)
                    IF (LocateInList(iReach_Imp,iReachesUpNetwork) .GT. 0) THEN
                        CALL SetLastMessage('Exporting stream node for by-pass '//TRIM(IntToText(ID))//' must be upstream from the receiving node!',iFatal,ThisProcedure)
                        iStat = -1
                        RETURN
                    END IF
                END IF
            END IF
        END IF

        !Make sure there is only one bypass from a node
        DO indxBypass1=1,indxBypass-1
            IF (pBypass%iNode_Exp .EQ. Bypasses(indxBypass1)%iNode_Exp) THEN
                MessageArray(1) = 'There are multiple bypassses defined at stream node '//TRIM(IntToText(iNode_Exp_ID))//'.'
                MessageArray(2) = 'Only one bypass is allowed from a stream node!'
                CALL SetLastMessage(MessageArray(1:2),iFatal,ThisProcedure)
                iStat = -1
                RETURN
            END IF
        END DO

        !Compute and save the points in rating table for bypass flows
        IF (iNum .LT. 0) THEN
            IF (NPoints .LT. 2) THEN
                CALL SetLastMessage('There should be at least 2 rating table points for bypass '//TRIM(IntToText(ID))//'!',iFatal,ThisProcedure)
                iStat = -1
                RETURN
            END IF
            CALL AllocArray(Dummy2DRealArray,NPoints,2,ThisProcedure,iStat)  ;  IF (iStat .EQ. -1) RETURN
            CALL InFile%ReadData(Dummy2DRealArray,iStat)                     ;  IF (iStat .EQ. -1) RETURN
            CALL pBypass%RatingTable%New(NPoints,Dummy2DRealArray(:,1)*FactFlow,Dummy2DRealArray(:,2)*FactBypass,iStat)
            IF (iStat .EQ. -1) RETURN
        END IF

        !Destination region
        SELECT CASE (pBypass%iDestType)
            CASE (f_iFlowDest_Outside)
                !Do nothing
            CASE (f_iFlowDest_StrmNode)
                pBypass%iDestRegion = 0
            CASE (f_iFlowDest_Lake)
                CALL StrmLakeConnector%AddData(f_iBypassToLakeFlow,indxBypass,pBypass%iDest)
        END SELECT

        !Store export and import reach IDs for stream network processing
        iBypassOutReachIDs(indxBypass) = StrmReach_GetReachNumber(pBypass%iNode_Exp,Reaches)
        IF (iBypassOutReachIDs(indxBypass) .GT. 0) iBypassOutReachIDs(indxBypass) = Reaches(iBypassOutReachIDs(indxBypass))%ID
        IF (pBypass%iDestType .EQ. f_iFlowDest_StrmNode) THEN
            iBypassInReachIDs(indxBypass) = StrmReach_GetReachNumber(pBypass%iDest,Reaches)
            IF (iBypassInReachIDs(indxBypass) .GT. 0) iBypassInReachIDs(indxBypass) = Reaches(iBypassInReachIDs(indxBypass))%ID
        ELSE
            iBypassInReachIDs(indxBypass) = 0
        END IF

    END DO

    !Read the recharge zones
    iBypassIDs = Bypasses%ID
    CALL RechargeZone_New(NBypass,iBypassIDs,iElemIDs,'Bypass',InFile,Bypasses%Recharge,iStat)
    IF (iStat .EQ. -1) RETURN

    !Close file
    CALL InFile%Kill()

    !Re-order reaches based on bypasses
    CALL StrmReach_CompileReachNetwork(SIZE(Reaches),Reaches,iStat,iBypassOutReachIDs,iBypassInReachIDs)

    !Free memory
    DEALLOCATE (Dummy2DRealArray , iBypassIDs , iBypassOutReachIDs , iBypassInReachIDs , STAT=ErrorCode)

  END SUBROUTINE Bypass_New

END MODULE

max_stars_repo_name: giorgk/IWFM-SAFE

max_issues_repo_name: giorgk/IWFM-SAFE

max_forks_repo_name: giorgk/IWFM-SAFE


text:
      subroutine put42
      common /sysplt/ mmajx  ,mmajy  ,mminx  ,mminy  ,mxlab  ,mylab  ,
     1                mflg   ,mtype  ,mxa    ,mya    ,mxb    ,myb    ,
     2                mx     ,my     ,mtypex ,mtypey ,xxa    ,yya    ,
     3                xxb    ,yyb    ,xxc    ,yyc    ,xxd    ,yyd    ,
     4                xfactr ,yfactr ,xadd   ,yadd   ,xx     ,yy     ,
     5                mfmtx(3)       ,mfmty(3)       ,mumx   ,mumy   ,
     6                msizx  ,msizy  ,mxdec  ,mydec  ,mxor   ,mop(19),
     7                mname(19)      ,mxold  ,myold  ,mxmax  ,mymax  ,
     8                mxfac  ,myfac  ,modef  ,mf2er  ,mshftx ,mshfty ,
     9                mmgrx  ,mmgry  ,mmnrx  ,mmnry  ,mfrend ,mfrlst ,
     +                mcrout ,mpair1 ,mpair2 ,msblen ,mflcnt ,mjxmin ,
     1                mjymin ,mjxmax ,mjymax ,mnxsto ,mnysto ,mxxsto ,
     2                mxysto ,mprint ,msybuf(360)    ,mncpw  ,minst  ,
     3                mbufa  ,mbuflu ,mfwa(12)       ,mlwa(12)       ,
     4                mipair ,mbprs(16)      ,mbufl  ,munit  ,mbswap ,
     5                small
      mjxmax = max0(mx,mjxmax)
      mjymax = max0(my,mjymax)
      mjxmin = min0(mx,mjxmin)
      mjymin = min0(my,mjymin)
c
c test if increment instruction will work
c
      if (iabs(mx-mxold).gt.mxmax .or. iabs(my-myold).gt.mymax)
     1    go to 101
c
c construct increment instructions
c
      incx = (mx-mxold)/mxfac+160
      incy = (my-myold)/myfac+32+minst*128
c
c put instruction in buffer
c
      mbpair = ior(ishift(incx,8),incy)
      mipair = mipair+1
      mbprs(mipair) = mbpair
      if (mipair .ge. 16) call flushb
      mxold = mx
      myold = my
      return
  101 continue
c
c mx is first half of the instruction as it stands
c
      mbpair = mx
      mipair = mipair+1
      mbprs(mipair) = mbpair
      if (mipair .ge. 16) call flushb
c
c my needs only pen bit
c
      mbpair = my+ishift(minst,15)
      mipair = mipair+1
      mbprs(mipair) = mbpair
      if (mipair .ge. 16) call flushb
      mxold = mx
      myold = my
      return
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
c
c Module to declare dynamic arrays in multiple subroutines:
c
      module geostat

      real,allocatable    :: prop(:)
      integer,allocatable :: var(:,:,:),ntrain(:,:,:,:),nact(:,:,:,:),
     +                       ntry(:,:,:,:,:),num(:,:),ixl(:),iyl(:),
     +                       izl(:),nlag(:),icross(:),icat(:)

      real      VERSION,tol
      integer   test,report,nsim,nx,ny,nz,nr,np,nv,lout,idbg,ldbg,
     +          maxit,ndir,mlag,maxroc
      character datafl*512,trainfl*512,outfl*512,dbgfl*512

      end module
c
c
c
      program main
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C                                                                      %
C Copyright (C) 1996, The Board of Trustees of the Leland Stanford     %
C Junior University.  All rights reserved.                             %
C                                                                      %
C The programs in GSLIB are distributed in the hope that they will be  %
C useful, but WITHOUT ANY WARRANTY.  No author or distributor accepts  %
C responsibility to anyone for the consequences of using them or for   %
C whether they serve any particular purpose or work at all, unless he  %
C says so in writing.  Everyone is granted permission to copy, modify  %
C and redistribute the programs in GSLIB, but only under the condition %
C that this notice and the above copyright notice remain intact.       %
C                                                                      %
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c-----------------------------------------------------------------------
c
c              Post Process an Integer Coded Realization
c              *****************************************
c
c Input image(s) are post-processed with a steepest descent approach
c so that the two-point histogram for specified directions and lags are
c identified to that of a training image.
c
c
c
c INPUT/OUTPUT Variables:
c
c    datafl     file containing realizations to be post-processed
c    trainfl    file containing the training image
c    outfl      output file for post-processed realizations
c    idbg       debugging level (0-3)
c    report     after ``report'' loops the image is written to the
c               output file
c    dbgfl      output file for debugging information
c    maxit      maximum number of iterations over all nodes
c    tol        tolerance for the objective function (starts at 1)
c    nsim       number of simulations to post-process
c    nx,ny,nz   grid size for the input realizations AND training image
c    ndir       the number of directions to consider for the two-point
c               histogram control
c    ixl(),iyl(),izl()  the offsets for each direction
c
c
c AUTHOR: Clayton V. Deutsch                             DATE: 1989-1999
c-----------------------------------------------------------------------
      use      geostat
      VERSION = 2.905
c
c Read the Parameter File, training image, and initialize random image:
c
      call readparm
c
c LOOP over all the simulations:
c
      do isim=1,nsim
            call initgrid
            call tphist(1)
            call sim
            call tphist(1)
      end do
c
c Finished:
c
      close(lout)
      close(ldbg)
      write(*,9998) VERSION
 9998 format(/' ANNEAL Version: ',f5.3, ' Finished'/)
      stop
      end

      subroutine readparm
c-----------------------------------------------------------------------
c
c                  Initialization and Read Parameters
c                  **********************************
c
c The input parameters are read from a file name provided from standard
c input (a default name will be tried if none is keyed in by the user).
c
c
c
c-----------------------------------------------------------------------
      use       msflib
      use       geostat
c
c ACORN parameters:
c
      parameter   (KORDEI=12,MAXOP1=KORDEI+1,MAXINT=2**30)
      common /iaco/   ixv(MAXOP1)
      real*8    acorni

      character title*80,str*512
      logical   testfl
      lin  = 1
      lout = 2
      ldbg = 3
      linp = 10
c
c Note VERSION number:
c
      write(*,9999) VERSION
 9999 format(/' ANNEAL Version: ',f5.3/)
c
c Get the name of the parameter file - try the default name if no input:
c
      do i=1,512
            str(i:i) = ' '
      end do
      call getarg(1,str)
      if(str(1:1).eq.' ')then
            write(*,*) 'Which parameter file do you want to use?'
            read (*,'(a)') str
      end if
      if(str(1:1).eq.' ') str(1:20) = 'anneal.par          '
      inquire(file=str,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR - the parameter file does not exist,'
            write(*,*) '        check for the file and try again  '
            write(*,*)
            if(str(1:20).eq.'anneal.par          ') then
                  write(*,*) '        creating a blank parameter file'
                  call makepar
                  write(*,*)
            end if
            stop
      endif
      open(lin,file=str,status='OLD')
c
c Find Start of Parameters:
c
 1    read(lin,'(a4)',end=97) str(1:4)
      if(str(1:4).ne.'STAR') go to 1
c
c Read Input Parameters:
c
      read(lin,'(a512)',err=97) datafl
      call chknam(datafl,512)
      write(*,*) ' input images = ',datafl(1:40)

      read(lin,'(a512)',err=97) trainfl
      call chknam(trainfl,512)
      write(*,*) ' training image = ',trainfl(1:40)

      read(lin,*,err=97) idbg,report
      write(*,*) ' idbg, report = ',idbg,report

      read(lin,'(a512)',err=97) dbgfl
      call chknam(dbgfl,512)
      write(*,*) ' debug file = ',dbgfl(1:40)

      read(lin,'(a512)',err=97) outfl
      call chknam(outfl,512)
      write(*,*) ' output images = ',outfl(1:40)

      read(lin,*,err=97) nsim
      write(*,*) ' # of simulations = ',nsim

      read(lin,*,err=97) nx,ny,nz
      write(*,*) ' nx,ny,nz = ',nx,ny,nz

      read(lin,*,err=97) ixv(1)
      write(*,*) ' random number seed',ixv(1)

      read(lin,*,err=97) maxit,tol
      write(*,*) ' maxit,tol = ',maxit,tol

      read(lin,*,err=97) ndir
      write(*,*) ' ndir = ',ndir

      allocate (ixl(ndir),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
      allocate (iyl(ndir),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
      allocate (izl(ndir),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
      allocate (nlag(ndir),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
      allocate (icross(ndir),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if

      mlag = 0
      do i=1,ndir
            read(lin,*,err=97)  ixl(i),iyl(i),izl(i),
     +                              nlag(i)
            icross(i) = 1
            if(nlag(i).gt.mlag) mlag = nlag(i)
            write(*,120) i,ixl(i),iyl(i),izl(i),nlag(i),icross(i)
      end do
      close(lin)
 120  format(' Direction',i2,' ixl iyl izl: ',3i3,' nlag: ',i3,
     +       ' icross: ',i2)
c
c ACORN Initialization (if the seed is chosen small then the first few
c                       random numbers are always small - loop KORDEI)
c
      do i=2,KORDEI+1
            ixv(i) = 0
      end do
      do i=1,KORDEI*KORDEI
            randnu = real(acorni(idum))
      end do
c
c Open the debugging and output files:
c
      open(ldbg,file=dbgfl,status='UNKNOWN')
      open(lout,file=outfl,status='UNKNOWN')
      title = 'ANNEAL SIMULATIONS:                      '//
     +        '                                         '
      inquire(file=trainfl,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR file ',trainfl(1:40),' does not exist!'
            write(*,*) '  you need a training image'
            stop
      end if
c
c Allocate the needed memory:
c
      allocate (var(nx,ny,nz),stat = test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
c
c Read through training data file:
c
      ncat = 0
      open(lin,file=trainfl,status='UNKNOWN')
      read(lin,*,err=98)
      read(lin,*,err=98) nvari
      do i=1,nvari
            read(lin,*,err=98)
      end do
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        read(lin,*,err=98) tempvar
                        do i=1,ncat
                              if(tempvar.eq.var(i,1,1)) go to 71
                        end do
                        ncat = ncat + 1
                        var(i,1,1) = tempvar
 71                     continue
                  end do
            end do
      end do
      close(lin)
c
c Finish allocation:
c
      allocate (prop(ncat+1),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
      allocate (ntrain(ncat+1,ncat+1,ndir,mlag),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
      allocate (nact(ncat+1,ncat+1,ndir,mlag),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
      allocate (ntry(ncat,ncat+1,ncat+1,ndir,mlag),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
      allocate (num(ndir,mlag),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
      allocate (icat(ncat+1),stat=test)
      if(test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
c
c Initialize the grid:
c
      maxroc = ncat + 1
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        var(ix,iy,iz) = maxroc
                  end do
            end do
      end do
c
c Open the file with the input images and read off the header:
c
      open(linp,file=datafl,status='UNKNOWN')
      read(linp,*)
      read(linp,*) nvari
      do i=1,nvari
            read(linp,*)
      end do

c
c Read through training data file:
c
      ncat = 0
      open(lin,file=trainfl,status='UNKNOWN')
      read(lin,'(a60)',err=98) title(21:80)
      read(lin,*,err=98) nvari
      do i=1,nvari
            read(lin,*,err=98)
      end do
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        read(lin,*,err=98) tempvar
                        icode = int(tempvar+0.5)
                        do i=1,ncat
                              if(icode.eq.icat(i)) go to 72
                        end do
                        ncat = ncat + 1
                        icat(ncat)    = icode
 72                     var(ix,iy,iz) = icode
                  end do
            end do
      end do
      close(lin)
c
c Initialize the total proportion and the missing value code:
c
      do i = 1,maxroc
            prop(i) = 0.0
      end do
      tprop   = 0.0
      imiss   = 1
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        ind   = var(ix,iy,iz)
                        icode = imiss
                        do i=1,ncat
                              if(ind.eq.icat(i)) icode = i
                        end do
                        if(icode.ge.1.and.icode.le.maxroc) then
                              prop(icode) = prop(icode) + 1.0
                              tprop       = tprop + 1
                        endif
                        var(ix,iy,iz) = icode
                  end do
            end do
      end do
c
c Number of rock types, compute two point histogram of training image:
c
      nr   = ncat
      nxyz = nx * ny * nz
      write(*,100) nr,(i,icat(i),prop(i),i=1,nr)
 100  format('There are ',i2,' rock types:',/,10
     +        ('   ',i2,' code: ',i3,' there are : ',f10.0,/))
      tprop = 1.0 / max(tprop,1.0)
      oldpro = 0.0
      do i=1,nr
            prop(i) = oldpro + prop(i)*tprop
            oldpro  = prop(i)
      end do
      call tphist(0)
c
c Write a header on the output file and return:
c
      write(lout,101) title
 101  format(a80)
      write(lout,102) nvar,nx,ny,nz
 102  format(4(1x,i4))
      write(lout,103)
 103  format('simulation')

      return
 97   stop 'ERROR in parameter file!'
 98   stop 'ERROR in distribution file!'
      end

      subroutine initgrid
c-----------------------------------------------------------------------
c
c                      Initialization of Grid
c                      **********************
c
c
c
c-----------------------------------------------------------------------
      use      geostat
      data linp/10/
c
c Read in a new initial image:
c
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        read(linp,*) tempvar
                        k = int(tempvar+0.5)
                        do i=1,nr
                              if(k.eq.icat(i)) icode = i
                        end do
                        var(ix,iy,iz) = icode
                  end do
            end do
      end do
      return
      end

      subroutine sim
c-----------------------------------------------------------------------
c
c               MAIN ROUTINE TO PERFORM THE SIMULATION
c
c
c
c
c-----------------------------------------------------------------------
      use      geostat
      real*8   acorni
c
c ACORN parameters:
c
      parameter   (KORDEI=12,MAXOP1=KORDEI+1,MAXINT=2**30)
      common /iaco/   ixv(MAXOP1)
c
c Initialize the starting two point histogram, and objective funtion:
c
      current = object(0)
      objrsc  = 1.0/max(current,0.000001)
      current = objrsc*current
c
c The random path generation procedure of Srivastava:
c
      nxy  = nx * ny
      nxyz = nx * ny * nz
      modlus = 1
 1    modlus = modlus * 2
      if(modlus.le.nxyz) go to 1
c
c MAIN LOOP TO MAXIMUM NUMBER OF ITERATIONS:
c
      do 2 iloop=1,maxit
            write(*,100) iloop,current
 100        format(' Working on ',i2,' objective function: ',f12.9)
            if(current.lt.tol) go to 31
c
c MAIN LOOP OVER ALL THE NODES (unless we're done?):
c
            index = int(real(acorni(idum))*nxyz+1)
            do 3 in=1,nxyz
c
c Figure out the location of this node:
c
            iz = int((index-1)/nxy) + 1
            iy = int((index-(iz-1)*nxy-1)/nx) + 1
            ix = index - (iz-1)*nxy - (iy-1)*nx
c
c Compute the two point histogram and objective function for rock types:
c
            call update(ix,iy,iz)
            best  = current
            ibest = var(ix,iy,iz)
            do ir=1,nr
                  if(ir.ne.var(ix,iy,iz)) then
                        objtry = objrsc*object(ir)
                        if(objtry.lt.best) then
                              best  = objtry
                              ibest = ir
                        endif
                  endif
            end do
c
c Decide on the best and keep it:
c
            if(ibest.ne.var(ix,iy,iz)) then
                  var(ix,iy,iz) = ibest
                  current = best
                  do ir=1,maxroc
                     do jr=1,maxroc
                        do id=1,ndir
                           do il=1,nlag(id)
                            nact(ir,jr,id,il) = ntry(ibest,ir,jr,id,il)
                           end do
                        end do
                     end do
                  end do
            endif
c
c Get a new node to consider:
c
 20         index = mod(5*index+1,modlus)
            if(index.gt.nxyz.or.index.lt.1) go to 20
c
c END MAIN LOOP OVER NODES:
c
 3          continue
c
c Write this realization to the output file?
c
 31         if((int(iloop/report)*report).eq.iloop.or.
     +            iloop.eq.maxit.or.current.lt.tol) then
                  do iz=1,nz
                        do iy=1,ny
                              do ix=1,nx
                                    ind = var(ix,iy,iz)
                                    write(lout,'(i1)')  icat(ind)
                              end do
                        end do
                  end do
                  if(current.lt.tol) return
            endif
c
c END MAIN LOOP OVER MAXIMUM ITERATIONS:
c
 2    continue
      return
      end

      subroutine update(ix,iy,iz)
c-----------------------------------------------------------------------
c
c       Considering a node - decide on the best value to keep
c
c
c
c
c-----------------------------------------------------------------------
      use      geostat
c
      inow = var(ix,iy,iz)
c
c Initialize:
c
      do irtry=1,nr
            do ir=1,maxroc
               do jr=1,maxroc
                  do id=1,ndir
                     do il=1,nlag(id)
                        ntry(irtry,ir,jr,id,il) = nact(ir,jr,id,il)
                     end do
                  end do
               end do
            end do
      end do
c
c MAIN LOOP OVER ALL ROCK TYPES:
c
      do 10 ir=1,nr
            if(ir.eq.inow) go to 10
c
c Consider the change to all lags and directions:
c
            do 3 id=1,ndir
            do 4 il=1,nlag(id)

c
c      Update the positive lag:
c
            jx = ix + il*ixl(id)
            jy = iy + il*iyl(id)
            jz = iz + il*izl(id)
            if(jx.gt.nx.or.jy.gt.ny.or.jz.gt.nz) go to 40
            if(jx.lt. 1.or.jy.lt. 1.or.jz.lt. 1) go to 40
            j  = var(jx,jy,jz)
            ntry(ir,inow,j,id,il) = ntry(ir,inow,j,id,il) - 1
            ntry(ir,ir,  j,id,il) = ntry(ir,ir,  j,id,il) + 1
c
c      Update the negative lag:
c
 40         jx = ix - il*ixl(id)
            jy = iy - il*iyl(id)
            jz = iz - il*izl(id)
            if(jx.gt.nx.or.jy.gt.ny.or.jz.gt.nz) go to 4
            if(jx.lt. 1.or.jy.lt. 1.or.jz.lt. 1) go to 4
            j  = var(jx,jy,jz)
            ntry(ir,j,inow,id,il) = ntry(ir,j,inow,id,il) - 1
            ntry(ir,j,ir,  id,il) = ntry(ir,j,ir,  id,il) + 1

 4          continue
 3          continue
c
c END loop over rock types:
c
 10   continue
      return
      end

      subroutine tphist(iflag)
c-----------------------------------------------------------------------
c
c Compute a two point histogram of an array "var" that is nx by ny by nz
c (dimensioned for MAXX by MAXY by MAXZ), for "nr" rock types, "nlag"
c lags, "ndir" directions specified by integer arrays "ixl", "iyl",
c "izl".  Put the result in the "tp" array that is dimensioned MP by MV
c
c
c-----------------------------------------------------------------------
      use      geostat
c
c Initialize:
c
      do ir=1,maxroc
         do jr=1,maxroc
            do id=1,ndir
               do il=1,nlag(id)
                  if(iflag.eq.0) then
                        ntrain(ir,jr,id,il) = 0
                  else
                        nact  (ir,jr,id,il) = 0
                  endif
                  num(id,il) = 0
               end do
            end do
         end do
      end do
c
c Calculate the Experimental two point histogram:
c
      do 2 iz=1,nz
      do 2 iy=1,ny
      do 2 ix=1,nx
c
c Consider the first value and all directions and lags:
c
            i = var(ix,iy,iz)
            do 3 id=1,ndir
            do 3 il=1,nlag(id)
                  ii = ix + il*ixl(id)
                  jj = iy + il*iyl(id)
                  kk = iz + il*izl(id)
                  if(ii.gt.nx.or.jj.gt.ny.or.kk.gt.nz) go to 3
                  if(ii.lt. 1.or.jj.lt. 1.or.kk.lt. 1) go to 3
c
c If the second point is within the grid then keep the sample:
c
                  j = var(ii,jj,kk)
                  if(iflag.eq.0) then
                        ntrain(i,j,id,il) = ntrain(i,j,id,il)+1
                  else
                        nact  (i,j,id,il) = nact  (i,j,id,il)+1
                  endif
                  num(id,il) = num(id,il) + 1
 3          continue
 2    continue
c
c Debugging output:
c
      if(idbg.ge.3) then
            do ir=1,nr
            do jr=1,nr
            do id=1,ndir
            do il=1,nlag(id)
                  write(ldbg,100) ir,jr,id,il,ntrain(ir,jr,id,il),
     +                                        nact  (ir,jr,id,il)
            end do
            end do
            end do
            end do
 100        format('ir jr ',2i5,' id il ',2i5,' train act ',2i12)
      endif
c
c Return with the two point histogram:
c
      return
      end

      real function object(iflag)
c-----------------------------------------------------------------------
c
c Compute the objective function given a training two point histogram
c and an experimental two point histogram.
c
c
c-----------------------------------------------------------------------
      use      geostat
      object = 0.0
c
c Loop over all direction vectors and rock codes:
c
      do 1 id=1,ndir
      do 1 il=1,nlag(id)
            if(num(id,il).gt.0) then
                  do 2 ir=1,nr
                  do 2 jr=1,nr
                  if(icross(id).eq.0.and.ir.ne.jr) go to 2
                  if(iflag.eq.0) then
                        actual = nact(ir,jr,id,il)
                  else
                        actual = ntry(iflag,ir,jr,id,il)
                  endif
                  object = object + (actual-ntrain(ir,jr,id,il))
     +                                  * (actual-ntrain(ir,jr,id,il))
     +                                  /    (num(id,il)*num(id,il))
 2                continue
            endif
 1    continue
c
c Return with the objective function:
c
      return
      end

      subroutine makepar
c-----------------------------------------------------------------------
c
c                      Write a Parameter File
c                      **********************
c
c
c
c-----------------------------------------------------------------------
      lun = 99
      open(lun,file='anneal.par',status='UNKNOWN')
      write(lun,10)
 10   format('                  Parameters for ANNEAL',/,
     +       '                  *********************',/,/,
     +       'START OF PARAMETERS:')

      write(lun,11)
 11   format('sisim.out                         ',
     +       '-file with input image(s)')
      write(lun,12)
 12   format('ellipsim.out                      ',
     +       '-file with training image')
      write(lun,13)
 13   format('3   10                            ',
     +       '-debug level, reporting interval')
      write(lun,14)
 14   format('anneal.dbg                        ',
     +       '-file for debug output')
      write(lun,15)
 15   format('anneal.out                        ',
     +       '-file for output simulation')
      write(lun,16)
 16   format('1                                 ',
     +       '-number of realizations')
      write(lun,17)
 17   format(' 25   25  1                       ',
     +       '-nx, ny, nz')
      write(lun,18)
 18   format('69069                             ',
     +       '-random number seed')
      write(lun,19)
 19   format('5   0.000001                      ',
     +       '-maximum iterations, tolerance')
      write(lun,20)
 20   format('4                                 ',
     +       '-number of directions')
      write(lun,21)
 21   format('1  0  0  10                       ',
     +       '-ixl, iyl, izl, nlag')
      write(lun,22)
 22   format('0  1  0  10                       ',
     +       '-ixl, iyl, izl, nlag')
      write(lun,23)
 23   format('1  1  0   5                       ',
     +       '-ixl, iyl, izl, nlag')
      write(lun,24)
 24   format('1 -1  0   5                       ',
     +       '-ixl, iyl, izl, nlag')

      close(lun)
      return
      end

      double precision function acorni(idum)
c-----------------------------------------------------------------------
c
c Fortran implementation of ACORN random number generator of order less
c than or equal to 12 (higher orders can be obtained by increasing the
c parameter value MAXORD).
c
c
c NOTES: 1. The variable idum is a dummy variable. The common block
c           IACO is used to transfer data into the function.
c
c        2. Before the first call to ACORN the common block IACO must
c           be initialised by the user, as follows. The values of
c           variables in the common block must not subsequently be
c           changed by the user.
c
c             KORDEI - order of generator required ( must be =< MAXORD)
c
c             MAXINT - modulus for generator, must be chosen small
c                      enough that 2*MAXINT does not overflow
c
c             ixv(1) - seed for random number generator
c                      require 0 < ixv(1) < MAXINT
c
c             (ixv(I+1),I=1,KORDEI)
c                    - KORDEI initial values for generator
c                      require 0 =< ixv(I+1) < MAXINT
c
c        3. After initialisation, each call to ACORN generates a single
c           random number between 0 and 1.
c
c        4. An example of suitable values for parameters is
c
c             KORDEI   = 10
c             MAXINT   = 2**30
c             ixv(1)   = an odd integer in the (approximate) range
c                        (0.001 * MAXINT) to (0.999 * MAXINT)
c             ixv(I+1) = 0, I=1,KORDEI
c
c
c
c Author: R.S.Wikramaratna,                           Date: October 1990
c-----------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      parameter (KORDEI=12,MAXOP1=KORDEI+1,MAXINT=2**30)
      common/iaco/ ixv(MAXOP1)
      do i=1,KORDEI
            ixv(i+1)=(ixv(i+1)+ixv(i))
            if(ixv(i+1).ge.MAXINT) ixv(i+1)=ixv(i+1)-MAXINT
      end do
      acorni=dble(ixv(KORDEI+1))/MAXINT
      return
      end

ext: for


text: Identity Profile removed because Facebook STALKERS BE Pest Control BUGGIN

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE MFFTOM(C,IMS,IVS,IES,NM,NV,NE,INDEX,ITEMP)
*
*   PURPOSE :
*     THIS ROUTINE PERFORMS A REORDERING OF A VECTOR-OF-2 VECTORS
*     OF COMPLEX C[IMS,NM [IVS,NV [IES,NE]]], ACCORDING TO A
*     PERMUTATION INDEX "INDEX".
*     SEE REF.[1] FOR NOTATIONS, AND COMMENTS TO MFFTDM ABOVE.
*  ARGUMENTS
*     C : VECTOR-OF-2 VECTORS TO BE REORDERED
*     IMS,IVS,IES,NM,NV,NE: THESE ARGUMENTS DESCRIBE THE STRUCTURE OF
*     C, ACCORDING TO THE ABOVE DEFINITION;
*     INDEX: INTEGER ARRAY , CONTAINING THE PERMUTATION INDEX  ;
*            IT IS NM ELEMENTS LONG ; PREPARED BY MFFTP.
*     IWORK: INTEGER ARRAY, OF LENGTH AT LEAST NM, USED AS WORKSPACE;
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER INDEX(0:*),ITEMP(0:*)
      DOUBLE COMPLEX C(0:IMS-1,0:*),T

      I3LIM=(NV-1)*IVS
      JLIM=(NE-1)*IES
      DO 1 I=0,NM-1
1     ITEMP(I)=INDEX(I)

      DO 4 I=1,NM-3
2       IF(ITEMP(I).NE.I)THEN
          IDEST=ITEMP(I)
          DO 3 I3=0,I3LIM,IVS
          DO 3 J=I3,I3+JLIM,IES
            T=C(J,I)
            C(J,I)=C(J,IDEST)
            C(J,IDEST)=T
3         CONTINUE
          ITEMP(I)=ITEMP(IDEST)
          ITEMP(IDEST)=IDEST
        GOTO 2
        ENDIF
4     CONTINUE
      END

max_stars_repo_name: cea-trust-platform/trust-code

max_issues_repo_name: pledac/trust-code

max_forks_repo_name: pledac/trust-code


text:
C -*- Mode: Fortran; -*-
C
C  (C) 2003 by Argonne National Laboratory.
C      See COPYRIGHT in top-level directory.
C
C
C Test user-defined operations.  This tests a simple commutative operation
C
      subroutine uop( cin, cout, count, datatype )
      implicit none
      include 'mpif.h'
      integer cin(*), cout(*)
      integer count, datatype
      integer i

      if (datatype .ne. MPI_INTEGER) then
         print *, 'Invalid datatype (',datatype,') passed to user_op()'
         return
      endif

      do i=1, count
         cout(i) = cin(i) + cout(i)
      enddo
      end

      program main
      implicit none
      include 'mpif.h'
      external uop
      integer ierr, errs
      integer count, sumop, vin(65000), vout(65000), i, size
      integer comm

      errs = 0

      call mtest_init(ierr)
      call mpi_op_create( uop, .true., sumop, ierr )

      comm = MPI_COMM_WORLD
      call mpi_comm_size( comm, size, ierr )
      count = 1
      do while (count .lt. 65000)
         do i=1, count
            vin(i) = i
            vout(i) = -1
         enddo
         call mpi_allreduce( vin, vout, count, MPI_INTEGER, sumop,
     *                       comm, ierr )
C         Check that all results are correct
         do i=1, count
            if (vout(i) .ne. i * size) then
               errs = errs + 1
               if (errs .lt. 10) print *, "vout(",i,") = ", vout(i)
            endif
         enddo
         count = count + count
      enddo

      call mpi_op_free( sumop, ierr )

      call mtest_finalize(errs)
      call mpi_finalize(ierr)
      end


text:
module lfortran_intrinsic_string
implicit none

contains

integer elemental function len_trim(string) result(r)
character(len=*), intent(in) :: string
r = len(string)
if (r == 0) return
do while(string(r:r) == " ")
    r = r - 1
    if (r == 0) exit
end do
end function

function trim(x) result(r)
character(len=*),intent(in) :: x
character(len=len_trim(x)) :: r
! This does not work yet in LFortran:
!r = x(1:len(r))
! So we do this workaroud that works:
integer :: i
do i = 1, len(r)
    r(i:i) = x(i:i)
end do
end function

integer elemental function index(string_, substring_) result(idx)
character(len=*), intent(in) :: string_
character(len=*), intent(in) :: substring_
integer :: i, j, k, pos
logical :: found
found = .true.
idx = 0
i = 1
do while (i < len(string_) .and. found)
    k = 0
    j = 1
    do while (j < len(substring_) .and. found)
        pos = i + k
        if( string_(pos:pos) /= substring_(j:j) ) then
            found = .false.
        end if
        k = k + 1
        j = j + 1
    end do
    if( found ) then
        idx = i
        found = .false.
    else
        found = .true.
    end if
    i = i + 1
end do
end function

end module

max_stars_repo_name: lfortran/lfortran

max_issues_repo_name: lfortran/lfortran

max_forks_repo_name: lfortran/lfortran


text:
      SUBROUTINE eqwidtb(lmin,lmax,eqwidth)
*
*-----------------------------------------------------------------------
*
* eqwidtb computes equivalent width for eqwidt (in spherical geometry)
*      06/07-2000 BPz
*
*-----------------------------------------------------------------------
*
      INCLUDE 'spectrum.inc'
*
      INTEGER    MAXNL
      PARAMETER (MAXNL=200)
      INTEGER    NLBL
      PARAMETER (NLBL=100)
      parameter (imax=1000)
*
      real pos(imax),intens(imax),mum
      character*50 mcode
      CHARACTER*20  LELE(100)
      DIMENSION ION(100),SRXLAL(2*MAXNL),
     &          GFELOG(100),ABUND(100),ETAK(2*MAXNL),ETAR(NDP*2*MAXNL),
     &          ETAD(NDP),XC(NDP),CHIE(100),UTREC(103)
      real eqwidth
      doubleprecision XL1,XL2
      doubleprecision XLBEG,XLEND,DEL,XLB(100),XL(100),XLAL(2*MAXNL)
      COMMON/ATMOS/ T(NDP),PE(NDP),PG(NDP),XI(NDP),MUM(NDP),RO(NDP),
     & nnNTAU
*
* Special for spherical
*
      COMMON /CTRAN/X(NDP),S(NDP),BPLAN(NDP),XJ(NDP),HFLUX(NDP),XK(NDP)
     & ,FJ(NDP),SOURCE(NDP),TAUS(NDP),DTAUS(NDP),JTAU0,JTAU1,ISCAT
      COMMON /TAUC/TAU(NDP),DTAULN(NDP),NTAU
      COMMON /ROSSC/ROSS(NDP),CDUMM(NDP)
      COMMON /RHOC/RHO(NDP)
      COMMON /CSPHER/NCORE,DIFLOG,RADIUS,RR(NDP)
      COMMON /CSTYR/MIHAL  /CTAUM/TAUM
      COMMON /SPACE2/ SPACEDUM1(NDP*7+NDP*NRAYS*5+NRAYS*3+1),
     &                MMU(NDP),SPACEDUM2(NDP*2),PFEAU(NRAYS,NDP),
     &                XMU(NRAYS,NDP)
      COMMON /TRDBUG/IDEBUG
      common /limbdk/ pos,intens,totintens,tottrans
      LOGICAL DEBUG
      DIMENSION FLUXME(nlbl)
      DIMENSION Y1C(NRAYS),Y1L(NRAYS,100),XMUC(NRAYS)
      dimension fcfc(lpoint),y1cy1c(nrays,lpoint),
     & xmucxmuc(nrays,lpoint),mmucmmuc(lpoint),xlm(lpoint),
     & jlcont(lpoint)
*
      COMMON/CSURF/ HSURF,Y1(NRAYS)
CCC      COMMON/CANGLE/ NMY,XMY(6),XMY2(6),WMY(6)
      COMMON/PIECES/ XL1,XL2,DEL,EPS,NMX,NLBLDU,IINT,XMYC,IWEAK
*
* extension for large number of wavelengths and lines (monster II)
      doubleprecision xlambda
      common/large/ xlambda(lpoint),maxlam,ABSO(NDP,lpoint),
     & absos(ndp,lpoint),absocont(ndp,lpoint)
*
      character*256 filterfil
      character*80 filttitle
      logical limbdark,first,findtau1
      common/filter/limbdark,ifilt,filtlam(1000),filttrans(1000),
     &              filterfil,filttitle
      logical hydrovelo
      real velocity
      common/velo/velocity(ndp),hydrovelo
*
*
* NLBLDU is a dummy, real val. of NLBL is set as parameter
*
      DATA first/.true./
      data debug/.false./
      PI=3.141593
*
* Initiate mode of calculation
* IINT =1  : intensity at all mu-points
* IINT =0  : flux
* XL1      : wavelength where synthetic spectrum starts
* XL2      : wavelength where synthetic spectrum stops
* DEL      : wavelength step
* IWEAK =1 : weak line approximation for L/KAPPA le. EPS
* note  XL2.gt.XL1
*
      if (first) then
        iweak=0
        IF(IINT.EQ.0) WRITE(7,201) XL1,XL2,DEL
        IF(IWEAK.GT.0) WRITE(7,202) EPS
        first=.false.
      endif
*
* Read model atmosphere
*
      REWIND 14
      READ(14) MCODE,nlcont,xlm(1),BPLAN,XC,S,XI
      REWIND 14
      WRITE(7,203) MCODE(1:lenstr(mcode))
*
* Continuum calculations:
*
      do 1963 jc=1,nlcont
        READ(14) MCODE,idum,xlm(jc),BPLAN,XC,S,XI
        DO 9 K=1,NTAU
          X(K)=XC(K)
    9   CONTINUE
        if (debug) then
          print*,'eqwidtb: calling traneq for continuum',jc,xlm(jc)
          do k=1,ntau,5
            print*,(x(kk),s(kk),kk=k,k+4)
          enddo
        endif
        CALL TRANEQ
*
* Spherical fluxes
*
C
C FLUX TO PRINT
        HFLUX1C=4.*PI*HSURF*(RR(1)/RADIUS)**2
ccc        print*,hflux1c
        HFLUX2C=4.*PI*HFLUX(NTAU)*(RR(NTAU)/RADIUS)**2
        FCFC(JC)=HFLUX1C/PI
C        GFLUX1C=HFLUX1/PI*XL(J)**2/CLIGHT
C        GFLUX2C=HFLUX2/PI*XL(J)**2/CLIGHT
C        FFLUX1C=-2.5*ALOG10(AMAX1(GFLUX1,1E-20))
C        FFLUX2C=-2.5*ALOG10(AMAX1(GFLUX2,1E-20))
        DO 900 MM=1,MMU(1)
          Y1CY1C(MM,jc)=Y1(MM)
          XMUCXMUC(MM,jc)=XMU(MM,1)
 900    CONTINUE
        MMUCMMUC(jc)=MMU(1)
C
*
* End of spherical fluxes
*
        if (debug) print*,jc,xlm(jc),fcfc(jc)
        IF(IINT.LE.0) WRITE(7,204) fcFC(jc),xlm(jc)
CC      IF(IINT.GT.0) WRITE(7,205) Y1C
1963  continue
*
* initialize
      eqwidth=0.e0
      profold=0.e0
*
* Now solve transfer equation in the wavelength range around the line
*
      DO 39 j=lmin,lmax
        xlsingle=xlambda(j)
        DO 30 K=1,NTAU
* the continuum opacity is already included in abso
ccc          X(K)=XC(K)+ABSO(K,J)
          X(K)=ABSO(K,J)
          s(k)=absos(k,j)
          BPLAN(k)=BPL(T(k),xlsingle)
   30   CONTINUE
        CALL TRANEQ
*
* Spherical fluxes
*
C
C FLUX TO PRINT
        HFLUX1=4.*PI*HSURF*(RR(1)/RADIUS)**2
        HFLUX2=4.*PI*HFLUX(NTAU)*(RR(NTAU)/RADIUS)**2
ccc        FLUXME(J)=HFLUX1/PI
C        GFLUX1=HFLUX1/PI*XL(J)**2/CLIGHT
C        GFLUX2=HFLUX2/PI*XL(J)**2/CLIGHT
C        FFLUX1=-2.5*ALOG10(AMAX1(GFLUX1,1E-20))
C        FFLUX2=-2.5*ALOG10(AMAX1(GFLUX2,1E-20))
C
* interpolate continuum flux
        do 3691 jc=2,nlcont
          if(xlm(jc)-xlsingle.gt.0.) then
            jjc=jc-1
            goto 3692
          endif
3691    continue
3692    continue
        if (nlcont.gt.1) then
          jjc=min(jjc,nlcont-1)
          fc=(fcfc(jjc+1)-fcfc(jjc))/
     &       (xlm(jjc+1)-xlm(jjc))*(xlsingle-xlm(jjc)) + fcfc(jjc)
        else
          fc=fcfc(1)
        endif
        PRF=hflux1/pi/FC
        if (debug) print*,j,xlsingle,hflux1/pi,prf
        IF(IINT.GT.0) THEN
          DO 901 MM=1,MMU(1)
            Y1L(MM,J)=Y1(MM)
 901      CONTINUE
        ENDIF
*
* End of spherical fluxes
*
        PROF=1.-PRF
        eqwidth=eqwidth+(prof+profold)*del/2.
        profold=prof
   39 CONTINUE
*
ccc      print 1111,eqwidth*1000.
1111  format(' eqwidth: ',f9.1,' mA'/'warning! flux calculation only!')
      call clock
      RETURN
*
  100 FORMAT(4X,I1,6X,I3)
  207 FORMAT(' SPECTRUM CALCULATED FOR THE FOLLOWING ',I3,' LINES'
     &      /' ELEMENT LAMBDA XI(EV) LOG(GF) LOG(ABUND) LINE NO')
  208 FORMAT('   ',A2,I2,F9.2,F5.2,1X,F6.2,3X,F7.2,4X,I5)
  200 FORMAT(' INTENSITY SPECTRUM AT MU=',F6.2,' BETWEEN',F10.3,
     &      ' AND',F10.3,' A WITH STEP',F6.3,' A')
  201 FORMAT(' FLUX SPECTRUM BETWEEN',F10.3,' A AND ',F10.3,' A WITH',
     &       ' STEP',F6.3,' A')
  202 FORMAT(' ***WEAK LINE APPROX FOR L/KAPPA L.T.',F7.4)
  203 FORMAT(' MODEL ATMOSPHERE:',A)
  204 FORMAT(' CONTINUUM FLUX=',E12.4, ' AT ',f10.2,'AA')
  205 FORMAT(' CONTINUUM INTENSITY=',E12.4)
  206 FORMAT(1X,10F8.4)
 2341 FORMAT(1X,'Block number: ',I8)
*
      END

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee


text:
The Muslim Student Association at the University of California Davis is a Student Organizations student Spiritual Organizations spiritual organization united for the sake of God. The mission of the MSA is to create an avenue for Muslim students to meet, exchange ideas, debunk misconceptions, and work to understand and disseminate the true message of Islam as in accordance with the Quran and the life of the Prophet Muhammad. The club is open to members of all religions.

The MSA hosts various events throughout the school year including general body meetings, a fast for a day event geared towards the nonMuslim members of the community, lectures by prominent Muslim speakers, and more.

The Muslim Student Association is an active participant in the Interfaith Campus Council

 Media Coverage
  The California Aggie about http://theaggie.org/article/2010/04/05/muslimstudentscelebrateannualislamawarenessmonth Islam Awareness Month.
  The California Aggie about http://www.californiaaggie.com/news/2008/01/09/CampusNews/UcdMuslim.Student.Association.Sponsors.A.Holiday.Celebration3148011.shtml Eid Al Adha Banquet.

  There are multiple groups in town called the Muslim Student Association.  See Muslim Student Association the disambiguation page for others

20060317 20:46:25 nbsp Awesome, I was completely unaware that UCD MSA was intune with the wonders of wiki... Users/AroojAhmad

20070831 11:48:41 nbsp I love MSA !!! ~ YOUR SISTER !  Users/Amna


!
!  Copyright 2019-2020 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!

#include "config.h"

module hamiltonian
  implicit none
  integer,private,parameter :: Nd = 4

contains

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

SUBROUTINE hpsi(tpsi,htpsi,info,mg,V_local,system,stencil,srg,ppg,ttpsi)
  use structures
  use stencil_sub
  use nonlocal_potential
  use pseudo_pt_plusU_sub, only: pseudo_plusU, PLUS_U_ON
  use pseudo_pt_so_sub, only: pseudo_so, SPIN_ORBIT_ON
  use nondiagonal_so_sub, only: nondiagonal_so
  use sendrecv_grid, only: s_sendrecv_grid, update_overlap_real8, update_overlap_complex8
  use salmon_global, only: yn_want_communication_overlapping,yn_periodic,yn_jm,yn_symmetrized_stencil
  use timer
  use code_optimization, only: stencil_is_parallelized_by_omp
  use communication, only: comm_summation
  implicit none
  type(s_dft_system)      ,intent(in) :: system
  type(s_parallel_info),intent(in) :: info
  type(s_rgrid)  ,intent(in) :: mg
  type(s_scalar) ,intent(in) :: V_local(system%Nspin)
  type(s_stencil),intent(in) :: stencil
  type(s_sendrecv_grid),intent(inout) :: srg
  type(s_pp_grid),intent(in) :: ppg
  type(s_orbital)            :: tpsi,htpsi
  type(s_orbital),optional   :: ttpsi
  !
  integer :: nspin,ispin,io,ik,im,im_s,im_e,ik_s,ik_e,io_s,io_e,norb,ix,iy,iz
  real(8) :: k_nabt(Nd,3),k_lap0,kAc(3)
  logical :: if_kAc,if_singlescale
  logical :: is_enable_overlapping
  !real(8) :: tmp,tmp1

  call timer_begin(LOG_UHPSI_ALL)

  im_s = info%im_s
  im_e = info%im_e
  ik_s = info%ik_s
  ik_e = info%ik_e
  io_s = info%io_s
  io_e = info%io_e
  nspin = system%nspin
  norb = Nspin* info%numo * info%numk * info%numm

  if_kAc = (yn_periodic=='y')
  if_singlescale = allocated(system%Ac_micro%v)

  ! check: can we execute computation/communication overlapping
  if (if_singlescale) then
    is_enable_overlapping = (yn_want_communication_overlapping == 'y') .and. &
                            stencil%if_orthogonal .and. &
                            info%if_divide_rspace
  else
    is_enable_overlapping = (yn_want_communication_overlapping == 'y') .and. &
                            stencil%if_orthogonal .and. &
                            info%if_divide_rspace .and. &
                            (im_e - im_s) <= 0    .and. &
                            (ik_e - ik_s) <= 0
  end if

  if(allocated(tpsi%rwf)) then

  ! overlap region communication
    call timer_begin(LOG_UHPSI_UPDATE_OVERLAP)
    if(info%if_divide_rspace) then
      call update_overlap_real8(srg, mg, tpsi%rwf)
    end if
    call timer_end(LOG_UHPSI_UPDATE_OVERLAP)

  ! stencil
    call timer_begin(LOG_UHPSI_STENCIL)
    do im=im_s,im_e
    do ik=ik_s,ik_e
    do io=io_s,io_e
    do ispin=1,Nspin
      call dstencil(mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz &
                   ,tpsi%rwf(:,:,:,ispin,io,ik,im),htpsi%rwf(:,:,:,ispin,io,ik,im) &
                   ,V_local(ispin)%f,stencil%coef_lap0,stencil%coef_lap)
    end do
    end do
    end do
    end do
    call timer_end(LOG_UHPSI_STENCIL)

    ! nonlocal potential
    if ( SPIN_ORBIT_ON ) then
      ! pseudopotential
      if(yn_jm=='n') call dpseudo(tpsi,htpsi,info,Nspin,ppg)
    else
      ! pseudopotential
      if(yn_jm=='n') call dpseudo(tpsi,htpsi,info,Nspin,ppg)
    end if

    ! DFT+U
    if ( PLUS_U_ON ) then
      call pseudo_plusU(tpsi,htpsi,system,info,ppg)
    end if

  else

  ! overlap region communication
    call timer_begin(LOG_UHPSI_UPDATE_OVERLAP)
    if(info%if_divide_rspace .and. .not. is_enable_overlapping) then
      call update_overlap_complex8(srg, mg, tpsi%zwf)
    end if
    call timer_end(LOG_UHPSI_UPDATE_OVERLAP)

  ! stencil
    call timer_begin(LOG_UHPSI_STENCIL)
    if(stencil%if_orthogonal .and. .not.if_singlescale) then
    ! orthogonal lattice (general)

      if(stencil_is_parallelized_by_omp .or. is_enable_overlapping) then

        do im=im_s,im_e
        do ik=ik_s,ik_e
          if(if_kAc) then
            kAc(1:3) = system%vec_k(1:3,ik) + system%vec_Ac(1:3)
            k_lap0 = stencil%coef_lap0 + 0.5d0* sum(kAc(1:3)**2)
            k_nabt(:,1) = kAc(1) * stencil%coef_nab(:,1)
            k_nabt(:,2) = kAc(2) * stencil%coef_nab(:,2)
            k_nabt(:,3) = kAc(3) * stencil%coef_nab(:,3)
          else
            k_lap0 = stencil%coef_lap0
            k_nabt = 0d0
          end if
          if (is_enable_overlapping) then
            call zstencil_overlapped
          else
            do io=io_s,io_e
            do ispin=1,Nspin
              call zstencil(mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz &
                            ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                            ,V_local(ispin)%f,k_lap0,stencil%coef_lap,k_nabt)
            end do
            end do
          end if
        end do
        end do

      else
      ! OpenMP parallelization: k-point & orbital indices

!$omp parallel do collapse(4) default(none) &
!$omp          private(im,ik,io,ispin,kAc,k_lap0,k_nabt) &
!$omp          shared(im_s,im_e,ik_s,ik_e,io_s,io_e,nspin,if_kac,system,stencil,mg,tpsi,htpsi,V_local)
        do im=im_s,im_e
        do ik=ik_s,ik_e
        do io=io_s,io_e
        do ispin=1,Nspin
          if(if_kAc) then
            kAc(1:3) = system%vec_k(1:3,ik) + system%vec_Ac(1:3)
            k_lap0 = stencil%coef_lap0 + 0.5d0* sum(kAc(1:3)**2)
            k_nabt(:,1) = kAc(1) * stencil%coef_nab(:,1)
            k_nabt(:,2) = kAc(2) * stencil%coef_nab(:,2)
            k_nabt(:,3) = kAc(3) * stencil%coef_nab(:,3)
          else
            k_lap0 = stencil%coef_lap0
            k_nabt = 0d0
          end if
          call zstencil(mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz &
                            ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                            ,V_local(ispin)%f,k_lap0,stencil%coef_lap,k_nabt)
        end do
        end do
        end do
        end do
!$omp end parallel do

      end if

    else if(stencil%if_orthogonal .and. if_singlescale) then
    ! orthogonal lattice, single-scale Maxwell-TDDFT

      if(yn_symmetrized_stencil=='y') then

!$omp parallel do collapse(4) default(none) &
!$omp          private(im,ik,io,ispin) &
!$omp          shared(im_s,im_e,ik_s,ik_e,io_s,io_e,nspin,mg,tpsi,htpsi,V_local,system,stencil)
        do im=im_s,im_e
        do ik=ik_s,ik_e
        do io=io_s,io_e
        do ispin=1,Nspin
          call zstencil_microAc_symmetrized(mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz &
                        ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                        ,V_local(ispin)%f,system%Ac_micro%v,stencil%coef_lap0 &
                        ,stencil%coef_lap,stencil%coef_nab,system%vec_k(1:3,ik))
        end do
        end do
        end do
        end do
!$omp end parallel do

      else if(stencil_is_parallelized_by_omp .or. is_enable_overlapping) then
        ! OpenMP parallelization: rgrid

        if (is_enable_overlapping) then
          call zstencil_microac_overlapped
        else
          do im=im_s,im_e
          do ik=ik_s,ik_e
          do io=io_s,io_e
          do ispin=1,Nspin
            call zstencil_microAc(mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz &
                          ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                          ,V_local(ispin)%f,system%Ac_micro%v,system%div_Ac%f,stencil%coef_lap0 &
                          ,stencil%coef_lap,stencil%coef_nab,system%vec_k(1:3,ik))
          end do
          end do
          end do
          end do
        end if

      else
        ! OpenMP parallelization: k-point & orbital indices

!$omp parallel do collapse(4) default(none) &
!$omp          private(im,ik,io,ispin) &
!$omp          shared(im_s,im_e,ik_s,ik_e,io_s,io_e,nspin,mg,tpsi,htpsi,V_local,system,stencil)
        do im=im_s,im_e
        do ik=ik_s,ik_e
        do io=io_s,io_e
        do ispin=1,Nspin
          call zstencil_microAc(mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz &
                        ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                        ,V_local(ispin)%f,system%Ac_micro%v,system%div_Ac%f,stencil%coef_lap0 &
                        ,stencil%coef_lap,stencil%coef_nab,system%vec_k(1:3,ik))
        end do
        end do
        end do
        end do
!$omp end parallel do

      end if

    else if(.not.stencil%if_orthogonal) then
    ! non-orthogonal lattice

!$omp parallel do collapse(4) default(none) &
!$omp private(im,ik,io,ispin,kAc,k_lap0) &
!$omp shared(im_s,im_e,ik_s,ik_e,io_s,io_e,nspin,if_kac,system,stencil,mg,tpsi,htpsi,V_local)
      do im=im_s,im_e
      do ik=ik_s,ik_e
      do io=io_s,io_e
      do ispin=1,Nspin
        kAc = 0d0
        k_lap0 = 0d0
        if(if_kAc) then
          kAc(1:3) = system%vec_k(1:3,ik) + system%vec_Ac(1:3) ! Cartesian vector k+A/c
          k_lap0 = stencil%coef_lap0 + 0.5d0* sum(kAc(1:3)**2)
          kAc(1:3) = matmul(system%rmatrix_B,kAc) ! B* (k+A/c)
        end if
          call zstencil_nonorthogonal(mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz &
                                     ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                                     ,V_local(ispin)%f,k_lap0,stencil%coef_lap,stencil%coef_nab,kAc,stencil%coef_F)
      end do
      end do
      end do
      end do
!$omp end parallel do

    end if
    call timer_end(LOG_UHPSI_STENCIL)

  ! subtraction
    call timer_begin(LOG_UHPSI_SUBTRACTION)
    if(present(ttpsi)) then
      if(allocated(tpsi%rwf)) then
        !$omp parallel do collapse(6) default(none) &
        !$omp          private(im,ik,io,ispin,iz,iy,ix) &
        !$omp          shared(im_s,im_e,ik_s,ik_e,io_s,io_e,nspin,mg) &
        !$omp          shared(ttpsi,htpsi,V_local,tpsi)
        do im=im_s,im_e
        do ik=ik_s,ik_e
        do io=io_s,io_e
        do ispin=1,Nspin
          do iz=mg%is(3),mg%ie(3)
          do iy=mg%is(2),mg%ie(2)
          do ix=mg%is(1),mg%ie(1)
            ttpsi%rwf(ix,iy,iz,ispin,io,ik,im) = htpsi%rwf(ix,iy,iz,ispin,io,ik,im) &
                                               - V_local(ispin)%f(ix,iy,iz) * tpsi%rwf(ix,iy,iz,ispin,io,ik,im)
          end do
          end do
          end do
        end do
        end do
        end do
        end do
        !$omp end parallel do
      else
        !$omp parallel do collapse(6) default(none) &
        !$omp          private(im,ik,io,ispin,iz,iy,ix) &
        !$omp          shared(im_s,im_e,ik_s,ik_e,io_s,io_e,nspin,mg) &
        !$omp          shared(ttpsi,htpsi,V_local,tpsi)
        do im=im_s,im_e
        do ik=ik_s,ik_e
        do io=io_s,io_e
        do ispin=1,Nspin
          do iz=mg%is(3),mg%ie(3)
          do iy=mg%is(2),mg%ie(2)
          do ix=mg%is(1),mg%ie(1)
            ttpsi%zwf(ix,iy,iz,ispin,io,ik,im) = htpsi%zwf(ix,iy,iz,ispin,io,ik,im) &
                                               - V_local(ispin)%f(ix,iy,iz) * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
          end do
          end do
          end do
        end do
        end do
        end do
        end do
        !$omp end parallel do
      end if
    end if
    call timer_end(LOG_UHPSI_SUBTRACTION)

  ! nonlocal potential
    if(yn_jm=='n') then
      if ( SPIN_ORBIT_ON ) then
        call nondiagonal_so(tpsi,htpsi,info,mg)
        call pseudo_so(tpsi,htpsi,info,nspin,ppg,mg)
      else
      ! pseudopotential
        call zpseudo(tpsi,htpsi,info,nspin,ppg)
      end if
      if ( PLUS_U_ON ) then
        call pseudo_plusU(tpsi,htpsi,system,info,ppg)
      end if
    end if

  end if

  call timer_end(LOG_UHPSI_ALL)

  return
contains
  subroutine zstencil_overlapped
    use sendrecv_grid, only: srg_pack, srg_communication, srg_unpack, &
                             update_overlap_complex8
    use code_optimization, only: modx,mody,modz,optimized_stencil_is_callable
    use communication, only: comm_proc_null, comm_get_groupinfo
    implicit none
    integer :: igs(3),ige(3)
    integer,parameter :: nyblk=8, nzblk=8
    integer :: ibx,iby,ibz
    integer :: iplane,ibs(3),ibe(3)
    integer :: is(3),ie(3)
    logical :: is_divided(3)
    integer :: myrank,nprocs

    call comm_get_groupinfo(srg%icomm, myrank, nprocs)
    do iplane=1,3
      is_divided(iplane) = srg%neig(1,iplane) /= comm_proc_null .and. &
                           srg%neig(1,iplane) /= myrank
    end do

    is(:) = mg%is(:)
    ie(:) = mg%ie(:)
    do iplane=1,3
      if (is_divided(iplane)) then
        is(iplane) = is(iplane) + 4
        ie(iplane) = ie(iplane) - 4
      end if
    end do

! phase 1. pack halo region
    call timer_begin(LOG_UHPSI_OVL_PHASE1)
    call update_overlap_complex8(srg, mg, tpsi%zwf, srg_pack)
    call timer_end  (LOG_UHPSI_OVL_PHASE1)

! phase 2. halo communication and computation without halo region
    call timer_begin(LOG_UHPSI_OVL_PHASE2)
!$omp parallel default(none) &
!$omp          private(io,ispin,igs,ige,ibx,iby,ibz) &
!$omp          shared(is,ie,ik,im,io_s,io_e,nspin,mg,tpsi,htpsi,V_local,k_lap0,stencil,k_nabt,srg,modx,mody,modz) &
!$omp          shared(optimized_stencil_is_callable)

! halo communication by master thread (tid = 0)
!$omp master
    call timer_begin(LOG_UHPSI_OVL_PHASE2_COMM)
    call update_overlap_complex8(srg, mg, tpsi%zwf, srg_communication)
    call timer_end  (LOG_UHPSI_OVL_PHASE2_COMM)
!$omp end master

! A computation with multi-thread except master thread,
! but master thread can join this loop if the communication completed before computation done.
!$omp do collapse(4) schedule(dynamic,1)
    do io=io_s,io_e
    do ispin=1,Nspin
    do ibz=is(3),ie(3),nzblk
    do iby=is(2),ie(2),nyblk
      igs(3) = ibz ; ige(3) = min(ibz + nzblk - 1, ie(3))
      igs(2) = iby ; ige(2) = min(iby + nyblk - 1, ie(2))
      igs(1) = is(1) ; ige(1) = ie(1)
#ifdef USE_OPT_EXPLICIT_VECTORIZATION
      if (optimized_stencil_is_callable) then
        call zstencil_tuned_seq(mg%is_array,mg%ie_array,mg%is,mg%ie,modx,mody,modz,igs,ige &
                               ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                               ,V_local(ispin)%f,k_lap0,stencil%coef_lap,k_nabt)
      else
#endif
        call zstencil_typical_seq(mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz,igs,ige &
                                 ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                                 ,V_local(ispin)%f,k_lap0,stencil%coef_lap,k_nabt)
#ifdef USE_OPT_EXPLICIT_VECTORIZATION
      end if
#endif
    end do
    end do
    end do
    end do
!$omp end do
!$omp end parallel
    call timer_end  (LOG_UHPSI_OVL_PHASE2)

! phase 3. unpack halo region
    call timer_begin(LOG_UHPSI_OVL_PHASE3)
    call update_overlap_complex8(srg, mg, tpsi%zwf, srg_unpack)
    call timer_end  (LOG_UHPSI_OVL_PHASE3)

    is(:) = mg%is(:)
    ie(:) = mg%ie(:)
! phase 4. computation with halo region
    call timer_begin(LOG_UHPSI_OVL_PHASE4)
!$omp parallel default(none) &
!$omp          firstprivate(is,ie) &
!$omp          private(io,ispin,iplane,igs,ige,ibx,iby,ibz,ibs,ibe) &
!$omp          shared(is_divided,ik,im,io_s,io_e,nspin,mg,tpsi,htpsi,V_local,k_lap0,stencil,k_nabt,srg,modx,mody,modz) &
!$omp          shared(optimized_stencil_is_callable)
    do iplane=1,6

      if (.not. is_divided((iplane+1)/2)) cycle

      ibs(:) = is(:)
      ibe(:) = ie(:)
      select case(iplane)
        case(1) ! update X (up)
          ibs(1) = mg%ie(1) - 4 + 1
          ibe(1) = mg%ie(1)
        case(2) ! update X (down)
          ibs(1) = mg%is(1)
          ibe(1) = mg%is(1) + 4 - 1
        case(3) ! update Y (up)
          ibs(2) = mg%ie(2) - 4 + 1
          ibe(2) = mg%ie(2)
        case(4) ! update Y (down)
          ibs(2) = mg%is(2)
          ibe(2) = mg%is(2) + 4 - 1
        case(5) ! update Z (up)
          ibs(3) = mg%ie(3) - 4 + 1
          ibe(3) = mg%ie(3)
        case(6) ! update Z (down)
          ibs(3) = mg%is(3)
          ibe(3) = mg%is(3) + 4 - 1
      end select

      select case(iplane)
        case(2)
          is(1) = is(1) + 4
          ie(1) = ie(1) - 4
        case(4)
          is(2) = is(2) + 4
          ie(2) = ie(2) - 4
      end select

!$omp do collapse(4) schedule(dynamic,1)
      do io=io_s,io_e
      do ispin=1,Nspin
      do ibz=ibs(3),ibe(3),nzblk
      do iby=ibs(2),ibe(2),nyblk
        igs(3) = ibz ; ige(3) = min(ibz + nzblk - 1, ibe(3))
        igs(2) = iby ; ige(2) = min(iby + nyblk - 1, ibe(2))
        igs(1) = ibs(1) ; ige(1) = ibe(1)
#ifdef USE_OPT_EXPLICIT_VECTORIZATION
        if (optimized_stencil_is_callable) then
          call zstencil_tuned_seq(mg%is_array,mg%ie_array,mg%is,mg%ie,modx,mody,modz,igs,ige &
                                 ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                                 ,V_local(ispin)%f,k_lap0,stencil%coef_lap,k_nabt)
        else
#endif
          call zstencil_typical_seq(mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz,igs,ige &
                                   ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
                                   ,V_local(ispin)%f,k_lap0,stencil%coef_lap,k_nabt)
#ifdef USE_OPT_EXPLICIT_VECTORIZATION
        end if
#endif
      end do
      end do
      end do
      end do
!$omp end do nowait
    end do
!$omp end parallel
    call timer_end  (LOG_UHPSI_OVL_PHASE4)
  end subroutine zstencil_overlapped

  subroutine zstencil_microac_overlapped
    use sendrecv_grid, only: srg_pack, srg_communication, srg_unpack, &
                             update_overlap_complex8
    use communication, only: comm_proc_null, comm_get_groupinfo
    implicit none
    integer :: igs(3),ige(3)
    integer,parameter :: nyblk=8, nzblk=8
    integer :: ibx,iby,ibz
    integer :: iplane,ibs(3),ibe(3)
    integer :: is(3),ie(3)
    logical :: is_divided(3)
    integer :: myrank,nprocs

    call comm_get_groupinfo(srg%icomm, myrank, nprocs)
    do iplane=1,3
      is_divided(iplane) = srg%neig(1,iplane) /= comm_proc_null .and. &
                           srg%neig(1,iplane) /= myrank
    end do

    is(:) = mg%is(:)
    ie(:) = mg%ie(:)
    do iplane=1,3
      if (is_divided(iplane)) then
        is(iplane) = is(iplane) + 4
        ie(iplane) = ie(iplane) - 4
      end if
    end do

! phase 1. pack halo region
    call timer_begin(LOG_UHPSI_OVL_PHASE1)
    call update_overlap_complex8(srg, mg, tpsi%zwf, srg_pack)
    call timer_end  (LOG_UHPSI_OVL_PHASE1)

! phase 2. halo communication and computation without halo region
    call timer_begin(LOG_UHPSI_OVL_PHASE2)
!$omp parallel default(none) &
!$omp          private(ik,im,io,ispin,igs,ige,ibx,iby,ibz) &
!$omp          shared(is,ie,im_s,im_e,ik_s,ik_e,io_s,io_e,nspin,mg,tpsi,htpsi,V_local,k_lap0,stencil,system,srg)

! halo communication by master thread (tid = 0)
!$omp master
    call timer_begin(LOG_UHPSI_OVL_PHASE2_COMM)
    call update_overlap_complex8(srg, mg, tpsi%zwf, srg_communication)
    call timer_end  (LOG_UHPSI_OVL_PHASE2_COMM)
!$omp end master

! A computation with multi-thread except master thread,
! but master thread can join this loop if the communication completed before computation done.
    do im=im_s,im_e
    do ik=ik_s,ik_e
!$omp do collapse(4) schedule(dynamic,1)
    do io=io_s,io_e
    do ispin=1,Nspin
    do ibz=is(3),ie(3),nzblk
    do iby=is(2),ie(2),nyblk
      igs(3) = ibz ; ige(3) = min(ibz + nzblk - 1, ie(3))
      igs(2) = iby ; ige(2) = min(iby + nyblk - 1, ie(2))
      igs(1) = is(1) ; ige(1) = ie(1)
      call zstencil_microAc_typical_seq( &
              mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz,igs,ige &
             ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
             ,V_local(ispin)%f,system%Ac_micro%v,system%div_Ac%f,stencil%coef_lap0 &
             ,stencil%coef_lap,stencil%coef_nab,system%vec_k(1:3,ik))
    end do
    end do
    end do
    end do
!$omp end do nowait
    end do
    end do
!$omp end parallel
    call timer_end  (LOG_UHPSI_OVL_PHASE2)

! phase 3. unpack halo region
    call timer_begin(LOG_UHPSI_OVL_PHASE3)
    call update_overlap_complex8(srg, mg, tpsi%zwf, srg_unpack)
    call timer_end  (LOG_UHPSI_OVL_PHASE3)

    is(:) = mg%is(:)
    ie(:) = mg%ie(:)
! phase 4. computation with halo region
    call timer_begin(LOG_UHPSI_OVL_PHASE4)
!$omp parallel default(none) &
!$omp          firstprivate(is,ie) &
!$omp          private(im,ik,io,ispin,iplane,igs,ige,ibx,iby,ibz,ibs,ibe) &
!$omp          shared(is_divided,im_s,im_e,ik_s,ik_e,io_s,io_e,nspin,mg,tpsi,htpsi,V_local,k_lap0,stencil,system,srg)
    do iplane=1,6

      if (.not. is_divided((iplane+1)/2)) cycle

      ibs(:) = is(:)
      ibe(:) = ie(:)
      select case(iplane)
        case(1) ! update X (up)
          ibs(1) = mg%ie(1) - 4 + 1
          ibe(1) = mg%ie(1)
        case(2) ! update X (down)
          ibs(1) = mg%is(1)
          ibe(1) = mg%is(1) + 4 - 1
        case(3) ! update Y (up)
          ibs(2) = mg%ie(2) - 4 + 1
          ibe(2) = mg%ie(2)
        case(4) ! update Y (down)
          ibs(2) = mg%is(2)
          ibe(2) = mg%is(2) + 4 - 1
        case(5) ! update Z (up)
          ibs(3) = mg%ie(3) - 4 + 1
          ibe(3) = mg%ie(3)
        case(6) ! update Z (down)
          ibs(3) = mg%is(3)
          ibe(3) = mg%is(3) + 4 - 1
      end select

      select case(iplane)
        case(2)
          is(1) = is(1) + 4
          ie(1) = ie(1) - 4
        case(4)
          is(2) = is(2) + 4
          ie(2) = ie(2) - 4
      end select

      do im=im_s,im_e
      do ik=ik_s,ik_e
!$omp do collapse(4) schedule(dynamic,1)
      do io=io_s,io_e
      do ispin=1,Nspin
      do ibz=ibs(3),ibe(3),nzblk
      do iby=ibs(2),ibe(2),nyblk
        igs(3) = ibz ; ige(3) = min(ibz + nzblk - 1, ibe(3))
        igs(2) = iby ; ige(2) = min(iby + nyblk - 1, ibe(2))
        igs(1) = ibs(1) ; ige(1) = ibe(1)
        call zstencil_microAc_typical_seq( &
                mg%is_array,mg%ie_array,mg%is,mg%ie,mg%idx,mg%idy,mg%idz,igs,ige &
               ,tpsi%zwf(:,:,:,ispin,io,ik,im),htpsi%zwf(:,:,:,ispin,io,ik,im) &
               ,V_local(ispin)%f,system%Ac_micro%v,system%div_Ac%f,stencil%coef_lap0 &
               ,stencil%coef_lap,stencil%coef_nab,system%vec_k(1:3,ik))
      end do
      end do
      end do
      end do
!$omp end do nowait
      end do
      end do
    end do
!$omp end parallel
    call timer_end  (LOG_UHPSI_OVL_PHASE4)
  end subroutine zstencil_microac_overlapped
end subroutine hpsi

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

subroutine update_vlocal(mg,nspin,Vh,Vpsl,Vxc,Vlocal)
  use structures
  use timer
  implicit none
  type(s_rgrid), intent(in) :: mg
  integer       ,intent(in) :: nspin
  type(s_scalar),intent(in) :: Vh,Vpsl,Vxc(nspin)
  type(s_scalar)            :: Vlocal(nspin)
  !
  integer :: is,ix,iy,iz

  do is=1,nspin
!$omp parallel do collapse(2) private(ix,iy,iz)
    do iz=mg%is(3),mg%ie(3)
    do iy=mg%is(2),mg%ie(2)
    do ix=mg%is(1),mg%ie(1)
      Vlocal(is)%f(ix,iy,iz) = Vpsl%f(ix,iy,iz) + Vh%f(ix,iy,iz) + Vxc(is)%f(ix,iy,iz)
    end do
    end do
    end do
  end do

  return
end subroutine update_vlocal

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

subroutine update_kvector_nonlocalpt(ik_s,ik_e,system,ppg)
  use math_constants,only : zi
  use structures
  use update_kvector_so_sub, only: update_kvector_so, SPIN_ORBIT_ON
  use update_kvector_plusU_sub, only: update_kvector_plusU, PLUS_U_ON
  implicit none
  integer           ,intent(in) :: ik_s,ik_e !,n_max
  type(s_dft_system),intent(in) :: system
  type(s_pp_grid)               :: ppg
  !
  integer :: ilma,iatom,j,ik
  real(8) :: x,y,z
  complex(8) :: ekr
  real(8),allocatable :: kAc(:,:)

  allocate(kAc(3,ik_s:ik_e))
  do ik=ik_s,ik_e
    kAc(1:3,ik) = system%vec_k(1:3,ik) + system%vec_Ac(1:3)
  end do

  if ( SPIN_ORBIT_ON ) then
    call update_kvector_so( ppg, kAc, ik_s, ik_e )
  end if
  if ( PLUS_U_ON ) then
    call update_kvector_plusU( ppg, kAc, ik_s, ik_e )
  end if

  if(.not.allocated(ppg%zekr_uV)) allocate(ppg%zekr_uV(ppg%nps,ppg%nlma,ik_s:ik_e))

!$omp parallel do collapse(2) private(ik,ilma,iatom,j,x,y,z,ekr)
  do ik=ik_s,ik_e
    do ilma=1,ppg%nlma
      iatom = ppg%ia_tbl(ilma)
      do j=1,ppg%mps(iatom)
        x = ppg%rxyz(1,j,iatom)
        y = ppg%rxyz(2,j,iatom)
        z = ppg%rxyz(3,j,iatom)
        ekr = exp(zi*(kAc(1,ik)*x+kAc(2,ik)*y+kAc(3,ik)*z))
        ppg%zekr_uV(j,ilma,ik) = conjg(ekr) * ppg%uv(j,ilma)
      end do
    end do
  end do
!$omp end parallel do

  deallocate(kAc)
  return
end subroutine update_kvector_nonlocalpt

subroutine update_kvector_nonlocalpt_microAc(ik_s,ik_e,system,ppg)
  use math_constants,only : zi
  use structures
  use timer
!  use fdtd_coulomb_gauge, only: line_integral
  implicit none
  integer           ,intent(in) :: ik_s,ik_e !,n_max
  type(s_dft_system),intent(in) :: system
!  type(s_rgrid)     ,intent(in) :: lg
!  real(8)           ,intent(in) :: vec_Ac(3,lg%is(1):lg%ie(1),lg%is(2):lg%ie(2),lg%is(3):lg%ie(3))
  type(s_pp_grid)               :: ppg
  !
  integer :: ilma,iatom,j,ik,ix,iy,iz
  real(8) :: kAc(3),x,y,z ! ,Hgs(3),r1_r0(3),r1(3),r0(3),integral
  complex(8) :: ekr
!  integer,allocatable :: index(:)
!  real(8),allocatable :: A_lerp(:,:),line(:,:),wrk(:)
!  Hgs = system%Hgs
!  allocate(A_lerp(3,n_max),line(3,n_max),wrk(n_max),index(n_max))

  call timer_begin(LOG_SS_UPDATE_NONLOCALPT_MICROAC)

  if(.not.allocated(ppg%zekr_uV)) allocate(ppg%zekr_uV(ppg%nps,ppg%nlma,ik_s:ik_e))
  do ilma=1,ppg%nlma
    iatom = ppg%ia_tbl(ilma)
    do j=1,ppg%mps(iatom)

!    ! C. Pickard & F. Mauri, PRL 91, 196401 (2003).
!      ix = ppg%jxyz(1,j,iatom) - lg%is(1) ! lg%is:lg%ie --> 0:ng%num-1
!      iy = ppg%jxyz(2,j,iatom) - lg%is(2)
!      iz = ppg%jxyz(3,j,iatom) - lg%is(3)
!      r1(1) = dble(ix)*hgs(1)
!      r1(2) = dble(iy)*hgs(2)
!      r1(3) = dble(iz)*hgs(3)
!      r1_r0 = ppg%rxyz(1:3,j,iatom) - system%Rion(1:3,iatom)
!      r0 = r1 - r1_r0
!      ! path: r0 --> r1 = (ix*hgs(1),iy*hgs(2),iz*hgs(3))
!      call line_integral(integral,r0,vec_Ac,lg%num(1),lg%num(2),lg%num(3),ix,iy,iz,Hgs(1),Hgs(2),Hgs(3) &
!            ,A_lerp,line,wrk,n_max,index)

      ix = ppg%jxyz(1,j,iatom)
      iy = ppg%jxyz(2,j,iatom)
      iz = ppg%jxyz(3,j,iatom)
      x = ppg%rxyz(1,j,iatom)
      y = ppg%rxyz(2,j,iatom)
      z = ppg%rxyz(3,j,iatom)
      do ik=ik_s,ik_e

!        k = system%vec_k(:,ik)
!        ekr = exp(zi*( k(1)*x+k(2)*y+k(3)*z + integral ))

      ! approximation: vector potential is almost constant in typical cutoff radius of pseudopotentials
        kAc = system%vec_k(:,ik) + system%Ac_micro%v(:,ix,iy,iz)
        ekr = exp(zi*( kAc(1)*x+kAc(2)*y+kAc(3)*z ))
        ppg%zekr_uV(j,ilma,ik) = conjg(ekr) * ppg%uv(j,ilma)
      end do
    end do
  end do
!  deallocate(A_lerp,line,wrk,index)

  call timer_end(LOG_SS_UPDATE_NONLOCALPT_MICROAC)

  return
end subroutine update_kvector_nonlocalpt_microAc

end module hamiltonian


text:
! =============================================================================
!
!> \file
!> \brief  Contains the default values for the objects in the Loggers module
!> \author CMJ (XCP-3; LANL)
!
! =============================================================================

    !> \brief Defines the number of message types contained
    integer(int32), private, parameter:: numMessageLevels = 7_int32

    !> \brief Defines numerics for the logging filters/types
    enum, bind(c)
        enumerator:: QUIETLog     =  0_c_int
        enumerator:: CRITICALLog  = 10_c_int
        enumerator:: ERRORLog     = 20_c_int
        enumerator:: WARNINGLog   = 30_c_int
        enumerator:: INFOLog      = 40_c_int
        enumerator:: DEBUGLog     = 50_c_int
        enumerator:: UNSETLog     = 60_c_int
    end enum
    public:: QUIETLog, CRITICALLog, ERRORLog, WARNINGLog, INFOLog, &
        & DEBUGLog, UNSETLog

    !> @{
    !> Labels prepended to a message for different message types
    character(*), private, parameter:: preCritical = "(F)   Fatal: "
    character(*), private, parameter:: preError    = "(E)   Error: "
    character(*), private, parameter:: preWarning  = "(W) Warning: "
    character(*), private, parameter:: preInfo     = "       Info: "
    character(*), private, parameter:: preDebug    = "      Debug: "
    character(*), private, parameter:: preUnset    = ""
    !> }

    !> @{
    !> Formats used to print messages
    character(*), private, parameter:: stringFormat = "(A)"
    !> @}

    !> \brief Defines the default message level printed
    integer(c_int), private, parameter:: defaultLevel = UNSETLog

    !> \brief Defines the default stream messages are sent to
    integer(int32), private, parameter:: defaultStream = output_unit

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense