text:
    !-------------------------------------------------------------------------------------------------------------
    !
    !> \file    RemPureConPhase.f90
    !> \brief   Remove a pure condensed phase from the system.
    !> \author  M.H.A. Piro
    !> \date    Apr. 26, 2012
    !> \sa      CheckPureConPhaseRem.f90
    !> \sa      CheckPhaseChange.f90
    !
    !
    ! Revisions:
    ! ==========
    !
    !   Date            Programmer          Description of change
    !   ----            ----------          ---------------------
    !   03/31/2011      M.H.A. Piro         Original code
    !   10/25/2011      M.H.A. Piro         Clean up code: modules, fragment code into multiple subroutines.
    !   03/25/2012      M.H.A. Piro         Added the capability of checking the iteration history when there
    !                                       is more than one pure condensed phase that needs to be removed.
    !   04/26/2012      M.H.A. Piro         Implementing Gibbs energy Minimization algorithm and dOxygen.
    !   05/30/2012      M.H.A. Piro         Store the previous phase assemblage info in case if a pure con phase
    !                                        is initially removed and one (or more) additional phases are removed,
    !                                        but only the first phase is returned to the system.
    !
    !
    ! Purpose:
    ! ========
    !
    !> \details The purpose of this subroutine is to remove a pure condensed phase from the estimated phase
    !! assemblage.
    !
    !
    ! Pertinent variables:
    ! ====================
    !
    !> \param[in]   iPhaseChange    An integer scalar representing the index of the pure condensed phase to be
    !!                               removed from the system.
    !> \param[out]  lSwapLater      A logical scalar indicating whether this phase should be swapped later.
    !> \param[out]  lPhasePass      A logical scalar indicating whether the new phase assemblage has passed.
    !
    !-------------------------------------------------------------------------------------------------------------

subroutine RemPureConPhase(iPhaseChange,lSwapLater,lPhasePass)

    USE ModuleThermo
    USE ModuleGEMSolver

    implicit none

    integer::                        i, j, k, iPhaseChange, INFO, nConPhasesLast
    integer,dimension(nConPhases)::  iTempVec
    real(8)::                        dTemp
    real(8),dimension(nConPhases)::  dTempVec
    logical::                        lSwapLater, lPhasePass

    ! Initialize variables:
    lSwapLater = .FALSE.
    lPhasePass = .FALSE.

    ! Store the phase assemblage and number of moles in a temporary vector:
    iTempVec(1:nConPhases) = iAssemblage(1:nConPhases)
    dTempVec(1:nConPhases) = dMolesPhase(1:nConPhases)
    nConPhasesLast         = nConPhases

    ! If iPhaseChange was the last one added, give the system a chance to converge:
    if ((iConPhaseLast == iAssemblage(iPhaseChange)).AND.(iterGlobal - iterLastCon < iterStep)) return

    ! Remove the pure condensed phase corresponding to iPhaseChange:
    iConPhaseLast             = iAssemblage(iPhaseChange)
    iAssemblage(iPhaseChange) = iAssemblage(nConPhases)
    iAssemblage(nConPhases)   = 0
    dTemp                     = dMolesPhase(iPhaseChange)
    dMolesPhase(iPhaseChange) = dMolesPhase(nConPhases)
    dMolesPhase(nConPhases)   = 0D0
    nConPhases                = nConPhases - 1
    dMolesPhase               = dMolesPhase * 0.95D0

    ! Check that this phase change is acceptable:
    k = MAX(1, nConPhases)
    do i = 1, k

        call CheckPhaseChange(lPhasePass,INFO)

        if (INFO > nElements + nSolnPhases) then
            ! A pure condensed phase should be removed.
            j = INFO - nElements - nSolnPhases
            iAssemblage(j)          = iAssemblage(nConPhases)
            dMolesPhase(j)          = dMolesPhase(nConPhases)
            iAssemblage(nConPhases) = 0
            dMolesPhase(nConPhases) = 0D0
            nConPhases              = nConPhases - 1
            iterLast                = iterGlobal
        else
            ! This phase assemblage is appropriate.
            exit
        end if
    end do

    if (lPhasePass .EQV. .FALSE.) then
        ! The phase in question cannot be removed.  Add it back to the assemblage:
        nConPhases                = nConPhasesLast
        iAssemblage(1:nConPhases) = iTempVec(1:nConPhases)
        dMolesPhase(1:nConPhases) = dTempVec(1:nConPhases)
        lSwapLater                = .TRUE.
    else
        ! The new phase assemblage is acceptable.
        iterLastCon               = iterGlobal
        iterLast                  = iterGlobal
        dMolesPhase               = DABS(dMolesPhase)

        ! Compute the number of moles of all phases:
        call CompMolAllSolnPhases

    end if

    return

end subroutine RemPureConPhase

max_stars_repo_name: nrc-fuels/thermochimica

max_issues_repo_name: nrc-fuels/thermochimica

max_forks_repo_name: nrc-fuels/thermochimica


text:
C Copyright(C) 2009 Sandia Corporation. Under the terms of Contract
C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
C certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of Sandia Corporation nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE SCAPRT (NAMECO, NAME, IVAR, TIMES,
     &   VALMIN, NUMMIN, XYZMIN, ISTMIN, VALMAX, NUMMAX, XYZMAX, ISTMAX)
C=======================================================================

C   --*** SCAPRT *** (BLOT) Print variable min/max
C   --   Written by Amy Gilkey - revised 04/01/88
C   --
C   --SCAPRT prints the minimum and maximum values for the desired variable.
C   --
C   --Parameters:
C   --   NAMECO - IN - the coordinate names
C   --   NAME - IN - the variable name
C   --   IVAR - IN - the variable index
C   --   TIMES - IN - the database times
C   --   VALMIN, VALMAX - IN - the minimum and maximum value
C   --      (with selected element block and birth/death)
C   --   NUMMIN, NUMMAX - IN - the node or element number of the minimum and
C   --      maximum value (nodal and element variable only)
C   --   XYZMIN, XYZMAX - IN - the coordinates of NUMMIN, NUMMAX
C   --   ISTMIN, ISTMAX - IN - the step number of the minimum and maximum
C   --      value
C   --
C   --Common Variables:
C   --   Uses NDIM of /DBNUMS/

      include 'dbnums.blk'

      CHARACTER*(*) NAMECO(*)
      CHARACTER*(*) NAME
      REAL TIMES(*)
      REAL XYZMIN(3), XYZMAX(3)

      REAL RVAL(2), RXYZ(2,3), RTIM(2)
      CHARACTER*20 SVAL(0:2), SXYZ0(3), SXYZ(2,3), STIM(0:2)
      CHARACTER*4 STR4
      CHARACTER TYP

      IF ((ISTMIN .LE. 0) .OR. (VALMIN .EQ. VALMAX)) THEN
         CALL NUMSTR (1, 4, VALMIN, SVAL(1), LVAL)
         WRITE (*, 10060) NAME(:LENSTR(NAME)),
     &      ' does not vary - all values = ', SVAL(1)(:LVAL)

      ELSE
         CALL DBVTYP_BL (IVAR, TYP, IDUM)

         RVAL(1) = VALMIN
         RVAL(2) = VALMAX
         CALL NUMSTR (2, 4, RVAL, SVAL(1), LVAL)
         CALL CENLAB ('Value', 2, SVAL(1), LVAL, SVAL(0))

         IF ((TYP .EQ. 'N') .OR. (TYP .EQ. 'E')) THEN
            NNDIM = MAX (MIN (NDIM, 3), 2)
            DO 100 I = 1, MIN (NDIM, 3)
               RXYZ(1,I) = XYZMIN(I)
               RXYZ(2,I) = XYZMAX(I)
  100       CONTINUE
            CALL NUMSTR (NNDIM * 2, 3, RXYZ, SXYZ(1,1), LX)
            LY = LX
            LZ = LX
            CALL CENLAB (NAMECO(1), 2, SXYZ(1,1), LX, SXYZ0(1))
            CALL CENLAB (NAMECO(2), 2, SXYZ(1,2), LY, SXYZ0(2))
            IF (NDIM .GE. 3) THEN
               CALL CENLAB (NAMECO(3), 2, SXYZ(1,3), LZ, SXYZ0(3))
            END IF
         END IF

         RTIM(1) = TIMES(ISTMIN)
         RTIM(2) = TIMES(ISTMAX)
         CALL NUMSTR (2, 4, RTIM, STIM(1), LTIM)
         CALL CENLAB ('Time', 2, STIM(1), LTIM, STIM(0))

         IF ((TYP .EQ. 'N') .OR. (TYP .EQ. 'E')) THEN
            IF (TYP .EQ. 'N') THEN
               STR4 = 'Node'
            ELSE
               STR4 = 'Elem'
            END IF
            IF (NDIM .LE. 2) THEN
               WRITE (*, 10000) 'Range: ', SVAL(0)(:LVAL),
     &            SXYZ0(1)(:LX), SXYZ0(2)(:LY),
     &            STIM(0)(:LTIM), STR4, 'Step'
               WRITE (*, 10010, IOSTAT=IDUM) 'Minimum', SVAL(1)(:LVAL),
     &            SXYZ(1,1)(:LX), SXYZ(1,2)(:LY),
     &            STIM(1)(:LTIM), NUMMIN, ISTMIN
               WRITE (*, 10010, IOSTAT=IDUM) 'Maximum', SVAL(2)(:LVAL),
     &            SXYZ(2,1)(:LX), SXYZ(2,2)(:LY),
     &            STIM(2)(:LTIM), NUMMAX, ISTMAX
10000           FORMAT (4X, A7, 1X, A, 3X, 2 (1X, A),
     &            4X, A, 3X, A4, 3X, A4)
10010           FORMAT (4X, A7, 1X, A, 3X, 2 (1X, A), 4X, A, I7, I7)
            ELSE
               WRITE (*, 10020) 'Range: ', SVAL(0)(:LVAL),
     &            SXYZ0(1)(:LX), SXYZ0(2)(:LY), SXYZ0(3)(:LZ),
     &            STIM(0)(:LTIM), STR4, 'Step'
               WRITE (*, 10030, IOSTAT=IDUM) 'Minimum', SVAL(1)(:LVAL),
     &            SXYZ(1,1)(:LX), SXYZ(1,2)(:LY), SXYZ(1,3)(:LZ),
     &            STIM(1)(:LTIM), NUMMIN, ISTMIN
               WRITE (*, 10030, IOSTAT=IDUM) 'Maximum', SVAL(2)(:LVAL),
     &            SXYZ(2,1)(:LX), SXYZ(2,2)(:LY), SXYZ(2,3)(:LZ),
     &            STIM(2)(:LTIM), NUMMAX, ISTMAX
10020           FORMAT (4X, A7, 1X, A, 3X, 3 (1X, A),
     &            4X, A, 3X, A4, 3X, A4)
10030           FORMAT (4X, A7, 1X, A, 3X, 3 (1X, A), 4X, A, I7, I7)
            END IF
         ELSE
            WRITE (*, 10040) 'Range: ', SVAL(0)(:LVAL),
     &         STIM(0)(:LTIM), 'Step'
            WRITE (*, 10050, IOSTAT=IDUM) 'Minimum', SVAL(1)(:LVAL),
     &         STIM(1)(:LTIM), ISTMIN
            WRITE (*, 10050, IOSTAT=IDUM) 'Maximum', SVAL(2)(:LVAL),
     &         STIM(2)(:LTIM), ISTMAX
10040        FORMAT (4X, A7, 1X, A, 4X, A, 3X, A4)
10050        FORMAT (4X, A7, 1X, A, 4X, A, I7)
         END IF
      END IF

      RETURN
10060  FORMAT (' *** WARNING - ', 5A)
      END

ext: f

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
!     Post-proc of vgrid.in after gen_vgs.f90, and do vertical interpolation if # of levels is changed
!     Inputs: nlev.gr3 (edit from output of gen_vqs.f90), vgrid.in
!     Output: vgrid.new
!     pgf90 -O2 -mcmodel=medium  -Mbounds -Bstatic -o edit_vgrid edit_vgrid.f90
      allocatable :: sig0(:),sig(:),sigout(:)

      open(14,file='nlev.gr3',status='old')
      open(19,file='vgrid.in',status='old')
      open(20,file='vgrid.new',status='replace')
      read(14,*)
      read(14,*)ne,np
      read(19,*)
      read(19,*)nvrt
      write(20,*)1; write(20,*)nvrt
      allocate(sig0(nvrt),sig(nvrt),sigout(nvrt))
      do i=1,np
        read(14,*)j,x,y,tmp
        nlev=nint(tmp)
        read(19,*)j,kbp,sig(kbp:nvrt)

        nlev0=nvrt-kbp+1
        if(nlev0==nlev) then
          kbp1=kbp
          sigout(kbp:nvrt)=sig(kbp:nvrt)
        else !do interpolation
          print*, 'Altering node ',i
          sig0(kbp:nvrt)=(/(-1+real(k-kbp)/(nvrt-kbp),k=kbp,nvrt)/) !even spaced base

          kbp1=nvrt-nlev+1
          sigout(kbp1)=-1; sigout(nvrt)=0
          do k=kbp1+1,nvrt-1
            sig1=-1+real(k-kbp1)/(nvrt-kbp1)
            k0=0 !flag
            do kk=kbp+1,nvrt
              if(sig1>=sig0(kk-1).and.sig1<=sig0(kk)) then
                k0=kk; exit
              endif
            enddo !kk
            if(k0==0) stop 'cannot find a level'

            rat=(sig1-sig0(k0-1))/(sig0(k0)-sig0(k0-1))
            sigout(k)=sig(k0)*rat+sig(k0-1)*(1-rat)
          enddo !k

          !Check
          do k=kbp1+1,nvrt
            if(sigout(k)<=sigout(k-1)) then
              print*, 'Inverted levels:',i,k
              stop
            endif
          enddo !k
        endif !nlev0

        !Output
        !write(20,'(2(1x,i10),100000(1x,f14.7))')i,kbp1,(sigout(max(k,kbp1)),k=1,nvrt)
        write(20,'(2(1x,i10),100000(1x,f14.7))')i,kbp1,sigout(kbp1:nvrt)
      enddo !i=1,np

      stop
      end

max_stars_repo_name: rustychris/schism

max_issues_repo_name: rustychris/schism

max_forks_repo_name: rustychris/schism


text:
program TargetReduction_3D

    implicit none

    integer, dimension ( 3,3,3 ) ::  iV
    integer i,j,k

!   !$OMP  target enter data map ( from : iV)

     print*,'Starting'
    !$OMP  target MAP(from:iV)
    !$OMP  teams
    !$OMP  distribute private (i)
    do i = 1, 3
      !$omp parallel  private ( j, k )
      do j = 1, 3
        !$omp do
        do k = 1, 3
            iV(i,j,k) = i+j+k
        end do
        !$omp end do
      end do
      !$omp end parallel
    end do
    !$OMP  end distribute
    !$OMP END TEAMS
    !$OMP END TARGET
    print*, 'MaxSpeed'
    print*, iV
end program TargetReduction_3D

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


text:
      SUBROUTINE PARTN1
