text:
module ExtFEcomparer_parameters

! Include basic Feat-2 modules
  use fsystem
  use storage
  use genoutput
  use paramlist

  use element

  use ExtFEcomparer_typedefs

  implicit none

  private

  public :: ExtFEcomparer_get_parameters
  public :: ExtFEcomparer_init_parameters
  public :: ExtFEcomparer_getLogFiles
  public :: ExtFEcomparer_parseCmdlArguments
  public :: ExtFEcomparer_getCmdlMasterDat

contains

subroutine ExtFEcomparer_get_parameters(rparamList)

!<description>
  ! Reads in all DAT files into the parameter list rparamlist
!</description>

!<inputoutput>
  ! The parameter list where the values of the DAT files should be stored.
  ! The structure must have been initialised, the parameters are just added
  ! to the list.
  type(t_parlist), intent(inout) :: rparamList
!</inputoutput>

!</subroutine>

    logical :: bexists
    character(LEN=ExtFE_STRLEN) :: smaster

    ! Check if a command line parameter specifies the master.dat file.
    ! It always is the last argument
    smaster = ''
    call ExtFEcomparer_getCmdlMasterDat(smaster)

    ! Let's see
    if(smaster .eq. '') then
        smaster = './data/master.dat'
    end if

    ! Read the file "master.dat".
    ! If that does not exist, try to manually read files with parameters from a
    ! couple of files.
    inquire(file=smaster, exist=bexists)

    if (bexists) then
      ! Read the master file. That either one contains all parameters or
      ! contains references to subfiles with data.
      call parlst_readfromfile (rparamList, smaster)
    else
      ! Each "readfromfile" command adds the parameter of the specified file
      ! to the parameter list.
      call parlst_readfromfile (rparamList, "./data/generalsettings.dat")
      call parlst_readfromfile (rparamList, "./data/function_specific_settings.dat")
      call parlst_readfromfile (rparamList, "./data/postprocessing.dat")
      call parlst_readfromfile (rparamList, "./data/logfile_settings.dat")

    end if

    ! Now we parse the command-line parameters in the parameter list
    call ExtFEcomparer_parseCmdlArguments(rparamList)

    ! Thats it

end subroutine

subroutine ExtFEcomparer_init_parameters(rproblem, section)
!<description>
  ! Initialises the structure rproblem with data from the initialisation
  ! files.
!</description>

!<input>
  ! A string to define the section
  ! Must be "ExtFE-FIRST" or "ExtFE-SECOND" to specify
  ! which function we init.
  character(LEN=*), intent(in) :: section
!<input>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(inout) :: rproblem
!</inputoutput>

    integer :: i,ilvmax, tmpvalue, ielemPairId
    integer(I32) :: ielemtype

    ! Variable for a filename:
    character(LEN=ExtFE_STRLEN) :: sString
    character(LEN=ExtFE_STRLEN) :: sPRMFile, sTRIFile
    character(LEN=ExtFE_STRLEN) :: sVectorFile
    character(LEN=ExtFE_STRLEN) :: smessage

    ! ilvmax receives the level where our solution lives on

    call parlst_getvalue_int (rproblem%rparamlist,section,&
                              "NLMAX",ilvmax)
    rproblem%NLMAX=ilvmax

    ! find out the dimension of the problem
    call parlst_getvalue_int(rproblem%rparamlist, "ExtFE-DOMAININFO", &
                            "dim", tmpvalue)
    rproblem%iDimension = tmpvalue

    ! Get the path of the .prm and the .tri file from the parameter list.
    ! .prm is not there in the 1D-Case
    if (rproblem%iDimension .ne. ExtFE_NDIM1) then
        call parlst_getvalue_string (rproblem%rparamList,section,&
                                    "sParametrisation",sPRMFile,bdequote=.true.)
        rproblem%sPRMFile = sPRMFile
    end if

    ! Get the path of the mesh
    call parlst_getvalue_string (rproblem%rparamList,section,&
                                 "sMesh",sTRIFile,bdequote=.true.)
    rproblem%sTRIFile = sTRIFile

    ! Get the path of the vector file
    call parlst_getvalue_string (rproblem%rparamList,section,&
                                 "sVector",sVectorFile,bdequote=.true.)
    rproblem%sVectorFile = sVectorFile

    ! Get the format identifier telling us what we have
    ! in the vector file
    call parlst_getvalue_int(rproblem%rparamlist,section, &
                            "sFileFormat",tmpvalue)
    rproblem%vectorFileFormat = tmpvalue
    ! Written out by write_blockVectorHR or something else?
    call parlst_getvalue_int(rproblem%rparamlist,section, &
                            "sVectorType", tmpvalue)
    rproblem%vectorType = tmpvalue
    ! How many variables are in the vector?
    call parlst_getvalue_int(rproblem%rparamlist,section, &
                            "iNVAR", tmpvalue)
    rproblem%NVAR = tmpvalue

    ! get the element type. Standard is -1 so that it will crash when
    ! we try to recreate the discretisation and nothing is specified.
    ! however, there is one more option to specify the element - a string!
    ! Get the element setting: One element or an ID of a pair?

    call parlst_getvalue_int(rproblem%rparamlist,section, &
                            "iElementSetting",tmpvalue)
    rproblem%elementSetting = tmpvalue

    select case(rproblem%elementSetting)
        case(ExtFE_ElementPair)
            call output_line("Searching for a numerical identifier of the element pair")
            call parlst_getvalue_int(rproblem%rparamlist,section,"ielementType",&
                        ielemPairId)
            call output_line("Found a numerical identifier of the element pair")
            rproblem%iElemPairID = ielemPairId

        case(ExtFE_OneElement)
            call output_line("Searching for the name of an element")
            call parlst_getvalue_string(rproblem%rparamlist,section, &
                    "ElementName",sString,bdequote=.true.)
            ielemtype = elem_igetID(sString)
            call output_line("Found a name of an element")
            rproblem%ielemType = ielemtype

        case(ExtFE_ElementList)
            call output_line("Search for an element list")

            ! How long is the list we find?
            tmpvalue = parlst_querysubstrings(rproblem%rparamlist,&
                                    section,"ElementList")
            if(tmpvalue .ne. rproblem%NVAR) then
                write(smessage,*) 'Input error: List length &
                &does not match the number of variables for function ', section
                call output_line(smessage , &
                OU_CLASS_ERROR,OU_MODE_STD,"ExtFEcomparer_parameters")
                call sys_halt()
            end if

            ! Allocate the array for the element list with
            ! one element per variable
            allocate(rproblem%iElemList(rproblem%NVAR))
            do i=1,rproblem%NVAR
                call parlst_getvalue_string(rproblem%rparamlist, &
                    section, "ElementList", &
                            sString,sdefault="",isubstring=i,bdequote=.TRUE.)
                rproblem%iElemList(i) = elem_igetID(sString)
            end do
            call output_line("Found an element list")

        case default
            write(smessage,*) "Input error: choice not allowed for &
            &the selection of iElementSetting for function ", section
            call output_line(smessage , &
            OU_CLASS_ERROR,OU_MODE_STD,"ExtFEcomparer_parameters")
            call sys_halt()
    end select

end subroutine

