!======================================================================
!======================================================================
! Utility Program for Listing, Deleting, Relabeling Labeled Solutions
!                        in AUTO97 Data Files
!======================================================================
!======================================================================
!
PROGRAM AUTLAB
  USE UTILITY

  IMPLICIT NONE
  INTEGER, ALLOCATABLE :: LBR(:),LPT(:),LTY(:),LLB(:),LNL(:)
  INTEGER, ALLOCATABLE :: LFR(:),LTO(:),LBI(:)
  INTEGER, ALLOCATABLE :: BBR(:),BNB(:)
  TYPE(LABEL), POINTER :: LIST
  CHARACTER(2) CMD
  INTEGER NBR,NLB,NL,ILB

  OPEN(27,FILE='fort.27',STATUS='old',ACCESS='sequential')
  OPEN(28,FILE='fort.28',STATUS='old',ACCESS='sequential')

  NULLIFY(LIST)
  CALL RDFILE8(LIST,NLB,28)
  ALLOCATE(LBR(NLB),LPT(NLB),LTY(NLB),LLB(NLB),LNL(NLB),&
       LFR(NLB),LTO(NLB),LBI(NLB))
  CALL GTFILE8(LIST,NLB,LBR,LPT,LTY,LLB,LNL)
  CALL RDFILE7(LIST,NBR,NLB,LBI,27)
  ALLOCATE(BBR(NBR),BNB(NBR))
  CALL GTFILE7(LIST,NBR,BBR,BNB)
  CALL SYSTEM('clear')
  DO
     CALL RDCMD(NLB,CMD,NL,LFR,LTO)
     SELECT CASE(TRIM(CMD))
     CASE('H')
        CALL SYSTEM('clear')
        CALL HELP
     CASE('L')
        CALL LISTLB(NLB,LBR,LPT,LTY,LLB,LNL,NL,LFR,LTO,.TRUE.)
     CASE('LB')
        CALL LISTBR(NBR,BBR,BNB,NLB,LBI,LNL,NL,LFR,LTO)
     CASE('D')
        CALL DELETE('label',NLB,LLB,LNL,NL,LFR,LTO)
     CASE('DB')
        CALL DELETE('branch',NBR,BBR,BNB,NL,LFR,LTO)
        ! delete all labels on the branch
        DO ILB=1,NLB
           IF(LBI(ILB)/=0)THEN
              IF(BNB(LBI(ILB))==0)THEN
                 LNL(ILB)=0
              ENDIF
           ENDIF
        ENDDO
     CASE('R')
        CALL RELABEL('label',NLB,LLB,LNL,NL,LFR,LTO)
     CASE('RB')
        CALL RELABEL('branch',NBR,BBR,BNB,NL,LFR,LTO)
        DO ILB=1,NLB ! sync solutions
           IF(LBI(ILB)/=0)THEN
              LBR(ILB)=BNB(LBI(ILB))
           ENDIF
        ENDDO
     CASE('W')
        IF(CHCKLB(NLB,LNL))THEN
           WRITE(6,"(/,' Rewriting files ... ')")
           OPEN(37,FILE='fort.37',STATUS='unknown',ACCESS='sequential')
           OPEN(38,FILE='fort.38',STATUS='unknown',ACCESS='sequential')
           CALL WRFILE7(NBR,BNB,NLB,LLB,LNL)
           CALL WRFILE8(LBR,LNL)
           STOP
        ENDIF
     CASE('Q')
        WRITE(6,"(' Relabeling discontinued. Recover original files')")
        STOP
     END SELECT
  ENDDO

CONTAINS

!--------- -----
SUBROUTINE RDCMD(NLB,CMD,NL,LFR,LTO)

  INTEGER NLB,NL
  INTEGER LFR(NLB),LTO(NLB)
  CHARACTER(1) CMD1
  CHARACTER(2) CMD

  INTEGER, SAVE :: IFIRST
  CHARACTER(1) CHR
  CHARACTER(80) LINE
  INTEGER IP,i

  DO
     IF(IFIRST.NE.1)THEN
        WRITE(6,'(/A)',ADVANCE="no")' Enter Command ( h for Help) : '
        IFIRST=1
     ELSE
        WRITE(6,'(/A)',ADVANCE="no")' Enter Command : '
     ENDIF

     LINE=' '
     READ(5,'(A80)')LINE
     IP=1
     CMD1=GETCHR(LINE,IP)
     i=INDEX('hldrwq',CMD1)
     IF(i>=1)THEN
        CMD1='HLDRWQ'(i:i)
        EXIT
     ENDIF
     WRITE(6,'(/A)',ADVANCE="no")' Invalid Command '
     IFIRST=0
  ENDDO
  CMD=CMD1
  IF(CMD1=='L'.OR.CMD1=='D'.OR.CMD1=='R')THEN
     NL=0
     DO
        CHR=GETCHR(LINE,IP)
        IF(CHR=='b'.OR.CHR=='B')THEN
           CMD(2:2)='B'
           CYCLE
        ENDIF
        IF(INDEX('0123456789',CHR)==0)EXIT
        IP=IP-1
        NL=NL+1
        LFR(NL)=GETNUM(LINE,IP)
        CHR=GETCHR(LINE,IP)
        IF(CHR=='-')THEN
           LTO(NL)=GETNUM(LINE,IP)
        ELSE
           IP=IP-1
           LTO(NL)=LFR(NL)
        ENDIF
     ENDDO
  ENDIF

END SUBROUTINE RDCMD

!----------- -------- ------
CHARACTER(1) FUNCTION GETCHR(LINE,IP)
  !
  ! Gets next character from LINE that is not a ' ' or ','
  !
  CHARACTER(80), INTENT(IN) :: LINE
  INTEGER, INTENT(INOUT) :: IP

  INTEGER I

  I=VERIFY(LINE(IP:),' ,')
  IF(I>0)THEN
     IP=IP+I
     GETCHR=LINE(IP-1:IP-1)
  ELSE
     IP=81
     GETCHR=' '
  ENDIF

END FUNCTION GETCHR

!------ -------- ------
INTEGER FUNCTION GETNUM(LINE,IP)

  CHARACTER(80), INTENT(IN) :: LINE
  INTEGER, INTENT(INOUT) :: IP

  INTEGER I

  GETNUM=0
  I=SCAN(LINE(IP+1:),'-, ')
  IF(I>0)THEN
     READ(LINE(IP:IP+I-1),*)GETNUM
     IP=IP+I
  ELSE
     READ(LINE(IP:),*)GETNUM
     IP=81
  ENDIF

END FUNCTION GETNUM

!--------- ----
SUBROUTINE HELP

  WRITE(6,"(/A//A/A/A/A/A/A/A)") &
       ' Available commands : ', &
       '   l  :  list labels', &
       '   lb :  list branch numbers', &
       '   d  :  delete labels', &
       '   db :  delete branches', &
       '   r  :  relabel', &
       '   rb :  relabel branches', &
       '   w  :  rewrite files', &
       '   q  :  quit ', &
       '   h  :  help '
  WRITE(6,"(//A/A//A/A/A/A)") &
       ' The l, d, r, lb, db, and rb commands can be followed on the ', &
       ' same line by a list of labels, for example, ', &
       ' l 13           (list label 13)', &
       ' d 7 13         (delete labels 7 and 13)', &
       ' r 1 13 6-9     (relabel 1, 13, and 6 to 9)', &
       ' rb 1, 3-5, 7-9 (relabel branches 1, 3 to 5, and 7 to 9)'
  WRITE(6,"(//A//A/A/A/A/)") &
       ' If a list is not specified then the actions are', &
       ' l              (list all labels)', &
       ' d              (delete/confirm all labels)', &
       ' r              (automatic relabeling)', &
       ' rb             (automatic relabeling of branches)'

END SUBROUTINE HELP

!--------- ------
SUBROUTINE DELETE(S,NLB,LLB,LNL,NL,LFR,LTO)

  USE UTILITY
  CHARACTER(*) S
  INTEGER NLB,LLB(NLB),LNL(NLB),LFR(NLB),LTO(NLB),NL
  CHARACTER(1) CH1
  INTEGER I

  IF(NL.EQ.0)THEN
     DO I=1,NLB
        IF(LNL(I).GT.0)THEN
           WRITE(6,"(A,A,I5,A)",ADVANCE="no") &
                ' Delete (old) ',TRIM(S),LLB(I),' ? : '
           READ(5,'(A1)')CH1
           IF(CH1.EQ.'y'.OR.CH1.EQ.'Y')THEN
              LNL(I)=0
           ENDIF
        ENDIF
     ENDDO
  ELSE
     DO I=1,NLB
        IF(INLIST(NLB,LLB(I),NL,LFR,LTO))THEN
           LNL(I)=0
        ENDIF
     ENDDO
  ENDIF

END SUBROUTINE DELETE

!------ -------- ------
LOGICAL FUNCTION CHCKLB(NLB,LNL)

  INTEGER, INTENT(IN) :: NLB, LNL(NLB)
  INTEGER I, J

  CHCKLB=.TRUE.
  DO I=1,NLB
     IF(LNL(I).GT.0)THEN
        DO J=1,NLB
           IF(J.NE.I .AND. LNL(J).GT.0)THEN
              IF(LNL(J).EQ.LNL(I))THEN
                 WRITE(6,101)
                 CHCKLB=.FALSE.
                 RETURN
              ENDIF
           ENDIF
        ENDDO
     ENDIF
  ENDDO

101 FORMAT(' Duplicate labels: File not written')

  RETURN
END FUNCTION CHCKLB

! ---------- ------
  SUBROUTINE LISTBR(NBR,BBR,BNB,NLB,LBI,LNL,NB,BFR,BTO)
!
    USE UTILITY
    INTEGER NBR,BBR(NBR),BNB(NBR),NLB,LBI(NLB),LNL(NLB)
    INTEGER NB,BFR(NBR),BTO(NBR)
    LOGICAL FIRST
    INTEGER I, J
!
    IF(NBR.EQ.0)THEN
       WRITE(6,"(/,' Empty bifurcation diagram file')")
       RETURN
    ENDIF

    FIRST=.TRUE.
    J=1
    DO I=1,NBR
       IF(NB==0 .OR. INLIST(NLB,I,NB,BFR,BTO))THEN
          IF(FIRST)THEN
             WRITE(6,"(/,'  BR NEW  Labels in branch')")
             FIRST=.FALSE.
          ENDIF
          WRITE(6,"(I4)",ADVANCE="no")BBR(I)
          IF(I==BNB(I).OR.BNB(I)/=0)THEN
             WRITE(6,"(I4)",ADVANCE="no")BNB(I)
             IF(J<=NLB)THEN
                IF(LBI(J)==I)THEN
                   WRITE(6,"(I5)",ADVANCE="no")LNL(J)
                   J=J+1
                   DO WHILE(LBI(J)==I)
                      WRITE(6,"(A,I5)",ADVANCE="no")", ",LNL(J)
                      J=J+1
                   ENDDO
                ENDIF
             ENDIF
             WRITE(6,"()")
          ELSE
             WRITE(6,"('        DELETED')")
          ENDIF
       ENDIF
    ENDDO

  END SUBROUTINE LISTBR
!======================================================================

END PROGRAM AUTLAB
!======================================================================


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     This routine is the driver for computing the addresses and weights
!     for interpolating between two grids on a sphere.
!
!-----------------------------------------------------------------------
!
!     CVS:$Id: scrip.f,v 1.6 2001/08/21 21:06:44 pwjones Exp $
!
!     Copyright (c) 1997, 1998 the Regents of the University of
!       California.
!
!     This software and ancillary information (herein called software)
!     called SCRIP is made available under the terms described here.
!     The software has been approved for release with associated
!     LA-CC Number 98-45.
!
!     Unless otherwise indicated, this software has been authored
!     by an employee or employees of the University of California,
!     operator of the Los Alamos National Laboratory under Contract
!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
!     Government has rights to use, reproduce, and distribute this
!     software.  The public may copy and use this software without
!     charge, provided that this Notice and any statement of authorship
!     are reproduced on all copies.  Neither the Government nor the
!     University makes any warranty, express or implied, or assumes
!     any liability or responsibility for the use of this software.
!
!     If software is modified to produce derivative works, such modified
!     software should be clearly marked, so as not to confuse it with
!     the version available from Los Alamos National Laboratory.
!
!***********************************************************************

      program scrip

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

      use kinds_mod                  ! module defining data types
      use constants                  ! module for common constants
      use iounits                    ! I/O unit manager
      use timers                     ! CPU timers
      use grids                      ! module with grid information
      use remap_vars                 ! common remapping variables
      use remap_conservative         ! routines for conservative remap
      use remap_distance_weight      ! routines for dist-weight remap
      use remap_bilinear             ! routines for bilinear interp
      use remap_bicubic              ! routines for bicubic  interp
      use remap_write                ! routines for remap output

      implicit none