C
C     THIS IS THE DMAP MODULE PARTN WHICH PARTITIONS A MATRIX -A- INTO
C     FOUR PARTITIONS, SOME OR ALL OF WHICH MAY BE PURGED.
C
C
C                             **                  **
C                             *       I            *
C                             *  A11  I    A12     *
C          **   **            *       I            *
C          *     *            * ------+----------- *
C          *  A  *  BECOMES   *       I            *
C          *     *            *       I            *
C          **   **            *  A21  I    A22     *
C                             *       I            *
C                             **                  **
C
C
C     BASED ON ROW PARTITION MATRIX -RP- AND COLUMN PARTITION MATRIX
C     -CP-
C
C     DMAP SEQUENCE.
C
C     PARTN A,CP,RP/A11,A21,A12,A22/V,Y,SYM/V,Y,TYPE/V,Y,F11/V,Y,F21/
C                                   V,Y,F12/V,Y,F22 $
C
      IMPLICIT INTEGER (A-Z)
      EXTERNAL        RSHIFT,ANDF
      LOGICAL         CPHERE,RPHERE,CPNULL,RPNULL
      DIMENSION       MCB(7,4),BLOCK(80),SUBR(2),AIJ(4),BUFFS(5),
     1                HEAD(2),MCBA(7),REFUS(3)
      CHARACTER       UFM*23,UWM*25,UIM*29,SFM*25,SWM*27
      COMMON /XMSSG / UFM,UWM,UIM,SFM,SWM
      COMMON /MAHCIN/ MACHX
      COMMON /ZNTPKX/ ELEM(4),ROW,EOL
      COMMON /NAMES / RD,RDREW,WRT,WRTREW,CLSREW,CLS
      COMMON /SYSTEM/ SYSBUF,OUTPT,XXX(37),NBPW
      COMMON /PRTMRG/ CPSIZE,RPSIZE,CPONES,RPONES,CPNULL,RPNULL,CPHERE,
     1                RPHERE,ICP,NCP,IRP,NRP
      COMMON /ZZZZZZ/ Z(1)
      COMMON /BLANK / SYM,TYPE,FORM(4),CPCOL,RPCOL,IREQCL
      DATA    SUBR  / 4HPART,4HN1    /, A,CP,RP /101,102,103/,
     1        AIJ   / 201,202,203,204/
      DATA    NAFORM/ 4HFORM /, NATYPE/ 4HTYPE /, REFUS/2*3H   ,3HREF /
      DATA    EOR   / 1      /
C
      CORE = KORSZ(Z)
      BUFFS(1) = CORE - SYSBUF - 2
      DO 10 I = 2,5
      BUFFS(I) = BUFFS(I-1) - SYSBUF - 2
   10 CONTINUE
      CORE = BUFFS(5) - 1
      IF (CORE .LT. 10) CALL MESAGE (-8,0,SUBR)
C
C     OPEN MATRIX TO BE PARTITIONED.  IF PURGED RETURN IS MADE
C
      BUFF = BUFFS(5)
      CALL OPEN (*425,A,Z(BUFF),RDREW)
      CALL SKPREC (A,1)
      MCBA(1) = A
      CALL RDTRL (MCBA)
      INFORM = MCBA(4)
C
C     CALL TO PARTN2 WILL PROCESS -CP- AND -RP- INTO BIT STRINGS AND
C     DETERMINE SIZES OF THE PARTITIONS.
C
      BUFF = BUFFS(4)
      CALL PARTN2 (CP,RP,CORE,Z(BUFF))
C
C     IF RPSIZE OR CPSIZE ARE 0 THEY ARE SET EQUAL TO THE RESPECTIVE
C     SIZE OF A
C
      IF (CPSIZE .EQ. 0) CPSIZE = MCBA(2)
      IF (RPSIZE .EQ. 0) RPSIZE = MCBA(3)
C
C     MATRIX COMPATIBILITY CHECKS
C
      IF (RPSIZE.EQ.MCBA(3) .AND. CPSIZE.EQ.MCBA(2)) GO TO 40
      WRITE  (OUTPT,30) SWM,MCBA(3),MCBA(2),RPSIZE,CPSIZE
   30 FORMAT (A27,' 2166, MATRIX TO BE PARTITIONED IS OF SIZE',I10,
     1       ' ROWS BY',I10,' COLUMNS.', /5X,'ROW PARTITION SIZE IS',
     2       I10,' COLUMN PARTITION SIZE IS',I10,' (INCOMPATIBLE).')
C
C     PREPARE OUTPUT DATA BLOCKS AS REQUIRED.
C
   40 CPZERO = MCBA(2) - CPONES
      RPZERO = MCBA(3) - RPONES
C
C     CHECK OF TYPE PARAMETER
C
      NTYPE = MCBA(5)
      IF (NTYPE .EQ. TYPE) GO TO 60
      IF (TYPE  .EQ.    0) GO TO 54
      IF (TYPE.LT.0 .OR. TYPE.GT.4) GO TO 52
      WRITE  (OUTPT,50) SWM,NATYPE,TYPE,REFUS(1),SUBR,NTYPE
   50 FORMAT (A27,' 2163, REQUESTED VALUE OF ',A4,I10,2X,A3,
     1        'USED BY ',2A4,'. LOGICAL CHOICE IS',I10)
      NTYPE = TYPE
      GO TO 60
   52 WRITE (OUTPT,50) SWM,NATYPE,TYPE,REFUS(3),SUBR,NTYPE
   54 TYPE = NTYPE
C
   60 DO 140 I = 1,4
      FILE = AIJ(I)
      MCB(1,I) = 0
      COLS = CPZERO
      ROWS = RPZERO
      IF (I.EQ.3 .OR. I.EQ.4) COLS = CPONES
      IF (I.EQ.2 .OR. I.EQ.4) ROWS = RPONES
C
C     IF ROWS OR COLS EQUAL ZERO NOTHING IS WRITTEN ON THIS PARTITION
C
      IF (ROWS.EQ.0 .OR. COLS.EQ.0) GO TO 140
      BUFF = BUFFS(I)
      CALL OPEN (*140,FILE,Z(BUFF),WRTREW)
      CALL FNAME (FILE,HEAD)
      CALL WRITE (FILE,HEAD,2,EOR)
C
C     CHECK OF THE FORM PARAMETER
C
      NFORM = FORM(I)
      IF (NFORM.LT.1 .OR. NFORM.GT.8) GO TO 110
      GO TO (70,130,100,70,70,70,100,70), NFORM
C
C     FORM IMPLIES SQUARE
C
   70 IF (ROWS .EQ. COLS) GO TO 130
   80 WRITE  (OUTPT,90) SWM,HEAD,NFORM,ROWS,COLS
   90 FORMAT (A27,' 2168, THE FORM PARAMETER AS GIVEN TO THE PARTITION',
     1        'ING MODULE FOR SUB-PARTITION ',2A4, /5X,'IS INCONSISTANT'
     2,       ' WITH ITS SIZE.  FORM =',I9,' SIZE =',I9,' ROWS BY',I9,
     3        ' COLUMNS.')
      GO TO 130
C
C     DIAGONAL OR ROW MATRIX
C
  100 IF (COLS .EQ. 1) GO TO 130
      GO TO 80
C
C     NO FORM SPECIFIED THUS IT IS SQUARE IF ROWS = COLS OR RECTANGULAR
C     OTHERWISE.
C
  110 NFORM = 2
      IF (ROWS .EQ. COLS) NFORM = 1
      IF (SYM.LT.0 .AND. INFORM.EQ.6 .AND. NFORM.EQ.1 .AND.
     1   (I.EQ.1 .OR. I.EQ.4)) NFORM = 6
      IF (FORM(I) .EQ. 0) GO TO 128
      JJ = 1
      IF (FORM(I).LT.1 .OR. FORM(I).GT.8) JJ = 3
      WRITE (OUTPT,50) SWM,NAFORM,FORM(I),REFUS(JJ),SUBR,NFORM
      IF (JJ .NE. 3) NFORM = FORM(I)
  128 FORM(I) = NFORM
C
C     TRAILER INITIALIZATION.  BLDPKN WILL SET MCB(2) AND MCB(6) LATER.
C
  130 CALL MAKMCB (MCB(1,I),FILE,ROWS,NFORM,NTYPE)
  140 CONTINUE
C
C     ROW PARTITIONING BIT STRING IS AT THIS POINT CONVERTED TO A CORE
C     VECTOR ONE WORD PER BIT.  EACH WORD CONTAINS THE ROW NUMBER OF THE
C     PARTITION TO WHICH THE ELEMENT OF -A- IS TO BE MOVED TO.  IF THE
C     NUMBER IS NEGATIVE THE ELEMENT IS MOVED TO THE LOWER PARTITIONS
C     AND IF THE NUMBER IS POSITIVE THE ELEMENT IS MOVED TO THE UPPER
C     PARTITION
C
      IZ = NRP + 1
      NZ = IZ  + RPSIZE - 1
      IF (NZ+NBPW .GT. CORE) CALL MESAGE (-8,0,SUBR)
      IF (.NOT.RPNULL .AND. RPONES.NE.0) GO TO 160
      K  = 0
      DO 150 I = IZ,NZ
      K  = K + 1
      Z(I) = K
  150 CONTINUE
      GO TO 210
  160 JZ   = IZ - 1
      ZERO = 0
      ONES = 0
C
C     NOTE THIS LOGIC WORKS ON CRAY WITH 48 OF 64 BIT INTEGER WORD
C
      DO 200 I = IRP,NRP
      DO 190 J = 1,NBPW
      SHIFT = NBPW - J
      BIT   = RSHIFT(Z(I),SHIFT)
      JZ    = JZ + 1
      IF (ANDF(BIT,1)) 180,170,180
  170 ZERO  = ZERO + 1
      Z(JZ) = ZERO
      GO TO 190
  180 ONES  = ONES - 1
      Z(JZ) = ONES
  190 CONTINUE
  200 CONTINUE
C
C     LOOP ON ALL THE COLUMNS OF -A-.
C
  210 IZM1  = IZ - 1
      DO 400 I = 1,CPSIZE
      IF (CPNULL) GO TO 220
      IL1   = I - 1
      BITWD = IL1/NBPW + ICP
      SHIFT = NBPW - MOD(IL1,NBPW) - 1
      BIT   = RSHIFT(Z(BITWD),SHIFT)
      IF (ANDF(BIT,1)) 230,220,230
C
C     ZERO-S COLUMN (LEFT PARTITIONS A11 AND A21)
C
  220 IFILE  = 1
      IBLOCK = 1
      GO TO 240
C
C     ONE-S COLUMN (RIGHT PARTITIONS A12 AND A22)
C
  230 IFILE  = 3
      IBLOCK = 41
      GO TO 240
C
C     START COLUMNS OF THE 2 AIJ PARTITIONS.
C
  240 KFILE  = IFILE
      KBLOCK = IBLOCK
      M = 0
      DO 270 J = 1,2
      IF (MCB(1,KFILE)) 260,260,250
  250 CALL BLDPK (NTYPE,MCB(5,KFILE),MCB(1,KFILE),BLOCK(KBLOCK),1)
      M = 1
  260 KFILE  = KFILE + 1
      KBLOCK = KBLOCK + 20
  270 CONTINUE
      IF (M) 280,390,280
C
C     START THE I-TH COLUMN OF THE MATRIX BEING PARTITIONED -A-.
C
  280 CALL INTPK (*350,A,0,NTYPE,0)
C
C     LOOP ON NON-ZEROS OF THE COLUMN
C
  290 IF (EOL) 300,300,350
C
C     PICK UP A NON-ZERO ELEMENT
C
  300 CALL ZNTPKI
C
C     DETERMINE ROW POSITION AND FILE DESTINATION.
C
      L = IZM1 + ROW
      IF (Z(L)) 320,310,310
C
C     ZERO-S ROW PARTITION.
C
  310 JROW   = Z(L)
      KFILE  = IFILE
      KBLOCK = IBLOCK
      GO TO 330
C
C     ONE-S ROW PARTITION.
C
  320 JROW   = -Z(L)
      KFILE  = IFILE  + 1
      KBLOCK = IBLOCK + 20
C
C     OUTPUT THE ELEMENT.
C
  330 IF (MCB(1,KFILE)) 290,290,340
  340 CALL BLDPKI (ELEM,JROW,MCB(1,KFILE),BLOCK(KBLOCK))
      GO TO 290
C
C     COMPLETE COLUMNS OF THE 2 AIJ PARTITIONS BEING WORKED ON.
C
  350 KFILE  = IFILE
      KBLOCK = IBLOCK
      DO 380 J = 1,2
      IF (MCB(1,KFILE)) 370,370,360
  360 CALL BLDPKN (MCB(1,KFILE),BLOCK(KBLOCK),MCB(1,KFILE))
  370 KFILE  = KFILE  + 1
      KBLOCK = KBLOCK + 20
  380 CONTINUE
      GO TO 400
C
C     COLUMN NOT BEING OUTPUT TO ANY PARTITIONS AT ALL THUS SKIP IT.
C
  390 CALL SKPREC (A,1)
C
  400 CONTINUE
C
C     WRAP UP.
C
      CALL CLOSE (A,CLSREW)
      DO 420 I = 1,4
      IF (MCB(1,I)) 420,420,410
  410 CALL WRTTRL (MCB(1,I))
      CALL CLOSE (MCB(1,I),CLSREW)
  420 CONTINUE
  425 RETURN
      END


text:
	include 'VICMAIN_FOR'
	subroutine main44

c  NIMSR2IOF:   convert a NIMS radiance cube to IOF ( = pi * I / F(Solar) ),
c  using the solar flux in the label.

c   2mar98 -lwk- initial version
c  27mar98 -lwk- fix processing of "special values"
c  15nov00 -lwk- check NS against buffer size
c  14nov02 -lwk- added FORMAT to xlget call for unix port

	parameter (maxl = 10000)		! max. line length
	parameter (ntasks=3)		! # of valid TASKs for label
	character*3 org
	character*4 fmt
	character*132 pbuf
	character*8 tasks(ntasks)/ 'NIMSCMM2', 'VISIS2', 'NIMSCMM'/,
	1 task
	real*4 wav(408), sflux(408), obuf(maxl)
	real*8 xpi/3.1415926536/
	real*4 valid_min/-1.7014101e38/
	logical nosv, xvptst

	call xvmessage('NIMSR2IOF version 2019-11-05',' ')

c  open input cube and check format:
	call xvunit( inp, 'INP', 1, i,' ')
	call xvopen( inp, i, 'OPEN_ACT', 'SA', 'IO_ACT', 'SA',' ')
	call xvget( inp, i, 'NL', nl, 'NS', ns, 'NB', nb, 'ORG', org,
	1 'FORMAT', fmt,' ')
	if (fmt.ne.'REAL') call mabend(' ** input cube must be Real **')
	if (org.ne.'BSQ') call mabend(' ** only BSQ is supported **')
	if (ns.gt.maxl) call mabend(' ** NS exceeds buffer size **')

c  get the Solar Flux from the label:
	do i=1,ntasks
	  task = tasks(i)
	  call xlget( inp, 'HISTORY', 'SOLAR_F', sflux, ist, 'NELEMENT',
	1  -1, 'NRET', nb0, 'FORMAT', 'REAL', 'HIST', task, ' ')
	  if (ist.eq.1) go to 10
	enddo
	if (ist.ne.1) call mabend(' ** no Solar Flux in label **')
10	if (nb0.ne.nb) call mabend(' ** Solar Flux has wrong length **')

c  open the output file
	call xvunit( iout, 'OUT', 1, i,' ')
	call xvopen( iout, i, 'U_NL', nl, 'U_NS', ns, 'U_NB', nb, 'OP',
	1  'WRITE', 'O_FORMAT', 'REAL', 'U_FORMAT', 'REAL', 'OPEN_ACT', 'SA',
	2  'IO_ACT', 'SA',' ')

	nosv = xvptst('NOSPCVAL')

c  main processing:
	do ib = 1,nb
	  do il = 1,nl
	    call xvread( inp, obuf, i, ' ')
	    do is = 1,ns
	      if (obuf(is).lt.valid_min) then
		if (nosv) obuf(is) = 0.0
	      else
		obuf(is) = obuf(is)*xpi/sflux(ib)
	      endif
	    enddo
	    call xvwrit( iout, obuf, i, ' ')
	  enddo
	enddo

	call xvclose( inp, i,' ')
	call xvclose( iout, i,' ')

	return
	end

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


! Ensure argument -fbackslash works as expected.

!--------------------------
! FLANG DRIVER (flang-new)
!--------------------------
! RUN: %flang -E %s  2>&1 | FileCheck %s --check-prefix=ESCAPED
! RUN: %flang -E -fbackslash -fno-backslash %s  2>&1 | FileCheck %s --check-prefix=ESCAPED
! RUN: %flang -E -fbackslash %s  2>&1 | FileCheck %s --check-prefix=UNESCAPED