!<subroutine>
subroutine ExtFEcomparer_getLogFiles(slogfile,serrorfile,sbenchlogfile)

    !<description>
    ! Temporarily reads the output DAT file to get the names of the output
    ! files.
    !</description>

    !<output>
    ! Name of the message log file.
    character(LEN=*), intent(out) :: slogfile

    ! Name of the error log file.
    character(LEN=*), intent(out) :: serrorfile

    ! Name of the benchmark log file.
    character(LEN=*), intent(out) :: sbenchlogfile
    !</output>

    type(t_parlist) :: rparlist
    character(LEN=ExtFE_STRLEN) :: smaster
    logical :: bexists

    ! Init parameter list that accepts parameters for output files
    call parlst_init (rparlist)

    smaster = ''
    ! Check if a command line parameter specifies the master.dat file.
    call ExtFEcomparer_getCmdlMasterDat(smaster)
    if(smaster .eq. '') then
        smaster = './data/master.dat'
    else
        call output_line("Using the following master.dat: "//smaster)
    end if

    ! Read parameters that configure the output
    inquire(file=smaster, exist=bexists)

    if (bexists) then
        ! Read the master file. That either one contains all parameters or
        ! contains references to subfiles with data.
        call parlst_readfromfile (rparlist, smaster)
    else
        call parlst_readfromfile (rparlist, "./data/logfile_settings.dat")
    end if

    ! Now the real initialisation of the output including log file stuff!
    call parlst_getvalue_string (rparlist,"ExtFE-LOGFILESETTINGS",&
        "smsgLog",slogfile,"",bdequote=.true.)

    call parlst_getvalue_string (rparlist,"ExtFE-LOGFILESETTINGS",&
        "serrorLog",serrorfile,"",bdequote=.true.)

    call parlst_getvalue_string (rparlist,"ExtFE-LOGFILESETTINGS",&
        "sbenchLog",sbenchlogfile,"",bdequote=.true.)

    ! That temporary parameter list is not needed anymore.
    call parlst_done (rparlist)

end subroutine

!************************************************************************

!<subroutine>

  subroutine ExtFEcomparer_parseCmdlArguments(rparlist)

!<description>
    ! This subroutine parses the commandline arguments and modifies the
    ! parameter values in the global parameter list.
!</description>

!<inputoutput>
    ! parameter list
    type(t_parlist), intent(inout) :: rparlist
!</inputoutput>
!</subroutine>

    ! local variables
    type(t_parlstSection), pointer :: p_rsection
    character(LEN=PARLST_MLSECTION) :: ssection
    character(LEN=PARLST_MLNAME) :: sparameter,sentry
    character(LEN=PARLST_MLDATA) :: svalue
    character(LEN=SYS_STRLEN) :: soption
    character(LEN=PARLST_MLDATA), dimension(:), pointer :: Svalues
    integer :: iarg,narg,iformat,itoken1,itoken2,isubstring,nsubstrings,i

    iarg = 1; narg = sys_ncommandLineArgs()

    cmdarg: do while(iarg .le. narg)
      ! Retrieve next command line argument
      call sys_getcommandLineArg(iarg,soption,svalue,iformat)

      select case(iformat)
      case (0)
        ! Options without parameter values

      case (1,2)
        ! Options with given parameter values

        ! What option are we?
        if (soption(1:1) .eq. 'D') then

          ! Tokenize string: soption=D<section>.variable:<entry>
          itoken1 = scan(soption,'.')
          if (itoken1 .ne. 0) then
            ssection = trim(adjustl(soption(2:itoken1-1)))
          else
            ssection = ''
          end if

          itoken2 = scan(soption,':')
          if (itoken2 .ne. 0) then
            sparameter = trim(adjustl(soption(max(2,itoken1+1):itoken2-1)))
            sentry     = trim(adjustl(soption(max(2,itoken2+1):)))
            read(sentry, fmt='(I10)') isubstring
          else
            sparameter = trim(adjustl(soption(max(2,itoken1+1):)))
            sentry     = ''
            isubstring = 0
          end if

          ! Query/add section in parameter list
          call parlst_querysection(rparlist, trim(adjustl(ssection)), p_rsection)
          if (.not.associated(p_rsection)) then
            call parlst_addsection (rparlist, trim(adjustl(ssection)))
            call parlst_querysection(rparlist, trim(adjustl(ssection)), p_rsection)
          end if

          ! We need to consider several cases
          if (parlst_queryvalue(p_rsection, trim(adjustl(sparameter))) .ne. 0) then

            ! Parameter exists already in section
            if (isubstring .eq. 0) then
              ! Overwrite existing value
              call parlst_addvalue(p_rsection, trim(adjustl(sparameter)),&
                  trim(adjustl(svalue)))
            else
              ! Get number of existing substrings
              nsubstrings = parlst_querysubstrings(p_rsection, trim(adjustl(sparameter)))
              if (isubstring .lt. nsubstrings) then
                ! Overwrite existing substring
                call parlst_setvalue(p_rsection, trim(adjustl(sparameter)),&
                    trim(adjustl(svalue)), isubstring=isubstring)
              else
                ! Make a backup of existing substrings
                allocate(Svalues(0:nsubstrings))
                do i=0,nsubstrings
                  call parlst_getvalue_string(p_rsection, trim(adjustl(sparameter)),&
                      Svalues(i), isubstring=i)
                end do
                ! Increase the number of substrings
                call parlst_addvalue(p_rsection, trim(adjustl(sparameter)),&
                    trim(adjustl(svalue)), isubstring)
                ! Copy existing substrings
                do i=0,nsubstrings
                  call parlst_setvalue(p_rsection, trim(adjustl(sparameter)),&
                      Svalues(i), isubstring=i)
                end do
                ! Add new substring
                call parlst_setvalue(p_rsection, trim(adjustl(sparameter)),&
                    trim(adjustl(svalue)), isubstring=isubstring)
                deallocate(Svalues)
              end if
            end if
          else
            ! Add new value to parameter list
            if (isubstring .eq. 0) then
              call parlst_addvalue(p_rsection, trim(adjustl(sparameter)),&
                  trim(adjustl(svalue)))
            else
              call parlst_addvalue(p_rsection, trim(adjustl(sparameter)),&
                  trim(adjustl(svalue)), isubstring)
              call parlst_setvalue(p_rsection, trim(adjustl(sparameter)),&
                  trim(adjustl(svalue)), isubstring=isubstring)
            end if
          end if
        elseif (soption(1:7) .eq. 'smaster') then
            ! This is a master.dat We search for it in another routine
            ! don't do anything
        else
          call output_line('Invalid option: '//trim(adjustl(soption))//'!',&
              OU_CLASS_WARNING,OU_MODE_STD,'ExtFEcomparer_parseCmdLineArg')
        end if
      end select

      ! Proceed with next command line argument
      iarg=iarg+1
    end do cmdarg

  end subroutine

  subroutine ExtFEcomparer_getCmdlMasterDat(masterDat)

!<description>
    ! This subroutine parses the commandline arguments and modifies the
    ! parameter values in the global parameter list.
!</description>

!<inputoutput>
    ! path of the master.dat
    character(LEN=ExtFE_STRLEN),intent(inout) :: masterDat
!</inputoutput>
!</subroutine>

    ! local variables
    character(LEN=SYS_STRLEN) :: soption
    character(LEN=PARLST_MLDATA) :: svalue
    integer :: iarg,narg,iformat
    integer :: ipathlen
    iarg = 1; narg = sys_ncommandLineArgs()

    cmdarg: do while(iarg .le. narg)
      ! Retrieve next command line argument
      call sys_getcommandLineArg(iarg,soption,svalue,iformat)

      select case(iformat)
      case (0)
        ! Options without parameter values

      case (1,2)
        ! Options with given parameter values

        ! What option are we?
        if (soption(1:7) .eq. 'smaster') then
            ! no need to tokenize the string. soption=smaster:somepath
            ! We just get everything after the smaster:
            call output_line("Found a master.dat on the command line")
            !ipathlen = len(adjustl(trim(soption(9:len(soption)))))
            masterDat = svalue
            masterDat = trim(masterDat)
            return
        end if
      end select

      ! Proceed with next command line argument
      iarg=iarg+1

    end do cmdarg

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE LOADB0

! Preliminary reading of the Bulk Data to count several data sizes so that arrays may be allocated prior to the final reading
! of the Bulk Data.

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06, IN1
      USE SCONTR, ONLY                :  BD_ENTRY_LEN, BLNK_SUB_NAM, FATAL_ERR, LCMASS, LDOFG, LELE,                               &
                                         LEDAT, LFORCE, LCONM2, LCORD, LGRAV, LGRID, LGUSERIN, LLOADC, LLOADR,                     &
                                         LMATL, LMPC, LMPCADDC, LMPCADDR, LPBAR, LPBEAM, LPBUSH, LPCOMP, LPCOMP_PLIES, LPDAT,      &
                                         LPELAS, LPLOAD, LPMASS, LPROD, LPSHEL, LPSHEAR, LPSOLID, LPUSER1, LPUSERIN, LRFORCE,      &
                                         LRIGEL, LSEQ, LSLOAD, LSPC, LSPC1, LSPCADDC, LSPCADDR, LSUSERIN, LTDAT,                   &
                                         MEDAT_CBAR, MEDAT_CBEAM, MEDAT_CBUSH,                                                     &
                                         MEDAT_CELAS1, MEDAT_CELAS2, MEDAT_CELAS3, MEDAT_CELAS4,                                   &
                                         MEDAT_CQUAD, MEDAT_CROD, MEDAT_CSHEAR, MEDAT_CTRIA, MEDAT_CUSER1, MEDAT0_CUSERIN, MMPC,   &
                                         MPDAT_PLOAD2, MPDAT_PLOAD4, MEDAT_PLOTEL, MRBE3, MRSPLINE, MTDAT_TEMPRB, MTDAT_TEMPP1,    &
                                         NPBARL, NSPOINT, PROG_NAME
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  LOADB0_BEGEND
      USE MODEL_STUF, ONLY            :  GRDSET3, GRDSET7, GRDSET8
      USE PARAMS, ONLY                :  GRIDSEQ

      USE LOADB0_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME   = 'LOADB0'
      CHARACTER( 7*BYTE), PARAMETER   :: END_CARD    = 'ENDDATA'

      CHARACTER(LEN=BD_ENTRY_LEN)     :: CARD              ! 16 col field card (either CARD1 if small field or CARD1 + CARD2 if
!                                                            a large field card. This is output from subr FFIELD

      CHARACTER(LEN=BD_ENTRY_LEN)     :: CARD1             ! BD card (a small field card or the 1st half of a large field card)
      CHARACTER(LEN=BD_ENTRY_LEN)     :: CARD2             ! 2nd half of a large field card
      CHARACTER( 9*BYTE)              :: DECK_NAME   = 'BULK DATA'
      CHARACTER( 1*BYTE)              :: LARGE_FLD_INP     ! If 'Y', card is in large field format

      INTEGER(LONG)                   :: COMMENT_COL       ! Col on CARD where a comment begins (if one exists)
      INTEGER(LONG)                   :: DELTA_LEDAT       ! Delta number of words to add to LEDAT for an element (for elements
!                                                            where the delta cannot be discerned from the BD card name, such as
!                                                            CHEXA which can have 8 or 20 nodes
      INTEGER(LONG)                   :: DELTA_SLOAD       ! Delta number of words to add to LSLOAD for SLOAD's
      INTEGER(LONG)                   :: DELTA_SPOINT      ! Delta number of words to add to LGRID for SPOINT's
      INTEGER(LONG)                   :: I                 ! DO loop index
      INTEGER(LONG)                   :: ICNT              ! Card count
      INTEGER(LONG)                   :: IERR              ! Error indicatro from subr FFIELD
      INTEGER(LONG)                   :: IOCHK             ! IOSTAT error number when reading a file
      INTEGER(LONG)                   :: ILOAD             ! Number of loads defined on 1 logical B.D. LOAD card
      INTEGER(LONG)                   :: IMPC              ! Number of grid/comp/coeff triplets defined on 1 logical B.D. MPC  card
      INTEGER(LONG)                   :: IMPCADD           ! Number of MPC set ID's defined on 1 B.D. MPCADD card
      INTEGER(LONG)                   :: IRBE3             ! Number of grid/comp/coeff triplets defined on 1 logical B.D. RBE3 card
      INTEGER(LONG)                   :: IRSPLINE          ! Number of fields that can have a grid or comp num on an RSPLINE entry
      INTEGER(LONG)                   :: ISPCADD           ! Number of SPC set ID's defined on 1 B.D. SPCADD card
      INTEGER(LONG)                   :: IPLIES            ! Number of composite layers on 1 B.D. PCOMP card
      INTEGER(LONG)                   :: NG_USERIN         ! Number of grids found on USERIN elems (not incl SPOINT's)
      INTEGER(LONG)                   :: NS_USERIN         ! Number of SPOINT's found on USERIN elems
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = LOADB0_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
      CARD(1:) = ' '

! Set MRBE3, MRSPLINE so that we will make sure that they will be at least 1

      MRBE3    = 1
      MRSPLINE = 1

! Process Bulk Data cards in a loop that runs until either an ENDDATA card is found or when an error or EOF/EOR occurs

      ICNT = 0
      DO

         ICNT = ICNT + 1
         READ(IN1,101,IOSTAT=IOCHK) CARD1
         CARD(1:)  = CARD1(1:)

! Quit if EOF/EOR occurs.

         IF (IOCHK < 0) THEN
            WRITE(ERR,1011) END_CARD
            WRITE(F06,1011) END_CARD
            FATAL_ERR = FATAL_ERR + 1
            CALL OUTA_HERE ( 'Y' )
         ENDIF

! Check if error occurs.

         IF (IOCHK > 0) THEN
            WRITE(ERR,1010) DECK_NAME
            WRITE(F06,1010) DECK_NAME
            WRITE(F06,'(A)') CARD1
            FATAL_ERR = FATAL_ERR + 1
            CYCLE
         ENDIF

! Remove any comments within the CARD1 by deleting everything from $ on (after col 1)

         COMMENT_COL = 1
         DO I=2,BD_ENTRY_LEN
            IF (CARD1(I:I) == '$') THEN
               COMMENT_COL = I
               EXIT
            ENDIF
         ENDDO

         IF (COMMENT_COL > 1) THEN
            CARD1(COMMENT_COL:) = ' '
         ENDIF

! Determine if the card is large or small format

         LARGE_FLD_INP = 'N'
         DO I=1,8
            IF (CARD1(I:I) == '*') THEN
               LARGE_FLD_INP = 'Y'
            ENDIF
         ENDDO

! FFIELD converts free-field card to fixed field and left justifies data in fields 2-9 and outputs a 10 field, 16 col/field CARD1

         IF ((CARD1(1:1) /= '$')  .AND. (CARD1(1:) /= ' ')) THEN

            IF (LARGE_FLD_INP == 'N') THEN

               CALL FFIELD ( CARD1, IERR )
               CARD(1:) = CARD1(1:)

            ELSE

               READ(IN1,101,IOSTAT=IOCHK) CARD2            ! Read 2nd physical entry for a large field parent B.D. entry

               IF (IOCHK < 0) THEN
                  WRITE(ERR,1011) END_CARD
                  WRITE(F06,1011) END_CARD
                  FATAL_ERR = FATAL_ERR + 1
                  CALL OUTA_HERE ( 'Y' )
               ENDIF

               IF (IOCHK > 0) THEN
                  WRITE(ERR,1010) DECK_NAME
                  WRITE(F06,1010) DECK_NAME
                  WRITE(F06,'(A)') CARD2
                  FATAL_ERR = FATAL_ERR + 1
                  CYCLE
               ENDIF

               COMMENT_COL = 1
               DO I=2,BD_ENTRY_LEN
                  IF (CARD2(I:I) == '$') THEN
                     COMMENT_COL = I
                     EXIT
                  ENDIF
               ENDDO

               IF (COMMENT_COL > 1) THEN
                  CARD2(COMMENT_COL:) = ' '
               ENDIF

!              IF (CARD2(1:8) /= CARD1(73:80)) THEN
!                 BACKSPACE(IN1)
!                 CARD2(1:) = ' '
!                 CARD2(1:8) = CARD1(73:80)
!              ENDIF

               IF      (CARD2( 1: 8) == CARD1(73:80)) THEN
                  CONTINUE     !ICONT = 1
               ELSE IF ((CARD2( 1: 8) == '*       ') .AND. (CARD1(73:80) == '        ')) THEN
                  CONTINUE     !ICONT = 1
               ELSE IF ((CARD2( 1: 8) == '        ') .AND. (CARD1(73:80) == '*       ')) THEN
                  CONTINUE     !ICONT = 1
               ELSE
                  BACKSPACE(IN1)
                  CARD2(1:) = ' '
                  CARD2(1:8) = CARD1(73:80)
               ENDIF

               CALL FFIELD2 ( CARD1, CARD2, CARD, IERR )

            ENDIF

            IF (IERR /= 0) THEN
               CYCLE
            ENDIF

         ENDIF

! No errors, so process Bulk Data card. No need to check for imbedded blanks found when FFIELD was run - this will be
! checked when LOADB reads the bulk data

          IF      (CARD(1:5) == 'BAROR'   )  THEN
            CALL BD_BAROR0 ( CARD )

         ELSE IF (CARD(1:6) == 'BEAMOR'  )  THEN
            CALL BD_BEAMOR0( CARD )

         ELSE IF ((CARD(1:4) == 'CBAR'    ) .OR. (CARD(1:5) == 'CBEAM'   ))  THEN
            LELE  = LELE + 1
            IF (CARD(1:4) == 'CBAR'    ) THEN
               LEDAT = LEDAT + MEDAT_CBAR
            ELSE
               LEDAT = LEDAT + MEDAT_CBEAM
            ENDIF
            CALL BD_CBAR0 ( CARD, LARGE_FLD_INP )

         ELSE IF (CARD(1:5) == 'CBUSH'   )  THEN
            LELE  = LELE + 1
            LEDAT = LEDAT + MEDAT_CBUSH
            CALL BD_CBUSH0 ( CARD, LARGE_FLD_INP )

         ELSE IF (CARD(1:6) == 'CELAS1'  )  THEN
            LELE  = LELE + 1
            LEDAT = LEDAT + MEDAT_CELAS1

         ELSE IF (CARD(1:6) == 'CELAS2'  )  THEN
            LELE   = LELE + 1
            LEDAT  = LEDAT + MEDAT_CELAS2
            LPELAS = LPELAS + 1                            ! CELAS2 has props on conn entry and we create a PELAS for them

         ELSE IF (CARD(1:6) == 'CELAS3'  )  THEN
            LELE  = LELE + 1
            LEDAT = LEDAT + MEDAT_CELAS3

         ELSE IF (CARD(1:6) == 'CELAS4'  )  THEN
            LELE   = LELE + 1
            LEDAT  = LEDAT + MEDAT_CELAS4
            LPELAS = LPELAS + 1                            ! CELAS4 has props on conn entry and we create a PELAS for them

         ELSE IF (CARD(1:5) == 'CHEXA'   ) THEN
            LELE  = LELE + 1
            CALL BD_CHEXA0 ( CARD, LARGE_FLD_INP, DELTA_LEDAT )
            LEDAT = LEDAT + DELTA_LEDAT

         ELSE IF (CARD(1:6) == 'CMASS1'  )  THEN
            LCMASS = LCMASS + 1

         ELSE IF (CARD(1:6) == 'CMASS2'  )  THEN
            LCMASS = LCMASS + 1
            LPMASS = LPMASS + 1

         ELSE IF (CARD(1:6) == 'CMASS3'  )  THEN
            LCMASS = LCMASS + 1

         ELSE IF (CARD(1:6) == 'CMASS4'  )  THEN
            LCMASS = LCMASS + 1
            LPMASS = LPMASS + 1

         ELSE IF (CARD(1:5) == 'CONM2'   )  THEN
            LCONM2 = LCONM2 + 1

         ELSE IF (CARD(1:6) == 'CONROD'  )  THEN
            LELE  = LELE  + 1
            LPROD = LPROD + 1
            LEDAT = LEDAT + MEDAT_CROD

         ELSE IF((CARD(1:6) == 'CORD1C'  ) .OR. (CARD(1:6) == 'CORD1R'  ) .OR. (CARD(1:6) == 'CORD1S'  )) THEN
            LCORD = LCORD + 1
            IF (CARD(41:48) /= '        ') LCORD = LCORD + 1

         ELSE IF((CARD(1:6) == 'CORD2C'  ) .OR. (CARD(1:6) == 'CORD2R'  ) .OR. (CARD(1:6) == 'CORD2S'  )) THEN
            LCORD = LCORD + 1

         ELSE IF (CARD(1:6) == 'CPENTA'  ) THEN
            LELE  = LELE + 1
            CALL BD_CPENTA0 ( CARD, LARGE_FLD_INP, DELTA_LEDAT )
            LEDAT = LEDAT + DELTA_LEDAT

         ELSE IF (CARD(1:6) == 'CQUAD4'  ) THEN
            LELE  = LELE + 1
            LEDAT = LEDAT + MEDAT_CQUAD
            CALL BD_CQUAD0 ( CARD, LARGE_FLD_INP )

         ELSE IF (CARD(1:4) == 'CROD'    )  THEN
            LELE  = LELE + 1
            LEDAT = LEDAT + MEDAT_CROD

         ELSE IF (CARD(1:6) == 'CSHEAR'  )  THEN
            LELE  = LELE + 1
            LEDAT = LEDAT + MEDAT_CSHEAR

         ELSE IF (CARD(1:6) == 'CTETRA'  ) THEN
            LELE  = LELE + 1
            CALL BD_CTETRA0 ( CARD, LARGE_FLD_INP, DELTA_LEDAT )
            LEDAT = LEDAT + DELTA_LEDAT

         ELSE IF (CARD(1:6) == 'CTRIA3'  ) THEN
            LELE  = LELE + 1
            LEDAT = LEDAT + MEDAT_CTRIA
            CALL BD_CTRIA0 ( CARD, LARGE_FLD_INP )

         ELSE IF (CARD(1:6) == 'CUSER1'  )  THEN
            LELE  = LELE + 1
            LEDAT = LEDAT + MEDAT_CUSER1

         ELSE IF (CARD(1:7) == 'CUSERIN' )  THEN
            LELE  = LELE + 1
            CALL BD_CUSERIN0 ( CARD, NG_USERIN, NS_USERIN )
            IF (NG_USERIN > LGUSERIN) THEN
               LGUSERIN = NG_USERIN
            ENDIF
            IF (NS_USERIN > LSUSERIN) THEN
               LSUSERIN = NS_USERIN
            ENDIF
                                                           ! LEDAT has "+ 1" term since last record is NUM_BDY_DOF not in MEDAT0
            LEDAT = LEDAT + MEDAT0_CUSERIN + 2*NG_USERIN + NS_USERIN +  1

         ELSE IF (CARD(1:5) == 'DEBUG'   )  THEN
            CALL BD_DEBUG0 ( CARD )

         ELSE IF((CARD(1:5) == 'FORCE'   ) .OR. (CARD(1:6) == 'MOMENT'  )) THEN
            LFORCE = LFORCE +1

         ELSE IF (CARD(1:4) == 'GRAV'    )  THEN
            LGRAV = LGRAV + 1

         ELSE IF (CARD(1:6) == 'GRDSET'  )  THEN
            CALL BD_GRDSET0 ( CARD )

         ELSE IF (CARD(1:4) == 'GRID'    )  THEN
            LGRID = LGRID + 1
            LDOFG = LDOFG + 6

         ELSE IF (CARD(1:4) == 'LOAD'    )  THEN
            LLOADR = LLOADR + 1
            CALL BD_LOAD0 ( CARD, LARGE_FLD_INP, ILOAD )
            IF (ILOAD > LLOADC) THEN
               LLOADC = ILOAD
            ENDIF

         ELSE IF((CARD(1:4) == 'MAT1'    ) .OR. (CARD(1:4) == 'MAT2'    )  .OR.                                                    &
                 (CARD(1:4) == 'MAT8'    ) .OR. (CARD(1:4) == 'MAT9'    )) THEN
            LMATL = LMATL + 1

         ELSE IF((CARD(1:4) == 'MPC '    ) .OR. (CARD(1:4) == 'MPC*'    )) THEN
            LMPC = LMPC + 1
            CALL BD_MPC0 ( CARD, LARGE_FLD_INP, IMPC )
            IF (IMPC > MMPC) THEN
               MMPC = IMPC
            ENDIF

         ELSE IF (CARD(1:6) == 'MPCADD'  )  THEN
            LMPCADDR = LMPCADDR + 1
            CALL BD_MPCADD0 ( CARD, LARGE_FLD_INP, IMPCADD )
            IF (IMPCADD > LMPCADDC) THEN
               LMPCADDC = IMPCADD
            ENDIF

         ELSE IF (CARD(1:6) == 'PARAM '  )  THEN
            CALL BD_PARAM0 ( CARD )

         ELSE IF((CARD(1:5) == 'PBAR '   ) .OR. (CARD(1:5) == 'PBAR*'   ))  THEN
            LPBAR = LPBAR + 1

         ELSE IF (CARD(1:5) == 'PBARL'   )  THEN
            LPBAR  = LPBAR  + 1
            NPBARL = NPBARL + 1

         ELSE IF (CARD(1:5) == 'PBEAM'   )  THEN
            LPBEAM = LPBEAM + 1

         ELSE IF (CARD(1:5) == 'PBUSH'   )  THEN
            LPBUSH = LPBUSH + 1

         ELSE IF (CARD(1:5) == 'PCOMP'   )  THEN
            LPCOMP = LPCOMP + 1
            CALL BD_PCOMP0 ( CARD, LARGE_FLD_INP, IPLIES )
            IF (IPLIES > LPCOMP_PLIES) THEN
               LPCOMP_PLIES = IPLIES
            ENDIF

         ELSE IF (CARD(1:6) == 'PCOMP1'  )  THEN
            LPCOMP = LPCOMP + 1
            CALL BD_PCOMP10 ( CARD, LARGE_FLD_INP, IPLIES )
            IF (IPLIES > LPCOMP_PLIES) THEN
               LPCOMP_PLIES = IPLIES
            ENDIF

         ELSE IF (CARD(1:5) == 'PELAS'   )  THEN
            LPELAS = LPELAS + 1

         ELSE IF (CARD(1:6) == 'PLOAD2'  )  THEN
            LPDAT  = LPDAT  + MPDAT_PLOAD2
            LPLOAD = LPLOAD + 1

         ELSE IF (CARD(1:6) == 'PLOAD4'  )  THEN
            LPDAT  = LPDAT  + MPDAT_PLOAD4
            LPLOAD = LPLOAD + 1

         ELSE IF (CARD(1:6) == 'PLOTEL'  )  THEN
            LELE  = LELE + 1
            LEDAT = LEDAT + MEDAT_PLOTEL

         ELSE IF (CARD(1:5) == 'PMASS'   )  THEN
            LPMASS = LPMASS + 4

         ELSE IF (CARD(1:4) == 'PROD'    )  THEN
            LPROD = LPROD + 1

         ELSE IF (CARD(1:6) == 'PSHEAR'  )  THEN
            LPSHEAR = LPSHEAR + 1

         ELSE IF (CARD(1:6) == 'PSHELL'  )  THEN
            LPSHEL = LPSHEL + 1

         ELSE IF (CARD(1:6) == 'PSOLID'  )  THEN
            LPSOLID = LPSOLID + 1

         ELSE IF (CARD(1:6) == 'PUSER1'  )  THEN
            LPUSER1 = LPUSER1 + 1

         ELSE IF (CARD(1:7) == 'PUSERIN' )  THEN
            LPUSERIN = LPUSERIN + 1

         ELSE IF (CARD(1:4) == 'RBAR'    ) THEN
            LRIGEL  = LRIGEL + 1

         ELSE IF (CARD(1:4) == 'RBE1'    ) THEN
            LRIGEL  = LRIGEL + 1

         ELSE IF (CARD(1:4) == 'RBE2'    ) THEN
            LRIGEL  = LRIGEL + 1

         ELSE IF (CARD(1:4) == 'RBE3'    )  THEN
            LRIGEL  = LRIGEL + 1
            CALL BD_RBE30 ( CARD, LARGE_FLD_INP, IRBE3 )
            IF (IRBE3 > MRBE3) THEN
               MRBE3 = IRBE3
            ENDIF

         ELSE IF (CARD(1:6) == 'RFORCE'  )  THEN
            LRFORCE = LRFORCE + 1

         ELSE IF (CARD(1:7) == 'RSPLINE' )  THEN
            CALL BD_RSPLINE0 ( CARD, LARGE_FLD_INP, IRSPLINE )
            LRIGEL = LRIGEL + 1
            IF (IRSPLINE > MRSPLINE) THEN
               MRSPLINE = IRSPLINE
            ENDIF

         ELSE IF (CARD(1:5) == 'SLOAD'   )  THEN
            CALL BD_SLOAD0 ( CARD, DELTA_SLOAD )
            LSLOAD = LSLOAD + DELTA_SLOAD

         ELSE IF (CARD(1:5) == 'SEQGP'   )  THEN
            LSEQ = LSEQ + 4                                ! Conservative estimate. There can only be 4 entries per card

         ELSE IF((CARD(1:4) == 'SPC '    ) .OR. (CARD(1:4) == 'SPC*'    )) THEN
            LSPC = LSPC + 1

         ELSE IF (CARD(1:4) == 'SPC1'    )  THEN
            LSPC1 = LSPC1 + 1

         ELSE IF (CARD(1:6) == 'SPCADD'  )  THEN
            LSPCADDR = LSPCADDR + 1
            CALL BD_SPCADD0 ( CARD, LARGE_FLD_INP, ISPCADD )
            IF (ISPCADD > LSPCADDC) THEN
               LSPCADDC = ISPCADD
            ENDIF

         ELSE IF (CARD(1:6) == 'SPOINT'  )  THEN
            CALL BD_SPOINT0 ( CARD, DELTA_SPOINT )
            NSPOINT = NSPOINT + DELTA_SPOINT               ! DELTA_SPOINT = number of SPOINTS defined on this SPOINT Bulk Data entry
            LGRID   = LGRID   + DELTA_SPOINT               ! Each SPOINT counts as 1 in the number of grids
            LDOFG   = LDOFG   + DELTA_SPOINT

         ELSE IF (CARD(1:6) == 'TEMPRB'  )  THEN
            LTDAT = LTDAT + MTDAT_TEMPRB

         ELSE IF (CARD(1:6) == 'TEMPP1'  )  THEN
            LTDAT = LTDAT + MTDAT_TEMPP1

         ELSE IF (CARD(1:7) == 'ENDDATA' )  THEN
            EXIT

         ENDIF

      ENDDO

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
  101 FORMAT(A)

 1010 FORMAT(' *ERROR  1010: ERROR READING FOLLOWING ',A,' ENTRY. ENTRY IGNORED')

 1011 FORMAT(' *ERROR  1011: NO ',A10,' ENTRY FOUND BEFORE END OF FILE OR END OF RECORD IN INPUT FILE')

! **********************************************************************************************************************************

      END SUBROUTINE LOADB0

max_stars_repo_name: luzpaz/MYSTRAN

max_issues_repo_name: luzpaz/MYSTRAN

max_forks_repo_name: luzpaz/MYSTRAN


MODULE INIT_VARS
    use READ_DATA
    IMPLICIT NONE
    REAL*8 :: energy_exp
    INTEGER*2 :: rnd_spin,rnd_layer
    REAL :: rnd_accept,time_start
    INTEGER :: index_MC,index_MC_step
    INTEGER*1 :: index_AB,index_PT
    !WINDOWS
    INTEGER :: n_windows,index_window
    INTEGER*8 :: window_measures_q
    REAL*8, DIMENSION(:,:), ALLOCATABLE :: measure_windows_e,measure_windows_q,measure_windows_q2,measure_windows_q4,measure_windows_sp
    INTEGER, DIMENSION(:), ALLOCATABLE :: windows_limits
    REAL*8 :: q_aux

    CONTAINS
    SUBROUTINE DIMENSION_VARS()
        IMPLICIT NONE
        INTEGER :: i,j,num
        ALLOCATE(spin_array(n_spins,n_pt,2))
        ALLOCATE(g(n_pt))
        ALLOCATE(e_intra(n_pt,2),number_neighbours(n_spins))
        ALLOCATE(q_overlap(n_pt),q_l_overlap(n_pt))
        ALLOCATE(pntr(n_pt,2))
        ALLOCATE(histo_q(-n_spins:n_spins,n_pt))
        !Initialize pointer, gamma array and K1 K2
        DO i=0,n_pt-1
            g(i+1)=low_gamma+ dfloat(i)*((grt_gamma-low_gamma)/dfloat(n_pt-1))
            pntr(i+1,:)=i+1
        ENDDO
        !Conigurin the windows
        IF (MC_factor==1)THEN
            n_windows=3+2*(MC_steps_expoent-2)
        elseif(MC_factor==3)THEN
            n_windows=3+2*(MC_steps_expoent-2)+1
        else
            print*,'bad input for MC factor, 1 or 3'
            stop
        endif
        ALLOCATE(measure_windows_q(n_pt,n_windows), measure_windows_e(n_pt,n_windows), windows_limits(n_windows))
        ALLOCATE(measure_windows_q2(n_pt,n_windows),measure_windows_q4(n_pt,n_windows), measure_windows_sp(n_pt,n_windows))
        j=1
        num=10
        DO i=1,n_windows
            IF(j.le.n_windows)THEN
                windows_limits(j)=num
                j=j+1
                num=num*3
                IF(j.le.n_windows)THEN
                    windows_limits(j)=num
                    j=j+1
                    num=int(real(10*num)/real(3))
                END IF
            END IF
        END DO
        DO i=n_windows,2,-1
            windows_limits(i)=windows_limits(i)-windows_limits(i-1)
        END DO
        window_measures_q=int(0,8)
        measure_windows_e=0d0
        measure_windows_q=0d0
        measure_windows_q2=0d0
        measure_windows_q4=0d0
        measure_windows_sp=0d0

    END SUBROUTINE DIMENSION_VARS
    SUBROUTINE RANDOM_SPIN_ARRAY() !UPDATED
        IMPLICIT NONE
        INTEGER i,j,z
        REAL aux_rand_number,r1279
            DO i=1,n_pt
                DO j=1,n_spins
                        DO z=1,2
                            aux_rand_number=r1279()
                            IF(aux_rand_number.ge.0.5) THEN
                                spin_array(j,i,z)=1
                            ELSE
                                spin_array(j,i,z)=-1
                            END IF
                        END DO
                END DO
            END DO
    END SUBROUTINE
END MODULE INIT_VARS


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!
!     Cuts a convex part of the master surface with a slave surface
!     inserts new active edges into iactiveline for current triangle
!     calling sutherland-hodgeman algorithm
!
!     [in]       nopes		number of slave nodes for current slave surface
!     [in]       slavstraight	plane equations for mean slave plane
!     [in]       xn			mean slave normal
!     [in]       xns		slave normals in nodes of slave surface
!     [in]       xl2s		current positions of lsave nodes
!     [in]       xl2sp		projected positions of slave nodes
!     [in]       ipe	 	(i) pointer to ime for node i
!     [in]       ime     		... cataloging the edges with node i
!     [in,out]   iactiveline 	field storing the active master triange lines for the active line search
!     [in,out]   nactiveline 	number of active lines
!     [in]       nelemm		current master element
!     [in,out]   nintpoint		number of generated integration points
!     [in,out]   pslavsurf		field storing position xil, etal and weight for integration point on slave side
!     [in,out]   imastsurf		pointer into pmastsurf
!     [in,out]   pmastsurf		field storing position and weights for integration points on master side
!     [in,out]   xl2m		coordinates of master nodes for current master face
!     [in,out]   nnodelem		number of nodes in curretn master element
!     [in,out]   xl2m2		resorted coordinates of master nodes for current master face
!     [in]       nmp		number of master nodes for current master face
!     [in]       nodem		node numbers for current master face
!     [in,out]   gapmints		gap evaluated at the integration points
!     [in]       issurf		current slave surface index
!     [in,out]   areaslav		current covering of the slave surface (in the reference element)
!     [in]       debug		debug output parameter
!
      subroutine treatmasterface_mortar(
     &     nopes,slavstraight,xn,xns,xl2s,xl2sp,
     &     ipe,ime,iactiveline,nactiveline,
     &     nelemm,nintpoint,pslavsurf,
     &     imastsurf,pmastsurf,xl2m,nnodelem,xl2m2,nmp,
     &     nodem,gapmints,issurf,areaslav,debug,clearance,
     &     cl2s,cl2m,shrink,reltime)
!
!     Author: Saskia Sitzmann
!
      implicit none
!
      logical debug,shrink
!
      integer nvertex,nopes,ipe(*),ime(4,*),iactiveline(3,*),
     &     nactiveline,ifreeintersec,nmp,i,j,k,nintpoint,imastsurf(*),
     &     nnodelem,ijk,nodem(*),modf,nelemm,k_max,issurf,ii,jj,nipold
!
      real*8 pvertex(3,13),slavstraight(36),xn(3),
     &     xilm,etlm,xnl(3),clearance,p(3),dist,
     &     xl2s(3,*),p1(3),p2(3),pslavsurf(3,*),
     &     xil,etl,area,areax,areay,areaz,pmastsurf(2,*),
     &     xl2m(3,8),xl2m2(3,8),al,gapmints(*),err,xns(3,8),
     &     xl2sp(3,*),xl2mp(3,8),cgp(3),spm,spmc,
     &     pm(3),ph(3),reltime,ps(3),xit(3),etat(3),phc(3),
     &     areaslav,cl2s(3,8),cl2m(3,8),psc(3),pmc(3)
!
!
!
      data ijk /0/
      save ijk
!
      include "gauss.f"
!
      ifreeintersec=0
      err=1.d-6
      nvertex=0
      nipold=nintpoint
!
      if(debug) write(20,*) 'TT:slavsurf',issurf,' melem',nelemm
      if(debug) write(20,*) 'TT:xn',(xn(1:3))
!
!     Project master nodes to meanplane, needed for Sutherland-Hodgman
!
      do j=1,nmp
        al=-xn(1)*xl2m2(1,j)-xn(2)*
     &       xl2m2(2,j)-xn(3)*
     &       xl2m2(3,j)-slavstraight(nopes*4+4)
        do k=1,3
          xl2mp(k,j)=xl2m2(k,j)+al*xn(k)
        enddo
      enddo
!
      if(debug) then
        write(20,*) 'TT: xl2sp'
        do j=1,nopes
          write(20,*)(xl2sp(k,j),k=1,3)
        enddo
        write(20,*)'TT:xl2mp'
        do j=1,nmp
          write(20,*)(xl2mp(k,j),k=1,3)
        enddo
      endif
!
 111  format(3(e27.20))
!
!     call Sutherland-Hodgman Algo
!
      call sutherland_hodgman(nopes,xn,xl2sp,xl2mp,nodem,
     &     ipe,ime,iactiveline,nactiveline,
     &     ifreeintersec,nelemm,nmp,
     &     nvertex,pvertex)
!
!     do we have a degenerated triangle?
!
      if(debug)then
        write(20,*) 'nelemm',nelemm ,'p_new'
        do k=1,nvertex
          write(20,111)(pvertex(i,k),i=1,3)
        enddo
      endif
!
      do k=1,3
        cgp(k)=0.0
      enddo
      if(nvertex.lt.3) return
!
      if(nvertex.eq.3)then
        do k=1,3
          cgp(k)=pvertex(k,nvertex)
        enddo
        nvertex=nvertex-1
        k_max=1
      else
        do i=1,nvertex
          do k=1,3
            cgp(k)=cgp(k)+pvertex(k,i)/nvertex
          enddo
        enddo
        k_max=nvertex
      endif
!
      if(debug)then
        write(20,*) 'TT: nactiveline',nactiveline,'iactiveline'
        write(20,*) (iactiveline(1,k),k=1,nactiveline)
        write(20,*) 'cg' ,(cgp(k),k=1,3)
        write(20,*)'********************************************'
      endif
!
!     Project center point back on slave face
!
      call attachline(xl2s,cgp,nopes,xit(3),etat(3),xn,p,dist)
!
!     generating integration points on the slave surface S
!
      do k=1,k_max
!
!     storing the triangulation of the slave surfaces
!
        if(debug) then
          ijk=ijk+1
          write(40,100) ijk,(cgp(i),i=1,3)
          ijk=ijk+1
          write(40,100) ijk,(pvertex(i,modf(nvertex,k)),i=1,3)
          ijk=ijk+1
          write(40,100) ijk,(pvertex(i,modf(nvertex,k+1)),i=1,3)
          write(40,101) ijk-2,ijk-2,ijk-1
          write(40,101) ijk-1,ijk-1,ijk
          write(40,101) ijk,ijk,ijk-2
        endif
 100    format('PNT ',i10,'P',3(1x,e21.14))
 101    format('LINE ',i10,'L',i10,'P ',i10,'P')
!
!     Project back on slave surface
!
        call attachline(xl2s,pvertex(1:3,modf(nvertex,k)),
     &       nopes,xit(1),etat(1),xn,p,dist)
        call attachline(xl2s,pvertex(1:3,modf(nvertex,k+1)),
     &       nopes,xit(2),etat(2),xn,p,dist)
        p1(1)=xit(1)-xit(3)
        p1(2)=etat(1)-etat(3)
        p1(3)=0
        p2(1)=xit(2)-xit(3)
        p2(2)=etat(2)-etat(3)
        p2(3)=0
        areax=((p1(2)*p2(3))-(p2(2)*p1(3)))**2
        areay=(-(p1(1)*p2(3))+(p2(1)*p1(3)))**2
        areaz=((p1(1)*p2(2))-(p2(1)*p1(2)))**2
        area=dsqrt(areax+areay+areaz)/2.
        if(area.lt.1.e-4) cycle
        if((nopes.eq.4 .or. nopes.eq.8)
     &       .and.areaslav+area-4.0.gt.1.e-3
     &       .and.nactiveline.gt.0)then
          write(*,*)'TT: face',issurf,'loop in slavintmortar'
          write(*,*)'area',areaslav,'+',area,'.gt.4!',nactiveline
          nactiveline=0
          return
        endif
        if((nopes.eq.3.or.nopes.eq.6)
     &       .and.areaslav+area-0.5.gt.1.e-4
     &       .and.nactiveline.gt.0)then
          write(*,*)'TT: face',issurf,'loop in slavintmortar'
          write(*,*)'area',areaslav,'+',area,'.gt.0.5!'
          write(20,*)'TT: face',issurf,'loop in slavintmortar'
          write(20,*)'area',areaslav,'+',area,'.gt.0.5!'
          nactiveline=0
          return
        endif
        areaslav=areaslav+area
        if(debug)then
          write(20,106) k,area,areaslav
 106      format('tri',i10,' area',e15.8,' atot',e15.8)
          write(20,*) '# fuer itri',k,'werden 7intp gen.'
        endif
!
!     7 points scheme
!
        do i=1,7
          xil=xit(3)*gauss2d6(1,i)+
     &         xit(1)*gauss2d6(2,i)+
     &         xit(2)*(1-gauss2d6(1,i)-gauss2d6(2,i))

          etl=etat(3)*gauss2d6(1,i)+
     &         etat(1)*gauss2d6(2,i)+
     &         etat(2)*(1-gauss2d6(1,i)-gauss2d6(2,i))
!
          call evalshapefunc(xil,etl,xns,nopes,xnl)
          call evalshapefunc(xil,etl,xl2s,nopes,ps)
!
          nintpoint=nintpoint+1
!
!     projection of the master integration point onto the
!     master surface in order to get the local coordinates
!     own xn for every integration point?
!
          call attachline(xl2m,ps,nnodelem,xilm,etlm,xn,p,dist)
          call evalshapefunc(xilm,etlm,xl2m,nnodelem,pm)
!
          if(debug)then
            ijk=ijk+1
            write(40,100) ijk,(ps(jj),jj=1,3)
            ijk=ijk+1
            write(40,100) ijk,(pm(jj),jj=1,3)
            write(40,101) ijk-1,ijk-1,ijk
          endif
!
!     Calculation of the gap function at the integration point
!
          do ii=1,3
            ph(ii)=pm(ii)-ps(ii)
          enddo

          al=sqrt(ph(1)**2+ph(2)**2+ph(3)**2)
          spm=ph(1)*xn(1)
     &         +ph(2)*xn(2)
     &         +ph(3)*xn(3)
!
          if(.not.((clearance.gt.1.2357111316d0).and.
     &         (clearance.lt.1.2357111318d0))) then
!
!     assuming zero displacements at the beginning of the calculation
!     only valid for small tangential movement in the contact zone
!
            call evalshapefunc(xil,etl,cl2s,nopes,psc)
            call evalshapefunc(xilm,etlm,cl2m,nnodelem,pmc)
            do ii=1,3
              phc(ii)=pmc(ii)-psc(ii)
            enddo
            spmc=phc(1)*xn(1)
     &           +phc(2)*xn(2)
     &           +phc(3)*xn(3)
            spm=spm-spmc+clearance
          endif
          if(shrink)then
            spm=reltime*spm
          endif
          gapmints(nintpoint)=spm
          pslavsurf(1,nintpoint)=xil
          pslavsurf(2,nintpoint)=etl
!
!     weight add to 0.5 \hat{w}_p=A_triaref*w_p=0.5*w_p
!     division by 0.5 to get to original weights...
!
          pslavsurf(3,nintpoint)=area*weight2d6(i)/0.5
          pmastsurf(1,nintpoint)=xilm
          pmastsurf(2,nintpoint)=etlm
          imastsurf(nintpoint)=nelemm
          if(debug)then
            write(30,201) xil,etl,pslavsurf(3,nintpoint),spm,nelemm
          endif
 201      format('xil ',1x,e15.8,2x,'etl',2x,e15.8,2x,'w ',e15.8,2x,
     &         e15.8,2x,'M',i10)
        enddo
!
      enddo
c     write(20,*)'********************'
      return
      end

ext: f

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
c     Program to separate second USAF update into station files
c     AJ_Kettle, 29Mar2021

c     05Apr2021: newsci6: 1 file test run to yield 12943 stn files
c     06Apr2021: newsci6: 1 file from 2019 processed c/c=0.34/1.73min
c     06Apr2021 0321: newsci6: stopped at i=44 group 2019; 5h run time 14197 stns
c     08Apr2021 0311: newsci6: stopped at i=78 group 2019; 10h run time 14274 stns
c      runtimes 1: 5.15min; 23: 3.75min...; 76: 12.80min; 77: 13.75min
c     09Apr2021 0300: newsci6: stopped at i=244 group 2019; 18h run; 14585 stns
c      runtimes 241: 14.32min; 242: 12.23min; 12.77min
c     09Apr2021 2008: newsci6: end of run; 14725 stns
c      runtimes 280: 0.33/11.53min; 281: 0.34/12.65min

c************************************************************************
c2-------10--------20--------30--------40--------50--------60--------70--------80
      IMPLICIT NONE
c************************************************************************
c     Declare variables

      INTEGER             :: i,j,k,ii,jj,kk

      REAL                :: f_time_st,f_time_en

      REAL                :: f_deltime_cpu_s
      REAL                :: f_deltime_clock_s

      CHARACTER(LEN=8)    :: s_date_st
      CHARACTER(LEN=10)   :: s_time_st
      CHARACTER(LEN=5)    :: s_zone_st
      INTEGER             :: i_values_st(8)
      CHARACTER(LEN=8)    :: s_date_en
      CHARACTER(LEN=10)   :: s_time_en
      CHARACTER(LEN=5)    :: s_zone_en
      INTEGER             :: i_values_en(8)

      REAL                :: f_clock_st_s,f_clock_en_s
c*****
      CHARACTER(LEN=300)  :: s_directory_filelist
      CHARACTER(LEN=300)  :: s_directory_dataroot1
      CHARACTER(LEN=300)  :: s_directory_dataroot2
      CHARACTER(LEN=300)  :: s_directory_dataroot3

c     directories for intermediate and final steps
      CHARACTER(LEN=300)  :: s_directory_1tar
      CHARACTER(LEN=300)  :: s_directory_2untarpack
      CHARACTER(LEN=300)  :: s_directory_3unsort
      CHARACTER(LEN=300)  :: s_directory_3headerlist
      CHARACTER(LEN=300)  :: s_directory_3archive

      CHARACTER(LEN=300)  :: s_subdir

      CHARACTER(LEN=300)  :: s_name_filelist1
      CHARACTER(LEN=300)  :: s_name_filelist2
      CHARACTER(LEN=300)  :: s_name_filelist3
c*****
c     Standard ASCII characters
      INTEGER             :: l_ascii
      INTEGER             :: i_list_ascii(62)
      CHARACTER(LEN=1)    :: s_list_asciichar(62)
c*****
      INTEGER             :: l_tar1
      CHARACTER(LEN=300)  :: s_vec_path1(2000)
      CHARACTER(LEN=300)  :: s_vec_name1(2000)

      INTEGER             :: l_tar2
      CHARACTER(LEN=300)  :: s_vec_path2(2000)
      CHARACTER(LEN=300)  :: s_vec_name2(2000)

      INTEGER             :: l_tar3
      CHARACTER(LEN=300)  :: s_vec_path3(2000)
      CHARACTER(LEN=300)  :: s_vec_name3(2000)
c*****
      CHARACTER(LEN=300)  :: s_command
      INTEGER             :: io
c************************************************************************
      print*,'start up2'
c************************************************************************
c     Find start time
      CALL CPU_TIME(f_time_st)

c     Find date & time
      CALL DATE_AND_TIME(s_date_st,s_time_st,s_zone_st,i_values_st)
c************************************************************************

c************************************************************************
c     Directory of files

      s_directory_filelist=
     +  '/home/users/akettle/Work/P20210325_usaf_update2/'//
     +  'Datafile_listing/'
c     Note root is in GWS for update1 & scratch-pw for update2
      s_directory_dataroot1=
     + '/work/scratch-pw/akettle/P20210325_usaf_update2/Data_source/'//
     + '2019/'
      s_directory_dataroot2=
     + '/work/scratch-pw/akettle/P20210325_usaf_update2/Data_source/'//
     + '2020/'
      s_directory_dataroot3=
     + '/work/scratch-pw/akettle/P20210325_usaf_update2/Data_source/'//
     + '2021/'

c      s_directory_dataroot=
c     + '/gws/nopw/j04/c3s311a_lot2/data/incoming/usaf_update'

      s_name_filelist1='file2019.dat'
      s_name_filelist2='file2020.dat'
      s_name_filelist3='file2021.dat'
c************************************************************************
      CALL get_directory_paths20210330(
     + s_directory_1tar,s_directory_2untarpack,
     + s_directory_3unsort,s_directory_3headerlist,s_directory_3archive)
c************************************************************************
      GOTO 41

c     Remove all output files

      s_command=
     +  'rm '//TRIM(s_directory_1tar)//'*'
      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

      s_command=
     +  'rm '//TRIM(s_directory_2untarpack)//'*'
      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

      s_command=
     +  'rm '//TRIM(s_directory_3unsort)//'*.csv'
      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

      s_command=
     +  'rm '//TRIM(s_directory_3headerlist)//'*'
      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

      s_command=
     +  'rm '//TRIM(s_directory_3archive)//'*'
      print*,'s_command=',TRIM(s_command)
      CALL SYSTEM(s_command,io)

c      STOP 'just erased output directories; up2'

 41   CONTINUE
c************************************************************************
c     Get list of ASCII characters

      CALL get_ascii_list(l_ascii,i_list_ascii,s_list_asciichar)
c************************************************************************
c     Import file listing
      CALL import_file_listing20210401(
     +  s_directory_filelist,s_name_filelist1,
     +  s_directory_dataroot1,
     +  l_tar1,
     +  s_vec_path1,s_vec_name1)

      CALL import_file_listing20210401(
     +  s_directory_filelist,s_name_filelist2,
     +  s_directory_dataroot2,
     +  l_tar2,
     +  s_vec_path2,s_vec_name2)

      CALL import_file_listing20210401(
     +  s_directory_filelist,s_name_filelist3,
     +  s_directory_dataroot3,
     +  l_tar3,
     +  s_vec_path3,s_vec_name3)

      print*,'l_tar1,l_tar2...=',l_tar1,l_tar2,l_tar3
      print*,'s_vec_path1=',(TRIM(s_vec_path1(i)),i=1,1)
      print*,'s_vec_name1=',(TRIM(s_vec_name1(i)),i=1,1)
      print*,'s_vec_path2=',(TRIM(s_vec_path2(i)),i=1,1)
      print*,'s_vec_name2=',(TRIM(s_vec_name2(i)),i=1,1)
      print*,'s_vec_path3=',(TRIM(s_vec_path3(i)),i=1,1)
      print*,'s_vec_name3=',(TRIM(s_vec_name3(i)),i=1,1)

c      s_command=TRIM(s_vec_path1(1))//TRIM(s_vec_name1(1))

c      print*,'s_commandup='//TRIM(s_command)

c      STOP 'up2'
c************************************************************************
c     Process data files
      CALL process_files(l_tar1,s_vec_path1,s_vec_name1,
     +  s_directory_1tar,s_directory_2untarpack,
     +  s_directory_3headerlist,s_directory_3unsort,
     +  s_directory_3archive,
     +  l_ascii,i_list_ascii,s_list_asciichar)
c************************************************************************

c************************************************************************
c     Find end time
      CALL CPU_TIME(f_time_en)

c     Find date & time
      CALL DATE_AND_TIME(s_date_en,s_time_en,s_zone_en,i_values_en)
c************************************************************************
      print*,'s_date_en=',TRIM(s_date_en)
      print*,'s_time_en=',TRIM(s_time_en)
      print*,'i_values_st',(i_values_st(i),i=1,8)
      print*,'i_values_en',(i_values_en(i),i=1,8)

      f_clock_st_s=FLOAT(i_values_st(3))*24.0*60.0*60.0+
     +  FLOAT(i_values_st(5))*60.0*60.0+FLOAT(i_values_st(6))*60.0+
     +  FLOAT(i_values_st(7))
      f_clock_en_s=FLOAT(i_values_en(3))*24.0*60.0*60.0+
     +  FLOAT(i_values_en(5))*60.0*60.0+FLOAT(i_values_en(6))*60.0+
     +  FLOAT(i_values_en(7))

      print*,'f_clock_st_s,f_clock_en_s=',f_clock_st_s,f_clock_en_s

      f_deltime_clock_s=f_clock_en_s-f_clock_st_s
c************************************************************************
      f_deltime_cpu_s=f_time_en-f_time_st
      print*,'f_deltime_cpu_s,f_deltime_cpu_min=',
     +  f_deltime_cpu_s,f_deltime_cpu_s/60.0
      print*,'f_deltime_clock_s,f_deltime_clock_min=',
     +  f_deltime_clock_s,f_deltime_clock_s/60.0
c************************************************************************

      END

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


text:
	SUBROUTINE GENANM  ( iret )
C************************************************************************
C* GENANM								*
C* 									*
C* This subroutine ends an animation sequence.   			*
C* 									*
C* GENANM  ( IRET )							*
C*									*
C* Output parameters:							*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log									*
C* A. Chang/EAI	 	12/93						*
C* L. Williams/EAi       3/94   Renamed file from GEPIXM to GENANM and  *
C*                              removed blank comments from header      *
C************************************************************************
	INCLUDE		'FUNCCODE.PRM'
	INCLUDE		'ERROR.PRM'
C*
	INTEGER		isend (2)
C------------------------------------------------------------------------
C*	Load input parameters into buffer and write them to the mailbox.
C
	isend (1) = 2
	isend (2) = FENANM
C
	CALL GPUT  ( isend, 2, iret )
	IF  ( iret .ne. NORMAL )  RETURN
	CALL GGET  ( iret, 1, ierr )
	IF  ( ierr .ne. NORMAL )  iret = ierr
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
! subroutine cal_angles
SUBROUTINE cal_angles(lam_p,phip_p,tol_lam,orb_elements,spheroid, &
             smodel,track,num_tpts,timet,theta_p,azimuth,istat)

!     Cal_Angles is the main routine in the set
!     It finds the segment of the track model where a line from a point
!     hits the track line at an angle approx perpendicular to
!     the heading. It then refines the estimated crossing point
!     The angle from the track point to the given point
!     assumes a great circle relationship
!     Returned data are:
!          timet (time of track point -t0) secs
!          theta_p (look angle from the track to the given point) radians
!          azimuth (azimuth from track point to given point) radians
!          (The azimuth uses the BRDF convention)

!   * Re-written as an independent subroutine by JS, Aug 2014

!   Inputs:
!       lam_p
!       phip_p
!       tol_lam
!       orb_elements
!           1. Orbital inclination (degrees)
!           2. Semi_major radius (m)
!           3. Angular velocity (rad sec-1)
!       spheroid
!           1. Spheroid major axis
!           2. Inverse flattening
!           3. Eccentricity squared
!           4. Earth rotational angular velocity rad/sec
!       smodel
!           1. phi0
!           2. phi0_p
!           3. rho0
!           4. t0
!           5. lam0
!           6. gamm0
!           7. beta0
!           8. rotn0
!           9. hxy0
!           10. N0
!           11. H0
!           12. th_ratio0
!       track
!           1. t
!           2. rho
!           3. phi_p
!           4. lam
!           5. beta
!           6. hxy
!           7. mj
!           8. skew
!   Outputs
!       timet
!       theta_p
!       azimuth
!       istat

    use sys_variables, only : pi, d2r, r2d

    implicit none

    double precision, intent(in) :: lam_p, phip_p, tol_lam
    double precision, dimension(4), intent(in) :: spheroid
    double precision, dimension(3), intent(in) :: orb_elements
!   smodel(phi0,phi0_p,rho0,t0,lam0,gamm0,beta0,rotn0,hxy0,N0,H0,th_ratio0)
    double precision, dimension(12), intent(in) :: smodel

    integer, intent(in) :: num_tpts
!   track(t,rho,phi_p,lam,beta,hxy,mj,skew)
    double precision, dimension(num_tpts,8), intent(in) :: track
    double precision, dimension(num_tpts) :: phi_p, lam, beta
    double precision, dimension(num_tpts) :: hxy, mj

!f2py depend(num_tpts), track, phi_p, lam, beta, hxy, mj

    double precision orad
    double precision asph, e2
    double precision t0, hxy0

    double precision, intent(out) :: timet, theta_p, azimuth
    integer, intent(out) :: istat

    integer west, neg, hit, j, negp
    integer left, right
    double precision rho_p, t_p, lamcal, betacal, adiff
    double precision delta_mu, tbet_p, phi_pt, lamt, dmu_prev, mjt, lam_test
    double precision phip_test, rho_c, t_c
    double precision theta_e, rn_p, ratio, temp

!   smodel parameters
    t0 = smodel(4)
    hxy0 = smodel(9)

!   track parameters
    phi_p(:) = track(:,3)
    lam(:) = track(:,4)
    beta(:) = track(:,5)
    hxy(:) = track(:,6)
    mj(:) = track(:,7)

!   Satellite orbital parameters
!   semi_major radius (m)
    orad = orb_elements(2)

!   Spheroid parameters
!   Spheroid major axis
!   Eccentricity squared
    asph = spheroid(1) !6378137d0
    e2 = spheroid(3)

!   Initialise the return values
    istat = 0
    timet = 0.0
    theta_p = 0.0
    azimuth = 0.0

    temp = phip_p
    call q_cal(temp,orb_elements,spheroid,smodel,rho_p,t_p,lamcal, &
           betacal,istat)

!   find where the azimuth of the track crossing at the same phi_p is
!   relative to input and use it to decide what side of the line
!   the point is or if it is on the track
    adiff = lam_p-lamcal

!   if on the track easy - done and go back
    if (abs(adiff) .lt. tol_lam) then
        timet = t_p-t0
        theta_p = 0.0d0
        azimuth = 0.0d0
        go to 99
    endif

!   look to see if the point is west or east of the track
    if (adiff .lt. 0.0) then
        west = 1
    else
        west = 0
    endif

!   now look over the sample points on the track
!   to find a bounding interval (the time interval between two track points)
!   for the solution
!
!   method is to define a line through given point expressed as
!   L(lam,phi)=0 and evaluate L at the track points
!   slope base on the above East/West point on the
!   the track
    delta_mu = 0.0d0

!   tbet_p is tan(beta)*hy/hx
    tbet_p = hxy0*tan(betacal)

    do j=1,num_tpts
        hit = 0
        phi_pt = dble(phi_p(j))
        lamt = dble(lam(j))
        if (j .gt. 1) then
            dmu_prev = delta_mu
        endif

!   Equation L(lamt,phi_pt)
        delta_mu = (phi_pt-phip_p)+tbet_p*(lamt-lam_p)
        if (j .eq. 1) then
            if (delta_mu .lt. 0.0) then
                neg = 1
            else
                neg = 0
            endif
            dmu_prev = delta_mu
        else
            if (delta_mu .lt. 0.0) then
                negp = 1
            else
                negp = 0
            endif
            if (negp .ne. neg) then
                hit = 1
                left = j-1
                right = j
                goto 10
            endif
        endif
    enddo

!   if no segment found it is bad news! Something is wrong.
    print*, 'error in cal_angles - no interval found in track'
    istat = 99
    go to 99

10  continue

!   Now solve for the crossing of a line from the given point
!   to the track interval (assumed a line joining ends)
!   using the local heading and scaling information
    tbet_p = (hxy(left)*tan(beta(left))+hxy(right)* &
      tan(beta(right)))/2.0d0
    mjt = mj(left)

    phi_pt = dble(phi_p(left))
    lamt = dble(lam(left))

    lam_test = ((phip_p-phi_pt)+tbet_p*lam_p+mjt*lamt)/(tbet_p+mjt)
    phip_test = phi_pt+mjt*(lam_test-lamt)

!   get info at the solution
    call q_cal(phip_test,orb_elements,spheroid,smodel,rho_c,t_c, &
           lamt,betacal,istat)

!   calculate the return results
!   theta_e is angle at Earth centre
!   theta_p is the look angle
    phi_pt = phip_test
    theta_e = acos(sin(phip_p)*sin(phi_pt)+ &
      cos(phip_p)*cos(phi_pt)*cos(lam_p-lamt))

    RN_p = dble(asph)/sqrt(1.0d0-dble(e2)*sin(phi_pt)**2)
    ratio = (RN_p*cos(theta_e))/(orad-RN_p*cos(theta_e))
    theta_p = atan(ratio*tan(theta_e))
    if (west.gt.0) then
        azimuth = betacal+(dble(pi)/2.0d0)
    else
        azimuth = 3.0d0*dble(pi)/2.0d0+betacal
    endif
    timet = t_c-t0

!   return to caller
99  continue
    return

END SUBROUTINE cal_angles

max_stars_repo_name: Oceancolour-RG/wagl

max_issues_repo_name: Oceancolour-RG/wagl

max_forks_repo_name: Oceancolour-RG/wagl


text:
${PRJ_DIR}/hw/top_earlgrey/rtl/top_pkg.sv
${PRJ_DIR}/hw/ip/tlul/rtl/tlul_pkg.sv
${PRJ_DIR}/hw/ip/prim/rtl/prim_pkg.sv

${PRJ_DIR}/hw/ip/prim/rtl/prim_assert.sv
${PRJ_DIR}/hw/ip/prim/rtl/prim_flop_2sync.sv
${PRJ_DIR}/hw/ip/prim/rtl/prim_fifo_async.sv
${PRJ_DIR}/hw/ip/prim/rtl/prim_fifo_sync.sv
${PRJ_DIR}/hw/ip/prim/rtl/prim_clock_inverter.sv
${PRJ_DIR}/hw/ip/prim/abstract/prim_clock_mux2.sv
${PRJ_DIR}/hw/ip/prim_generic/rtl/prim_generic_clock_mux2.sv

${PRJ_DIR}/hw/vendor/pulp_riscv_dbg/src/dm_pkg.sv
${PRJ_DIR}/hw/vendor/pulp_riscv_dbg/debug_rom/debug_rom.sv
${PRJ_DIR}/hw/vendor/pulp_riscv_dbg/src/dm_sba.sv
${PRJ_DIR}/hw/vendor/pulp_riscv_dbg/src/dm_csrs.sv
${PRJ_DIR}/hw/vendor/pulp_riscv_dbg/src/dm_mem.sv
${PRJ_DIR}/hw/vendor/pulp_riscv_dbg/src/dmi_cdc.sv
${PRJ_DIR}/hw/vendor/pulp_riscv_dbg/src/dmi_jtag.sv
${PRJ_DIR}/hw/vendor/pulp_riscv_dbg/src/dmi_jtag_tap.sv
${PRJ_DIR}/hw/ip/tlul/rtl/tlul_adapter_sram.sv
${PRJ_DIR}/hw/ip/tlul/rtl/tlul_err.sv


        !COMPILER-GENERATED INTERFACE MODULE: Thu Oct 25 10:21:43 2018
        MODULE READ_VREF__genmod
          INTERFACE
            SUBROUTINE READ_VREF(AR,MD)
              CHARACTER(LEN=8) :: AR
              CHARACTER(LEN=8) :: MD
            END SUBROUTINE READ_VREF
          END INTERFACE
        END MODULE READ_VREF__genmod


text:
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 2015 Adrian Martinez Vargas                            %
!                                                                      %
! This software may be modified and distributed under the terms        %
! of the MIT license.  See the LICENSE.txt file for details.           %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 1996, The Board of Trustees of the Leland Stanford     %
! Junior University.  All rights reserved.                             %
!                                                                      %
! The programs in GSLIB are distributed in the hope that they will be  %
! useful, but WITHOUT ANY WARRANTY.  No author or distributor accepts  %
! responsibility to anyone for the consequences of using them or for   %
! whether they serve any particular purpose or work at all, unless he  %
! says so in writing.  Everyone is granted permission to copy, modify  %
! and redistribute the programs in GSLIB, but only under the condition %
! that this notice and the above copyright notice remain intact.       %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-----------------------------------------------------------------------
!
! The functions and subroutines below were modified from the
! version 2.0 of the gslib code written in fortran 77
!
! The objective is to add functionality to GSLIB and to link
! this code with python using f2py. It uses f90 code convention.
! The arrays are dynamic and externally declared (from python)
!
! the code was converted from Fortran 77 to Fortran 90 using F2F.pl
!
! for more information please refer to:
! - gslib77 source code: http://www.statios.com/software/gslib77_ls.tar.gz
! - GSLIB: Geostatistical Software Library and User's Guide. Second edition
!   by Clayton V. Deutsch, Andre G. Journel, 1997.
! - F2PY Users Guide: http://docs.scipy.org/doc/numpy-dev/f2py/
! - F2F https://bitbucket.org/lemonlab/f2f
!-----------------------------------------------------------------------

Module Commons

    real :: UNEST    ! asign this externally at init , otherwise will be zero by default

End Module Commons

subroutine set_unest(unest_)

    use Commons

    real, intent(in) :: unest_

    UNEST = unest_

end subroutine set_unest

subroutine get_unest(unest_)

    use Commons

    real, intent(out) :: unest_

    unest_ = UNEST

    if (isnan(UNEST)) write(*,*) 'UNEST is a NaN'

    if (UNEST==anan) write(*,*) 'UNEST == NaN'

    return

end subroutine get_unest

subroutine postik( &
                  iout,outpar, &  ! output option, output parameter
                  nccut, ccut1, & ! number of thresholds, the thresholds
                  ivol,ivtyp,varred, & ! volume support?, type, varred
                  zmin,zmax, & ! minimum and maximum Z value
                  ltail,ltpar, & ! lower tail: option, parameter
                  middle,mpar, & ! middle    : option, parameter
                  utail,utpar, & ! upper tail: option, parameter
                  maxdis, & ! maximum discretization
                  vr, wt, nc, & ! global distribution
                  p, na, &  ! input data
                  out1,out2,out3, &  ! output data
                  error) ! to control the error

    !-----------------------------------------------------------------------

    !                   Post Process IK Distributions
    !                   *****************************

    ! Reads IK-generated distributions and post processes them
    ! according to user specifications.

    !      - affine or indirect lognormal volume support correction
    !      - E-type mean
    !      - probability of exceeding specified threshold and the mean
    !           value above threshold.
    !      - compute a z-threshold value for a specified CDF value.

    ! See the text for detailed discussion of parameters to extrapolate the
    ! discretely coded cdf values.

    ! INPUT/OUTPUT Parameters:

    !   distin           the input distributions (output from IK3D)
    !   outfl            the output file for E-type,....
    !   iout,outpar         =1 E-type,
    !                       =2 prob and grade > and <= outpar;
    !                       =3 Z percentile corresponding to outpar,
    !                       =4 conditional variance.
    !   nccut            the number of cutoffs
    !   ccut(i)          cutoffs
    !   ivol             =1, the consider volume support corrections
    !   ivtyp            =1, affine; =2 indirect lognormal
    !   varred           variance adjustment factor, between 0 and 1
    !   datafl           the global data distribution
    !   ivr,iwt,tmin     parameters to read in global data distribution
    !   varred           variance reduction parameter
    !   zmin,zmax        minimum and maximum Z values
    !   ltail,ltpar      option and parameter to handle values in lower tail
    !   middle,mpar      option and parameter to handle values in the middle
    !   utail,utpar      option and parameter to handle values in upper tail
    !   maxdis           discretization to compute E-type and mean > cutoff

    !-----------------------------------------------------------------------

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !    iout == 1 out1 -> 'mean', out2 -> man, out3 -> man
    !    iout == 2 out1 -> 'prob > cutoff', out2 -> 'mean > cutoff', out3 -> 'mean < cutoff'
    !    iout == 3 out1 -> 'Z value corresponding to CDF = value', out2 -> man, out3 -> man
    !    iout == 4 out1 -> 'Conditional Variance', out2 -> man, out3 -> man
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use commons

    ! external variables (in)
    real, intent(in) :: outpar, varred, zmin,zmax,ltpar,mpar,utpar
    integer, intent(in) :: iout, nccut, ivol,ivtyp, ltail,middle,utail,maxdis,nc, na
    real, intent(in), dimension(nccut) :: ccut1
    real, intent(in), dimension(nc) :: vr, wt
    real, intent(in), dimension(na,nccut) :: p

    ! external variables (out)
    real, intent(out), dimension(na) :: out1,out2,out3
    integer, intent(out) :: error

    ! local variables
    real :: EPSLON, outpar_, varred_
    real, dimension(nc) :: cdf,cut
    real, dimension(nccut) :: ccdf,ccdf1,ccdf2,ccut
    logical ::   testfl

    EPSLON=1.0e-6
    error = 0
    outpar_ = outpar
    varred_ = varred

    ! write(*,*) ' output option & par = ',iout,outpar_
    if(iout /= 1 .AND. iout /= 2 .AND. iout /= 3 .AND. iout /= 4) then
        write(*,*) ' ERROR: invalid output option ',iout
        error = 1
        return
    end if
    if(iout == 3) then
        if(outpar_ < 0.0) then
            write (*,*) 'Invalid p-value for iout=3'
            error = 2
            return ! stop 'Invalid p-value for iout=3'
        end if
        if(outpar_ > 1.0) then
            write (*,*) 'Invalid p-value for iout=3'
            error = 3
            return ! 'Invalid p-value for iout=3'
        end if
        outpar_ = min(max(outpar_,EPSLON),(1.0-EPSLON))
    end if

    ! write(*,*) ' number of cutoffs = ',nccut
    ! write(*,*) ' cutoffs = ',(ccut1(i),i=1,nccut)
    ! write(*,*) ' volume variance = ',ivol,ivtyp,varred_
    if(varred_ < 0.0 .OR. varred_ > 1.0) then
        write(*,*) ' ERROR: invalid variance reduction ',varred_
        write(*,*) '        must be between 0 and 1'
        error = 3
        return
    end if

    !write(*,*) ' minimum and maximum = ',zmin,zmax
    if(iout == 2 .AND. outpar_ < zmin) then
        write(*,*) 'Invalid z-value for iout=2'
        error = 4
    end if
    if(iout == 2 .AND. outpar_ > zmax) then
        write(*,*) 'Invalid z-value for iout=2'
        error = 5
    end if

    ! write(*,*) ' ltail, ltpar = ',ltail,ltpar
    ! write(*,*) ' middle, mpar = ',middle,mpar
    ! write(*,*) ' utail, utpar = ',utail, utpar
    ! write(*,*) ' discretization = ',maxdis

    ! Do we have a global distribution and do we need it anyway?

    ncut = 0
    if(ltail /= 3 .AND. middle /= 3 .AND. utail /= 3) testfl = .FALSE.

    ! Read in the global cdf ("cut" and "cdf" arrays):

    if(testfl) then
        do l=1, nc
            tcdf = 0.0
            gmean = 0.0
            ncut  = 0
            if(wt(l) <= 0.0 .or. wt(l) == UNEST .or. isnan(wt(l))) then
                write(*,*) ' ERROR: negative weights at index ', l
                error = 6
                return
            endif

            if(vr(l) == UNEST .or. isnan(vr(l)) ) then
                write(*,*) ' ERROR: nans in globa CDF at index ', l
                error = 6
                return
            endif

            ncut = ncut + 1
            cut(ncut) = vr(l)
            gmean     = gmean + vr(l)
            cdf(ncut) = wt(l)
            tcdf = tcdf + cdf(ncut)

            if(tcdf <= 0) then
                write (*,*) ' total global CDF <= 0'
                return
            end if

            gmean = gmean / tcdf

            ! Turn the (possibly weighted) distribution into a cdf that is useful:

            call sortem(1,ncut,cut,1,cdf,c,d,e,f,g,h)
            oldcp = 0.0
            cp    = 0.0
            tcdf  = 1.0 / tcdf
            do i=1,ncut
                cp     = cp + cdf(i) * tcdf
                cdf(i) =(cp + oldcp) * 0.5
                oldcp  = cp
            end do

            ! Get median and write some info to the screen:

            call locate(cdf,ncut,1,ncut,0.5,j)
            gmedian = powint(cdf(j),cdf(j+1),cut(j),cut(j+1),0.5,1.)
            ! write(*,*) 'Global cdf from file: ',datafl
            ! write(*,*) '   number of data: ',ncut
            ! write(*,*) '   global mean:    ',gmean
            ! write(*,*) '   global median:  ',gmedian
            ! write(*,*)
        end do
    endif

    ! If we are doing the affine correction then compute the square root
    ! of the variance ``adjustment'' factor:

    if(ivol == 1 .AND. ivtyp == 1) varred_ = sqrt(max(varred_,0.0))

    ! BIG LOOP reading in each distribution in turn:

    nproc = 0
    procm = 0.0

    do l=1, na
        do i=1,nccut
            ccdf(i) = p(l,i)
        end do

        ! initialize output
        out1(l) = UNEST
        out2(l) = UNEST
        out3(l) = UNEST

        ! Check for missing values:

        if(ccdf(nccut) < -0.1) cycle  ! skip this iteration and do next iteration
        if(ccdf(nccut) == UNEST) cycle  ! skip this iteration and do next iteration
        if(isnan(ccdf(nccut))) cycle

        ! Reinstate "ccut1" because volume support may have changed "ccut":

        do i=1,nccut
            ccut(i) = ccut1(i)
        end do

        ! Correct Order Relations (the distributions coming from IK3D have
        ! already been corrected):

        do i=1,nccut
            ccdf(i)=max(min(1.0,ccdf(i)),0.0)
        end do
        ccdf1(1) = ccdf(1)
        do i=2,nccut
            ccdf1(i) = max(ccdf1(i-1),ccdf(i))
        end do
        ccdf2(nccut) = ccdf(nccut)
        do i=nccut-1,1,-1
            ccdf2(i) = min(ccdf2(i+1),ccdf(i))
        end do
        do i=1,nccut
            ccdf(i) = 0.5*(ccdf1(i)+ccdf2(i))
        end do

        ! Volume support correction with parameters (ivol,ivtyp,varred_):

        if(ivol == 1) then
            dis    = 1.0 / real(maxdis)
            cdfval = -0.5*dis
            etype  = 0.0
            ecv    = 0.0
            do i=1,maxdis
                cdfval  = cdfval + dis
                zval    = -1.0
                call beyond(1,nccut,ccut,ccdf,ncut,cut,cdf,zmin, &
                zmax,ltail,ltpar,middle,mpar,utail, &
                utpar,zval,cdfval,ierr)
                if(ierr /= 0) write(*,*) 'ERROR: ',ierr,' continuing1'
                etype = etype + zval
                ecv   = ecv   + zval*zval
            end do
            etype = etype / real(maxdis)
            ecv   = sqrt(max((ecv/real(maxdis)-etype*etype),0.0)) &
            / max(etype,EPSLON)
            if(etype == 0.0) then
                write(*,*) 'NO support correction with 0 mean'
                etype = EPSLON
            endif

            ! Affine Correction:

            if(ivtyp == 1) then
                do i=1,nccut
                    ccut(i) = etype+varred_*(ccut(i)-etype)
                end do
            else

                ! Indirect lognormal:  1. Compute parameters "a" and "b"
                !                      2. Correct quantiles
                !                      3. recompute etype
                !                      4. Recorrect quantiles to correct mean

                !            Parameters "a" and "b":

                b = sqrt(max((alog(varred_*ecv*ecv+1) &
                /alog(ecv*ecv+1)),0.0))
                a = (etype/sqrt(max((varred_*ecv*ecv+1),0.0))) * &
                (sqrt(max((ecv*ecv+1),0.0))/etype)**b

                !            Correct quantiles:

                do i=1,nccut
                    ccut(i) = a*ccut(i)**b
                end do

                !            New etype:

                cdfval = -0.5*dis
                enew  = 0.0
                do i=1,maxdis
                    cdfval  = cdfval + dis
                    zval    = -1.0
                    call beyond(1,nccut,ccut,ccdf,ncut,cut,cdf,zmin, &
                    zmax,ltail,ltpar,middle,mpar,utail, &
                    utpar,zval,cdfval,ierr)
                    if(ierr /= 0) then
                        write(*,*) 'ERROR: ',ierr,' continuing2'
                    endif
                    enew = enew + zval
                end do
                enew = enew / real(maxdis)

                !            Recorrect "quantiles":

                do i=1,nccut
                    ccut(i) = (etype/enew)*ccut(i)
                end do
            endif
        endif

        ! Compute mean of local distribution if E-type (iout=1), or need the
        ! mean above (below) a cutoff, or if performing volume support:

        if(iout <= 2 .OR. iout == 4) then
            dis    = 1.0 / real(maxdis)
            cdfval = -0.5*dis
            eabove = 0.0
            nabove = 0
            ebelow = 0.0
            nbelow = 0
            etype  = 0.0
            ecv    = 0.0

            do i=1,maxdis
                cdfval  = cdfval + dis
                zval    = -1.0
                call beyond(1,nccut,ccut,ccdf,ncut,cut,cdf,zmin, &
                zmax,ltail,ltpar,middle,mpar,utail, &
                utpar,zval,cdfval,ierr)
                if(ierr /= 0) then
                    write(*,*) 'ERROR: ',ierr,' continuing3'
                endif
                etype = etype + zval
                ecv   = ecv   + zval*zval
                if(zval <= outpar_) then
                    nbelow = nbelow + 1
                    ebelow = ebelow + zval
                else
                    nabove = nabove + 1
                    eabove = eabove + zval
                end if
            end do

            ! e-type and conditional variance:

            etype = etype / real(maxdis)
            ecv   = ecv/real(maxdis)-etype*etype
        endif

        ! Write out E-type?

        if(iout == 1) out1(l) = etype

        ! Do we need probability and mean value above a threshold?

        if(iout == 2) then

        !      Get probability:

            cdfval  = -1.0
            call beyond(1,nccut,ccut,ccdf,ncut,cut,cdf,zmin, &
            zmax,ltail,ltpar,middle,mpar,utail, &
            utpar,outpar_,cdfval,ierr)
            if(ierr /= 0) then
                write(*,*) 'ERROR: ',ierr,' continuing4'
            endif
            prob   = 1.0 - cdfval
            eabove = eabove / real(max(1,nabove))
            ebelow = ebelow / real(max(1,nbelow))
            out1(l) = prob
            out2(l) = eabove
            out3(l) = ebelow

            if (nabove < 1 ) then
                out1(l) = UNEST
                out2(l) = UNEST
                out3(l) = UNEST
            end if

        endif

        ! Do we need the "Z" value corresponding to a particular CDF?

        if(iout == 3) then
            zval    = -1.0
            call beyond(1,nccut,ccut,ccdf,ncut,cut,cdf,zmin, &
            zmax,ltail,ltpar,middle,mpar,utail, &
            utpar,zval,outpar_,ierr)
            if(ierr /= 0) then
                write(*,*) 'ERROR: ',ierr,' continuing5'
            endif
            out1(l) = zval
        endif

        ! Write out conditional variance:

        if(iout == 4) out1(l) = ecv

        ! Return for another:

        nproc = nproc + 1
        procm = procm + etype

    end do

    ! Finished:

    procm = procm / max(real(nproc),1.0)
    ! write(*,*)
    ! write(*,*) 'Number of distributions ',nproc
    ! if(procm /= 0.0) write(*,*) 'Overall mean:           ',procm
    ! if(iout == 1) write(*,*) 'Local Means (E-type):   ',outfl
    ! if(iout == 2) write(*,*) 'Prob and mean > cutoff: ',outfl
    ! if(iout == 3) write(*,*) 'Z values for outpar in: ',outfl
    ! write(*,*)

    ! Finished:

    ! write (*,*) ' POSTIK Finished'

    return

end subroutine postik

max_stars_repo_name: cmrajan/pygslib

max_issues_repo_name: kaufmanno/pygslib

max_forks_repo_name: kaufmanno/pygslib


!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 3D boundary conditon and guard cell
!!

!! X direction, X-Z surface
!!
subroutine bc_val_3d_l(sxyz, ng, p, val)
implicit none
    integer :: ng, sy, sz
    integer, dimension(3) :: sxyz
    real, dimension(1-ng:sxyz(1)-ng, 1:sxyz(2), 1:sxyz(3)) :: p
    real :: val

    sy = sxyz(2)
    sz = sxyz(3)
    p(1-ng:0, 1:sy, 1:sz) = val
end subroutine bc_val_3d_l

subroutine bc_val_3d_r(sxyz, ng, p, val)
implicit none
    integer :: ng, sy, sz, nx
    integer, dimension(3) :: sxyz
    real, dimension(1-ng:sxyz(1)-ng, 1:sxyz(2), 1:sxyz(3)) :: p
    real :: val

    nx = sxyz(1) - 2*ng
    sy = sxyz(2)
    sz = sxyz(3)
    p(nx+1:nx+ng, 1:sy, 1:sz) = val
end subroutine bc_val_3d_r

subroutine bc_der_3d_l(sxyz, ng, p)
implicit none
    integer :: ng, sy, sz
    integer, dimension(3) :: sxyz
    real, dimension(1-ng:sxyz(1)-ng, 1:sxyz(2), 1:sxyz(3)) :: p

    sy = sxyz(2)
    sz = sxyz(3)
    p(1-ng:0, 1:sy, 1:sz) = p(1:1, 1:sy, 1:sz);
end subroutine bc_der_3d_l

subroutine bc_der_3d_r(sxyz, ng, p)
implicit none
    integer :: ng, sy, sz, nx
    integer, dimension(3) :: sxyz
    real, dimension(1-ng:sxyz(1)-ng, 1:sxyz(2), 1:sxyz(3)) :: p

    nx = sxyz(1) - 2*ng
    sy = sxyz(2)
    sz = sxyz(3)
    p(nx+1:nx+ng, 1:sy, 1:sz) = p(nx:nx, 1:sy, 1:sz)
end subroutine bc_der_3d_r

!! 3D's guard cell of X direction is a cuboid (ng:Y:Z), Y-Z cross section
subroutine packgc_3d_l(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sy, sz
    integer, dimension(3) :: sxyz
    real, dimension(1-ng:sxyz(1)-ng, 1:sxyz(2), 1:sxyz(3)) :: p
    real, dimension(1:ng,            1:sxyz(2), 1:sxyz(3)) :: gdata

    sy = sxyz(2)
    sz = sxyz(3)
    gdata(1:ng, 1:sy, 1:sz) = p(1:ng, 1:sy, 1:sz)
end subroutine packgc_3d_l

subroutine unpackgc_3d_l(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sy, sz
    integer, dimension(3) :: sxyz
    real, dimension(1-ng:sxyz(1)-ng, 1:sxyz(2), 1:sxyz(3)) :: p
    real, dimension(1:ng,            1:sxyz(2), 1:sxyz(3)) :: gdata

    sy = sxyz(2)
    sz = sxyz(3)
    p(1-ng:0, 1:sy, 1:sz) = gdata(1:ng, 1:sy, 1:sz)
end subroutine unpackgc_3d_l

subroutine packgc_3d_r(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sy, nx, sz
    integer, dimension(3) :: sxyz
    real, dimension(1-ng:sxyz(1)-ng, 1:sxyz(2), 1:sxyz(3)) :: p
    real, dimension(1:ng,            1:sxyz(2), 1:sxyz(3)) :: gdata

    nx = sxyz(1) - 2*ng
    sy = sxyz(2)
    sz = sxyz(3)
    gdata(1:ng, 1:sy, 1:sz) = p(nx-ng+1:nx, 1:sy, 1:sz)
end subroutine packgc_3d_r

subroutine unpackgc_3d_r(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sy, sz, nx
    integer, dimension(3) :: sxyz
    real, dimension(1-ng:sxyz(1)-ng, 1:sxyz(2), 1:sxyz(3)) :: p
    real, dimension(1:ng,            1:sxyz(2), 1:sxyz(3)) :: gdata

    nx = sxyz(1) - 2*ng
    sy = sxyz(2)
    sz = sxyz(3)
    p(nx+1:nx+ng, 1:sy, 1:sz) = gdata(1:ng, 1:sy, 1:sz)
end subroutine unpackgc_3d_r

!!
!! Y direction, X-Z surface
!!

subroutine bc_val_3d_f(sxyz, ng, p, val)
implicit none
    integer :: ng, sx, sz
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1-ng:sxyz(2)-ng, 1:sxyz(3)) :: p
    real :: val

    sx = sxyz(1)
    sz = sxyz(3)
    p(1:sx, 1-ng:0, 1:sz) = val
end subroutine bc_val_3d_f

subroutine bc_val_3d_b(sxyz, ng, p, val)
implicit none
    integer :: ng, sx, sz, ny
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1-ng:sxyz(2)-ng, 1:sxyz(3)) :: p
    real :: val

    ny = sxyz(2) - 2*ng
    sx = sxyz(1)
    sz = sxyz(3)
    p(1:sx, ny+1:ny+ng, 1:sz) = val
end subroutine bc_val_3d_b

subroutine bc_der_3d_f(sxyz, ng, p)
implicit none
    integer :: ng, sx, sz
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1-ng:sxyz(2)-ng, 1:sxyz(3)) :: p
    real :: val

    sx = sxyz(1)
    sz = sxyz(3)
    p(1:sx, 1-ng:0, 1:sz) = p(1:sx, 1:1, 1:sz)
end subroutine bc_der_3d_f

subroutine bc_der_3d_b(sxyz, ng, p)
implicit none
    integer :: ng, sx, sz, ny
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1-ng:sxyz(2)-ng, 1:sxyz(3)) :: p
    real :: val

    ny = sxyz(2) - 2*ng
    sx = sxyz(1)
    sz = sxyz(3)
    p(1:sx, ny+1:ny+ng, 1:sz) = p(1:sx, ny:ny, 1:sz)
end subroutine bc_der_3d_b

subroutine packgc_3d_f(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sx, sz
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1-ng:sxyz(2)-ng, 1:sxyz(3)) :: p
    real, dimension(1:sxyz(1), 1:ng,            1:sxyz(3)) :: gdata

    sx = sxyz(1)
    sz = sxyz(3)
    gdata(1:sx, 1:ng, 1:sz) = p(1:sx, 1:ng, 1:sz)
end subroutine packgc_3d_f

subroutine unpackgc_3d_b(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sx, sz, ny
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1-ng:sxyz(2)-ng, 1:sxyz(3)) :: p
    real, dimension(1:sxyz(1), 1:ng,            1:sxyz(3)) :: gdata

    ny = sxyz(2) - 2*ng
    sx = sxyz(1)
    sz = sxyz(3)
    p(1:sx, ny+1:ny+ng, 1:sz) = gdata(1:sx, 1:ng, 1:sz)
end subroutine unpackgc_3d_b

subroutine packgc_3d_b(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sx, sz, ny
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1-ng:sxyz(2)-ng, 1:sxyz(3)) :: p
    real, dimension(1:sxyz(1), 1:ng,            1:sxyz(3)) :: gdata

    ny = sxyz(2) - 2*ng
    sx = sxyz(1)
    sz = sxyz(3)
    gdata(1:sx, 1:ng, 1:sz) = p(1:sx, ny-ng+1::ny, 1:sz)
end subroutine packgc_3d_b

subroutine unpackgc_3d_f(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sx, sz
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1-ng:sxyz(2)-ng, 1:sxyz(3)) :: p
    real, dimension(1:sxyz(1), 1:ng,            1:sxyz(3)) :: gdata

    sx = sxyz(1)
    sz = sxyz(3)
    p(1:sx, 1-ng:0, 1:sz) = gdata(1:sx, 1:ng, 1:sz)
end subroutine unpackgc_3d_f

!!
!! Z direction, X-Y surface
!!
subroutine bc_val_3d_d(sxyz, ng, p, val)
implicit none
    integer :: ng, sx, sy
    integer, dimension(3) :: sxyz
    real, dimension( 1:sxyz(1), 1:sxyz(2), 1-ng:sxyz(3)-ng ) :: p
    real :: val

    sx = sxyz(1)
    sy = sxyz(2)
    p(1:sx, 1:sy, 1-ng:0) = val
end subroutine bc_val_3d_d

subroutine bc_val_3d_u(sxyz, ng, p, val)
implicit none
    integer :: ng, sx, sy, nz
    integer, dimension(3) :: sxyz
    real, dimension( 1:sxyz(1), 1:sxyz(2), 1-ng:sxyz(3)-ng ) :: p
    real :: val

    sx = sxyz(1)
    sy = sxyz(2)
    nz = sxyz(3) - 2*ng
    p(1:sx, 1:sy, nz+1:nz+ng) = val
end subroutine bc_val_3d_u

subroutine bc_der_3d_d(sxyz, ng, p)
implicit none
    integer :: ng, sx, sy
    integer, dimension(3) :: sxyz
    real, dimension( 1:sxyz(1), 1:sxyz(2), 1-ng:sxyz(3)-ng ) :: p

    sx = sxyz(1)
    sy = sxyz(2)
    p(1:sx, 1:sy, 1-ng:0) = p(1:sx, 1:sy, 1:1)
end subroutine bc_der_3d_d

subroutine bc_der_3d_u(sxyz, ng, p)
implicit none
    integer :: ng, sx, sy, nz
    integer, dimension(3) :: sxyz
    real, dimension( 1:sxyz(1), 1:sxyz(2), 1-ng:sxyz(3)-ng ) :: p

    sx = sxyz(1)
    sy = sxyz(2)
    nz = sxyz(3) - 2*ng
    p(1:sx, 1:sy, nz+1:nz+ng) = p(1:sx, 1:sy, nz:nz)
end subroutine bc_der_3d_u

subroutine packgc_3d_d(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sx, sy
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1:sxyz(2), 1-ng:sxyz(3)-ng ) :: p
    real, dimension(1:sxyz(1), 1:sxyz(2),            1:ng ) :: gdata

    sx = sxyz(1)
    sy = sxyz(2)
    gdata(1:sx, 1:sy, 1:ng) = p(1:sx, 1:sy, 1:ng)
end subroutine packgc_3d_d

subroutine unpackgc_3d_d(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sx, sy
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1:sxyz(2), 1-ng:sxyz(3)-ng ) :: p
    real, dimension(1:sxyz(1), 1:sxyz(2),            1:ng ) :: gdata

    sx = sxyz(1)
    sy = sxyz(2)
    p(1:sx, 1:sy, 1-ng:0) = gdata(1:sx, 1:sy, 1:ng)
end subroutine unpackgc_3d_d

subroutine packgc_3d_u(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sx, sy, nz
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1:sxyz(2), 1-ng:sxyz(3)-ng ) :: p
    real, dimension(1:sxyz(1), 1:sxyz(2),            1:ng ) :: gdata

    nz = sxyz(3) - 2*ng
    sx = sxyz(1)
    sy = sxyz(2)
    gdata(1:sx, 1:sy, 1:ng) = p(1:sx, 1:sy, nz-ng+1:nz)
end subroutine packgc_3d_u

subroutine unpackgc_3d_u(sxyz, ng, p, gdata)
implicit none
    integer :: ng, sx, sy, nz
    integer, dimension(3) :: sxyz
    real, dimension(1:sxyz(1), 1:sxyz(2), 1-ng:sxyz(3)-ng ) :: p
    real, dimension(1:sxyz(1), 1:sxyz(2),            1:ng ) :: gdata

    nz = sxyz(3) - 2*ng
    sx = sxyz(1)
    sy = sxyz(2)
    p(1:sx, 1:sy, nz+1:nz+ng) = gdata(1:sx, 1:sy, 1:ng)
end subroutine unpackgc_3d_u

!! pull out the inner data of one grid
subroutine pack_3d(nxyz, ng, p, gdata)
implicit none
    integer :: ng, nx, ny, nz
    integer, dimension(3) :: nxyz
    real, dimension( 1-ng:nxyz(1)+ng, 1-ng:nxyz(2)+ng, 1-ng:nxyz(3)+ng ) :: p
    real, dimension( 1:nxyz(1), 1:nxyz(2), 1:nxyz(3) ) :: gdata

    nx = nxyz(1)
    ny = nxyz(2)
    nz = nxyz(3)
    gdata(1:nx, 1:ny, 1:nz) = p(1:nx, 1:ny, 1:nz)
end subroutine pack_3d

subroutine update_soln_3d(nxyz, ng, soln, delta)
implicit none
    integer, dimension(3) :: nxyz
    integer ::  ng, i, ix, j, jy, k, kz
    real, dimension( 1-ng:nxyz(1)+ng, 1-ng:nxyz(2)+ng, 1-ng:nxyz(3)+ng ) :: soln, delta

    ix = nxyz(1)
    jy = nxyz(2)
    kz = nxyz(3)
    do k=1, kz
    do j=1, jy
    do i=1, ix
        soln(i,j,k) = soln(i,j,k) + delta(i,j,k)
    end do
    end do
    end do
end subroutine update_soln_3d