!-----------------------------------------------------------------------
!
!     input namelist variables
!
!-----------------------------------------------------------------------

      character (char_len) ::
     &           grid1_file,   ! filename of grid file containing grid1
     &           grid2_file,   ! filename of grid file containing grid2
     &           interp_file1, ! filename for output remap data (map1)
     &           interp_file2, ! filename for output remap data (map2)
     &           map1_name,    ! name for mapping from grid1 to grid2
     &           map2_name,    ! name for mapping from grid2 to grid1
     &           map_method,   ! choice for mapping method
     &           normalize_opt,! option for normalizing weights
     &           output_opt    ! option for output conventions

      integer (kind=int_kind) ::
     &           nmap          ! number of mappings to compute (1 or 2)

      namelist /remap_inputs/ grid1_file, grid2_file,
     &                        interp_file1, interp_file2,
     &                        map1_name, map2_name, num_maps,
     &                        luse_grid1_area, luse_grid2_area,
     &                        map_method, normalize_opt, output_opt,
     &                        restrict_type, num_srch_bins

!-----------------------------------------------------------------------
!
!     local variables
!
!-----------------------------------------------------------------------

      integer (kind=int_kind) :: n,     ! dummy counter
     &                           iunit  ! unit number for namelist file

!-----------------------------------------------------------------------
!
!     initialize timers
!
!-----------------------------------------------------------------------

      call timers_init
      do n=1,max_timers
        call timer_clear(n)
      end do

!-----------------------------------------------------------------------
!
!     read input namelist
!
!-----------------------------------------------------------------------

      grid1_file    = 'unknown'
      grid2_file    = 'unknown'
      interp_file1  = 'unknown'
      interp_file2  = 'unknown'
      map1_name     = 'unknown'
      map2_name     = 'unknown'
      luse_grid1_area = .false.
      luse_grid2_area = .false.
      num_maps      = 2
      map_type      = 1
      normalize_opt = 'fracarea'
      output_opt    = 'scrip'
      restrict_type = 'latitude'
      num_srch_bins = 900

      call get_unit(iunit)
      open(iunit, file='scrip_in', status='old', form='formatted')
      read(iunit, nml=remap_inputs)
      call release_unit(iunit)

      select case(map_method)
      case ('conservative')
        map_type = map_type_conserv
        luse_grid_centers = .false.
      case ('bilinear')
        map_type = map_type_bilinear
        luse_grid_centers = .true.
      case ('bicubic')
        map_type = map_type_bicubic
        luse_grid_centers = .true.
      case ('distwgt')
        map_type = map_type_distwgt
        luse_grid_centers = .true.
      case default
        stop 'unknown mapping method'
      end select

      select case(normalize_opt(1:4))
      case ('none')
        norm_opt = norm_opt_none
      case ('frac')
        norm_opt = norm_opt_frcarea
      case ('dest')
        norm_opt = norm_opt_dstarea
      case default
        stop 'unknown normalization option'
      end select

!-----------------------------------------------------------------------
!
!     initialize grid information for both grids
!
!-----------------------------------------------------------------------

      call grid_init(grid1_file, grid2_file)

      write(stdout, *) ' Computing remappings between: ',grid1_name
      write(stdout, *) '                          and  ',grid2_name

!-----------------------------------------------------------------------
!
!     initialize some remapping variables.
!
!-----------------------------------------------------------------------

      call init_remap_vars

!-----------------------------------------------------------------------
!
!     call appropriate interpolation setup routine based on type of
!     remapping requested.
!
!-----------------------------------------------------------------------

      select case(map_type)
      case(map_type_conserv)
        call remap_conserv
      case(map_type_bilinear)
        call remap_bilin
      case(map_type_distwgt)
        call remap_distwgt
      case(map_type_bicubic)
        call remap_bicub
      case default
        stop 'Invalid Map Type'
      end select

!-----------------------------------------------------------------------
!
!     reduce size of remapping arrays and then write remapping info
!     to a file.
!
!-----------------------------------------------------------------------

      if (num_links_map1 /= max_links_map1) then
        call resize_remap_vars(1, num_links_map1-max_links_map1)
      endif
      if ((num_maps > 1) .and. (num_links_map2 /= max_links_map2)) then
        call resize_remap_vars(2, num_links_map2-max_links_map2)
      endif

      call write_remap(map1_name, map2_name,
     &                 interp_file1, interp_file2, output_opt)

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

      end program scrip

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

ext: f


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

C=======================================================================
      SUBROUTINE PRNPS (OPTION, NOUT, NUMNPS, LISNPS, LNPSNL,
     &   IDNPS, NNNPS, IXNNPS, LTNNPS, FACNPS, NSNAME, MAPND)
C=======================================================================

C   --*** PRNPS *** (BLOT) Display database node set
C   --   Written by Amy Gilkey - revised 01/18/88
C   --
C   --PRNPS displays the node sets.
C   --
C   --Parameters:
C   --   OPTION - IN - '*' to print all, else print options:
C   --      ' ' - set summary (number of nodes, etc)
C   --      'N' - nodes in set
C   --      'F' - distribution factors for set
C   --   NOUT - IN - the output file, <=0 for standard
C   --   NUMNPS - IN - the number of node sets
C   --   LISNPS - IN - the indices of the selected node sets
C   --   LNPSNL - IN - the number of nodes for all sets
C   --   IDNPS - IN - the node set ID for each set
C   --   NNNPS - IN - the number of nodes for each set
C   --   IXNNPS - IN - the index of the first node for each set
C   --   LTNNPS - IN - the nodes for all sets
C   --   FACNPS - IN - the distribution factors for all sets

      CHARACTER*(*) OPTION
      INTEGER LISNPS(0:*)
      INTEGER IDNPS(*)
      INTEGER NNNPS(*)
      INTEGER IXNNPS(*)
      INTEGER LTNNPS(*)
      REAL FACNPS(*)
      CHARACTER*(*) NSNAME(*)
      INTEGER MAPND(*)

      LOGICAL ISABRT
      LOGICAL DONOD, DOFAC
      CHARACTER*20 STRA, STRB

      DONOD = ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'N') .GT. 0))
      DOFAC = ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'F') .GT. 0))

      IF (NOUT .GT. 0) THEN
         IF (DONOD .AND. DOFAC) THEN
            WRITE (NOUT, 10020) 'Node List and Distribution Factors'
         ELSE IF (DONOD) THEN
            WRITE (NOUT, 10020) 'Node List'
         ELSE IF (DOFAC) THEN
            WRITE (NOUT, 10020) 'Distribution Factors'
         ELSE
            WRITE (NOUT, 10020)
         END IF
      END IF

      IF (NOUT .GT. 0) THEN
         WRITE (NOUT, *)
      ELSE
         WRITE (*, *)
      END IF

      WRITE (STRA, 10000, IOSTAT=IDUM) NUMNPS
10000  FORMAT ('(#', I4, ')')
      CALL PCKSTR (1, STRA)
      LSTRA = LENSTR (STRA)
      WRITE (STRB, 10010, IOSTAT=IDUM) LNPSNL
10010  FORMAT ('(index=', I9, ')')
      CALL PCKSTR (1, STRB)
      LSTRB = LENSTR (STRB)

      DO 100 IX = 1, LISNPS(0)
         IF (ISABRT ()) RETURN
         INPS = LISNPS(IX)
         WRITE (STRA, 10000, IOSTAT=IDUM) INPS
         CALL PCKSTR (1, STRA)
         WRITE (STRB, 10010, IOSTAT=IDUM) IXNNPS(INPS)
         CALL PCKSTR (1, STRB)
         IF (NOUT .GT. 0) THEN
            WRITE (NOUT, 10030, IOSTAT=IDUM)
     &         IDNPS(INPS), STRA(:LSTRA),
     &         NNNPS(INPS), STRB(:LSTRB),
     $         NSNAME(INPS)(:LENSTR(NSNAME(INPS)))
         ELSE
            WRITE (*, 10030, IOSTAT=IDUM)
     &         IDNPS(INPS), STRA(:LSTRA),
     &         NNNPS(INPS), STRB(:LSTRB),
     $         NSNAME(INPS)(:LENSTR(NSNAME(INPS)))
         END IF

         IF (DONOD .AND. (NNNPS(INPS) .GT. 0)) THEN
            IS = IXNNPS(INPS)
            IE = IS + NNNPS(INPS) - 1
            IF (NOUT .GT. 0) THEN
               WRITE (NOUT, 10040, IOSTAT=IDUM)
     &          (MAPND(LTNNPS(I)), I=IS,IE)
            ELSE
               WRITE (*, 10040, IOSTAT=IDUM)
     &          (MAPND(LTNNPS(I)), I=IS,IE)
            END IF
         END IF

         IF (DOFAC .AND. (NNNPS(INPS) .GT. 0)) THEN
            IS = IXNNPS(INPS)
            IE = IS + NNNPS(INPS) - 1
            IF (NOUT .GT. 0) THEN
               WRITE (NOUT, 10050, IOSTAT=IDUM)
     &            (FACNPS(I), I=IS,IE)
            ELSE
               WRITE (*, 10050, IOSTAT=IDUM)
     &            (FACNPS(I), I=IS,IE)
            END IF
         END IF
  100 CONTINUE

      RETURN

10020  FORMAT (/, 1X, 'Node Sets (Global Node Ids)', :, ' - ', A)
10030  FORMAT (1X, 'Set', I9, 1X, A, ':',
     &   I9, ' nodes', 1X, A,' name = "',A,'"')
10040  FORMAT ((1X, 8I11))
10050  FORMAT ((1X, 6 (2X, E11.4)))
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


program drive

  use regression
  implicit none

  integer, parameter :: nr = 1, np = 1000
  real, dimension(nr) :: x, y, a, b, eps
  type (reg_type), dimension(nr) :: r
  integer :: i, j
  real :: avg

  do
     write (*,*) 'x, y'
     read (*,*,err=100) x(1), y(1)
     call regress (r, x, y)
     write (*,*) r%a, r%b
  end do

100 continue

  write (*,*) 'x, y'
  read (*,*,err=100) x(1), y(1)
  call regress (r, x, y, a, b)

  write (*,*) 'Final values:'
  write (*,*) 'a(1) = ',a(1),'     b(1) = ',b(1)
!  write (*,*) 'a(2) = ',a(2),'     b(2) = ',b(2)

  x(1) = 10.
!  x(2) = 10.

  write (*,*)
  write (*,*) 'x(1) = ',x(1)!,' x(2) = ',x(2)
  write (*,*) 'y= ',yp(r, x)

  y(1) = 36.
!  y(2) = -3.

  write (*,*)
  write (*,*) 'y(1) = ',y(1)!,' y(2) = ',y(2)
  write (*,*) 'x= ',xp(r, y)

end program drive


text:
program demo_load
  use stdlib_io_npy, only: load_npy
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), allocatable :: array(:, :, :)

  call load_npy("array.npy", array)

  print '(*(g0,1x))', "shape:", shape(array)
  ! shape: 5 2 6
end program demo_load

max_stars_repo_name: awvwgk/stdlib-docs

max_issues_repo_name: awvwgk/stdlib-docs

max_forks_repo_name: awvwgk/stdlib-docs