!-----------------------------------------
! FRONTEND FLANG DRIVER (flang-new -fc1)
!-----------------------------------------
! RUN: %flang_fc1 -E %s  2>&1 | FileCheck %s --check-prefix=ESCAPED
! RUN: %flang_fc1 -E -fbackslash -fno-backslash %s  2>&1 | FileCheck %s --check-prefix=ESCAPED
! RUN: %flang_fc1 -E -fbackslash %s  2>&1 | FileCheck %s --check-prefix=UNESCAPED

!-----------------------------------------
! EXPECTED OUTPUT FOR ESCAPED BACKSLASHES
!-----------------------------------------
! ESCAPED:program backslash
! ESCAPED-NEXT:New\\nline
! ESCAPED-NOT:New\nline

!-------------------------------------------
! EXPECTED OUTPUT FOR UNESCAPED BACKSLASHES
!-------------------------------------------
! UNESCAPED:program backslash
! UNESCAPED-NEXT:New\nline
! UNESCAPED-NOT:New\\nline

program Backslash
    print *, 'New\nline'
end


text:
subroutine makevoro2c(voro,npt,V0d,nptref,voro2c,V0d2c,npt2c)
implicit none
include 'params.h'

real, intent(in) :: voro(malay,4)
real, intent(in) :: V0d(maxlay,2)
integer, intent(in) :: nptref,npt

integer, intent(out) :: npt2c
real, intent(out) :: voro2c(maxlay,4), V0d2c(maxlay)

integer num,j,i,order(npt+nptref-2),ind
real pos(npt+nptref-2),maxx,minn,temp(4)

num = nptref+npt-2

pos=[voro(2:npt,1), V0d(2:nptref,1)]

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

ind=1
!!!dismin=99999999999999
!print *, npt
do i=1,num
	!write(*,*)'ordre',i
	maxx=99999999
	if (i==1) then
		minn=-9999999
	else
		minn=pos(order(i-1))
	endif

	!write(*,*)'min','max',minn,maxx
	do j=1,num
		if ((minn<pos(j)).and.(pos(j)<=maxx)) then
	!		write(*,*)j,minn,voro(j,1),maxx
			ind=j
			maxx=pos(j)
	!		write(*,*)ind,maxx
		endif

! ! 		if ((abs(voro(j,1)-voro(i,1))<dismin).and.(i.ne.j)) then
! ! 		dismin=abs(voro(j,1)-voro(i,1))
! ! 		endif

	enddo
	order(i)=ind
 	!write(*,*)ind
 	!write(*,*)
!   	write(*,*)

enddo

!5555555555555555555555555555555555555555555555555555555
voro2c(1,1)=0
voro2c(1,2)=V0d(1,2)*(1+voro(1,2))
voro2c(1,3)=voro(1,3)
voro2c(1,4)=voro(1,4) !!CHANGE!!
V0d2c(1)=V0d(1,2)

temp(1)=V0d(1,2)
temp(2)=voro(1,2)
temp(3)=voro(1,3)
temp(4)=voro(1,4)  !!CHANGE!!

do i=1,num
      !write(*,*)i,num
      if (order(i)<npt) then ! it is a dv/v

		ind=order(i)+1
		!write(*,*)'1',ind
		if ((ind<2).or.(ind>npt)) then
			stop "stop2"
		endif

		voro2c(i+1,1)=voro(ind,1)
		voro2c(i+1,2)=temp(1)*(1+voro(ind,2))
		voro2c(i+1,3)=voro(ind,3)
		voro2c(i+1,4)=voro(ind,4)  !!CHANGE!!
		V0d2c(i+1)=temp(1)

		temp(2)=voro(ind,2)
		temp(3)=voro(ind,3)
		temp(4)=voro(ind,4) !!CHANGE!!

      elseif ((order(i)>=npt).and.(order(i)<=num)) then ! it is a V0ref

		ind=order(i)-npt+2
		!write(*,*)ind,order(i),npt
		if ((ind<2).or.(ind>nptref)) then
			  stop "stop1"
		endif

		voro2c(i+1,1)=V0d(ind,1)
		voro2c(i+1,2)=V0d(ind,2)*(1+temp(2))
		voro2c(i+1,3)=temp(3)
		voro2c(i+1,4)=temp(4) !!CHANGE!!
		V0d2c(i+1)=V0d(ind,2)

		temp(1)=V0d(ind,2)
      else

      write(*,*) 'stop2'
      stop
      endif

enddo

npt2c=num+1
! ! ! ! ! do i=1,npt2c-1
! ! ! ! !       if (voro2c(i,1)==voro2c(i+1,1)) then
! ! ! ! !
! ! ! ! !
! ! ! ! ! 	    write(*,*)'----Y--Y--Y--Y--Y--Y-------------'
! ! ! ! ! 	    write(*,*)nptref
! ! ! ! ! 	    do j=1,nptref
! ! ! ! ! 	    write(*,*)V0d(j,1),V0d(j,2)
! ! ! ! ! 	    enddo
! ! ! ! ! 	    write(*,*)'---------------------------'
! ! ! ! ! 	    write(*,*)'npt',npt
! ! ! ! ! 	    do j=1,npt
! ! ! ! ! 	    write(*,*)voro(j,1),voro(j,2)
! ! ! ! ! 	    enddo
! ! ! ! !
! ! ! ! ! ! ! 	    write(*,*)
! ! ! ! ! ! ! 	    write(*,*)'---------------------------'
! ! ! ! ! ! ! 	    do j=1,npt2c
! ! ! ! ! ! !             write(*,*)voro2c(j,1),voro2c(j,2),V0d2c(j)
! ! ! ! ! ! !             enddo
! ! ! ! ! ! !             write(*,*)'---------------------------'
! ! ! ! ! ! ! 	    stop "same depth"
! ! ! ! !       endif
! ! ! ! ! enddo

! ! ! if (voro2c(12,2)>50) then
! ! ! write(*,*)'----Y--Y--Y--Y--Y--Y-------------'
! ! ! 	    write(*,*)nptref
! ! ! 	    do j=1,nptref
! ! ! 	    write(*,*)V0d(j,1),V0d(j,2)
! ! ! 	    enddo
! ! ! 	    write(*,*)'---------------------------'
! ! ! 	    write(*,*)'npt',npt
! ! ! 	    do j=1,npt
! ! ! 	    write(*,*)voro(j,1),voro(j,2)
! ! ! 	    enddo
! ! !
! ! ! 	    write(*,*)'---------------------------'
! ! ! 	    do j=1,npt2c
! ! !             write(*,*)voro2c(j,1),voro2c(j,2),V0d2c(j)
! ! !             enddo
! ! ! endif

return
end

max_stars_repo_name: AntonBiryukovUofC/RayTracerFortran

max_issues_repo_name: AntonBiryukovUofC/RayTracerFortran

max_forks_repo_name: AntonBiryukovUofC/RayTracerFortran


text:
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! HND X
! HND X   TurboGAP
! HND X
! HND X   TurboGAP is copyright (c) 2019-2021, Miguel A. Caro and others
! HND X
! HND X   TurboGAP is published and distributed under the
! HND X      Academic Software License v1.0 (ASL)
! HND X
! HND X   This file, neighbors.f90, is copyright (c) 2019-2021, Miguel A. Caro
! HND X
! HND X   TurboGAP is distributed in the hope that it will be useful for non-commercial
! HND X   academic research, but WITHOUT ANY WARRANTY; without even the implied
! HND X   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! HND X   ASL for more details.
! HND X
! HND X   You should have received a copy of the ASL along with this program
! HND X   (e.g. in a LICENSE.md file); if not, you can write to the original
! HND X   licensor, Miguel Caro (mcaroba@gmail.com). The ASL is also published at
! HND X   http://github.com/gabor1/ASL
! HND X
! HND X   When using this software, please cite the following reference:
! HND X
! HND X   Miguel A. Caro. Phys. Rev. B 100, 024112 (2019)
! HND X
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

module neighbors

  use soap_turbo_functions

  contains

!**************************************************************************
!
! This subroutine returns the distance between ri and rj under certain
! boundary conditions.
!
  subroutine get_distance(posi, posj, a, b, c, PBC, dist, d, i_shift)

    implicit none

    real*8, intent(in) :: posi(1:3), posj(1:3), a(1:3), b(1:3), c(1:3)
    integer, intent(out) :: i_shift(1:3)
    logical, intent(in) :: PBC(1:3)
    real*8, intent(out) :: d
    real*8, intent(out) :: dist(1:3)
    real*8 :: d2, L(1:3), d_tol = 1.d-6
    real*8 :: mat(1:3,1:3), md, indices_real(1:3), res, res_opt, dist_opt(1:3), dist_temp(1:3)
    real*8, save :: a0(1:3) = 0.d0, b0(1:3) = 0.d0, c0(1:3) = 0.d0, mat_inv(1:3,1:3) = 0.d0
    integer :: i, j, k, indices(1:3)
    logical :: lattice_check_a(1:3), lattice_check_b(1:3), lattice_check_c(1:3)

    if( dabs(a(2)) < d_tol .and. dabs(a(3)) < d_tol .and. &
        dabs(b(1)) < d_tol .and. dabs(b(3)) < d_tol .and. &
        dabs(c(1)) < d_tol .and. dabs(c(2)) < d_tol )then
      dist = posj - posi
      i_shift = floor([dist(1)/a(1), dist(2)/b(2), dist(3)/c(3)])
!     Fast solution for orthorhombic cells
      L = (/ a(1), b(2), c(3) /)
      d2 = 0.d0
      do i = 1, 3
        if( PBC(i) )then
          dist(i) = modulo(posj(i) - posi(i), L(i))
          if( dist(i) > L(i)/2.d0 )then
            dist(i) = dist(i) - L(i)
          end if
        else
          dist(i) = posj(i) - posi(i)
        end if
        d2 = d2 + dist(i)**2
      end do
      d = dsqrt(d2)
    else if( all( PBC ) )then
!     Slow solution for other unit cells
      lattice_check_a = ( a /= a0 )
      lattice_check_b = ( b /= b0 )
      lattice_check_c = ( c /= c0 )
      if( any(lattice_check_a) .or. any(lattice_check_b) .or. any(lattice_check_c) )then
        a0 = a
        b0 = b
        c0 = c
!       We construct our matrix to get the MIC only if the lattice vectors have changed
        mat(1,1) = dot_product(a, a)
        mat(1,2) = dot_product(a, b)
        mat(1,3) = dot_product(a, c)
        mat(2,1) = mat(1,2)
        mat(2,2) = dot_product(b, b)
        mat(2,3) = dot_product(b, c)
        mat(3,1) = mat(1,3)
        mat(3,2) = mat(2,3)
        mat(3,3) = dot_product(c, c)
!       We compute the inverse of this matrix analytically
        md = -mat(1,3)**2*mat(2,2) + 2.d0*mat(1,2)*mat(1,3)*mat(2,3) - mat(1,1)*mat(2,3)**2 &
             - mat(1,2)**2*mat(3,3) + mat(1,1)*mat(2,2)*mat(3,3)
        mat_inv(1,1) = mat(2,2)*mat(3,3) - mat(2,3)**2
        mat_inv(1,2) = mat(1,3)*mat(2,3) - mat(1,2)*mat(3,3)
        mat_inv(1,3) = mat(1,2)*mat(2,3) - mat(1,3)*mat(2,2)
        mat_inv(2,1) = mat_inv(1,2)
        mat_inv(2,2) = mat(1,1)*mat(3,3) - mat(1,3)**2
        mat_inv(2,3) = mat(1,2)*mat(1,3) - mat(1,1)*mat(2,3)
        mat_inv(3,1) = mat_inv(1,3)
        mat_inv(3,2) = mat_inv(2,3)
        mat_inv(3,3) = mat(1,1)*mat(2,2) - mat(1,2)**2
        mat_inv = mat_inv / md
      end if
      dist = posj - posi
      indices_real = -1.d0 * (/ dot_product(dist, a), dot_product(dist, b), dot_product(dist,c) /)
      indices_real = matmul(mat_inv, indices_real)
      i_shift = floor(-indices_real)
!     Closest integer solution
      indices(1:3) = nint(indices_real(1:3))
!     We bruteforce the integer solution among the 27 points surrounding the real solution
      res_opt = 1.d10
      do i = indices(1)-1, indices(1)+1
        do j = indices(2)-1, indices(2)+1
          do k = indices(3)-1, indices(3)+1
            dist_temp(1:3) = dist(1:3) + dfloat(i)*a(1:3) + dfloat(j)*b(1:3) + dfloat(k)*c(1:3)
            res = dot_product(dist_temp, dist_temp)
            if( res < res_opt )then
              res_opt = res
              dist_opt = dist_temp
            end if
          end do
        end do
      end do
      dist = dist_opt
      d = dsqrt( dot_product(dist,dist) )
    else
      write(*,*) "Sorry, non-orthorhombic unit cells only work in combination with full PBC"
    end if

  return
  end subroutine get_distance
!**************************************************************************

!**************************************************************************
!
! This subroutine returns the number of primitive unit cells required to
! construct a supercell whose unit cell places are at least 2*rcut apart.
! This construct allows us to search for all the neighbors within the given
! cutoff. indices(1:3) tell the user how many repetitions are required to
! construct a unit cell with the properties outlined above. (1,1,1) means
! that the primitive unit cell is already enough.
!
  subroutine number_of_unit_cells_for_given_cutoff(a, b, c, rcut, PBC, indices)

    implicit none

    real*8, intent(in) :: a(1:3), b(1:3), c(1:3), rcut
    logical, intent(in) :: PBC(1:3)
    integer, intent(out) :: indices(1:3)
    real*8 :: axb(1:3), bxc(1:3), cxa(1:3), indices_real(1:3)
    real*8 :: mat(1:3,1:3), mat_inv(1:3,1:3), md
    integer :: i

    axb = cross_product(a, b)
    axb = axb / dsqrt( dot_product(axb, axb) )
    cxa = cross_product(c, a)
    cxa = cxa / dsqrt( dot_product(cxa, cxa) )
    bxc = cross_product(b, c)
    bxc = bxc / dsqrt( dot_product(bxc, bxc) )

!   We construct our matrix to get the MIC only if the lattice vectors have changed
    mat(1,1) = dot_product(axb, a)
    mat(1,2) = dot_product(axb, b)
    mat(1,3) = dot_product(axb, c)
    mat(2,1) = dot_product(cxa, a)
    mat(2,2) = dot_product(cxa, b)
    mat(2,3) = dot_product(cxa, c)
    mat(3,1) = dot_product(bxc, a)
    mat(3,2) = dot_product(bxc, b)
    mat(3,3) = dot_product(bxc, c)
!   We compute the inverse of this matrix analytically
    md = -mat(1,3)*mat(2,2)*mat(3,1) + mat(1,2)*mat(2,3)*mat(3,1) + mat(1,3)*mat(2,1)*mat(3,2) &
         -mat(1,1)*mat(2,3)*mat(3,2) - mat(1,2)*mat(2,1)*mat(3,3) + mat(1,1)*mat(2,2)*mat(3,3)
    mat_inv(1,1) = mat(2,2)*mat(3,3) - mat(2,3)*mat(3,2)
    mat_inv(1,2) = mat(1,3)*mat(3,2) - mat(1,2)*mat(3,3)
    mat_inv(1,3) = mat(1,2)*mat(2,3) - mat(1,3)*mat(2,2)
    mat_inv(2,1) = mat(2,3)*mat(3,1) - mat(2,1)*mat(3,3)
    mat_inv(2,2) = mat(1,1)*mat(3,3) - mat(1,3)*mat(3,1)
    mat_inv(2,3) = mat(1,3)*mat(2,1) - mat(1,1)*mat(2,3)
    mat_inv(3,1) = mat(2,1)*mat(3,2) - mat(2,2)*mat(3,1)
    mat_inv(3,2) = mat(1,2)*mat(3,1) - mat(1,1)*mat(3,2)
    mat_inv(3,3) = mat(1,1)*mat(2,2) - mat(1,2)*mat(2,1)
    mat_inv = mat_inv / md

    indices_real = 2.d0 * (/ rcut, rcut, rcut /)
    indices_real = matmul(mat_inv, indices_real)

    indices_real = dabs(indices_real)
    do i = 1, 3
      indices(i) = int( indices_real(i) )
      if( indices_real(i) > dfloat(indices(i)) )then
        indices(i) = indices(i) + 1
      end if
    end do

  end subroutine