text:
C    5 SEP 94    ...AS...    (CRI) MSTP S/W CONVERSION (VICAR PORTING)
C         *******    SUBROUTINE DC   *******
C   10 OCT 88    ...SP...    CHANGED TO IGNORE FORMAT PARAMETER BECAUSE VICAR2
C                            USES ONLY THE FORMAT IN LABEL.
C   11 MAY 84    ...HBD...   CONVERT FOR VAX and to vicar2 format
C   22 OCT 79    ...JAM...   INITIAL RELEASE
C
      INCLUDE 'VICMAIN_FOR'
C
      SUBROUTINE MAIN44
      IMPLICIT NONE
C
C          LOCAL VARIABLES
C  STATUS -- OPEN AND CLOSE STATUS INDICATOR
C  INUNIT -- CONTAINS UNIT NUMBERS FOR INPUT DATASETS
C  NOLINE -- CONTAINS THE NUMBER OF LINES FOR EACH INPUT DATASET
C  NOSAMP -- CONTAINS THE NUMBER OF SAMPLES FOR EACH INPUT DATASET
C  REFIN  -- STARTING BYTE OF EACH PIX IN INPUT BUFFER (INBUF in DCB)
C
      INTEGER NOLINE(10), NOSAMP(10), REFIN(10), INUNIT(10)
      INTEGER STATUS, NUMIPT, NUMOPT, I, INDEX
      CHARACTER*10 LBUF(10)
C
      CHARACTER*25 MSG1

      COMMON /C1/ NOSAMP,NOLINE,REFIN,NUMIPT,NUMOPT,INUNIT

      CALL IFMESSAGE('DC version 5-SEP-94')
      CALL XVEACTION('SA',' ')
C--- Zero out buffers
      MSG1(1:25)='NUMBER OF INPUTS NOT GT 1'
      DO I = 1, 10
         REFIN(I) = 0
         NOLINE(I) = 0
         NOSAMP(I) = 0
      END DO
      CALL XVP('INP',LBUF,NUMIPT)
      CALL XVP('OUT',LBUF,NUMOPT)
      IF (NUMIPT .LE. 1) THEN
         CALL XVMESSAGE(MSG1,' ')
         CALL ABEND
      ENDIF
      INDEX = 1
C
C--- Open input datasets
      DO I = 1,NUMIPT
         REFIN(I) = INDEX
         CALL XVUNIT(INUNIT(I),'INP',I,STATUS,' ')
	 CALL XVOPEN(INUNIT(I),STATUS,'U_FORMAT','HALF',' ')
	 CALL XVGET(INUNIT(I),STATUS,'NL',NOLINE(I),'NS',NOSAMP(I),' ')
         INDEX = NOSAMP(I)+INDEX
      END DO

      CALL DCB

      RETURN
      END

      SUBROUTINE DCB
C
      IMPLICIT NONE
C
      INTEGER*2 INBUF(10000)
      INTEGER*2 OBUF1(50000),OBUF2(50000),OBUF3(50000),OBUF4(50000)
      INTEGER NOLINE(10),NOSAMP(10),REFIN(10),NUMIPT,NUMOPT,INUNIT(10)
      INTEGER OUTUNIT(4),CNT,DEF, STATUS, I, J, LINE, IPIXSIZ
      INTEGER SLO,SSO,NLO,NSO,MAXDN,OSAMP,ISAMP,DN,NLI,NSI
      INTEGER INT,IREC,ISIGN,ISLOPE,XPNTST
      REAL*4 EXPO(10),ACONST,CONST,DELTA,DENOM,DENOM2,SLOPE,SCALE
      REAL*4 X,X2,XY,Y,RMAXDN,RINT,YPOINT,PDN,RDN,TDELTA, XPOINT
      CHARACTER*4 FORMAT(2)
      LOGICAL ISWTCH, NEGA
      LOGICAL XVPTST

      COMMON /C1/ NOSAMP,NOLINE,REFIN,NUMIPT,NUMOPT,INUNIT
      DATA FORMAT/'WORD','BYTE'/
C
C--- Get SIZE parameters
      CALL XVSIZE(SLO,SSO,NLO,NSO,NLI,NSI)
      IREC = SLO

C--- Test for halfword and negative DC picture
      CALL XVGET(INUNIT(1), STATUS, 'PIX_SIZE', IPIXSIZ,' ')  !BYTES PER PIXEL.
      IF ( IPIXSIZ .NE. 1 .AND. IPIXSIZ .NE. 2 )
     .     CALL MABEND('ERROR:INP. IMAGE IS NOT BYTE OR HALF FORMAT')

      ISWTCH = IPIXSIZ .EQ. 2
      ISIGN = 1
      NEGA = XVPTST('NEGATIVE')
      IF (NEGA) ISIGN = -1

C--- Determine MAXDN value
      CALL XVPARM('MAXDN',MAXDN,CNT,DEF,1)
      IF (ISWTCH .AND. DEF .NE. 0) MAXDN = 32767
      RMAXDN = MAXDN
      CONST = ISIGN / 2.0

C--- Get EXPOSURE times
      CALL XVP('EXPOSURE',EXPO,CNT)
      IF (CNT .EQ. 0) THEN
	 CALL XVMESSAGE('DC - Exposure values required.',' ')
	 CALL ABEND
      END IF
      DO I = 1, NUMIPT
	 X = X + EXPO(I)
	 X2 = X2 + EXPO(I) * EXPO(I)
      END DO
      DENOM = NUMIPT * X2 - X * X

C--- Get scaling factor
      CALL XVP('SCALE',SCALE,CNT)
      DENOM2 = DENOM / SCALE
      CALL XVP('X',XPOINT,XPNTST)
      CALL XVP('CONS',ACONST,CNT)

C--- Open output data sets
      DO I = 1, NUMOPT
	 CALL XVUNIT(OUTUNIT(I),'OUT',I,STATUS,' ')
	 CALL XVOPEN(OUTUNIT(I),STATUS,'OP','WRITE',' ')
      END DO
      DO 400 LINE = 1, NLO
         DO J = 1, NUMIPT
	    CALL XVREAD(INUNIT(J),INBUF(REFIN(J)),STATUS,'SAMP',SSO,
     &		       'NSAMPS',NSO,'LINE',IREC,' ')
         END DO
         IREC = IREC + 1
         OSAMP = 0

C--- This part does linear fit (DN = SLOPE * EXPO + INT)
         DO 200 ISAMP = 1, NSO
            OSAMP = OSAMP + 1
            Y = 0.0
            XY = 0.0
            DO I = 1, NUMIPT
               DN = INBUF(REFIN(I)+ISAMP-1)
               Y = Y + DN
               XY = XY + DN * EXPO(I)
            END DO
            RINT = FLOAT(ISIGN) * ((Y * X2 - X * XY) / DENOM)
            INT = RINT + CONST + ACONST
            IF (ISIGN .EQ. -1) INT = INT + 1
            IF (INT .LT. 1) INT = 0
            IF (INT .GT. MAXDN) INT=MAXDN
            OBUF1(OSAMP) = INT

C--- Calculate SLOPE picture
            IF (2 .LE. NUMOPT) THEN
               SLOPE = ((NUMIPT * XY - X * Y) / DENOM2)
               ISLOPE = SLOPE + 0.5
               IF (ISLOPE .LT. 1) ISLOPE = 0
               IF (ISLOPE .GT. MAXDN) ISLOPE = MAXDN
               OBUF2(OSAMP) = ISLOPE
            ELSE
               GOTO 200
            END IF

C--- Find predicted Y value
            IF (XPNTST .EQ. 0) GO TO 250
            IF (3 .LE. NUMOPT) THEN
               YPOINT = SLOPE * XPOINT + RINT
               IF (YPOINT .LT. 1.0) YPOINT=0.0
               IF (YPOINT .GT. RMAXDN) YPOINT = RMAXDN
               OBUF3(OSAMP) = YPOINT + .5
            ELSE
               GOTO 200
            END IF

C--- Find point which deviates the most from fitted line
  250       IF (4 .LE. NUMOPT) THEN
               DELTA = 0.0
               DO I = 1, NUMIPT                  ! Calculate predicted DN (PDN)
                  PDN = SLOPE * EXPO(I) + RINT
                  RDN = INBUF(REFIN(I)+ISAMP-1)  ! RDN is actual DN
                  TDELTA = ABS(PDN-RDN)
                  IF (TDELTA .GT. DELTA) DELTA = ABS(RDN-PDN)
               END DO
               IF (XPNTST .EQ. 0) THEN
		  OBUF3(OSAMP) = DELTA + ACONST
               ELSE
		  OBUF4(OSAMP) = DELTA + ACONST
   	       END IF
            END IF
  200    CONTINUE

C--- Write out buffers to files
         CALL XVWRIT(OUTUNIT(1),OBUF1,STATUS,' ')
         IF (NUMOPT .EQ. 2) THEN
             CALL XVWRIT(OUTUNIT(2),OBUF2,STATUS,' ')
         ELSE IF (NUMOPT .EQ. 3) THEN
             CALL XVWRIT(OUTUNIT(2),OBUF2,STATUS,' ')
             CALL XVWRIT(OUTUNIT(3),OBUF3,STATUS,' ')
         ELSE IF (NUMOPT .EQ. 4) THEN
             CALL XVWRIT(OUTUNIT(2),OBUF2,STATUS,' ')
             CALL XVWRIT(OUTUNIT(3),OBUF3,STATUS,' ')
             CALL XVWRIT(OUTUNIT(4),OBUF4,STATUS,' ')
         END IF

C--- Close all files
  400 CONTINUE
      DO I = 1, NUMIPT
        CALL XVCLOSE(INUNIT(I),STATUS,' ')
      END DO
      DO I = 1, NUMOPT
         CALL XVCLOSE(OUTUNIT(I),STATUS,' ')
      END DO
      RETURN
      END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
program loop
    do i = 1, 4
        do j = 1, i
            write(*, fmt ='(a)', advance='no') "* "
        end do
        print *
    end do
end program loop

max_stars_repo_name: Zeocode/Hello-World

max_issues_repo_name: Zeocode/Hello-World

max_forks_repo_name: Zeocode/Hello-World


text:
program sto_mg_screen_gf

! Screened Green's function calculation
! The results for D->oo are to be compared against [1].

! [1] Doll, J. D., & Reinhardt, W. P. (1972). Many-Body Green’s Functions for
! Finite, Nonuniform Systems: Applications to Closed Shell Atoms. The Journal
! of Chemical Physics, 57(3), 1169–1184. doi:10.1063/1.1678374

use types, only: dp
use sto, only: stoints2, get_basis2, slater_sto_screen, sto_V_screen
use utils, only: assert
use radialscf, only: doscf, kinetic_energy, slater2int22, &
    get_basis, radialC2C, &
    radiallam2lam
use hfprint, only: printall, printlam
use mbpt, only: transform_int2, mbpt2, mbpt3, mbpt4, transform_int22
use gf2, only: find_poles, find_pole_diag, total_energy, plot_poles
use sorting, only: argsort
use scf, only: ijkl2intindex, ijkl2intindex2, create_intindex_sym4
implicit none

integer, allocatable :: nl(:, :), nbfl(:)
real(dp), allocatable :: zl(:, :), focc(:, :)

real(dp), allocatable :: S(:, :, :), T(:, :, :), V(:, :, :), slater(:, :)
integer :: n, Z, m, Nscf, Lmax, ndof
real(dp) :: alpha, Etot, tolE, tolP, Ekin
real(dp), allocatable :: H(:, :, :), P_(:, :, :), C(:, :, :), lam(:, :)

integer, allocatable :: nlist(:), llist(:), mlist(:)
real(dp), allocatable :: int2(:), moint2(:), Ctot(:, :), lamtot(:), Egreen(:)
integer, allocatable :: intindex(:, :, :, :)
integer :: Nelec

real(dp) :: D
integer, allocatable :: idx(:)
integer :: i, it, it_N, u
real(dp) :: Dmin, Dmax, it_a, it_b, t1, t2

Lmax = 1
allocate(nbfl(0:Lmax))
nbfl = [10, 7]
allocate(nl(maxval(nbfl), 0:Lmax), zl(maxval(nbfl), 0:Lmax), focc(3, 0:Lmax))
focc = 0
focc(:3, 0) = [2, 2, 2]
focc(:1, 1) = [6]
nl(:nbfl(0), 0) = [2, 1, 2, 1, 1, 2, 1, 1, 2, 1]
zl(:nbfl(0), 0) = [35.988252_dp, 19.193253_dp, 16.093419_dp, 10.773377_dp, &
                7.853579_dp, 6.023976_dp, 2.905152_dp, 1.694174_dp, &
                0.750076_dp, 0.711410_dp]
nl(:nbfl(1), 1) = [3, 2, 3, 2, 2, 2, 2]
zl(:nbfl(1), 1) = [32.077221_dp, 12.495413_dp, 9.644163_dp, 5.042248_dp, &
                3.070882_dp, 2.086080_dp, 0.711410_dp]

Z = 12
tolE = 1e-10_dp
tolP = 1e-4_dp
alpha = 0.6_dp
Nscf = 100

n = maxval(nbfl)
ndof = sum(nbfl)
Nelec = round(sum(focc))
print *, "total  DOFs =", ndof
allocate(S(n, n, 0:Lmax), T(n, n, 0:Lmax), V(n, n, 0:Lmax))
m = ndof*(ndof+1)/2
allocate(slater(m*(m+1)/2, 0:2*Lmax))
call stoints2(Z, nbfl, nl, zl, S, T, V, slater)
allocate(P_(n, n, 0:Lmax), C(n, n, 0:Lmax), H(n, n, 0:Lmax), lam(n, 0:Lmax))
call get_basis(nbfl, nlist, llist, mlist)
allocate(Ctot(size(nlist), size(nlist)), lamtot(size(nlist)), idx(size(nlist)))
allocate(Egreen(Nelec))
m = size(nlist)**2 * (size(nlist)**2 + 3) / 4
allocate(int2(m), moint2(m))
allocate(intindex(size(nlist), size(nlist), size(nlist), size(nlist)))
call create_intindex_sym4(intindex)

open(newunit=u, file="mg_eigs_D.txt", status="replace")
close(u)
Dmin = 1._dp
Dmax = 1e5
it_N = 100
it_a = Dmin * (Dmin/Dmax)**(1._dp/(it_N-1))
it_b = log(Dmax/Dmin)/(it_N-1)
do it = it_N, 1, -1
    D = it_a*exp(it_b*it)
    print *, "D =", D
    call sto_V_screen(Z, nbfl, nl, zl, V, D)
    call slater_sto_screen(nbfl, nl, zl, slater, D)
    H = T + V
    print *, "SCF cycle:"
    call doscf(nbfl, H, slater, S, focc, Nscf, tolE, tolP, alpha, C, P_, lam, Etot)
    Ekin = kinetic_energy(nbfl, P_, T)
    !call printall(nbfl, nl, zl, lam, C, Ekin, Etot)
    call printlam(nbfl, lam, Ekin, Etot)
    call radiallam2lam(nlist, llist, lam, lamtot)
    idx = argsort(lamtot)
    nlist = nlist(idx)
    llist = llist(idx)
    mlist = mlist(idx)
    lamtot = lamtot(idx)
    call slater2int22(nbfl, nlist, llist, mlist, slater, int2)
    call radialC2C(nlist, llist, mlist, C, Ctot)
    call cpu_time(t1)
    moint2 = transform_int22(int2, intindex, Ctot)
    call cpu_time(t2)
    print *, t2-t1
!    stop "OK"
    print *, "Green's function calculation:"
    print *, "i         lam(i)            dE         lam(i)+dE"
    Egreen = 0
    do i = 1, Nelec
        Egreen(i) = find_pole_diag(i, moint2, intindex, lamtot, Nelec/2, &
            200, 1e-10_dp)
        print "(i4, f15.6, f15.6, f15.6)", i, lamtot(i), Egreen(i)-lamtot(i), &
            Egreen(i)
    end do

    open(newunit=u, file="mg_eigs_D.txt", position="append", status="old")
    write(u, "(100(es23.16, ' '))") D, Etot, lamtot(:Nelec), Egreen(:Nelec)
    close(u)

end do

contains

integer function round(x) result(r)
! Rounds "x" to the nearest integer
real(dp), intent(in) :: x
r = int(x + 0.5_dp)
end function

end program

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
      SUBROUTINE RSSORT(IMIN,NXTLEN,NXTLST)
*
      INCLUDE 'common6.h'
      INTEGER NXTLEN,NXTLST(NMAX)
*
      TMIN = TIMENW(IFIRST)
      IMIN = IFIRST
*
      DO 3 J = IFIRST+1, NTOT
         IF(TIMENW(J).LT.TMIN)THEN
         TMIN = TIMENW(J)
         IMIN = J
         END IF
  3   CONTINUE
*
      DO 5 J = IFIRST, NTOT
         IF(DABS(TIMENW(J)-TMIN).LT.DTK(40)) THEN
            NXTLEN = NXTLEN + 1
            NXTLST(NXTLEN) = J
         END IF
  5   CONTINUE
*
      RETURN
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
       subroutine mixing  !actuall diffusivity

       use geometry_module!, only: mp_, np_, hzp
       use grid_module!, only: mtime,dx,dy,dz
       use signals_module!, only:  bit5k, bit2k
       use viscosity_module!, only:  axmin, aymin, prandtl, azzmin, ss_v, ss, chi_v, chi_h

       use tracer_module!  , only:   kxt,kyt,kzt,t,t_bkg
       include 'array_size.inc'
       real tprime(0:nx+1,0:ny+1,0:nz+1)

          integer i,j,k

          kxt = 0.0d0
          kyt = 0.0d0
          kzt = 0.0d0

!         write(6,*) 'ss'

!$OMP PARALLEL DO

          do k = 1, nz    !-1   ! set to nz-1 so that artifically large horizontal dispersion won't occur
                           ! at seafloor just be cause of large vertical shear there; alternative is to split smag formula
                           ! into horizontal and vertical production terms
                           ! or change bottom boundary so partial slip and not so much vertical shear
                           ! or modulate vertical shear term in a way that boudnary layer doesn't make large bottom
                           ! horiztonal mixing
          do j = 0, ny+1
          do i = 0, nx
                kxt(i,j,k) = bit5k * (( .5d0 * ( ss(i,j,k) + ss(i+1,j,k) ) ) / prandtl + axmin / prandtl)
          enddo;enddo;enddo

!$OMP END PARALLEL DO

!$OMP PARALLEL DO

          do k = 1, nz!-1
          do j = 0, ny
          do i = 0, nx+1
                kyt(i,j,k) = bit5k * (( .5d0 * ( ss(i,j,k) + ss(i,j+1,k) ) ) / prandtl + aymin / prandtl)
          enddo;enddo;enddo

!$OMP END PARALLEL DO

                kzt = 0.0
!$OMP PARALLEL DO

          do k = 1, nz-1
          do j = 0, ny+1
          do i = 0, nx+1
                kzt(i,j,k) = bit2k * ((0.5d0 * ( ss_v(i,j,k) + ss_v(i,j,k+1) ) ) / prandtl
     &                  +  azzmin(k) / prandtl)
           enddo;enddo;enddo

!$OMP END PARALLEL DO

!$OMP PARALLEL DO

            do k = 0, nz+1
            do j = 0, ny+1
            do i = 0, nx+1
               tprime(i,j,k) = t(i,j,k)- t_bkg(k)
            enddo;enddo;enddo

!$OMP END PARALLEL DO

               chi_v = 0; chi_h = 0

!$OMP PARALLEL DO
! Aug 2014 corrected the calculation of chi with /2dz /2dx /2dy since difference is over 2 grid cells
            do k = 1, nz
            do j = 1, ny
            do i = 1, nx
               chi_v(i,j,k) = kzt(i,j,k)* ( (tprime(i,j,k-1) - tprime(i,j,k+1))/ (hzp(k)*2.0d0* dz))**2
            enddo;enddo;enddo
!$OMP END PARALLEL DO
!$OMP PARALLEL DO
            do k = 1, nz
            do j = 1, ny
            do i = 1, nx
               chi_h(i,j,k) = kxt(i,j,k)* ( (tprime(i+1,j,k) - tprime(i-1,j,k))* mp_(i,j)/ (2.0d0*dx) )**2
     &                      + kyt(i,j,k)* ( (tprime(i,j+1,k) - tprime(i,j-1,k))* np_(i,j)/ (2.0d0*dy) )**2
            enddo;enddo;enddo

!$OMP END PARALLEL DO

      end subroutine mixing

max_stars_repo_name: diioriod/PlumeModel-Hotcross

max_issues_repo_name: diioriod/PlumeModel-Hotcross

max_forks_repo_name: diioriod/PlumeModel-Hotcross


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 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.
!
      subroutine mafillsmas(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
     &  xboun,nboun,
     &  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
     &  nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
     &  ad,au,bb,nactdof,icol,jq,irow,neq,nzl,nmethod,
     &  ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,
     &  nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
     &  t0,t1,ithermal,prestr,
     &  iprestr,vold,iperturb,sti,nzs,stx,adb,aub,iexpl,plicon,
     &  nplicon,plkcon,nplkcon,xstiff,npmat_,dtime,
     &  matname,mi,ncmat_,mass,stiffness,buckling,rhsi,intscheme,
     &  physcon,shcon,nshcon,cocon,ncocon,ttime,time,istep,iinc,
     &  coriolis,ibody,xloadold,reltime,veold,springarea,nstate_,
     &  xstateini,xstate,thicke,
     &  integerglob,doubleglob,tieset,istartset,iendset,
     &  ialset,ntie,nasym,pslavsurf,pmastsurf,mortar,clearini,ielprop,
     &  prop,ne0,kscale,iponoel,inoel,network)
!
!     filling the stiffness matrix in spare matrix format (sm)
!     asymmetric contributions
!
      implicit none
!
      logical mass(2),stiffness,buckling,rhsi,coriolis
!
      character*8 lakon(*)
      character*20 sideload(*)
      character*80 matname(*)
      character*81 tieset(3,*)
!
      integer kon(*),nodeboun(*),ndirboun(*),ipompc(*),nodempc(3,*),
     &  nodeforc(2,*),ndirforc(*),nelemload(2,*),icol(*),jq(*),ikmpc(*),
     &  ilmpc(*),ikboun(*),ilboun(*),mi(*),integerglob(*),ist,mpc1,
     &  nactdof(0:mi(2),*),irow(*),istartset(*),iendset(*),
     &  nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),index,
     &  ielorien(mi(3),*),ialset(*),ntie,ne0,nstate_,index1,index2,
     &  ipkon(*),intscheme,ncocon(2,*),nshcon(*),ipobody(2,*),nbody,
     &  ibody(3,*),nk,ne,nboun,nmpc,nforc,nload,neq,nzl,nmethod,
     &  ithermal(2),iprestr,iperturb(*),nzs(3),i,j,k,l,m,idist,jj,
     &  ll,jdof1,jdof2,node1,node2,id,i0,id1,id2,idof1,idof2,nasym,
     &  ntmat_,indexe,nope,norien,iexpl,ncmat_,istep,iinc,mpc2,
     &  nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,ist1,ist2,
     &  mortar,ielprop(*),kscale,iponoel(*),inoel(2,*),network
!
      real*8 co(3,*),xboun(*),coefmpc(*),xforc(*),xload(2,*),p1(3),
     &  p2(3),ad(*),au(*),bodyf(3),bb(*),xloadold(2,*),value,
     &  t0(*),t1(*),prestr(6,mi(1),*),vold(0:mi(2),*),s(60,60),
     &  ff(60),
     &  sti(6,mi(1),*),sm(60,60),stx(6,mi(1),*),adb(*),aub(*),
     &  elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),reltime,
     &  alcon(0:6,ntmat_,*),physcon(*),cocon(0:6,ntmat_,*),
     &  shcon(0:3,ntmat_,*),alzero(*),orab(7,*),xbody(7,*),cgr(4,*),
     &  xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
     &  springarea(2,*),thicke(mi(3),*),clearini(3,9,*),prop(*),
     &  plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
     &  xstiff(27,mi(1),*),veold(0:mi(2),*),doubleglob(*),
     &  om,dtime,ttime,time,pslavsurf(3,*),pmastsurf(6,*)
!
      i0=0
!
!     storing the symmetric matrix in asymmetric format
!
      do i=1,nzs(2)
         au(nzs(3)+i)=au(i)
      enddo