!**************************************************************************

!**************************************************************************
! This subroutine reads in the XYZ file and builds the lists of neighbors, the spherical
! coordinates, etc.
!
  subroutine build_neighbors_list(positions, a_box, b_box, c_box, do_timing, &
                                  species_supercell, rcut_max, n_atom_pairs, rjs, &
                                  thetas, phis, xyz, n_neigh, neighbors_list, neighbor_species, &
                                  n_sites, indices, rebuild_neighbors_list, do_list, rank )

    implicit none

!   Input variables
    real*8, intent(in) :: rcut_max, positions(:,:)
!    integer, intent(in) :: species_multiplicity(:), n_species
    integer, intent(in) :: species_supercell(:), indices(1:3), n_sites, rank
    logical, intent(in) :: do_timing, rebuild_neighbors_list, do_list(:)

!   Output variables
    integer, intent(out) :: n_atom_pairs
!    logical, allocatable, intent(out) :: mask_species(:,:)

!   In and out variables
    real*8, allocatable, intent(inout) :: rjs(:), thetas(:), phis(:), xyz(:,:)
    integer, allocatable, intent(inout) :: neighbor_species(:)
    real*8, intent(inout) :: a_box(1:3), b_box(1:3), c_box(1:3)
    integer, allocatable, intent(inout) :: neighbors_list(:), n_neigh(:)
!    integer, allocatable, intent(inout) :: species_supercell(:,:)

!   Internal variables
    real*8 :: time1, time2, dist(1:3), d, neigh_time, time3, tol, d_tol = 1.d-6
    integer, allocatable :: neighbors_list_temp(:), head(:), this_list(:)
    integer :: n_neigh_max, i, j, n_sites_supercell, &
               k, k2, i2, j2, i3, j3, k3, mx, my, mz, i_shift(1:3)
    logical :: is_box_square, is_box_small
    logical, save :: print_cutoff_warning = .true., print_shape_warning = .true.

    if( do_timing )then
      call cpu_time(time1)
      time3 = time1
    end if

    n_sites_supercell = size(positions,2)

!   Very inefficiently build neighbor lists. I should write some routine that performs                              <-- FIX THIS
!   overlapping domain decomposition to make this more efficient
!    if( a_box(2) == 0.d0 .and. a_box(3) == 0.d0 .and. b_box(1) == 0.d0 .and. &
!        b_box(3) == 0.d0 .and. c_box(1) == 0.d0 .and. c_box(2) == 0.d0 )then
!        .and. size(positions,2) == n_sites )then
!   This assumes that if the non-diagonal components of the lattice vectors are zero,
!   it is due to numerical noise, and thus the box is square
    if( dabs(a_box(2)) < d_tol .and. dabs(a_box(3)) < d_tol .and. &
        dabs(b_box(1)) < d_tol .and. dabs(b_box(3)) < d_tol .and. &
        dabs(c_box(1)) < d_tol .and. dabs(c_box(2)) < d_tol )then
      a_box(2) = 0.d0
      a_box(3) = 0.d0
      b_box(1) = 0.d0
      b_box(3) = 0.d0
      c_box(1) = 0.d0
      c_box(2) = 0.d0
      is_box_square = .true.
    else
      is_box_square = .false.
      if( rank == 0 .and. print_shape_warning )then
        print_shape_warning = .false.
        write(*,*)'                                       |'
        write(*,*)'WARNING: your simulation box is not    |  <-- WARNING'
        write(*,*)'orthorhombic; this will lead to slower |'
        write(*,*)'code execution. This warning will be   |'
        write(*,*)'printed only once. If you have more    |'
        write(*,*)'than one structure in your XYZ file,   |'
        write(*,*)'you may also have several instances of |'
        write(*,*)'non-orthorhombic cells.                |'
        write(*,*)'                                       |'
      end if
    end if
!
!   This is also inefficient <---------------------------------------- FIX THIS
    is_box_small = .false.
    if( any(indices > 1) )then
      is_box_small = .true.
      if( rank == 0 .and. print_cutoff_warning )then
        print_cutoff_warning = .false.
        write(*,*)'                                       |'
        write(*,*)'WARNING: your simulation box is smaller|  <-- WARNING'
        write(*,*)'than a neighbors cutoff sphere; this   |'
        write(*,*)'will lead to slow code execution due to|'
        write(*,*)'inefficient neighbor list builds. This |'
        write(*,*)'warning will be printed only once. If  |'
        write(*,*)'you have more than one structure in    |'
        write(*,*)'your XYZ file, you may also have       |'
        write(*,*)'several instances of non-orthorhombic  |'
        write(*,*)'cells.                                 |'
        write(*,*)'                                       |'
      end if
    end if
!
!   This is an initial guess for the maximum number of neighbors that we expect within a cutoff
!   It is used for the purpose of memory allocation. At the moment we are making this big, but
!   we should find a way to increase this value automatically whenever a bigger array is needed <---- fix
!   A way to fix it is, whenever an atom has more neighbors than n_neigh_max, to have store all
!   the neighbors lists into a temporary array, increase the size of the lists array, and then
!   store the lists back to it
    if( rebuild_neighbors_list )then
      n_neigh_max = 100
      allocate( neighbors_list(1:n_neigh_max*n_sites) )
      neighbors_list = 0
      allocate( n_neigh(1:n_sites) )
      n_neigh = 0
      n_atom_pairs = 0
    end if
!   We have an efficient algorithm for square boxes and inefficient for non-square boxes (sorry!)
!   Another requirement is that the minimum unit cell length is at least twice the cutoff
!
!   Tolerance in Angstrom for the ratio of positions to lattice vector. We add this because otherwise
!   atoms right at the periodic boundary can become problematic
    tol = 1.d-10
    if( is_box_square .and. (.not. is_box_small) .and. rebuild_neighbors_list )then
      mx = int( a_box(1) / rcut_max )
      my = int( b_box(2) / rcut_max )
      mz = int( c_box(3) / rcut_max )
      allocate( head(1:mx*my*mz) )
      head = 0
      allocate( this_list(1:n_sites) )
      do i = 1, n_sites
        call get_distance( [a_box(1)/2.d0, b_box(2)/2.d0, c_box(3)/2.d0], positions(1:3, i), &
                           a_box(1:3), b_box(1:3), c_box(1:3), (/ .true., .true., .true. /), dist, d, i_shift)
!       This is the position within the supercell, we must make sure it really is within the supercell
        dist = dist + [a_box(1)/2.d0, b_box(2)/2.d0, c_box(3)/2.d0]
        j = 1 + modulo(int( dist(1) / (a_box(1)+tol) * mx ), mx) &
              + modulo(int( dist(2) / (b_box(2)+tol) * my ), my) * mx &
              + modulo(int( dist(3) / (c_box(3)+tol) * mz ), mz) * my * mx
        this_list(i) = head(j)
        head(j) = i
      end do
      do i = 1, n_sites
        if( do_list(i) )then
!         We always count atom i as its own neighbor. This is useful when building the derivatives
          n_neigh(i) = n_neigh(i) + 1
          n_atom_pairs = n_atom_pairs + 1
          if( n_atom_pairs > n_neigh_max*n_sites )then
            allocate( neighbors_list_temp( 1:(n_neigh_max+10)*n_sites ) )
            neighbors_list_temp(1:n_atom_pairs-1) = neighbors_list(1:n_atom_pairs-1)
            deallocate( neighbors_list )
            allocate( neighbors_list( 1:(n_neigh_max+10)*n_sites ) )
            neighbors_list = neighbors_list_temp
            deallocate( neighbors_list_temp )
            n_neigh_max = n_neigh_max+10
          end if
          neighbors_list(n_atom_pairs) = i
!         Cell coordinates for this atom
          call get_distance( [a_box(1)/2.d0, b_box(2)/2.d0, c_box(3)/2.d0], positions(1:3, i), &
                             a_box(1:3), b_box(1:3), c_box(1:3), (/ .true., .true., .true. /), dist, d, i_shift)
          dist = dist + [a_box(1)/2.d0, b_box(2)/2.d0, c_box(3)/2.d0]
          i2 = 1 + int( dist(1) / (a_box(1)+tol) * mx )
          j2 = 1 + int( dist(2) / (b_box(2)+tol) * my )
          k2 = 1 + int( dist(3) / (c_box(3)+tol) * mz )
!         Look for other atoms in this and neighboring cells
          do k3 = k2-1, k2+1
            if( mz == 1 .and. k3 /= 1 )cycle
            if( mz == 2 .and. k2 == 1 .and. k3 == 0 )cycle
            if( mz == 2 .and. k2 == 2 .and. k3 == 3 )cycle
            do j3 = j2-1, j2+1
              if( my == 1 .and. j3 /= 1 )cycle
              if( my == 2 .and. j2 == 1 .and. j3 == 0 )cycle
              if( my == 2 .and. j2 == 2 .and. j3 == 3 )cycle
              do i3 = i2-1, i2+1
                if( mx == 1 .and. i3 /= 1 )cycle
                if( mx == 2 .and. i2 == 1 .and. i3 == 0 )cycle
                if( mx == 2 .and. i2 == 2 .and. i3 == 3 )cycle
                j = 1 + modulo(i3-1,mx) + modulo(j3-1,my)*mx + modulo(k3-1,mz)*mx*my
                k = head(j)
                do while( k /= 0 )
                  if( k /= i )then
                    call get_distance(positions(1:3, i), positions(1:3, k), a_box(1:3), b_box(1:3), &
                                      c_box(1:3), (/ .true., .true., .true. /), dist, d, i_shift)
                    if( d < rcut_max )then
                      n_neigh(i) = n_neigh(i) + 1
                      n_atom_pairs = n_atom_pairs + 1
                      if( n_atom_pairs > n_neigh_max*n_sites )then
                        allocate( neighbors_list_temp( 1:(n_neigh_max+10)*n_sites ) )
                        neighbors_list_temp(1:n_atom_pairs-1) = neighbors_list(1:n_atom_pairs-1)
                        deallocate( neighbors_list )
                        allocate( neighbors_list( 1:(n_neigh_max+10)*n_sites ) )
                        neighbors_list = neighbors_list_temp
                        deallocate( neighbors_list_temp )
                        n_neigh_max = n_neigh_max+10
                      end if
                      neighbors_list(n_atom_pairs) = k
                    end if
                  end if
                  k = this_list(k)
                end do
              end do
            end do
          end do
        end if
      end do
      deallocate(head, this_list)
!   Very inefficient algorithm for non-square boxes
    else if( rebuild_neighbors_list )then
      do i = 1, n_sites
        if( do_list(i) )then
!         We always count atom i as its own neighbor. This is useful when building the derivatives
          n_neigh(i) = n_neigh(i) + 1
          n_atom_pairs = n_atom_pairs + 1
          if( n_atom_pairs > n_neigh_max*n_sites )then
              allocate( neighbors_list_temp( 1:(n_neigh_max+10)*n_sites ) )
              neighbors_list_temp(1:n_atom_pairs-1) = neighbors_list(1:n_atom_pairs-1)
              deallocate( neighbors_list )
              allocate( neighbors_list( 1:(n_neigh_max+10)*n_sites ) )
              neighbors_list = neighbors_list_temp
              deallocate( neighbors_list_temp )
              n_neigh_max = n_neigh_max+10
          end if
          neighbors_list(n_atom_pairs) = i
          do j = 1, n_sites_supercell
            if( j /= i )then
              call get_distance(positions(1:3, i), positions(1:3, j), a_box(1:3), b_box(1:3), &
                                c_box(1:3), (/ .true., .true., .true. /), dist, d, i_shift)
              if( d < rcut_max )then
                n_neigh(i) = n_neigh(i) + 1
                n_atom_pairs = n_atom_pairs + 1
                if( n_atom_pairs > n_neigh_max*n_sites )then
                  allocate( neighbors_list_temp( 1:(n_neigh_max+10)*n_sites ) )
                  neighbors_list_temp(1:n_atom_pairs-1) = neighbors_list(1:n_atom_pairs-1)
                  deallocate( neighbors_list )
                  allocate( neighbors_list( 1:(n_neigh_max+10)*n_sites ) )
                  neighbors_list = neighbors_list_temp
                  deallocate( neighbors_list_temp )
                  n_neigh_max = n_neigh_max+10
                end if
!                j2 = mod(j-1, n_sites) + 1
!                neighbors_list(n_atom_pairs) = j2
                neighbors_list(n_atom_pairs) = j
              end if
            end if
          end do
        end if
      end do
    end if

    if( rebuild_neighbors_list )then
      allocate( neighbors_list_temp(1:n_atom_pairs) )
      neighbors_list_temp = neighbors_list(1:n_atom_pairs)
      deallocate( neighbors_list )
      allocate( neighbors_list(1:n_atom_pairs) )
      neighbors_list = neighbors_list_temp
      deallocate( neighbors_list_temp )
    end if

    if( do_timing )then
      call cpu_time(time2)
      neigh_time = time2 - time1
      time1 = time2
    end if

!   NOTE on performance: looping over interactions I could have chosen to calculate each interaction only once,
!   i.e., (i,j) = (j,i), because if j is i's neighbor, the reverse is also true. This would reduce the
!   calculation load by half. However, this is a design feature, since the code is easier to parallelize if
!   each atom is accompanied by its full list of neighbors. It also prevents ackward memory access.
!   In any case, I may want to rethink this in the future if I want to further gain extra performance (at the
!   expense of complicating the code, that is).
!
    n_atom_pairs = 0
    do i = 1, n_sites
      n_atom_pairs = n_atom_pairs + n_neigh(i)
    end do
!    allocate( mask_species(1:n_atom_pairs, 1:n_species) )
!    mask_species = .false.
    if( rebuild_neighbors_list )then
      allocate( rjs(1:n_atom_pairs) )
      allocate( xyz(1:3, 1:n_atom_pairs) )
      allocate( thetas(1:n_atom_pairs) )
      allocate( phis(1:n_atom_pairs) )
      allocate( neighbor_species(1:n_atom_pairs) )
    end if
    k2 = 0
    do i = 1, n_sites
      if( do_list(i) )then
        do k = 1, n_neigh(i)
          k2 = k2 + 1
          j = neighbors_list(k2)
          if( k == 1 )then
            rjs(k2) = 0.d0
            xyz(1:3, k2) = (/ 0.d0, 0.d0, 0.d0 /)
            thetas(k2) = 0.d0
            phis(k2) = 0.d0
            neighbor_species(k2) = species_supercell(i)
!            do i2 = 1, species_multiplicity(i)
!              mask_species(k2, species_supercell(i2, j)) = .true.
!            end do
          else
            call get_distance(positions(1:3,i), positions(1:3,j), a_box(1:3), b_box(1:3), &
                              c_box(1:3), (/ .true., .true., .true. /), dist, d, i_shift)
            rjs(k2) = d
            xyz(1:3, k2) = dist
!           Avoid numerical artifacts
            if( dabs(dist(3)) >= d )then
              if( dist(3) > 0.d0 )then
                thetas(k2) = 0.d0
              else
                thetas(k2) = dacos(-1.d0)
              end if
            else
              thetas(k2) = dacos( dist(3) / d )
            end if
            phis(k2) = datan2( dist(2), dist(1) )
            neighbor_species(k2) = species_supercell(j)
!            do i2 = 1, species_multiplicity(i)
!              mask_species(k2, species_supercell(i2, j)) = .true.
!            end do
          end if
        end do
      else
        k2 = k2 + n_neigh(i)
      end if
    end do

    if( do_timing )then
      call cpu_time(time2)
      write(*,*)'                                       |'
      write(*,*)'Atoms timings (build):                 |'
      write(*,*)'                                       |'
      write(*,'(A, F9.3, A)') '  *) Neighbors build: ', neigh_time, ' seconds |'
      write(*,'(A, F7.3, A)') '  *) Spherical coords.: ', time2-time1, ' seconds |'
      write(*,'(A, F19.3, A)') '  *) Total: ', time2-time3, ' seconds |'
      write(*,*)'                                       |'
      write(*,*)'.......................................|'
    end if

  end subroutine
!**************************************************************************

!**************************************************************************
  subroutine get_number_of_atom_pairs( n_neigh, rjs, rcut, l_max, n_max, max_Gbytes_per_process, &
                                       i_beg_list, i_end_list, j_beg_list, j_end_list )

    implicit none

!   Input variables
    real*8, intent(in) :: rjs(:), rcut, max_Gbytes_per_process
    integer, intent(in) :: n_neigh(:), l_max, n_max

!   Output variables
    integer, allocatable, intent(out) :: i_beg_list(:), i_end_list(:), j_beg_list(:), j_end_list(:)

!   Internal variables
    real*8 :: estimated_memory_in_Gbytes, mem_ratio, pairs_per_chunk
    integer :: n_sites, n_atom_pairs, k_max, n_chunks, i, j, k, k2, i_chunk, n_atom_pairs_in

    n_sites = size(n_neigh)
    n_atom_pairs = size(rjs)

    k = 0
    n_atom_pairs_in = 0
    do i = 1, n_sites
      do j = 1, n_neigh(i)
        k = k + 1
        if( rjs(k) < rcut )then
          n_atom_pairs_in = n_atom_pairs_in + 1
        end if
      end do
    end do

    k_max = 1 + l_max*(l_max+1)/2 + l_max
!   This is a conservative estimate of the maximum memory that this run will need
    estimated_memory_in_Gbytes = dfloat(n_atom_pairs_in) * dfloat(n_max) * dfloat(k_max) * 150.d0 / 1024.d0**3
    mem_ratio = estimated_memory_in_Gbytes/max_Gbytes_per_process
    n_chunks = ceiling(mem_ratio)
    if( n_chunks > n_sites )then
      n_chunks = n_sites
    end if

    pairs_per_chunk = float(n_atom_pairs_in) / float(n_chunks)

    allocate( i_beg_list(1:n_chunks) )
    allocate( i_end_list(1:n_chunks) )
    allocate( j_beg_list(1:n_chunks) )
    allocate( j_end_list(1:n_chunks) )

    i_beg_list(1) = 1
    j_beg_list(1) = 1
    i_end_list(n_chunks) = n_sites
    j_end_list(n_chunks) = n_atom_pairs

    if( n_chunks == 1 )then
      return
    end if

    k = 0
    k2 = 0
    i_chunk = 1
    do i = 1, n_sites
      do j = 1, n_neigh(i)
        k = k + 1
        if( rjs(k) < rcut )then
          k2 = k2 + 1
        end if
      end do
      if( k2 >= int( float(i_chunk)*pairs_per_chunk ) )then
        i_end_list(i_chunk) = i
        j_end_list(i_chunk) = k
        i_chunk = i_chunk + 1
        i_beg_list(i_chunk) = i+1
        j_beg_list(i_chunk) = k+1
        if( i_chunk == n_chunks )then
          exit
        end if
      end if
    end do
    return

  end subroutine
!**************************************************************************

end module neighbors

max_stars_repo_name: richardjana/turbogap

max_issues_repo_name: richardjana/turbogap

max_forks_repo_name: richardjana/turbogap


text:
The University Club Conference Center, often referred to as the http://uclub.ucdavis.edu/Uclub/ U Club, is a Campus Buildings building on the UC Davis campus just west of Wyatt Pavilion.  It features two primary rooms: the Clubroom and the Lounge (2436 and 1680 square feet, respectively).  The center is used for conferences and meetings, and is also the site of many Art Department guest lectures and Department of Theatre and Dance choreography/movement classes and rehearsals.

The space features a functional kitchen, which is used primarily by catering staff now. However, the club used to serve daily weekday lunches to mostly faculty and guests in the 70s and 80s, when it was aptly named the Faculty Club. As lunchtime faculty diners dwindled, the university changed the clubs name and opened its doors to the entire campus. Revenues continued to decline, until the club stopped serving daily lunches and focused more on conferences and events in the late 90s/early 00s.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module mb_data
implicit none
type::blk1
integer(4) ni,nj,nk
integer(4) num_i0,num_idim,num_j0,num_jdim,num_k0,num_kdim
real(8),allocatable,dimension(:,:,:)::x,y,z
integer(4),allocatable,dimension(:,:,:)::ibank
end type blk1

type::bc1
integer(4) key,blk,is,ie,js,je,ks,ke,face,side
character(20) typ,group,nam
integer(4) typ_num
integer(4) bcg
end type bc1

type::if1
integer(4) key,blk1,is,ie,js,je,ks,ke,blk2,pair
integer(4) face,side
end type if1

type::ifp1
integer(4) blk1,is1,ie1,js1,je1,ks1,ke1
integer(4) blk2,is2,ie2,js2,je2,ks2,ke2
integer(4) ii1,jj1,kk1,ss1
integer(4) ii2,jj2,kk2,ss2
integer(4) i1,j1,k1
integer(4) i2,j2,k2
real(8) x1(4),y1(4),z1(4)
real(8) x2(4),y2(4),z2(4)
integer(4) is3,ie3,js3,je3,ks3,ke3
integer(4) jj3,kk3
integer(4) link,dir1,dir2
end type ifp1

end module mb_data

program cfdrc_to_cfl3d
use mb_data
implicit none
namelist/cfdrcfile/file,pfgd,ibc
namelist/nondimension/cref,bref,sref,xmc,ymc,zmc
namelist/flowcondition/xmach,alpha,beta,pinf,tinf
namelist/solvercontrol/i2d,mglevg,irest,ncyc,dtmax,iflagts,fmax,ifullns,cprec,ifds,epsa_r,rkap0,idiag,iflim,ivisc,edvislim,iturbord,nwrest
namelist/output/isurfg,surfg
namelist/bc2002/bc2002name,pratio
namelist/controlsurface/csname
character(20) csname
!integer(4) ncs
character(20) bc2002name
real(8) pratio
character(20) file,pfgd,ibc
character(30) file1,file2
real(8) cref,bref,sref,xmc,ymc,zmc
real(8) xmach,alpha,beta,pinf,tinf
integer(4) i2d,mglevg,irest,ncyc(3),iflagts,ifullns,ifds,idiag,iflim,ivisc,iturbord
real(8) dtmax,fmax,cprec,epsa_r,rkap0,edvislim
real(8) muref,tref,tsut,csut,muinf,cinf,deninf
integer(4) isurfg
character(20) surfg
real(8),parameter::pi=3.141592653589793238d0
real(8) tinf_dr
real(8) reue
integer(4) ialph,ihstry
real(8) dt,cfltau
integer(4) iunst
integer(4) ngrid,nplot3d,nprint,nwrest,ichk,ntstep,ita
type(blk1),allocatable,target,dimension(:)::b
type(blk1),pointer::b1,b2
type(bc1),allocatable,target,dimension(:)::c
type(bc1),pointer::c1,c2
type(if1),allocatable,target,dimension(:)::f
type(if1),pointer::f1,f2
type(ifp1),allocatable,target,dimension(:)::p
type(ifp1),pointer::p1
character(20) ccc,ccc1,ccc2,ccc3
integer(4) num_blk,num_bc,num_if,num_ifp
integer(4) num_bcg
integer(4) i,j,k,ii,jj,kk
integer(4) d(3,3)
real(8) dist1,dist2,dmin
character(20),allocatable::bc_name(:)
integer(4),allocatable::ftyp(:)
integer(4) iflag
integer(4) ncg,iem,iadvance
character(3) iforce
integer(4) ilamlo,ilamhi,jlamlo,jlamhi,klamlo,klamhi
integer(4) inewg,igridc
integer(4) iovrlp
integer(4) ndata
real(8) twratio,cq
integer(4) mseq,mgflag,iconsf,mtt,ngam
integer(4) issc,issr
real(8) epsssc,epsssr
integer(4) nemgl,nitfo
integer(4) mit
integer(4) ninter
integer(4) iptype
integer(4) imovie
integer(4) ncs
integer(4) mglevg_max,mglevg_allow,imgflag,mglevg_allow_1,mglevg_allow_2,mglevg_allow_3
integer(4) nii,njj
integer(4),allocatable::bcg_num(:)
integer(4),allocatable::bcg_pair(:)
integer(4),allocatable::ipatch(:)
integer(4) ipatchflag
integer(4) ioflag,itrace
integer(4) iifit,llmit,iitmax,mmcxie,mmceta,iic0,iiorph
integer(4) iis,iie,jjs,jje

imgflag=0
mglevg_allow=1
mglevg_allow_1=1
mglevg_allow_2=1
mglevg_allow_3=1
mglevg_max=3

ipatchflag=0

open(501,file='cfdrc_to_cfl3d-v1.inp')
open(601,file='cfl3d.inp')
open(602,file='cfl3d.grd',form='unformatted')

read(501,NML=cfdrcfile)

file1=trim(file)//'.'//trim(pfgd)
file2=trim(file)//'.'//trim(ibc)
write(*,*) 'INPUT CFDRC GRID FILES:'
write(*,*) file1,file2
open(502,file=file1)
open(503,file=file2)

cref=1d0
bref=-100d0
sref=-100d0
xmc=0d0
ymc=0d0
zmc=0d0

read(501,NML=nondimension)

if(bref<0d0) then
bref=cref

endif

if(sref<0d0) then
sref=pi*cref**2/4d0
endif
xmach=2d0;alpha=0d0;beta=0d0;pinf=101325d0;tinf=300d0
read(501,NML=flowcondition)
i2d=0;irest=0;ifullns=0;cprec=0d0;ifds=1
epsa_r=0.3d0;rkap0=0.3333d0;idiag=0;iflim=4;ivisc=-5;edvislim=100000d0;iturbord=1
read(501,NML=solvercontrol)
isurfg=0;surfg='none'
read(501,NML=output)
read(501,NML=bc2002)
read(501,NML=controlsurface)

muref=1.716d-5
tref=273.15d0
tsut=110.4d0 !
csut=muref*(tref+tsut)/tref**1.5 !1.458d-6
muinf=csut*tinf**1.5/(tinf+tsut)
cinf=dsqrt(1.4d0*286.9d0*tinf)
deninf=pinf/286.9d0/tinf
reue=deninf*xmach*cinf/muinf/1d6
tinf_dr=tinf*9d0/5d0

ialph=0
ihstry=2

dt=dtmax/fmax
iunst=0
cfltau=5d0

!!!read ibc
read(503,*) ccc
read(503,*) ccc
read(503,*) num_if
if(num_if/=0)then
allocate(f(num_if))
num_ifp=num_if/2
allocate(p(num_ifp))

do i=1,num_if
f1=>f(i)
read(503,*) f1%key,f1%blk1,f1%is,f1%ie,f1%js,f1%je,f1%ks,f1%ke,f1%blk2
f1%pair=-1

if(f1%is==f1%ie) then
f1%face=1
if(f1%is>1)then
f1%side=2
else
f1%side=1
endif
elseif(f1%js==f1%je) then
f1%face=2
if(f1%js>1)then
f1%side=2
else
f1%side=1
endif
elseif(f1%ks==f1%ke) then
f1%face=3
if(f1%ks>1)then
f1%side=2
else
f1%side=1
endif
else
write(*,*) 'IBC interface error!'
stop
endif
enddo
num_blk=f1%blk1 !!!
endif

read(503,*) num_bc
if(num_bc/=0)then
allocate(c(num_bc),ftyp(num_blk))
allocate(bc_name(num_bc))
allocate(bcg_num(num_bc))
allocate(bcg_pair(num_bc))
allocate(ipatch(num_bc))

do i=1,num_blk
ftyp(i)=0
enddo

num_bcg=0

do i=1,num_bc
c1=>c(i)
c1%bcg=0
read(503,*) c1%key,c1%blk,c1%is,c1%ie,c1%js,c1%je,c1%ks,c1%ke,c1%typ,c1%group,c1%nam
if(c1%typ=='Symmetry') then
c1%typ_num=1001
elseif(c1%typ=='Inlet') then
c1%typ_num=1000
elseif(c1%typ=='Outlet') then
c1%typ_num=1002
elseif(c1%typ=='Wall') then
c1%typ_num=2004
ftyp(c1%blk)=1
c1%bcg=1
elseif(c1%typ=='R_wall') then
c1%typ_num=1013 !!! singular
elseif(c1%typ=='Thin_wall') then
c1%typ_num=0 !!! patched interface
else
write(*,*) 'IBC bc error!'
write(*,*) 'Undefined bc type:'
write(*,*) c1%typ
stop
endif

!!! bc2002
if(c1%nam==bc2002name) then
c1%typ_num=2002
endif

if(c1%is==c1%ie) then
c1%face=1
if(c1%is>1)then
c1%side=2
else
c1%side=1
endif
elseif(c1%js==c1%je) then
c1%face=2
if(c1%js>1)then
c1%side=2
else
c1%side=1
endif
elseif(c1%ks==c1%ke) then
c1%face=3
if(c1%ks>1)then
c1%side=2
else
c1%side=1
endif
else
write(*,*) 'IBC bc error!'
stop
endif

iflag=0
if(num_bcg==0) then
bc_name(1)=c1%nam
c1%bcg=1
else

do j=1,num_bcg

if(bc_name(j)==c1%nam) then
iflag=1
c1%bcg=j
endif

enddo

endif

if(iflag/=1) then
num_bcg=num_bcg+1
bc_name(num_bcg)=c1%nam
c1%bcg=num_bcg
endif

enddo

do j=1,num_bcg
bcg_num(j)=0
bcg_pair(j)=0
ipatch(j)=0

do i=1,num_bc
c1=>c(i)
if(c1%bcg==j) then
bcg_num(j)=bcg_num(j)+1
if(c1%typ=='Thin_wall') then
ipatch(j)=1
endif
endif

enddo

if(ipatch(j)==1) then

do k=1,num_bcg

if(k/=j)then

if(bc_name(j)(3:20)==bc_name(k)(3:20)) then
bcg_pair(j)=k
endif

endif

enddo

endif

enddo

ninter=0
do j=1,num_bcg
if(ipatch(j)==1) then

do i=1,num_bc

c1=>c(i)
if(c1%bcg==j) then
ninter=ninter+1

do k=1,num_bc
c2=>c(k)

if(c2%bcg==bcg_pair(j)) then

endif
enddo

endif
enddo

endif
enddo

do j=1,num_bcg
if(ipatch(j)==1) then
ipatchflag=1
endif
enddo

endif

if(i2d /= 1) then

open(603,file='plot3d.fvbnd')

!!! gen plot3d.fvbnd
write(*,*) 'BC GROUPS:',num_bcg
write(603,*) 'FVBND 1 4'
do i=1,num_bcg
write(*,*) i,bc_name(i)
write(603,*) bc_name(i)
enddo
write(603,*) 'BOUNDARIES'
do i=1,num_bcg
do j=1,num_bc
c1=>c(j)
if(c1%bcg==i)then
write(603,'(8I5,A5)') i,c1%blk,c1%is,c1%ie,c1%js,c1%je,c1%ks,c1%ke,'F 0'
endif
enddo
enddo

write(*,*) '3D case'
write(*,*) 'Generate plot3d.fvbnd'

else

write(*,*) '2D case'

endif

!!!read pfgd!!! gen cfl3d.grd
read(502,*) num_blk
write(602) num_blk !!!cfl3d.grd
allocate(b(num_blk))
do i=1,num_blk
b1=>b(i)
read(502,*) b1%ni,b1%nj,b1%nk
allocate(b1%x(b1%ni,b1%nj,b1%nk),b1%y(b1%ni,b1%nj,b1%nk),b1%z(b1%ni,b1%nj,b1%nk))
allocate(b1%ibank(b1%ni,b1%nj,b1%nk))
b1%num_i0=0
b1%num_idim=0
b1%num_j0=0
b1%num_jdim=0
b1%num_k0=0
b1%num_kdim=0
enddo

write(602) (b(i)%ni,b(i)%nj,b(i)%nk,i=1,num_blk)