!
!     dynamic applications
!
      if((nmethod.eq.2).or.(nmethod.eq.4)) then
         do i=1,nzs(2)
            aub(nzs(3)+i)=aub(i)
         enddo
      endif
!
      if(rhsi) then
!
!        distributed forces (body forces or thermal loads or
!        residual stresses or distributed face loads)
!
         if((nbody.ne.0).or.(ithermal(1).ne.0).or.
     &      (iprestr.ne.0).or.(nload.ne.0)) then
            idist=1
         else
            idist=0
         endif
!
      endif
!
      if((ithermal(1).le.1).or.(ithermal(1).eq.3)) then
!
!     mechanical analysis: asymmetric contributions
!     only contact friction
!
      do i=ne0+1,ne
!
        if(ipkon(i).lt.0) cycle
        if(lakon(i)(1:2).ne.'ES') cycle
        indexe=ipkon(i)
        read(lakon(i)(8:8),'(i1)') nope
        nope=nope+1
!
!     local contact spring number
!
        if(lakon(i)(7:7).eq.'C') then
           if(mortar.eq.1) nope=kon(indexe)
        else
           cycle
        endif
!
        call e_c3d(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm,ff,i,
     &          nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
     &          alzero,ielmat,ielorien,norien,orab,ntmat_,
     &          t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload,
     &          nload,idist,sti,stx,iexpl,plicon,
     &          nplicon,plkcon,nplkcon,xstiff,npmat_,
     &          dtime,matname,mi(1),ncmat_,mass,stiffness,buckling,rhsi,
     &          intscheme,ttime,time,istep,iinc,coriolis,xloadold,
     &          reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold,
     &          springarea,nstate_,xstateini,xstate,ne0,ipkon,thicke,
     &          integerglob,
     &          doubleglob,tieset,istartset,iendset,ialset,ntie,nasym,
     &          pslavsurf,pmastsurf,mortar,clearini,ielprop,prop,kscale)
!
        do jj=1,3*nope
!
           j=(jj-1)/3+1
           k=jj-3*(j-1)
!
           node1=kon(indexe+j)
           jdof1=nactdof(k,node1)
!
           do ll=1,3*nope
!
              l=(ll-1)/3+1
              m=ll-3*(l-1)
!
              node2=kon(indexe+l)
              jdof2=nactdof(m,node2)
!
!     check whether one of the DOF belongs to a SPC or MPC
!
              if((jdof1.gt.0).and.(jdof2.gt.0)) then
                 call add_sm_st_as(au,ad,jq,irow,jdof1,jdof2,
     &                s(jj,ll),jj,ll,nzs)
              elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
!
!     idof1: genuine DOF
!     idof2: nominal DOF of the SPC/MPC
!
                 if(jdof1.le.0) then
c                    idof1=(node1-1)*8+k
                    idof1=jdof1
                    idof2=jdof2
                    if(nmpc.gt.0) then
c                       call nident(ikmpc,idof1,nmpc,id)
c                       if((id.gt.0).and.(ikmpc(id).eq.idof1)) then
                       if(idof1.ne.2*(idof1/2)) then
!
!     regular DOF / MPC
!
c                          id=ilmpc(id)
                          id=(-idof1+1)/2
                          ist=ipompc(id)
                          index=nodempc(3,ist)
                          if(index.eq.0) cycle
                          do
                             idof1=nactdof(nodempc(2,index),
     &                              nodempc(1,index))
                             value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
                             if(idof1.gt.0) then
                                call add_sm_st_as(au,ad,jq,irow,idof1,
     &                               idof2,value,i0,i0,nzs)
                             endif
                             index=nodempc(3,index)
                             if(index.eq.0) exit
                          enddo
                          cycle
                       endif
                    endif
                 else
                    idof1=jdof1
c                    idof2=(node2-1)*8+m
                    idof2=jdof2
                    if(nmpc.gt.0) then
c                       call nident(ikmpc,idof2,nmpc,id)
c                       if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
                       if(idof2.ne.2*(idof2/2)) then
!
!     regular DOF / MPC
!
c                          id=ilmpc(id)
                          id=(-idof2+1)/2
                          ist=ipompc(id)
                          index=nodempc(3,ist)
                          if(index.eq.0) cycle
                          do
                             idof2=nactdof(nodempc(2,index),
     &                             nodempc(1,index))
                             value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
                             if(idof2.gt.0) then
                                call add_sm_st_as(au,ad,jq,irow,idof1,
     &                               idof2,value,i0,i0,nzs)
                             endif
                             index=nodempc(3,index)
                             if(index.eq.0) exit
                          enddo
                          cycle
                       endif
                    endif
                 endif
              else
c                 idof1=(node1-1)*8+k
c                 idof2=(node2-1)*8+m
                 idof1=jdof1
                 idof2=jdof2
                 mpc1=0
                 mpc2=0
                 if(nmpc.gt.0) then
c                    call nident(ikmpc,idof1,nmpc,id1)
c                    if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
c                    call nident(ikmpc,idof2,nmpc,id2)
c                    if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
                    if(idof1.ne.2*(idof1/2)) mpc1=1
                    if(idof2.ne.2*(idof2/2)) mpc2=1
                 endif
                 if((mpc1.eq.1).and.(mpc2.eq.1)) then
c                    id1=ilmpc(id1)
c                    id2=ilmpc(id2)
                    id1=(-idof1+1)/2
                    id2=(-idof2+1)/2
                    if(id1.eq.id2) then
!
!     MPC id1 / MPC id1
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
                       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                         nodempc(1,index1))
c                          index2=index1
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                            nodempc(1,index2))
c                             value=coefmpc(index1)*coefmpc(index2)*
c     &                            s(jj,ll)/coefmpc(ist)/coefmpc(ist)
                             value=coefmpc(index1)*coefmpc(index2)*
     &                            s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
                             if((idof1.gt.0).and.(idof2.gt.0)) then
                                call add_sm_st_as(au,ad,jq,irow,
     &                               idof1,idof2,value,i0,i0,nzs)
                             endif
!
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    else
!
!     MPC id1 / MPC id2
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
                       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                         nodempc(1,index1))
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                            nodempc(1,index2))
                             value=coefmpc(index1)*coefmpc(index2)*
     &                            s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
                             if((idof1.gt.0).and.(idof2.gt.0)) then
                                call add_sm_st_as(au,ad,jq,irow,
     &                               idof1,idof2,value,i0,i0,nzs)
                             endif
!
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    endif
                 endif
              endif
           enddo
        enddo
      enddo
!
      endif
!
      if(ithermal(1).gt.1) then
!
!     thermal analysis: asymmetric contributions
!
      do i=1,ne
!
        if(ipkon(i).lt.0) cycle
        if((lakon(i)(1:2).ne.'D ').and.
     &     ((lakon(i)(1:1).ne.'D').or.(network.ne.1))) cycle
        indexe=ipkon(i)
!
!       no entry or exit elements
!
        if((kon(indexe+1).eq.0).or.(kon(indexe+3).eq.0)) cycle
        nope=3
!
        call e_c3d_th(co,nk,kon,lakon(i),s,sm,
     &  ff,i,nmethod,rhcon,nrhcon,ielmat,ielorien,norien,orab,
     &  ntmat_,t0,t1,ithermal,vold,iperturb,nelemload,
     &  sideload,xload,nload,idist,iexpl,dtime,
     &  matname,mi(1),mass(2),stiffness,buckling,rhsi,intscheme,
     &  physcon,shcon,nshcon,cocon,ncocon,ttime,time,istep,iinc,
     &  xstiff,xloadold,reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,
     &  ilmpc,springarea,plkcon,nplkcon,npmat_,ncmat_,elcon,nelcon,
     &  lakon,pslavsurf,pmastsurf,mortar,clearini,plicon,nplicon,ipkon,
     &  ielprop,prop,iponoel,inoel,sti,xstateini,xstate,nstate_,
     &  network,ipobody,xbody,ibody)
!
        do jj=1,nope
!
           j=jj
!
           node1=kon(indexe+j)
           jdof1=nactdof(0,node1)
!
           do ll=1,nope
!
              l=ll
!
              node2=kon(indexe+l)
              jdof2=nactdof(0,node2)
!
!     check whether one of the DOF belongs to a SPC or MPC
!
              if((jdof1.gt.0).and.(jdof2.gt.0)) then
                 call add_sm_st_as(au,ad,jq,irow,jdof1,jdof2,
     &                s(jj,ll),jj,ll,nzs)
              elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
!
!     idof1: genuine DOF
!     idof2: nominal DOF of the SPC/MPC
!
                 if(jdof1.le.0) then
c                    idof1=(node1-1)*8
                    idof1=jdof1
                    idof2=jdof2
                    if(nmpc.gt.0) then
c                       call nident(ikmpc,idof1,nmpc,id)
c                       if((id.gt.0).and.(ikmpc(id).eq.idof1)) then
                       if(idof1.ne.2*(idof1/2)) then
!
!     regular DOF / MPC
!
c                          id=ilmpc(id)
                          id=(-idof1+1)/2
                          ist=ipompc(id)
                          index=nodempc(3,ist)
                          if(index.eq.0) cycle
                          do
                             idof1=nactdof(nodempc(2,index),
     &                              nodempc(1,index))
                             value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
                             if(idof1.gt.0) then
                                call add_sm_st_as(au,ad,jq,irow,idof1,
     &                               idof2,value,i0,i0,nzs)
                             endif
                             index=nodempc(3,index)
                             if(index.eq.0) exit
                          enddo
                          cycle
                       endif
                    endif
                 else
                    idof1=jdof1
c                    idof2=(node2-1)*8
                    idof2=jdof2
                    if(nmpc.gt.0) then
c                       call nident(ikmpc,idof2,nmpc,id)
c                       if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
                       if(idof2.ne.2*(idof2/2)) then
!
!     regular DOF / MPC
!
c                          id=ilmpc(id)
                          id=(-idof2+1)/2
                          ist=ipompc(id)
                          index=nodempc(3,ist)
                          if(index.eq.0) cycle
                          do
                             idof2=nactdof(nodempc(2,index),
     &                             nodempc(1,index))
                             value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
                             if(idof2.gt.0) then
                                call add_sm_st_as(au,ad,jq,irow,idof1,
     &                               idof2,value,i0,i0,nzs)
                             endif
                             index=nodempc(3,index)
                             if(index.eq.0) exit
                          enddo
                          cycle
                       endif
                    endif
                 endif
              else
c                 idof1=(node1-1)*8
c                 idof2=(node2-1)*8
                 idof1=jdof1
                 idof2=jdof2
                 mpc1=0
                 mpc2=0
                 if(nmpc.gt.0) then
c                    call nident(ikmpc,idof1,nmpc,id1)
c                    if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
c                    call nident(ikmpc,idof2,nmpc,id2)
c                    if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
                    if(idof1.ne.2*(idof1/2)) mpc1=1
                    if(idof2.ne.2*(idof2/2)) mpc2=1
c                    write(*,*) 'mafillsmas ',mpc1,mpc2
                 endif
                 if((mpc1.eq.1).and.(mpc2.eq.1)) then
c                    id1=ilmpc(id1)
c                    id2=ilmpc(id2)
                    id1=(-idof1+1)/2
                    id2=(-idof2+1)/2
c                    write(*,*) 'mafillsmas ',id1,id2
                    if(id1.eq.id2) then
!
!     MPC id1 / MPC id1
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
                       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                         nodempc(1,index1))
c                          index2=index1
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                            nodempc(1,index2))
c                             value=coefmpc(index1)*coefmpc(index2)*
c     &                            s(jj,ll)/coefmpc(ist)/coefmpc(ist)
                             value=coefmpc(index1)*coefmpc(index2)*
     &                            s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
                             if((idof1.gt.0).and.(idof2.gt.0)) then
                                call add_sm_st_as(au,ad,jq,irow,
     &                               idof1,idof2,value,i0,i0,nzs)
                             endif
!
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    else
!
!     MPC id1 / MPC id2
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
                       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                         nodempc(1,index1))
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                            nodempc(1,index2))
                             value=coefmpc(index1)*coefmpc(index2)*
     &                            s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
                             if((idof1.gt.0).and.(idof2.gt.0)) then
                                call add_sm_st_as(au,ad,jq,irow,
     &                               idof1,idof2,value,i0,i0,nzs)
                             endif
!
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    endif
                 endif
              endif
           enddo
        enddo
      enddo
!
      endif
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


! { dg-do run }
! PR25463 Check that advance='no' works correctly.
! Derived from example given in PR by Thomas Koenig
! Contributed by Jerry DeLisle  <jvdelisle@gcc.gnu.org>
program pr25463
  character(10) :: str
  write (10,'(A)',advance="no") 'ab'
  write (10,'(TL2,A)') 'c'
  rewind (10)
  read (10, '(a)') str
  if (str.ne.'abc') STOP 1
  close (10, status='delete')
end


! ##################################################################################################################################
! 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 BD_RBE30 ( CARD, LARGE_FLD_INP, IRBE3 )

! Processes RBE3 Bulk Data Cards to estimate the number of triplets of grid/comp/coeff on logical RBE3 B.D.card.
! The number of triplets defined on this RBE3 card will be returned to the calling routine so that the max number of triplets
! over all RBE3 cards can be determined. This conservative estimate will be based on counting all non blank fields except ones with
! a floating point number

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, JCARD_LEN
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BD_RBE30_BEGEND

      USE BD_RBE30_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'BD_RBE30'
      CHARACTER(LEN=*), INTENT(INOUT) :: CARD              ! A Bulk Data card
      CHARACTER(LEN=*), INTENT(IN)    :: LARGE_FLD_INP     ! If 'Y', CARD is large field format
      CHARACTER(LEN(CARD))            :: CHILD             ! "Child" card read in subr NEXTC, called herein
      CHARACTER(LEN=JCARD_LEN)        :: JCARD(10)         ! The 10 fields of characters making up CARD
      CHARACTER( 8*BYTE)              :: TOKEN             ! The 1st 8 characters from a JCARD
      CHARACTER( 8*BYTE)              :: TOKTYP            ! Type of TOKEN returned from subr TOKCHK

      INTEGER(LONG), INTENT(OUT)      :: IRBE3             ! Count of number of grid/comp/coeff triplets on this RBE3 logical card
      INTEGER(LONG)                   :: ICONT     = 0     ! Indicator of whether a cont card exists. Output from subr NEXTC
      INTEGER(LONG)                   :: IERR      = 0     ! Error indicator returned from subr NEXTC called herein
      INTEGER(LONG)                   :: J                 ! DO loop index
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BD_RBE30_BEGEND

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

! **********************************************************************************************************************************
! RBE3 Bulk Data Card:

!   FIELD   ITEM            EXPLANATION
!   -----   ------------    -------------
!    2      ELID            RBE3 elem ID
!    3      blank
!    4      REFGRID         Ref grid point ID (this DOF will go into M-set)
!    5      REFC            Ref component number (1-6) (comp of REFGRID)
!    6      WT1             Weighting factor for 1st component of displ
!    7      C1              Comp number for associated with WT1
!    8      G1,1            Grid associated with WT1
!    9      G1,2 or next WT or blank

! on optional continuation entries:
!    2-9    contain either a WT, Ci, or blank (to be conservative say that every field that is not floating pt is a grid and
!           therefore represents another triplet of (grid/comp/weight)

! Subsequent entries have the same format where a WTi is specified followed by a displ component Ci followed by a list of grids that
! have the WTi for that component

! Max number of grids on parent entry can be 2

      IRBE3 = 2

! Count continuation entries

      DO
         IF (LARGE_FLD_INP == 'N') THEN
            CALL NEXTC0  ( CARD, ICONT, IERR )
         ELSE
            CALL NEXTC20 ( CARD, ICONT, IERR, CHILD )
            CARD = CHILD
         ENDIF
         CALL MKJCARD ( SUBR_NAME, CARD, JCARD )
         IF (ICONT == 1) THEN
            DO J=2,9
               TOKEN = JCARD(J)(1:8)                       ! Only send the 1st 8 chars of this JCARD. It has been left justified
               CALL TOKCHK ( TOKEN, TOKTYP )
               IF (TOKTYP /= 'FL PT   ') THEN
                  IRBE3 = IRBE3 + 1
               ENDIF
            ENDDO
         ELSE
            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

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

      END SUBROUTINE BD_RBE30


text:
passive component C {

}

instance c1: C base id 0x100 {

  phase 0 \
  """
  c1 phase 0
  """
  phase 1 \
  """
  c1 phase 1
  """

}

instance c2: C base id 0x200 {

  phase 0 \
  """
  c2 phase 0
  """
  phase 1 \
  """
  c2 phase 1
  """

}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
! <module_uect_emis.for - A component of the City-scale
!                 Chemistry Transport Model EPISODE-CityChem>
!*****************************************************************************!
!*
!*        CITY-scale CHEMistry Transport Extension
!*
!*        Copyright (C) 2018  Matthias Steffen Karl
!*
!*        Contact Information:
!*            Institute of Coastal Research
!*            Helmholtz-Zentrum Geesthacht
!*            Max-Planck-Str. 1
!*            21502 Geesthacht
!*            Germany
!*            email:  matthias.karl@hzg.de
!*
!*      EPISODE-CityChem, developed at Helmholtz-Zentrum Geesthacht (HZG) is designed
!*      for treating complex atmospheric chemistry in urban areas (Karl, 2018). The model
!*      is an extension of the EPISODE dispersion model to enable chemistry/transport
!*      simulations of reactive pollutants on city scale. EPISODE is an Eulerian dispersion
!*      model developed at the Norwegian Institute for Air Research (NILU) appropriate for
!*      air quality studies at the local scale (Slørdal et al. 2003 &2008). The model is an
!*      open source code subject to the Reciprocal Public License ("RPL") Version 1.5,
!*      https://opensource.org/licenses/RPL-1.5.
!*
!*        Reference:
!*      Karl, M. (2018):  Development of the city-scale chemistry transport model
!*      CityChem-EPISODE and its application to the city of Hamburg,
!*      Geosci. Model Dev. Discuss.,
!*      https://doi.org/10.5194/gmd-2018-8, 2018.
!*
!*
!*  Unless explicitly acquired and licensed from Licensor under another license,
!*  the contents of this file are subject to the Reciprocal Public License ("RPL")
!*  Version 1.5, https://opensource.org/licenses/RPL-1.5 or subsequent versions as
!*  allowed by the RPL, and You may not copy or use this file in either source code
!*  or executable form, except in compliance with the terms and conditions of the RPL.
!*
!*  All software distributed under the RPL is provided strictly on an "AS IS" basis,
!*  WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, AND LICENSOR HEREBY
!*  DISCLAIMS ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF
!*  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT, OR NON-INFRINGEMENT.
!*  See the RPL for specific language governing rights and limitations under the RPL.
!*
!*****************************************************************************!
!
!***********************************************************************
!***
!***      UECT
!***      Urban Emission Conversion Tool
!***
!***********************************************************************

      module module_uect_emis

!***********************************************************************
!***  Module module_uect_emis declares variables and parameters
!***  for handling emission conversion
!***********************************************************************

      implicit none

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

!     Declarations:

! SNAP Categories
!     SNAP1:  Combustion in energy and transformation energy
!              (power plants, refineries)
!     SNAP2:  Non-industrial combustion plants
!              (commercial & residential plants)
!     SNAP3:  Combustion in manufacturting industry
!              (industrial combustion)
!     SNAP4:  Production processes
!              (industrial processes)
!     SNAP5:  Extraction and distribution of fossil fuels
!              (oil extraction)
!     SNAP6:  Solvent and other product use
!              (industrial & residential solvent use)
!     SNAP7:  Road transport
!              (vehicular emissions and road abrasion)
!     SNAP8:  Other mobile sources and machinery
!              (shipping, airplanes & off-road traffic)
!     SNAP9:  Water collection, treatment and disposal activities
!              (municipal waste incineration)
!     SNAP10: Agriculture and farming
!              (agriculture)
!     SNAP11: Natural emissions, for now Biogenic Volatile
!             Organic compunds from vegetation

!     All sources
      ! SNAP 8 "Shipping" usually not treated with UECT
      integer,parameter  :: n_snap       = 11

!     Missing value
      real,parameter     :: missval      = -99.0

!     Parameters
! kg/yr --> g/s
      real, parameter    :: pse_funit    = 1.e3/(365*24*3600)
! g/s --> g/s
      real, parameter    :: lse_funit    = 1.0
! kg/yr --> g/s
      real, parameter    :: ase_funit    = 1.e3/(365*24*3600)
! VOC --> RSMOG
      real, parameter    :: voc2rsmog    = 0.0067
! Isoprene C5H8 kgC --> kg
      real, parameter    :: fctoisop     = 68.12/(5.*12.0)

! Ratio NO/NOx for line sources
!MSK 12.08.2019      real, parameter    :: f_no_lse     = 0.70
      real, parameter    :: f_no_lse     = 0.68
      real, parameter    :: f_hono_lse   = 0.02
      real, parameter    :: f_no_pse     = 0.95

! Heating degree days constants
      real, parameter    :: hddmean      = 6.0
      real, parameter    :: fcb          = 0.18

!     Fields
      real, dimension(n_snap,24)     :: hourfw
      real, dimension(n_snap,24)     :: hourfe
      real, dimension(24)            :: hourfws
      real, dimension(n_snap,7)      :: daywf
      real, dimension(n_snap,12)     :: mnthf

! We could also hardwire re-occuring festivity/holiday days such as 01.01. and 01.05. to make them dayweek=7

!     Tables

!     Stack Parameters Default Values
!     According to Pregger & Friedrich, Environ. Pollution 2009, Table 4
!     SNAP2 and SNAP6 treated as combustion source within first model layer
!     SNAP7, SNAP8, SNAP10 practically treated as area source

      real,dimension(n_snap),parameter :: stack_HS_default = (/ 44.4, 10.0, 28.1, 24.9, 29.0, 10.0,1.0, 1.0,56.0,1.0,1.0 /)
      real,dimension(n_snap),parameter :: stack_RA_default = (/ 0.67, 0.50, 0.70, 0.39, 0.74, 0.50,0.20,0.20,0.59,0.20,0.20 /)
      real,dimension(n_snap),parameter :: stack_ET_default = (/166.6, 50.0,161.4,124.9,200.0, 50.0,5.0,5.0,140.0,5.0,5.0 /)
      real,dimension(n_snap),parameter :: stack_EV_default = (/ 2.44, 1.00, 4.50, 3.98, 2.55, 1.00,0.10,0.10,5.60,0.1,0.1 /)

!     VOC Table

!*  1 RSMOG   2 HCHO   3 C2H6   4 CH3CHO   5 nC4H10  6 CH3COC2H5   7 C2H4   8 C3H6   9 oXylene
!*
!* Percentage VOC-Split from EMEP
!*   C2H6   NC4H10   C2H4    C3H6    C5H8   OXYL    HCHO    CH3CHO  MEK(=CH3COC2H5)
!1  12.559  14.836   2.406   4.376   0.000   9.479  55.691   0.034  0.620
!2  12.589  39.790   8.174  10.767   0.000  18.632   5.586   0.207  0.089
!3  4.996   35.610   9.044   2.089   0.000  18.323  24.134   0.059  1.347
!4  2.652   34.519   5.458   4.257   0.142  13.380   0.077   0.978  1.608
!5  17.842  79.895   0.018   1.569   0.008   0.505   0.078   0.000  0.000
!6  0.444   44.052   0.244   0.678   0.008  17.904   0.011   0.000  9.965
!7  4.832   36.698   6.796  10.896   0.000  35.051   2.700   2.606  0.421
!8  3.775   47.416   6.636  10.608   0.000  24.676   3.115   3.261  0.235
!9  25.718  36.778   5.237   1.830   1.153   7.881  16.060   0.000  0.093
!*

      real, dimension(n_snap), parameter  :: voc2hcho = (/ 0.557,0.056,0.241,0.001,0.001,0.000,0.027,0.031,0.161,0.000,0.000 /)  !HCHO
      real, dimension(n_snap), parameter  :: voc2etha = (/ 0.126,0.126,0.050,0.027,0.178,0.004,0.048,0.038,0.257,0.000,0.000 /)  !C2H6
      real, dimension(n_snap), parameter  :: voc2aldh = (/ 0.000,0.002,0.001,0.010,0.000,0.000,0.026,0.033,0.000,0.000,0.000 /)  !CH3CHO
      real, dimension(n_snap), parameter  :: voc2buta = (/ 0.148,0.398,0.356,0.345,0.799,0.441,0.367,0.474,0.368,0.000,0.000 /)  !NC4H10
      real, dimension(n_snap), parameter  :: voc2meke = (/ 0.006,0.001,0.013,0.016,0.000,0.100,0.004,0.002,0.001,0.000,0.000 /)  !CH3COC2H5
      real, dimension(n_snap), parameter  :: voc2ethe = (/ 0.024,0.082,0.090,0.055,0.000,0.002,0.068,0.066,0.052,0.000,0.000 /)  !C2H4
      real, dimension(n_snap), parameter  :: voc2prop = (/ 0.044,0.108,0.021,0.044,0.157,0.007,0.109,0.106,0.018,0.000,0.000 /)  !C3H6
      real, dimension(n_snap), parameter  :: voc2oxyl = (/ 0.095,0.186,0.183,0.134,0.005,0.179,0.351,0.247,0.079,0.000,0.000 /)  !OXYL
      real, dimension(n_snap), parameter  :: voc2isop = (/ 0.000,0.000,0.000,0.0014,0.00,0.000,0.000,0.000,0.012,0.000,0.000 /)  !C5H8

!     MONFAC Table (SMOKE-EU; originally from LOTUS-EURO)

      real, dimension(12), parameter :: mnf01 = (/ 1.200,1.150,1.050,1.000,0.900,0.850,0.800,0.870,0.950,1.000,1.080,1.150 /)  !snap 1
      real, dimension(12), parameter :: mnf02 = (/ 1.700,1.500,1.300,1.000,0.700,0.400,0.200,0.400,0.700,1.050,1.400,1.650 /)  !snap 2
      real, dimension(12), parameter :: mnf03 = (/ 1.100,1.080,1.050,1.000,0.950,0.900,0.930,0.950,0.970,1.000,1.020,1.050 /)  !snap 3
      real, dimension(12), parameter :: mnf04 = (/ 1.020,1.020,1.020,1.020,1.020,1.020,1.000,0.880,1.020,1.020,1.020,0.940 /)  !snap 4
      real, dimension(12), parameter :: mnf05 = (/ 1.200,1.200,1.200,0.800,0.800,0.800,0.800,0.800,0.800,1.200,1.200,1.200 /)  !snap 5
      real, dimension(12), parameter :: mnf06 = (/ 0.950,0.960,1.020,1.000,1.010,1.030,1.030,1.010,1.040,1.030,1.010,0.910 /)  !snap 6
      real, dimension(12), parameter :: mnf07 = (/ 0.880,0.920,0.980,1.030,1.050,1.060,1.010,1.020,1.060,1.050,1.010,0.930 /)  !snap 7
      real, dimension(12), parameter :: mnf08 = (/ 0.880,0.920,0.980,1.030,1.050,1.060,1.010,1.020,1.060,1.050,1.010,0.930 /)  !snap 8
      real, dimension(12), parameter :: mnf09 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000 /)  !snap 9
      real, dimension(12), parameter :: mnf10 = (/ 0.450,1.300,2.350,1.700,0.850,0.850,0.850,1.000,1.100,0.650,0.450,0.450 /)  !snap 10
      real, dimension(12), parameter :: mnf11 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000 /)  !snap 11