do ii=1,num_blk
b1=>b(ii)
read(502,*) (((b1%x(i,j,k),i=1,b1%ni),j=1,b1%nj),k=1,b1%nk),(((b1%y(i,j,k),i=1,b1%ni),j=1,b1%nj),k=1,b1%nk),(((b1%z(i,j,k),i=1,b1%ni),j=1,b1%nj),k=1,b1%nk),(((b1%ibank(i,j,k),i=1,b1%ni),j=1,b1%nj),k=1,b1%nk)
write(602) (((b1%x(i,j,k),i=1,b1%ni),j=1,b1%nj),k=1,b1%nk),(((b1%y(i,j,k),i=1,b1%ni),j=1,b1%nj),k=1,b1%nk),(((b1%z(i,j,k),i=1,b1%ni),j=1,b1%nj),k=1,b1%nk)
enddo
!!!

ncs=0
if(num_bc /= 0) then
do i=1,num_bc
c1=>c(i)
b1=>b(c1%blk)

!!! control surface
if(c1%nam==csname) then
ncs=ncs+1
endif

if(c1%face==1) then

if(c1%side==1) then
b1%num_i0=b1%num_i0+1
else
b1%num_idim=b1%num_idim+1
endif

elseif(c1%face==2) then

if(c1%side==1) then
b1%num_j0=b1%num_j0+1
else
b1%num_jdim=b1%num_jdim+1
endif

else

if(c1%side==1) then
b1%num_k0=b1%num_k0+1
else
b1%num_kdim=b1%num_kdim+1
endif

endif

enddo
endif
!!!
if(num_if /= 0) then
do i=1,num_if
f1=>f(i)
b1=>b(f1%blk1)
!!! control surface
!if(f1%nam==csname) then
!ncs=ncs+1
!endif

if(f1%face==1) then

if(f1%side==1) then
b1%num_i0=b1%num_i0+1
else
b1%num_idim=b1%num_idim+1
endif

elseif(f1%face==2) then

if(f1%side==1) then
b1%num_j0=b1%num_j0+1
else
b1%num_jdim=b1%num_jdim+1
endif

else

if(f1%side==1) then
b1%num_k0=b1%num_k0+1
else
b1%num_kdim=b1%num_kdim+1
endif

endif

enddo
endif

!write(*,*) ncs

!!!
ngrid=-num_blk
nplot3d=num_blk
nprint=0
ichk=0
ntstep=1
ita=-2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(mglevg>mglevg_max)then
mglevg=mglevg_max
!write(*,*) 'Auto change mglevg due to mglevg_max=3!'
endif

if(i2d==0) then !!! MGL 3D

mglevg_allow=1

!mg_block check
if(mglevg /= 1) then
do j=2,mglevg_max
do i=1,num_blk
b1=>b(i)

if(mod(b1%ni-1,2**(j-1))/=0)then
imgflag=1
endif
if(mod(b1%nj-1,2**(j-1))/=0)then
imgflag=1
endif
if(mod(b1%nk-1,2**(j-1))/=0)then
imgflag=1
endif

enddo
if(imgflag/=1)then
mglevg_allow=j
else
imgflag=0
endif
enddo

!write(*,*) 'mg block check, allowable mg level',mglevg_allow

if(mglevg>mglevg_allow)then
mglevg=mglevg_allow
!write(*,*) 'Auto change mglevg due to block dimensions!'
endif

endif

mglevg_allow=1

!mg_bc check
if(num_bc /= 0) then
if(mglevg /= 1) then
do j=2,mglevg_max
do i=1,num_bc
c1=>c(i)
if(c1%face==1) then
nii=c1%je-c1%js+1
njj=c1%ke-c1%ks+1
elseif(c1%face==2) then
nii=c1%ke-c1%ks+1
njj=c1%ie-c1%is+1
else
nii=c1%ie-c1%is+1
njj=c1%je-c1%js+1
endif

if(mod(nii-1,2**(j-1))/=0)then
imgflag=1
endif

if(mod(njj-1,2**(j-1))/=0)then
imgflag=1
endif

enddo

if(imgflag/=1)then
mglevg_allow=j
else
imgflag=0
endif

enddo

!write(*,*) 'mg bc check, allowable mg level',mglevg_allow

if(mglevg>mglevg_allow)then
mglevg=mglevg_allow
!write(*,*) 'Auto change mglevg due to bc dimensions!'
endif
endif
endif

mglevg_allow=1

if(num_if /= 0) then
if(mglevg /= 1) then
do j=2,mglevg_max
do i=1,num_if
f1=>f(i)
if(f1%face==1) then
nii=f1%je-f1%js+1
njj=f1%ke-f1%ks+1
elseif(f1%face==2) then
nii=f1%ke-f1%ks+1
njj=f1%ie-f1%is+1
else
nii=f1%ie-f1%is+1
njj=f1%je-f1%js+1
endif
!if(mod(nii-1,(j-1)*2)/=0)then
if(mod(nii-1,2**(j-1))/=0)then
imgflag=1
endif
!if(mod(njj-1,(j-1)*2)/=0)then
if(mod(njj-1,2**(j-1))/=0)then
imgflag=1
endif
!if(nii<5+(j-2)*4) then
!imgflag=1
!endif
!if(njj<5+(j-2)*4) then
!imgflag=1
!endif
enddo
if(imgflag/=1)then
mglevg_allow=j
else
imgflag=0
endif
enddo
!write(*,*) 'mg interface check, allowable mg level',mglevg_allow
if(mglevg>mglevg_allow)then
mglevg=mglevg_allow
!write(*,*) 'Auto change mglevg due to interface dimensions!'
endif
endif
endif

else !!! 2D case

mglevg_allow=1

!mg_block check
if(mglevg /= 1) then
do j=2,mglevg_max
do i=1,num_blk
b1=>b(i)
if(mod(b1%nj-1,2**(j-1))/=0)then
imgflag=1
endif
if(mod(b1%nk-1,2**(j-1))/=0)then
imgflag=1
endif

enddo
if(imgflag/=1)then
mglevg_allow=j
else
imgflag=0
endif

enddo

!write(*,*) 'mg block check, allowable mg level',mglevg_allow

if(mglevg>mglevg_allow)then
mglevg=mglevg_allow
write(*,*) 'Auto change mglevg due to block dimensions!'
endif

endif

mglevg_allow=1

!mg_bc check
if(num_bc /= 0) then
if(mglevg /= 1) then
do j=2,mglevg_max
do i=1,num_bc
c1=>c(i)

if(c1%face==2) then
nii=c1%ke-c1%ks+1
njj=c1%ie-c1%is+1

if(mod(nii-1,2**(j-1))/=0)then
imgflag=1
endif

elseif(c1%face==3) then
nii=c1%ie-c1%is+1
njj=c1%je-c1%js+1

if(mod(njj-1,2**(j-1))/=0)then
imgflag=1
endif

endif

enddo

if(imgflag/=1)then
mglevg_allow=j
else
imgflag=0
endif

enddo

!write(*,*) 'mg bc check, allowable mg level',mglevg_allow

if(mglevg>mglevg_allow)then
mglevg=mglevg_allow
write(*,*) 'Auto change mglevg due to bc dimensions!'
endif
endif
endif

mglevg_allow=1

if(num_if /= 0) then
if(mglevg /= 1) then
do j=2,mglevg_max
do i=1,num_if
f1=>f(i)

if(f1%face==2) then
nii=f1%ke-f1%ks+1
njj=f1%ie-f1%is+1

if(mod(nii-1,2**(j-1))/=0)then
imgflag=1
endif

elseif(f1%face==3) then
nii=f1%ie-f1%is+1
njj=f1%je-f1%js+1

if(mod(njj-1,2**(j-1))/=0)then
imgflag=1
endif

endif

enddo

if(imgflag/=1)then
mglevg_allow=j
else
imgflag=0
endif
enddo
!write(*,*) 'mg interface check, allowable mg level',mglevg_allow
if(mglevg>mglevg_allow)then
mglevg=mglevg_allow
write(*,*) 'Auto change mglevg due to interface dimensions!'
endif
endif
endif

!write(*,*) mglevg_allow_1,mglevg_allow_2,mglevg_allow_3

endif !!!
!write(*,*) 'Current mglevg=',mglevg

write(*,*) 'Current mglevg=',mglevg

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!cfl3d.inp

write(601,'(A)') 'I/O FILES'
write(601,'(A)') 'cfl3d.grd'
write(601,'(A)') 'plot3dg.bin'
write(601,'(A)') 'plot3dq.bin'
write(601,'(A)') 'cfl3d.out'
write(601,'(A)') 'cfl3d.res'
write(601,'(A)') 'cfl3d.turres'
write(601,'(A)') 'cfl3d.blomax'
write(601,'(A)') 'cfl3d.out15'
write(601,'(A)') 'cfl3d.prout'
write(601,'(A)') 'cfl3d.out20'
write(601,'(A)') 'ovrlp.bin'
write(601,'(A)') 'patch.bin'
write(601,'(A)') 'restart.bin'
write(601,'(A)') '>--begin keyword-driven input section'
write(601,'(A,I3)') 'ifullns',ifullns
write(601,'(A,F5.1)') 'cprec',cprec
write(601,'(A,F5.1)') 'epsa_r',epsa_r
write(601,'(A,F10.1)') 'edvislim',edvislim
write(601,'(A,I3)') 'iturbord',iturbord
write(601,'(A)') '<--end keyword-driven input section'
write(601,'(A)') 'CFDRC_FILES:'//trim(file1)//'&'//trim(file2)
write(601,'(7A10)') 'XMACH','ALPHA','BETA','REUE,MIL','TINF,DR','IALPH','IHSTRY'
write(601,'(5F10.4,2I10)') xmach,alpha,beta,reue,tinf_dr,ialph,ihstry
write(601,'(6A10)') 'SREF','CREF','BREF','XMC','YMC','ZMC'
write(601,'(6F10.6)') sref,cref,bref,xmc,ymc,zmc
write(601,'(6A10)') 'DT','IREST','IFLAGTS','FMAX','IUNST','CFLTAU'
write(601,'(F10.4,2I10,F10.4,I10,F10.4)') dt,irest,iflagts,fmax,iunst,cfltau
write(601,'(8A10)') 'NGRID','NPLOT3D','NPRINT','NWREST','ICHK','I2D','NTSTEP','ITA'
write(601,'(8I10)') ngrid,nplot3d,nprint,nwrest,ichk,i2d,ntstep,ita

write(601,'(7A10)') 'NCG','IEM','IADVANCE','IFORCE','IVISC(I)','IVISC(J)','IVISC(K)'

ncg=mglevg-1

iem=0
iadvance=0
do i=1,num_blk
if(ftyp(i)==1) then
iforce='333'
else
iforce='000'
endif
write(601,'(3I10,A10,3I10)') ncg,iem,iadvance,iforce,ivisc,ivisc,ivisc
enddo

write(601,'(3A10)') 'IDIM','JDIM','KDIM'
do i=1,num_blk
b1=>b(i)
write(601,'(3I10)') b1%ni,b1%nj,b1%nk
enddo

write(601,'(6A10)') 'ILAMLO','ILAMHI','JLAMLO','JLAMHI','KLAMLO','KLAMHI'
ilamlo=0
ilamhi=0
jlamlo=0
jlamhi=0
klamlo=0
klamhi=0
do i=1,num_blk
write(601,'(6I10)') ilamlo,ilamhi,jlamlo,jlamhi,klamlo,klamhi
enddo

write(601,'(8A10)') 'INEWG','IGRIDC','IS','JS','KS','IE','JE','KE'
inewg=0
igridc=0
do i=1,num_blk
write(601,'(8I10)') inewg,igridc,inewg,igridc,inewg,igridc,inewg,igridc
enddo

write(601,'(6A10)') 'IDIAG(I)','IDIAG(J)','IDIAG(K)','IFLIM(I)','IFLIM(J)','IFLIM(K)'
do i=1,num_blk
write(601,'(6I10)') idiag,idiag,idiag,iflim,iflim,iflim
enddo

write(601,'(6A10)') 'IFDS(I)','IFDS(J)','IFDS(K)','RKAP0(I)','RKAP0(J)','RKAP0(K)'
do i=1,num_blk
write(601,'(3I10,3F10.4)') ifds,ifds,ifds,rkap0,rkap0,rkap0
enddo

write(601,'(8A10)') 'GRID','NBCI0','NBCIDIM','NBCJ0','NBCJDIM','NBCK0','NBCKDIM','IOVRLP'
iovrlp=0
do i=1,num_blk
b1=>b(i)
write(601,'(8I10)') i,b1%num_i0,b1%num_idim,b1%num_j0,b1%num_jdim,b1%num_k0,b1%num_kdim,iovrlp
enddo

twratio=0d0
cq=0d0

write(601,'(8A10)') 'I0:   GRID','SEGMENT','BCTYPE','JSTA','JEND','KSTA','KEND','NDATA'

do i=1,num_blk
k=1
if(num_if /= 0) then
do j=1,num_if
f1=>f(j)
if(f1%blk1==i) then
if(f1%face==1) then
if(f1%side==1) then
ndata=0
write(601,'(8I10)') i,k,ndata,f1%js,f1%je,f1%ks,f1%ke,ndata
k=k+1
endif
endif
endif

enddo
endif

if(num_bc /= 0) then

do j=1,num_bc
c1=>c(j)
if(c1%blk==i) then
if(c1%face==1) then
if(c1%side==1) then

if(c1%typ_num==2004)then
ndata=2

if(isurfg==0)then

kk=1

elseif(isurfg<0)then

if(c1%nam==surfg)then
kk=-1
else
kk=1
endif

else

if(c1%nam==surfg)then
kk=1
else
kk=-1
endif

endif

write(601,'(8I10)') i,kk*k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
write(601,'(2A10)') 'tw/tinf','cq'
write(601,'(2F10.4)') twratio,cq

elseif(c1%typ_num==2002)then
ndata=1
write(601,'(8I10)') i,kk*k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
write(601,'(A10)') 'p/pinf'
write(601,'(F10.4)') pratio

else
ndata=0
write(601,'(8I10)') i,k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
endif
k=k+1
endif
endif
endif

enddo

endif

enddo

write(601,'(8A10)') 'IDIM: GRID','SEGMENT','BCTYPE','JSTA','JEND','KSTA','KEND','NDATA'

do i=1,num_blk
k=1
if(num_if /= 0) then
do j=1,num_if
f1=>f(j)
if(f1%blk1==i) then
if(f1%face==1) then
if(f1%side==2) then
ndata=0
write(601,'(8I10)') i,k,ndata,f1%js,f1%je,f1%ks,f1%ke,ndata
k=k+1
endif
endif
endif

enddo
endif

if(num_bc /= 0) then

do j=1,num_bc
c1=>c(j)
if(c1%blk==i) then
if(c1%face==1) then
if(c1%side==2) then

if(c1%typ_num==2004)then
ndata=2

if(isurfg==0)then

kk=1

elseif(isurfg<0)then

if(c1%nam==surfg)then
kk=-1
else
kk=1
endif

else

if(c1%nam==surfg)then
kk=1
else
kk=-1
endif

endif

write(601,'(8I10)') i,kk*k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
write(601,'(2A10)') 'tw/tinf','cq'
write(601,'(2F10.4)') twratio,cq

elseif(c1%typ_num==2002)then
ndata=1
write(601,'(8I10)') i,kk*k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
write(601,'(A10)') 'p/pinf'
write(601,'(F10.4)') pratio

else
ndata=0
write(601,'(8I10)') i,k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
endif
k=k+1
endif
endif
endif

enddo

endif

enddo

write(601,'(8A10)') 'J0:   GRID','SEGMENT','BCTYPE','ISTA','IEND','KSTA','KEND','NDATA'

do i=1,num_blk
k=1
if(num_if /= 0) then
do j=1,num_if
f1=>f(j)
if(f1%blk1==i) then
if(f1%face==2) then
if(f1%side==1) then
ndata=0
write(601,'(8I10)') i,k,ndata,f1%is,f1%ie,f1%ks,f1%ke,ndata
k=k+1
endif
endif
endif

enddo
endif

if(num_bc /= 0) then

do j=1,num_bc
c1=>c(j)
if(c1%blk==i) then
if(c1%face==2) then
if(c1%side==1) then

if(c1%typ_num==2004)then
ndata=2

if(isurfg==0)then

kk=1

elseif(isurfg<0)then

if(c1%nam==surfg)then
kk=-1
else
kk=1
endif

else

if(c1%nam==surfg)then
kk=1
else
kk=-1
endif

endif

write(601,'(8I10)') i,kk*k,c1%typ_num,c1%is,c1%ie,c1%ks,c1%ke,ndata
write(601,'(2A10)') 'tw/tinf','cq'
write(601,'(2F10.4)') twratio,cq

elseif(c1%typ_num==2002)then
ndata=1
write(601,'(8I10)') i,kk*k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
write(601,'(A10)') 'p/pinf'
write(601,'(F10.4)') pratio

else
ndata=0
write(601,'(8I10)') i,k,c1%typ_num,c1%is,c1%ie,c1%ks,c1%ke,ndata
endif
k=k+1
endif
endif
endif

enddo

endif

enddo
write(601,'(8A10)') 'JDIM: GRID','SEGMENT','BCTYPE','ISTA','IEND','KSTA','KEND','NDATA'

do i=1,num_blk
k=1
if(num_if /= 0) then
do j=1,num_if
f1=>f(j)
if(f1%blk1==i) then
if(f1%face==2) then
if(f1%side==2) then
ndata=0
write(601,'(8I10)') i,k,ndata,f1%is,f1%ie,f1%ks,f1%ke,ndata
k=k+1
endif
endif
endif

enddo
endif

if(num_bc /= 0) then

do j=1,num_bc
c1=>c(j)
if(c1%blk==i) then
if(c1%face==2) then
if(c1%side==2) then

if(c1%typ_num==2004)then
ndata=2

if(isurfg==0)then

kk=1

elseif(isurfg<0)then

if(c1%nam==surfg)then
kk=-1
else
kk=1
endif

else

if(c1%nam==surfg)then
kk=1
else
kk=-1
endif

endif

write(601,'(8I10)') i,kk*k,c1%typ_num,c1%is,c1%ie,c1%ks,c1%ke,ndata
write(601,'(2A10)') 'tw/tinf','cq'
write(601,'(2F10.4)') twratio,cq

elseif(c1%typ_num==2002)then
ndata=1
write(601,'(8I10)') i,kk*k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
write(601,'(A10)') 'p/pinf'
write(601,'(F10.4)') pratio

else
ndata=0
write(601,'(8I10)') i,k,c1%typ_num,c1%is,c1%ie,c1%ks,c1%ke,ndata
endif
k=k+1
endif
endif
endif

enddo

endif

enddo

write(601,'(8A10)') 'K0:   GRID','SEGMENT','BCTYPE','ISTA','IEND','JSTA','JEND','NDATA'

do i=1,num_blk
k=1
if(num_if /= 0) then
do j=1,num_if
f1=>f(j)
if(f1%blk1==i) then
if(f1%face==3) then
if(f1%side==1) then
ndata=0
write(601,'(8I10)') i,k,ndata,f1%is,f1%ie,f1%js,f1%je,ndata
k=k+1
endif
endif
endif

enddo
endif

if(num_bc /= 0) then

do j=1,num_bc
c1=>c(j)
if(c1%blk==i) then
if(c1%face==3) then
if(c1%side==1) then

if(c1%typ_num==2004)then
ndata=2

if(isurfg==0)then

kk=1

elseif(isurfg<0)then

if(c1%nam==surfg)then
kk=-1
else
kk=1
endif

else

if(c1%nam==surfg)then
kk=1
else
kk=-1
endif

endif

write(601,'(8I10)') i,kk*k,c1%typ_num,c1%is,c1%ie,c1%js,c1%je,ndata
write(601,'(2A10)') 'tw/tinf','cq'
write(601,'(2F10.4)') twratio,cq

elseif(c1%typ_num==2002)then
ndata=1
write(601,'(8I10)') i,kk*k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
write(601,'(A10)') 'p/pinf'
write(601,'(F10.4)') pratio

else
ndata=0
write(601,'(8I10)') i,k,c1%typ_num,c1%is,c1%ie,c1%js,c1%je,ndata
endif
k=k+1
endif
endif
endif

enddo

endif

enddo

write(601,'(8A10)') 'KDIM: GRID','SEGMENT','BCTYPE','ISTA','IEND','JSTA','JEND','NDATA'

do i=1,num_blk
k=1
if(num_if /= 0) then
do j=1,num_if
f1=>f(j)
if(f1%blk1==i) then
if(f1%face==3) then
if(f1%side==2) then
ndata=0
write(601,'(8I10)') i,k,ndata,f1%is,f1%ie,f1%js,f1%je,ndata
k=k+1
endif
endif
endif

enddo
endif

if(num_bc /= 0) then

do j=1,num_bc
c1=>c(j)
if(c1%blk==i) then
if(c1%face==3) then
if(c1%side==2) then

if(c1%typ_num==2004)then
ndata=2

if(isurfg==0)then

kk=1

elseif(isurfg<0)then

if(c1%nam==surfg)then
kk=-1
else
kk=1
endif

else

if(c1%nam==surfg)then
kk=1
else
kk=-1
endif

endif

write(601,'(8I10)') i,kk*k,c1%typ_num,c1%is,c1%ie,c1%js,c1%je,ndata
write(601,'(2A10)') 'tw/tinf','cq'
write(601,'(2F10.4)') twratio,cq

elseif(c1%typ_num==2002)then
ndata=1
write(601,'(8I10)') i,kk*k,c1%typ_num,c1%js,c1%je,c1%ks,c1%ke,ndata
write(601,'(A10)') 'p/pinf'
write(601,'(F10.4)') pratio

else
ndata=0
write(601,'(8I10)') i,k,c1%typ_num,c1%is,c1%ie,c1%js,c1%je,ndata
endif
k=k+1
endif
endif
endif

enddo

endif

enddo

write(601,'(5A10)') 'MSEQ','MGFLAG','ICONSF','MTT','NGAM'
if(mglevg==1) then

mgflag=0
else

mgflag=1
endif
mseq=mglevg

iconsf=1
mtt=0
ngam=2
write(601,'(5I10)') mseq,mgflag,iconsf,mtt,ngam

write(601,'(8A10)') 'ISSC','EPSSSC(1)','EPSSSC(2)','EPSSSC(3)','ISSR','EPSSSR(1)','EPSSSR(2)','EPSSSR(3)'

issc=0
issr=0
epsssc=0.3333d0
epsssr=0.3333d0

write(601,'(I10,3F10.4,I10,3F10.4)') issc,epsssc,epsssc,epsssc,issr,epsssr,epsssr,epsssr

write(601,'(4A10)') 'NCYC','MGLEVG','NEMGL','NITFO'

nemgl=0
nitfo=0
do i=1,mglevg
write(601,'(4I10)') ncyc(i),i,nemgl,nitfo
enddo

write(601,'(8A10)') 'MIT1','MIT2','MIT3','MIT4','MIT5','MIT6','MIT7','MIT8'

mit=1
do i=1,mglevg
write(601,'(8I10)') mit,mit,mit,mit,mit,mit,mit,mit
enddo

write(601,'(A21)') '1-1 BLOCKING DATA:'
write(601,'(A10)') 'NBLI'
write(601,'(I10)') num_ifp

write(601,'(10A8)') 'NUMBER','GRID :','ISTA','JSTA','KSTA','IEND','JEND','KEND','ISVA1','ISVA2'

!!!
k=0
do i=1,num_if
f1=>f(i)
do j=1,num_if
f2=>f(j)

if(i /= j) then
if(f1%pair<0) then
if(f2%pair<0) then
if(f1%key==f2%key)then
f1%pair=1
f2%pair=1
k=k+1
p1=>p(k)
p1%blk1=f1%blk1
p1%blk2=f2%blk1

p1%ii1=f1%face
p1%ii2=f2%face

p1%ss1=f1%side
p1%ss2=f2%side

p1%is1=f1%is
p1%is2=f2%is
p1%ie1=f1%ie
p1%ie2=f2%ie

p1%js1=f1%js
p1%js2=f2%js
p1%je1=f1%je
p1%je2=f2%je

p1%ks1=f1%ks
p1%ks2=f2%ks
p1%ke1=f1%ke
p1%ke2=f2%ke

endif
endif
endif
endif

!write(*,*) k

enddo
enddo

do j=1,3
do i=1,3
d(i,j)=0
enddo
enddo
d(1,1)=1
d(2,2)=1
d(3,3)=1

do i=1,num_ifp
p1=>p(i)
p1%jj1=1*d(p1%ii1,3)+2*d(p1%ii1,1)+3*d(p1%ii1,2)
p1%kk1=1*d(p1%ii1,2)+2*d(p1%ii1,3)+3*d(p1%ii1,1)
p1%jj2=1*d(p1%ii2,3)+2*d(p1%ii2,1)+3*d(p1%ii2,2)
p1%kk2=1*d(p1%ii2,2)+2*d(p1%ii2,3)+3*d(p1%ii2,1)
enddo

!pause

do i=1,num_ifp
p1=>p(i)
b1=>b(p1%blk1)
b2=>b(p1%blk2)

p1%i1=d(p1%ii1,1)*(d(p1%ss1,1)*p1%is1+d(p1%ss1,2)*p1%ie1)+d(p1%jj1,1)*p1%is1+d(p1%kk1,1)*p1%is1
p1%j1=d(p1%ii1,2)*(d(p1%ss1,1)*p1%js1+d(p1%ss1,2)*p1%je1)+d(p1%jj1,2)*p1%js1+d(p1%kk1,2)*p1%js1
p1%k1=d(p1%ii1,3)*(d(p1%ss1,1)*p1%ks1+d(p1%ss1,2)*p1%ke1)+d(p1%jj1,3)*p1%ks1+d(p1%kk1,3)*p1%ks1

p1%x1(1)=b1%x(p1%i1,p1%j1,p1%k1)
p1%y1(1)=b1%y(p1%i1,p1%j1,p1%k1)
p1%z1(1)=b1%z(p1%i1,p1%j1,p1%k1)

p1%i1=d(p1%ii1,1)*(d(p1%ss1,1)*p1%is1+d(p1%ss1,2)*p1%ie1)+d(p1%jj1,1)*p1%ie1+d(p1%kk1,1)*p1%is1
p1%j1=d(p1%ii1,2)*(d(p1%ss1,1)*p1%js1+d(p1%ss1,2)*p1%je1)+d(p1%jj1,2)*p1%je1+d(p1%kk1,2)*p1%js1
p1%k1=d(p1%ii1,3)*(d(p1%ss1,1)*p1%ks1+d(p1%ss1,2)*p1%ke1)+d(p1%jj1,3)*p1%ke1+d(p1%kk1,3)*p1%ks1

p1%x1(2)=b1%x(p1%i1,p1%j1,p1%k1)
p1%y1(2)=b1%y(p1%i1,p1%j1,p1%k1)
p1%z1(2)=b1%z(p1%i1,p1%j1,p1%k1)

p1%i1=d(p1%ii1,1)*(d(p1%ss1,1)*p1%is1+d(p1%ss1,2)*p1%ie1)+d(p1%jj1,1)*p1%is1+d(p1%kk1,1)*p1%ie1
p1%j1=d(p1%ii1,2)*(d(p1%ss1,1)*p1%js1+d(p1%ss1,2)*p1%je1)+d(p1%jj1,2)*p1%js1+d(p1%kk1,2)*p1%je1
p1%k1=d(p1%ii1,3)*(d(p1%ss1,1)*p1%ks1+d(p1%ss1,2)*p1%ke1)+d(p1%jj1,3)*p1%ks1+d(p1%kk1,3)*p1%ke1

p1%x1(3)=b1%x(p1%i1,p1%j1,p1%k1)
p1%y1(3)=b1%y(p1%i1,p1%j1,p1%k1)
p1%z1(3)=b1%z(p1%i1,p1%j1,p1%k1)

p1%i1=d(p1%ii1,1)*(d(p1%ss1,1)*p1%is1+d(p1%ss1,2)*p1%ie1)+d(p1%jj1,1)*p1%ie1+d(p1%kk1,1)*p1%ie1
p1%j1=d(p1%ii1,2)*(d(p1%ss1,1)*p1%js1+d(p1%ss1,2)*p1%je1)+d(p1%jj1,2)*p1%je1+d(p1%kk1,2)*p1%je1
p1%k1=d(p1%ii1,3)*(d(p1%ss1,1)*p1%ks1+d(p1%ss1,2)*p1%ke1)+d(p1%jj1,3)*p1%ke1+d(p1%kk1,3)*p1%ke1

p1%x1(4)=b1%x(p1%i1,p1%j1,p1%k1)
p1%y1(4)=b1%y(p1%i1,p1%j1,p1%k1)
p1%z1(4)=b1%z(p1%i1,p1%j1,p1%k1)

!pause

p1%i2=d(p1%ii2,1)*(d(p1%ss2,1)*p1%is2+d(p1%ss2,2)*p1%ie2)+d(p1%jj2,1)*p1%is2+d(p1%kk2,1)*p1%is2
p1%j2=d(p1%ii2,2)*(d(p1%ss2,1)*p1%js2+d(p1%ss2,2)*p1%je2)+d(p1%jj2,2)*p1%js2+d(p1%kk2,2)*p1%js2
p1%k2=d(p1%ii2,3)*(d(p1%ss2,1)*p1%ks2+d(p1%ss2,2)*p1%ke2)+d(p1%jj2,3)*p1%ks2+d(p1%kk2,3)*p1%ks2

p1%x2(1)=b2%x(p1%i2,p1%j2,p1%k2)
p1%y2(1)=b2%y(p1%i2,p1%j2,p1%k2)
p1%z2(1)=b2%z(p1%i2,p1%j2,p1%k2)

p1%i2=d(p1%ii2,1)*(d(p1%ss2,1)*p1%is2+d(p1%ss2,2)*p1%ie2)+d(p1%jj2,1)*p1%ie2+d(p1%kk2,1)*p1%is2
p1%j2=d(p1%ii2,2)*(d(p1%ss2,1)*p1%js2+d(p1%ss2,2)*p1%je2)+d(p1%jj2,2)*p1%je2+d(p1%kk2,2)*p1%js2
p1%k2=d(p1%ii2,3)*(d(p1%ss2,1)*p1%ks2+d(p1%ss2,2)*p1%ke2)+d(p1%jj2,3)*p1%ke2+d(p1%kk2,3)*p1%ks2

p1%x2(2)=b2%x(p1%i2,p1%j2,p1%k2)
p1%y2(2)=b2%y(p1%i2,p1%j2,p1%k2)
p1%z2(2)=b2%z(p1%i2,p1%j2,p1%k2)

p1%i2=d(p1%ii2,1)*(d(p1%ss2,1)*p1%is2+d(p1%ss2,2)*p1%ie2)+d(p1%jj2,1)*p1%is2+d(p1%kk2,1)*p1%ie2
p1%j2=d(p1%ii2,2)*(d(p1%ss2,1)*p1%js2+d(p1%ss2,2)*p1%je2)+d(p1%jj2,2)*p1%js2+d(p1%kk2,2)*p1%je2
p1%k2=d(p1%ii2,3)*(d(p1%ss2,1)*p1%ks2+d(p1%ss2,2)*p1%ke2)+d(p1%jj2,3)*p1%ks2+d(p1%kk2,3)*p1%ke2

p1%x2(3)=b2%x(p1%i2,p1%j2,p1%k2)
p1%y2(3)=b2%y(p1%i2,p1%j2,p1%k2)
p1%z2(3)=b2%z(p1%i2,p1%j2,p1%k2)

p1%i2=d(p1%ii2,1)*(d(p1%ss2,1)*p1%is2+d(p1%ss2,2)*p1%ie2)+d(p1%jj2,1)*p1%ie2+d(p1%kk2,1)*p1%ie2
p1%j2=d(p1%ii2,2)*(d(p1%ss2,1)*p1%js2+d(p1%ss2,2)*p1%je2)+d(p1%jj2,2)*p1%je2+d(p1%kk2,2)*p1%je2
p1%k2=d(p1%ii2,3)*(d(p1%ss2,1)*p1%ks2+d(p1%ss2,2)*p1%ke2)+d(p1%jj2,3)*p1%ke2+d(p1%kk2,3)*p1%ke2