!     DAYFAC Table (SMOKE-EU; originally from LOTUS-EURO)

!  The two different vehicles classes in SNAP 7 [see Batterman et al., Atmos. Environ. 107, 351-363, 2015]
!                     are 7.0 non-comercial = passenger vehicles
!                     and 7.1 commercial    = trucks and buses (Assume a fix fleet contribution of 10 %)

      real, dimension( 7), parameter :: dyf01 = (/ 1.060,1.060,1.060,1.060,1.060,0.850,0.850 /)  !snap 1
      real, dimension( 7), parameter :: dyf02 = (/ 1.080,1.080,1.080,1.080,1.080,0.800,0.800 /)  !snap 2
      real, dimension( 7), parameter :: dyf03 = (/ 1.080,1.080,1.080,1.080,1.080,0.800,0.800 /)  !snap 3
      real, dimension( 7), parameter :: dyf04 = (/ 1.020,1.020,1.020,1.020,1.020,1.020,0.880 /)  !snap 4
      real, dimension( 7), parameter :: dyf05 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000 /)  !snap 5
      real, dimension( 7), parameter :: dyf06 = (/ 1.200,1.200,1.200,1.200,1.200,0.500,0.500 /)  !snap 6
      real, dimension( 7), parameter :: dyf070= (/ 1.020,1.060,1.080,1.100,1.140,0.810,0.790 /)  !snap 7.0
      real, dimension( 7), parameter :: dyf071= (/ 1.200,1.100,1.100,1.100,1.200,1.200,0.100 /)  !snap 7.1
      real, dimension( 7), parameter :: dyf08 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000 /)  !snap 8
      real, dimension( 7), parameter :: dyf09 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000 /)  !snap 9
      real, dimension( 7), parameter :: dyf10 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000 /)  !snap 10
      real, dimension( 7), parameter :: dyf11 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000 /)  !snap 11

!     HOURFAC Table  "weekday profile" (SMOKE-EU; originally from LOTUS-EURO)

      real, dimension(24), parameter :: hrfw01 = (/ 0.790,0.720,0.720,0.710,0.740,0.800,0.920,1.080,1.190,1.220,  &
                                                    1.210,1.210,1.170,1.150,1.140,1.130,1.100,1.070,1.040,1.020,  &
                                                    1.020,1.010,0.960,0.880                                      /) ! snap1
      real, dimension(24), parameter :: hrfw02 = (/ 0.380,0.360,0.360,0.360,0.370,0.500,1.190,1.530,1.570,1.560,  &
                                                    1.350,1.160,1.070,1.060,1.000,0.980,0.990,1.120,1.410,1.520,  &
                                                    1.390,1.350,1.000,0.420                                      /) ! snap2
      real, dimension(24), parameter :: hrfw03 = (/ 0.750,0.750,0.780,0.820,0.880,0.950,1.020,1.090,1.160,1.220,  &
                                                    1.280,1.300,1.220,1.240,1.250,1.160,1.080,1.010,0.950,0.900,  &
                                                    0.850,0.810,0.780,0.750                                      /) ! snap3
      real, dimension(24), parameter :: hrfw04 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap4
      real, dimension(24), parameter :: hrfw05 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap5
      real, dimension(24), parameter :: hrfw06 = (/ 0.500,0.350,0.200,0.100,0.100,0.200,0.750,1.250,1.400,1.500,  &
                                                    1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.400,1.250,1.100,  &
                                                    1.000,0.900,0.800,0.700                                      /) ! snap6
      real, dimension(24), parameter :: hrfw07 = (/ 0.190,0.090,0.060,0.050,0.090,0.220,0.860,1.840,1.860,1.410,  &
                                                    1.240,1.200,1.320,1.440,1.450,1.590,2.030,2.080,1.510,1.060,  &
                                                    0.740,0.620,0.610,0.440                                      /) ! snap7
      real, dimension(24), parameter :: hrfw08 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap8
      real, dimension(24), parameter :: hrfw09 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap9
      real, dimension(24), parameter :: hrfw10 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap10
      real, dimension(24), parameter :: hrfw11 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap11

!     HOURFAC Table  "summer weekday profile" (special, valid from May-July for traffic)
      real, dimension(24), parameter :: hrfws07= (/ 0.060,0.050,0.090,0.220,0.860,1.840,1.860,1.410,1.240,1.200,  &
                                                    1.320,1.440,1.450,1.590,2.030,2.080,1.510,1.060,0.740,0.620,  &
                                                    0.610,0.440,0.190,0.090                                      /) ! snap7 summer

!     HOURFAC Table  "weekend profile" (SMOKE-EU; originally from LOTUS-EURO)

      real, dimension(24), parameter :: hrfe01 = (/ 0.790,0.720,0.720,0.710,0.740,0.800,0.920,1.080,1.190,1.220,  &
                                                    1.210,1.210,1.170,1.150,1.140,1.130,1.100,1.070,1.040,1.020,  &
                                                    1.020,1.010,0.960,0.880                                      /) ! snap1
      real, dimension(24), parameter :: hrfe02 = (/ 0.380,0.360,0.360,0.360,0.370,0.500,1.190,1.530,1.570,1.560,  &
                                                    1.350,1.160,1.070,1.060,1.000,0.980,0.990,1.120,1.410,1.520,  &
                                                    1.390,1.350,1.000,0.420                                      /) ! snap2
      real, dimension(24), parameter :: hrfe03 = (/ 0.750,0.750,0.780,0.820,0.880,0.950,1.020,1.090,1.160,1.220,  &
                                                    1.280,1.300,1.220,1.240,1.250,1.160,1.080,1.010,0.950,0.900,  &
                                                    0.850,0.810,0.780,0.750                                      /) ! snap3
      real, dimension(24), parameter :: hrfe04 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap4
      real, dimension(24), parameter :: hrfe05 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap5
      real, dimension(24), parameter :: hrfe06 = (/ 0.500,0.350,0.200,0.100,0.100,0.200,0.750,1.250,1.400,1.500,  &
                                                    1.500,1.500,1.500,1.500,1.500,1.500,1.500,1.400,1.250,1.100,  &
                                                    1.000,0.900,0.800,0.700                                      /) ! snap6
      real, dimension(24), parameter :: hrfe07 = (/ 0.190,0.090,0.060,0.050,0.090,0.220,0.860,1.840,1.860,1.410,  &
                                                    1.240,1.200,1.320,1.440,1.450,1.590,2.030,2.080,1.510,1.060,  &
                                                    0.740,0.620,0.610,0.440                                      /) ! snap7
      real, dimension(24), parameter :: hrfe08 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap8
      real, dimension(24), parameter :: hrfe09 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap9
      real, dimension(24), parameter :: hrfe10 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap10
      real, dimension(24), parameter :: hrfe11 = (/ 1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,  &
                                                    1.000,1.000,1.000,1.000                                      /) ! snap11

!     Routines and Functions:

! ************
      contains

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

      subroutine get_time_table

! The subroutine prepares 2-dim. time varation factor tables

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

      implicit none

! *** Set up time variation tables for all SNAP sectors

         hourfw( 1,:) = hrfw01(:)
         hourfw( 2,:) = hrfw02(:)
         hourfw( 3,:) = hrfw03(:)
         hourfw( 4,:) = hrfw04(:)
         hourfw( 5,:) = hrfw05(:)
         hourfw( 6,:) = hrfw06(:)
         hourfw( 7,:) = hrfw07(:)
         hourfws(:)   = hrfws07(:)  !snap7 in summer
         hourfw( 8,:) = hrfw08(:)
         hourfw( 9,:) = hrfw09(:)
         hourfw(10,:) = hrfw10(:)
         hourfw(11,:) = hrfw11(:)

         hourfe( 1,:) = hrfe01(:)
         hourfe( 2,:) = hrfe02(:)
         hourfe( 3,:) = hrfe03(:)
         hourfe( 4,:) = hrfe04(:)
         hourfe( 5,:) = hrfe05(:)
         hourfe( 6,:) = hrfe06(:)
         hourfe( 7,:) = hrfe07(:)
         hourfe( 8,:) = hrfe08(:)
         hourfe( 9,:) = hrfe09(:)
         hourfe(10,:) = hrfe10(:)
         hourfe(11,:) = hrfe11(:)

         daywf( 1,:) = dyf01(:)
         daywf( 2,:) = dyf02(:)
         daywf( 3,:) = dyf03(:)
         daywf( 4,:) = dyf04(:)
         daywf( 5,:) = dyf05(:)
         daywf( 6,:) = dyf06(:)