p1%x2(4)=b2%x(p1%i2,p1%j2,p1%k2)
p1%y2(4)=b2%y(p1%i2,p1%j2,p1%k2)
p1%z2(4)=b2%z(p1%i2,p1%j2,p1%k2)

!pause

dmin=1d10
do j=1,4
dist1=(p1%x2(1)-p1%x1(j))**2+(p1%y2(1)-p1%y1(j))**2+(p1%z2(1)-p1%z1(j))**2
if(dist1<dmin)then
dmin=dist1
ii=j
endif
enddo

!pause

if(ii==1)then

dist1=(p1%x2(2)-p1%x1(2))**2+(p1%y2(2)-p1%y1(2))**2+(p1%z2(2)-p1%z1(2))**2
dist2=(p1%x2(2)-p1%x1(3))**2+(p1%y2(2)-p1%y1(3))**2+(p1%z2(2)-p1%z1(3))**2
if(dist1<dist2)then
p1%link=1
p1%dir1=1
p1%dir2=1
else
p1%link=2
p1%dir1=1
p1%dir2=1
endif

elseif(ii==2)then

dist1=(p1%x2(2)-p1%x1(1))**2+(p1%y2(2)-p1%y1(1))**2+(p1%z2(2)-p1%z1(1))**2
dist2=(p1%x2(2)-p1%x1(4))**2+(p1%y2(2)-p1%y1(4))**2+(p1%z2(2)-p1%z1(4))**2
if(dist1<dist2)then
p1%link=1
p1%dir1=2
p1%dir2=1
else
p1%link=2
p1%dir1=1
p1%dir2=2
endif

elseif(ii==3)then

dist1=(p1%x2(2)-p1%x1(4))**2+(p1%y2(2)-p1%y1(4))**2+(p1%z2(2)-p1%z1(4))**2
dist2=(p1%x2(2)-p1%x1(1))**2+(p1%y2(2)-p1%y1(1))**2+(p1%z2(2)-p1%z1(1))**2
if(dist1<dist2)then
p1%link=1
p1%dir1=1
p1%dir2=2
else
p1%link=2
p1%dir1=2
p1%dir2=1
endif

elseif(ii==4)then

dist1=(p1%x2(2)-p1%x1(3))**2+(p1%y2(2)-p1%y1(3))**2+(p1%z2(2)-p1%z1(3))**2
dist2=(p1%x2(2)-p1%x1(2))**2+(p1%y2(2)-p1%y1(2))**2+(p1%z2(2)-p1%z1(2))**2
if(dist1<dist2)then
p1%link=1
p1%dir1=2
p1%dir2=2
else
p1%link=2
p1%dir1=2
p1%dir2=2
endif

endif
write(601,'(10I8)') i,p1%blk1,p1%is1,p1%js1,p1%ks1,p1%ie1,p1%je1,p1%ke1,p1%jj1,p1%kk1

p1%jj3=p1%jj2
p1%kk3=p1%kk2
p1%is3=p1%is2
p1%ie3=p1%ie2
p1%js3=p1%js2
p1%je3=p1%je2
p1%ks3=p1%ks2
p1%ke3=p1%ke2

if(p1%link /= 1) then
p1%jj3=p1%kk2
p1%kk3=p1%jj2
endif

if(p1%dir1 /= 1) then
if(p1%ii2 == 1) then
p1%js3=p1%je2
p1%je3=p1%js2
elseif(p1%ii2 == 2) then
p1%ks3=p1%ke2
p1%ke3=p1%ks2
elseif(p1%ii2 == 3) then
p1%is3=p1%ie2
p1%ie3=p1%is2
endif
endif

if(p1%dir2 /= 1) then
if(p1%ii2 == 1) then
p1%ks3=p1%ke2
p1%ke3=p1%ks2
elseif(p1%ii2 == 2) then
p1%is3=p1%ie2
p1%ie3=p1%is2
elseif(p1%ii2 == 3) then
p1%js3=p1%je2
p1%je3=p1%js2
endif
endif

enddo

write(601,'(10A8)') 'NUMBER','GRID :','ISTA','JSTA','KSTA','IEND','JEND','KEND','ISVA1','ISVA2'

do i=1,num_ifp
p1=>p(i)
write(601,'(10I8)') i,p1%blk2,p1%is3,p1%js3,p1%ks3,p1%ie3,p1%je3,p1%ke3,p1%jj3,p1%kk3
enddo

write(601,'(A21)') 'PATCH SURFACE DATA:'

write(601,'(A10)') 'NINTER'

if(ipatchflag==1)then
open(605,file='ronnie.inp')
write(*,*) 'Generate ronnie.inp'
write(605,'(A)') 'I/O FILES'
write(605,'(A)') 'cfl3d.grd'
write(605,'(A)') 'ronnie.out'
write(605,'(A)') 'patch.bin'
write(605,'(A)') 'IOFLAG ITRACE'
ioflag=1
itrace=-1
write(605,'(2I5)') ioflag,itrace
write(605,'(A)') 'CFDRC_FILES:'//trim(pfgd)//'&'//trim(ibc)
write(605,'(A10)') 'NGRID'
write(605,'(I10)') ngrid
write(605,'(5A10)') 'NCG','IEM','IDIM','JDIM','KDIM'
do i=1,num_blk
b1=>b(i)
write(605,'(5I10)') ncg,iem,b1%ni,b1%nj,b1%nk
enddo
write(605,'(A10)') 'NINTER'
write(605,'(I10)') ninter
write(605,'(8A10)') 'INT','IIFIT','LLMIT','IITMAX','MMCXIE','MMCETA','IIC0','IIORPH'

iifit=1
llmit=1
iitmax=100
mmcxie=0
mmceta=0
iic0=0
iiorph=0
do i=1,ninter
write(605,'(8I10)') i,iifit,llmit,iitmax,mmcxie,mmceta,iic0,iiorph
enddo

write(605,'(A)') 'INT    TO   XIE1  XIE2   ETA1   ETA2   NFB (one per int)'
write(605,'(A)')  'FROM   XIE1  XIE2   ETA1   ETA2       (repeat nfb times for each int)'

ninter=0
do j=1,num_bcg
if(ipatch(j)==1) then

do i=1,num_bc

c1=>c(i)
if(c1%bcg==j) then
ninter=ninter+1
!write(*,*) 'TO',i,c1%blk*100+c1%face*10+c1%side

if(c1%face==1) then
iis=c1%js
iie=c1%je
jjs=c1%ks
jje=c1%ke
elseif(c1%face==2) then
iis=c1%ks
iie=c1%ke
jjs=c1%is
jje=c1%ie
else
iis=c1%js
iie=c1%je
jjs=c1%is
jje=c1%ie
endif
!bcg_num(bcg_pair(j))
!bcg_pair(j)
write(605,'(7I10)') ninter,c1%blk*100+c1%face*10+c1%side,iis,iie,jjs,jje,bcg_num(bcg_pair(j))
do k=1,num_bc
c2=>c(k)

if(c2%bcg==bcg_pair(j)) then
!write(*,*) 'FROM',k,c2%blk*100+c2%face*10+c2%side
if(c2%face==1) then
iis=c2%js
iie=c2%je
jjs=c2%ks
jje=c2%ke
elseif(c2%face==2) then
iis=c2%ks
iie=c2%ke
jjs=c2%is
jje=c2%ie
else
iis=c2%js
iie=c2%je
jjs=c2%is
jje=c2%ie
endif
write(605,'(I20,4I10)') c2%blk*100+c2%face*10+c2%side,iis,iie,jjs,jje
endif
enddo

endif
enddo

endif
enddo

ninter=-1
else
ninter=0
endif

write(601,'(I10)') ninter

write(601,'(A16)') 'PLOT3D OUTPUT:'

write(601,'(11A7)') 'GRID','IPTYPE','ISTART','IEND','IINC','JSTART','JEND','JINC','KSTART','KEND','KINC'

iptype=0

do i=1,num_blk
write(601,'(11I7)') i,iptype,iptype,iptype,iptype,iptype,iptype,iptype,iptype,iptype,iptype
enddo

write(601,'(A10)') 'IMOVIE'

imovie=0
write(601,'(I10)') imovie

write(601,'(A16)') 'PRINT OUT:'

write(601,'(11A7)') 'GRID','IPTYPE','ISTART','IEND','IINC','JSTART','JEND','JINC','KSTART','KEND','KINC'

write(601,'(A21)') 'CONTROL SURFACE:'

write(601,'(A10)') 'NCS'

!ncs=0

write(601,'(I10)') ncs
write(601,'(9A7)') 'GRID','ISTART','IEND','JSTART','JEND','KSTART','KEND','IWALL','INORM'

if(ncs/=0)then

do i=1,num_bc
c1=>c(i)
!b1=>b(c1%blk)

!!! control surface
if(c1%nam==csname) then

write(601,'(9I8)') c1%blk,c1%is,c1%ie,c1%js,c1%je,c1%ks,c1%ke,0,1

endif

enddo

endif

endprogram cfdrc_to_cfl3d


program main
	include 'mpif.h'
	parameter(msize=500)
	real(kind=8), allocatable :: a(:,:),b(:,:),c(:,:)
	real(kind=8) :: dinit, dflop, dbyte
	integer nWatch
	allocate (a(msize,msize), b(msize,msize), c(msize,msize), stat=istat)
	if (istat.ne.0) then
		write(*,*) "*** Allocate() failed."
		stop
	endif

	call MPI_Init(ierr )
	call MPI_Comm_rank( MPI_COMM_WORLD, myid, ierr )
	call MPI_Comm_size( MPI_COMM_WORLD, ncpus, ierr )
	write(6,'(a,i3,a)') "fortran <main> started process:", myid

	n=msize
	nWatch=4
	call f_pm_initialize (nWatch)

	call f_pm_start ("Initial-section")
	call subinit (msize,n,a,b,c)
	call f_pm_stop ("Initial-section")
!cx
    Max_power_knob=5
    do i=0,Max_power_knob-1
    call f_pm_getpowerknob (i,kval)
    if(myid.eq.0) write(*,'(a,2i6)') "preset value of power knob:", i, kval
    enddo
    do i=0,Max_power_knob-1
    kval=1;
    if(i.eq.0) kval=2200;
    if(i.eq.1) kval=5;
    call f_pm_setpowerknob (i,kval)
    if(myid.eq.0) write(*,'(a,2i6)') "preset value of power knob:", i, kval
    enddo
!cx
	call f_pm_start ("Kernel-Fast")
	call submtxm (msize,n,dflop,a,b,c)
	call f_pm_stop ("Kernel-Fast")
	call f_pm_report ("")

	call MPI_Finalize( ierr )
end program

subroutine subinit (msize,n,a,b,c)
	real(kind=8) :: a(msize,msize), b(msize,msize), c(msize,msize)
	write(6,'(a,i3,a)') " starting  <subinit>"
	!$omp parallel
	!$omp do
	do j=1, n
	do i=1, n
	a(i,j)=sin(real(i)/real(n))
	b(i,j)=cos(real(i)/real(n))
	end do
	end do
	!$omp end do
	!$omp end parallel
	write(6,'(a,i3,a)') " returns   <subinit>"
	return
end

subroutine submtxm (msize,n,dflop,a,b,c)
	real(kind=8) :: a(msize,msize), b(msize,msize), c(msize,msize)
	real(kind=8) :: dflop
	real(kind=8) :: x
	write(6,'(a,i3,a)') " starting  <submtxm>"
	dflop=2.0*dble(n)**3
	!$omp parallel
	!$omp do private(x)
	do j=1, n
	do i=1, n
	x=0
	do k=1, n
	x=x+a(i,k)*b(k,j)
	end do
	c(i,j) = x
	end do
	end do
	!$omp end do
	!$omp end parallel
	write(6,'(a,i3,a)') " returns   <submtxm>"
	return
end


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine make_set

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     This function allocates space for a set of cells and fills the set
c     such that communication between cells on different nodes is only
c     nearest neighbor
c---------------------------------------------------------------------

      include 'header.h'
      include 'mpinpb.h'

      integer p, i, j, c, dir, size, excess, ierr,ierrcode

c---------------------------------------------------------------------
c     compute square root; add small number to allow for roundoff
c     (note: this is computed in setup_mpi.f also, but prefer to do
c     it twice because of some include file problems).
c---------------------------------------------------------------------
      ncells = dint(dsqrt(dble(no_nodes) + 0.00001d0))

c---------------------------------------------------------------------
c     this makes coding easier
c---------------------------------------------------------------------
      p = ncells

c---------------------------------------------------------------------
c     determine the location of the cell at the bottom of the 3D
c     array of cells
c---------------------------------------------------------------------
      cell_coord(1,1) = mod(node,p)
      cell_coord(2,1) = node/p
      cell_coord(3,1) = 0

c---------------------------------------------------------------------
c     set the cell_coords for cells in the rest of the z-layers;
c     this comes down to a simple linear numbering in the z-direct-
c     ion, and to the doubly-cyclic numbering in the other dirs
c---------------------------------------------------------------------
      do c=2, p
         cell_coord(1,c) = mod(cell_coord(1,c-1)+1,p)
         cell_coord(2,c) = mod(cell_coord(2,c-1)-1+p,p)
         cell_coord(3,c) = c-1
      end do

c---------------------------------------------------------------------
c     offset all the coordinates by 1 to adjust for Fortran arrays
c---------------------------------------------------------------------
      do dir = 1, 3
         do c = 1, p
            cell_coord(dir,c) = cell_coord(dir,c) + 1
         end do
      end do

c---------------------------------------------------------------------
c     slice(dir,n) contains the sequence number of the cell that is in
c     coordinate plane n in the dir direction
c---------------------------------------------------------------------
      do dir = 1, 3
         do c = 1, p
            slice(dir,cell_coord(dir,c)) = c
         end do
      end do

c---------------------------------------------------------------------
c     fill the predecessor and successor entries, using the indices
c     of the bottom cells (they are the same at each level of k
c     anyway) acting as if full periodicity pertains; note that p is
c     added to those arguments to the mod functions that might
c     otherwise return wrong values when using the modulo function
c---------------------------------------------------------------------
      i = cell_coord(1,1)-1
      j = cell_coord(2,1)-1

      predecessor(1) = mod(i-1+p,p) + p*j
      predecessor(2) = i + p*mod(j-1+p,p)
      predecessor(3) = mod(i+1,p) + p*mod(j-1+p,p)
      successor(1)   = mod(i+1,p) + p*j
      successor(2)   = i + p*mod(j+1,p)
      successor(3)   = mod(i-1+p,p) + p*mod(j+1,p)

c---------------------------------------------------------------------
c     now compute the sizes of the cells
c---------------------------------------------------------------------
      do dir= 1, 3
c---------------------------------------------------------------------
c     set cell_coord range for each direction
c---------------------------------------------------------------------
         size   = grid_points(dir)/p
         excess = mod(grid_points(dir),p)
         do c=1, ncells
            if (cell_coord(dir,c) .le. excess) then
               cell_size(dir,c) = size+1
               cell_low(dir,c) = (cell_coord(dir,c)-1)*(size+1)
               cell_high(dir,c) = cell_low(dir,c)+size
            else
               cell_size(dir,c) = size
               cell_low(dir,c)  = excess*(size+1)+
     >              (cell_coord(dir,c)-excess-1)*size
               cell_high(dir,c) = cell_low(dir,c)+size-1
            endif
            if (cell_size(dir, c) .le. 2) then
               write(*,50)
 50            format(' Error: Cell size too small. Min size is 3')
               ierrcode = 1
               call MPI_Abort(mpi_comm_world,ierrcode,ierr)
               stop
            endif
         end do
      end do

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

max_stars_repo_name: jonarbo/KUBE

max_issues_repo_name: jonarbo/KUBE

max_forks_repo_name: jonarbo/KUBE