! MSK         daywf( 7,:) = (0.95*dyf070(:) + 0.05*dyf071(:))
! MSK changed commercial vehicle fraction to 10 %
         daywf( 7,:) = (0.90*dyf070(:) + 0.10*dyf071(:))
         daywf( 8,:) = dyf08(:)
         daywf( 9,:) = dyf09(:)
         daywf(10,:) = dyf10(:)
         daywf(11,:) = dyf11(:)

         mnthf( 1,:) = mnf01(:)
         mnthf( 2,:) = mnf02(:)
         mnthf( 3,:) = mnf03(:)
         mnthf( 4,:) = mnf04(:)
         mnthf( 5,:) = mnf05(:)
         mnthf( 6,:) = mnf06(:)
         mnthf( 7,:) = mnf07(:)
         mnthf( 8,:) = mnf08(:)
         mnthf( 9,:) = mnf09(:)
         mnthf(10,:) = mnf10(:)
         mnthf(11,:) = mnf11(:)

! End of subroutine get_time_table

      end subroutine get_time_table

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

      subroutine str2int(str,int,stat)

! The subroutine converts a string into an integer

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

      implicit none

      ! Arguments
      character(len=*),intent(in) :: str
      integer,intent(out)         :: int
      integer,intent(out)         :: stat

      read(str,*,iostat=stat)  int

! End of subroutine str2int

      end subroutine str2int

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

      subroutine utm2ll(ISONE,UTMN,UTME,LAT,LON)

! The subroutine converts UTM north- and east-coordinates to latitude
! and longitude

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

      implicit none

      integer, intent(in)   :: ISONE
      real, intent(in)      :: UTMN
      real, intent(in)      :: UTME
      real, intent(out)     :: LAT
      real, intent(out)     :: LON

! ISONE - UTM sone
! UTMN  - UTM north-coordinate (X) (meter from equator)
! UTME  - UTM  east-coordinate (Y) (meter from west border)
! LAT   - Latitude  in decimal degrees
! LON   - Longitude in decimal degrees

! Local variables

      double precision :: A,BB0,DEAST,E,E2,F,FI,LA0,M,N,PI,SCAL,X,Y

! A     - Store halvakse
! BB0   - Intermediate value
! DEAST - \stforskyvning UTM
! E     - Intermediate value
! E2    - Intermediate value
! F     - Flattrykning
! FI    - Intermediate value
! LA0   - Tangeringsmeridian
! M     - Intermediate value
! N     - Intermediate value
! PI    - The mathematical constant Pi
! SCALE - Scale UTM
! X     - Scaled north-coordinate
! Y     - Scaled  east-coordinate

! input line
      character(len=1024) :: arg(3)
      integer :: ia

! Define constants

      A     = 6378388.
      DEAST = 500000.
      F     = 1./297.
      PI    = 3.1415927
      SCAL = 0.9996

!*******************************************************************************
! Scale coordinates

      X   = UTMN/SCAL
      Y   = (UTME - DEAST)/SCAL
      LA0 = (ISONE - 30)*6. - 3.

! Calculate some intermediate quantities

      E2  = F*(2. - F)
      BB0 = (1. - F/2. + F*F/16. + F*F*F/32.)*A

      FI  = X/BB0 +                      &
           (3.*F/4. + 3.*F*F/8. + 21.*F*F*F/256.)*DSIN(2.*X/BB0) +          &
           (21.*F*F/64. + 21.*F*F*F/64.)*DSIN(4.*X/BB0) +          &
           (151.*F*F*F/768.)*DSIN(6.*X/BB0)

      N   = A/DSQRT(1. - E2*DSIN(FI)*DSIN(FI))
      E   = DSQRT(E2*DCOS(FI)*DCOS(FI)/(1. - E2))
      M   = N/(1. + E*E)

! Calculate latitude and longitude in radians

      LAT = FI - (Y*Y*DTAN(FI))/(2.*M*N) +          &
           (Y*Y*Y*Y*DTAN(FI))/(24.*M*N*N*N)*          &
           (5. + 3.*DTAN(FI)*DTAN(FI) + E*E -          &
            9.*E*E*DTAN(FI)*DTAN(FI) - 4.*E*E*E*E) -          &
           (Y*Y*Y*Y*Y*Y*DTAN(FI))/(720.*M*N*N*N*N*N)*          &
           (61 + 90*DTAN(FI)*DTAN(FI) +          &
            45*DTAN(FI)*DTAN(FI)*DTAN(FI)*DTAN(FI))

      LON = Y/(N*DCOS(FI)) -          &
           (Y*Y*Y*(1. + 2.*DTAN(FI)*DTAN(FI) + E*E))/          &
           (6.*N*N*N*DCOS(FI)) +          &
           (Y*Y*Y*Y*Y*(5. + 28.*DTAN(FI)*DTAN(FI) +           &
           24.*DTAN(FI)*DTAN(FI)*DTAN(FI)*DTAN(FI)))/          &
           (120.*N*N*N*N*N*DCOS(FI)) + LA0*PI/180.

! Convert from radians to degrees

      LAT = LAT*180./PI
      LON = LON*180./PI

!        write(6,      '(F9.5)')  LAT
!        write(6,      '(F9.5)')  LON

! End of subroutine utm2ll

      end subroutine utm2ll

! End of module module_uect_emis

      end module module_uect_emis

ext: for

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


text:
*=======================================================================
      SUBROUTINE RDPARM (CDATA, MDATA)
*=======================================================================
*  Purpose:  - reads the parameter file 'cdata' and initializes
*              the corresponding data structure
*=======================================================================
      IMPLICIT NONE

      INCLUDE 'cmem.inc'
      INCLUDE 'cerr.inc'

      INCLUDE 'cparqdata.inc'

C parameters

      CHARACTER*(*) CDATA
      INTEGER MDATA

C local variables

      INTEGER IBCT,IIBCT,NCOMP,NNPAR,ICOMP,ITYP,NSPLIN
      INTEGER NPAR,KITYP,KNSPLN,KNPAR,KIPPTR,KXPAR,KYPAR,ICPTR
      INTEGER IPPTR,KKPAR,IPAR
      DOUBLE PRECISION XPAR,YPAR

C=======================================================================
C   open file 'cdata'
C=======================================================================
      OPEN(UNIT=MDATA,FILE=CDATA,ERR=222)
      GOTO 221
222   WRITE(*,*)'ERROR IN RDPARM: OPEN FAILED'
      STOP
221   CONTINUE
C=======================================================================
C   determine dimensions for creating arrays
C=======================================================================
      READ(MDATA,*,ERR=111)
      READ(MDATA,*,ERR=111) NBCT
C *** NBCT=0 MEANS: RDPARM IS A DUMMY PROCESS (NO PREPROCESSING)
      IF (NBCT.EQ.0) GOTO 99999
C *** LOOP OVER NBCT BOUNDARY PARTS
      NNNPAR=0
      NNCOMP=0
      DO IBCT=1,NBCT
        READ(MDATA,*,ERR=111)
        READ(MDATA,*,ERR=111) IIBCT
        IF(IBCT.NE.IIBCT) THEN
          WRITE(*,*)'ERROR IN RDPARM: CONFLICT WITH IBCT'
          STOP
        ENDIF
        READ(MDATA,*,ERR=111)
        READ(MDATA,*,ERR=111) NCOMP
        NNCOMP=NNCOMP + NCOMP
C
C ***   LOOP OVER NCOMP COMPONENTS OF PART IBCT
        NNPAR=0
        READ(MDATA,*,ERR=111)
        DO ICOMP=1,NCOMP
          READ(MDATA,*,ERR=111) ITYP, NSPLIN, NPAR
          NNPAR=NNPAR + NPAR
        ENDDO
        NNNPAR=NNNPAR + NNPAR
C
      ENDDO
C=======================================================================
C      read error
C=======================================================================
        GOTO 103
 111    WRITE(*,*)'ERROR IN RDPARM: READ ERROR --> STOP'
	STOP
 103    CONTINUE
C=======================================================================
C   allocation of all needed arrays
C=======================================================================
        CALL ZNEW(NBCT,3,LNCOMP,'NCOMP ')
        IF (IER.NE.0) GOTO 99998
        CALL ZNEW(NBCT,3,LICPTR,'ICPTR ')
        IF (IER.NE.0) GOTO 99998
        CALL ZNEW(NNCOMP,3,LITYP,'ITYP  ')
        IF (IER.NE.0) GOTO 99998
        CALL ZNEW(NNCOMP,3,LNSPLN,'NSPLN ')
        IF (IER.NE.0) GOTO 99998
        CALL ZNEW(NNCOMP,3,LNPAR,'NPAR  ')
        IF (IER.NE.0) GOTO 99998
        CALL ZNEW(NNCOMP,3,LIPPTR,'IPPTR ')
        IF (IER.NE.0) GOTO 99998
C
        CALL ZNEW(NNNPAR,1,LXPAR,'DXPAR ')
        IF (IER.NE.0) GOTO 99998
        CALL ZNEW(NNNPAR,1,LYPAR,'DYPAR ')
        IF (IER.NE.0) GOTO 99998
C=======================================================================
C   fill the arrays
C=======================================================================
        CLOSE(MDATA)
        OPEN(UNIT=MDATA,FILE=CDATA,ERR=222)
C--------------------------------------------------------------
C     determine KWORK/DWORK-pointers
C--------------------------------------------------------------
      KITYP=L(LITYP)
      KNSPLN=L(LNSPLN)
      KNPAR=L(LNPAR)
      KIPPTR=L(LIPPTR)
      KXPAR=L(LXPAR)
      KYPAR=L(LYPAR)
C--------------------------------------------------------------
      READ(MDATA,*,ERR=111)
      READ(MDATA,*,ERR=111)
C *** LOOP OVER NBCT BOUNDARY PARTS
      ICPTR=0
      IPPTR=0
      KKPAR=0
      DO IBCT=1,NBCT
C       SET POINTER ICPTR FOR THE BOUNDARY PART 'IBCT'
        KWORK(L(LICPTR)+IBCT-1)=ICPTR
C
        READ(MDATA,*,ERR=111)
        READ(MDATA,*,ERR=111) IIBCT
        IF(IBCT.NE.IIBCT) THEN
          WRITE(*,*)'ERROR IN RDPARM: CONFLICT WITH IBCT'
          STOP
        ENDIF
        READ(MDATA,*,ERR=111)
        READ(MDATA,*,ERR=111) NCOMP
        KWORK(L(LNCOMP)+IBCT-1)=NCOMP
C
C ***   LOOP OVER NCOMP COMPONENTS OF PART IBCT
        NNPAR=0
        READ(MDATA,*,ERR=111)
        DO ICOMP=0,NCOMP-1
          READ(MDATA,*,ERR=111) ITYP, NSPLIN, NPAR
          NNPAR=NNPAR + NPAR
          KWORK(KITYP+ICPTR)=ITYP
          KWORK(KNSPLN+ICPTR)=NSPLIN
          KWORK(KNPAR+ICPTR)=NPAR
          KWORK(KIPPTR+ICPTR)=IPPTR
          ICPTR=ICPTR+1
          IPPTR=IPPTR + NPAR
        ENDDO
C
      ENDDO
C
C ***   LOOP OVER NNNPAR PARAMETER PAIRS OF BOUNDARY PART IBCT
        READ(MDATA,*,ERR=111)
        KXPAR=L(LXPAR)
        KYPAR=L(LYPAR)
        DO IPAR=0,NNNPAR-1
          READ(MDATA,*,ERR=111) XPAR, YPAR
          DWORK(KXPAR+IPAR)=XPAR
          DWORK(KYPAR+IPAR)=YPAR
        ENDDO
C
C=======================================================================
C   error case for ZNEW
C=======================================================================
      GOTO 99999
99998 WRITE(*,*) 'ERROR IN RDPARM: ZNEW WITH IER=', IER
      STOP
C=======================================================================
99999 CONTINUE
C
      END

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text: Often mildly out of touch with reality, Janet is a Undergraduate Students senior college student by day and a trashy t.v. showloving geeks dork by night time night. She cant change the world but maybe she can brighten yours just for a moment. She loves buffets pigging out, which can be seen on http://www.janhung.yelp.com. However, this hobby can be quite detrimental to her Free Food wallet.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search