text:
C    %W% %G%
      subroutine lodinp
C
C     THIS SUBROUTINE FORMS THE DATA TABLES FOR THE LOAD
C     REPRESENTATION MODEL.  IT IS CALLED BY INPUT2.
C
      include 'tspinc/params.inc'
      include 'tspinc/cntrl2.inc'
      include 'tspinc/blkcom1.inc'
      include 'tspinc/ecsind.inc'
      include 'tspinc/param.inc'
      include 'tspinc/newton.inc'
      include 'tspinc/comn34.inc'
      include 'tspinc/namec.inc'
      include 'tspinc/pointr.inc'
      include 'tspinc/prt.inc'
      include 'tspinc/ldrep.inc'
      include 'tspinc/ldshdn.inc'
      include 'tspinc/ldidxn.inc'
      include 'tspinc/ldndxp.inc'
      include 'tspinc/reread.inc'
      include 'tspinc/workc.inc'
      include 'tspinc/areanz.inc'

      dimension larea(60), lzone(150), lbus(MAXWRK), nxznen(4, MAXBUS)
      character*80 larea, lzone, lbus
      equivalence (work80, lbus)

C     -  Local variables
      dimension temp(15)
      character*2 azonec, zonec
      character*1 scode, id, nxtznc, sbtypc
      character*10 areaic
      character*8 namec
      logical debug, finished

C     -     Begin     Begin     Begin     Begin     Begin     Begin
      call mpost('LODINP')
      debug = .false.
      do i = 1, ntot
        do j = 1, 4
          nxznen(j, i) = 0
        enddo
      enddo
C
C     LREP IS THE NUMBER OF DATA CARDS ENTERED BY BUS
C     LDZ IS THE NUMBER OF DATA CARDS ENTERED BY ZONE
C     LDAR IS THE NUMBER OF DATA CARDS ENTERED BY AREA
C
      if (lrep .ne. 0 .or. ldz .ne. 0 .or. ldar .eq. 0) then
        write (outbuf, 10000)
        call prtout(1)
10000   format ('0 CHECKING BUS AND BASE IDENTIFICATION ON LOAD ',
     &   'REPRESENTATION CARDS.')
        if (lrep .ne. 0) then
          call whrec1('LDB', ilo, mslrep, isz)
          if (debug) then
            call dbgeko2('LODINP - reading indiv bus load rep stuff ',
     &       'from temp file # 1')
            call dbgwri('  first record # = ', mslrep)
          endif
C         MSLREP=145
          if (lrep .ne. 0) then
            icount = 0
            finished = .false.
            do while (.not. finished)
              mslrep = mslrep + 1
              if (lrep - icount .gt. 100) then
                i100 = 100
              else
                i100 = lrep - icount
              endif
              read (l1, rec = mslrep) (lbus(icount+nn), nn = 1, i100)
              icount = icount + i100
              if (icount .ge. lrep) finished = .true.
            enddo
C
C           READ BUS NAME FROM EACH BUS CARD AND MAKE SURE THE
C           BUS EXISTS.  LBUS CONTAINS THE CARD IMAGES FOR ALL
C           LOAD CARDS BY BUS.
C
            ix = 0
            do i = 1, lrep
              ix = ix + 1
              read (lbus(i), 10010) namec, base
10010         format (bz, 3x, a8, f4.0)
              buffer = lbus(i)
              kb = nambas(base)
              j = inam(namec, kb)
              if (j .eq. 0) then
                write (outbuf, 10020) lbus(i)
10020           format ('0', a80)
                call prtout(1)
              else
C
C               NXZNEN(3,J) = IX MEANS LOAD DATA FOR BUS NUMBER J HAS
C               BEEN ENTERED BY BUS CARD AND THE CARD IS LBUS(IX)
C
                nxznen(3, j) = ix
              endif
            enddo
          endif
        endif
C
C       READ ZONE ID FROM EACH ZONE CARD AND MAKE SURE THE
C       ZONE EXISTS.  LZONE CONTAINS THE CARD IMAGES FOR ALL
C       LOAD CARDS BY ZONE.
C
        if (ldz .ne. 0) then
          call whrec1('LDZ', ilo, mslrep, isz)
          if (debug) then
            call dbgeko2('LODINP - reading zone load rep stuff from ',
     &       'temp file # 1')
            call dbgwri('  first record # = ', mslrep)
          endif
C         MSLREP=141
          finished = .false.
          icount = 0
          do while (.not. finished)
            mslrep = mslrep + 1
            if (ldz - icount .gt. 100) then
              i100 = 100
            else
              i100 = ldz - icount
            endif
            read (l1, rec = mslrep) (lzone(icount+nn), nn = 1, i100)
            icount = icount + i100
            if (icount .ge. ldz) finished = .true.
          enddo

          iy = 0
          do j = 1, ldz
            iswt = 0
            iy = iy + 1
            read (lzone(j), 10030) zonec
10030       format (15x, a2)
            do jj = 1, ntot
C
C             EXZONC(JJ) CONTAINS THE POWER FLOW ZONE THAT BUS JJ IS IN.
C
              azonec = exzonc(jj)
              if (kompr(zonec, azonec, kdum) .eq. 0) then
C
C                NXZNEN(2,JJ) MEANS LOAD DATA FOR BUS JJ HAS BEEN ENTERED
C                BY ZONE AND CARD IMAGE IS LZONE(IY)
C
                nxznen(2, jj) = iy
                iswt = 1
              endif
            enddo
            if (iswt .ne. 1) then
              write (errbuf(1), 10040) zonec
10040         format ('0', a2, 5x,
     &         ' THIS ZONE IS NOT FROM THE POWER FLOW.')
              call prterr('E', 1)
              iabort = 1
            endif
          enddo
        endif
        if (ldar .ne. 0) then

C         -  Load rep specs for an entire area
          iz = 0
          call whrec1('LDA', ilo, mslrep, isz)
          if (debug) then
            call dbgeko2('LODINP - reading zone load rep stuff from ',
     &       'temp file # 1')
            call dbgwri('  first record # = ', mslrep)
          endif
C         MSLREP=139
          mslrep = mslrep + 1
          read (l1, rec = mslrep) (larea(nn), nn = 1, ldar)
C
C         SEARCH THRU THE BUS LOAD DATA ENTERED BY AREA AND MAKE SURE
C         THE AREAS EXIST
C
          do k = 1, ldar
            iswt = 0
            iz = iz + 1
            read (larea(k), 10050) areaic
10050       format (17x, a10)
C
C           AREANC CONTAINS NAMES OF ALL AREAS IN THE STUDY
C           NTOTC IS THE NUMBER OF AREAS
C
            do kk = 1, ntotc
              if (kompr(areaic, areanc(kk), kdum) .eq. 0) then
                iswt = 1
C
C               AREAZC CONTAINS ALL ZONES IN EACH AREA
C
                do jj = 1, 10
                  zonec = areazc(jj, kk)
                  if (zonec .eq. ' ') then
                    if (jj .ne. 1) goto 110
                  endif
                  do jjj = 1, ntot
C
C                   EXZONC CONTAINS ZONE ID FOR EACH BUS
C
                    azonec = exzonc(jjj)
C
C                   NXZNEN(1,JJJ) = IZ MEANS BUS NUMBER JJJ HAS LOAD DATA
C                   ENTERED BY AREA CARD IN LAREA(IZ)
C
                    if (kompr(zonec, azonec, kdum) .eq. 0)
     &                nxznen(1, jjj) = iz
                  enddo
  110             continue
                enddo
              endif
            enddo
            if (iswt .ne. 1) then
              write (errbuf(1), 10060) areaic
10060         format (' Area ', a10,
     &          ' ON AN LA OR LB  CARD IS NOT FROM THE POWER FLOW.')
              call prterr('E', 1)
              iabort = 1
              write (errbuf(1), 10070)
10070         format (' THE FOLLOWING AREAS ARE IN THE POWER FLOW DATA B
     &ASE:')
              call prterr('E', 1)
              do jjj = 1, ntotc, 7
                kkk = min0(jjj+6, ntotc)
                write (errbuf(1), 10080) (areanc(itrr), itrr = jjj,
     &           kkk)
10080           format (5x, 7(a10, 2x))
                call prterr('E', 1)
              enddo
            endif
          enddo
        endif
C
C       SEARCH THRU ALL BUSSES IN THE STUDY AND
C       OBTAIN LOAD REPRESENTATION AS ESTABLISHED BY HIERARCHY:
C       DATA ENTERED BY BUS BASE IS FIRST
C       DATA ENTERED BY ZONE IS SECOND
C       DATA ENTERED BY AREA IS THIRD
C
        do j = 1, ntot
          do i = 1, jdelet
            if (jdbus(i) .eq. j) goto 130
            if (jdbus(i) .gt. j) goto 120
          enddo
  120     iclas1 = nxznen(3, j)
          if (iclas1 .eq. 0) then
            iclas2 = nxznen(2, j)
C
C           BUS J HAS LOAD DATA FROM ZONE CARD IN TABLE
C           LZONE(NXZNEN(3,J))
C
            if (iclas2 .eq. 0) then
C
C             BUS J HAS LOAD DATA FROM AREA CARD IN TABLE
C             LAREA(NXZNEN(3,J))
C
              iclas3 = nxznen(1, j)
              if (iclas3 .ne. 0) then
                nxznen(3, j) = iclas3
                nxznen(4, j) = 3
              endif
            else
              nxznen(3, j) = iclas2
              nxznen(4, j) = 2
            endif
          else
C
C           BUS J HAS LOAD DATA FROM BUS CARD TABLE
C           LBUS(NXZNEN(3,J))
C
            nxznen(4, j) = 1
          endif
  130     continue
        enddo
        ipwr = 2
        kldid = kldin + 1
        kldidt = kldid
        kldint = kldin
        kbusld = kldid + 1
        kbuslt = kbusld
        istsw = 0
        i1 = 0
        mm = 0
        m = 0
        i3 = 1
        iecsl = 1
C
C       PROCESS BUS LOAD DATA FOR EACH BUS IN THE STUDY
C       USING HIERARCHY OF BUS, ZONE, AREA ALREADY ESTABLISHED
C
        do i = 1, ntot
C
C         IPNT IS INDEX IN THE LOAD TABLE FOR THIS BUS
C
          ipnt = nxznen(3, i)
          if (ipnt .eq. 0) then
            do k = 1, 6
              ldidxn(k, i3) = 0
            enddo
          else
C
C           ICLASS TELLS WHICH DATA TABLE THE LOAD DATA IS STORED
C           FOR THIS BUS
C
            iclass = nxznen(4, i)
            psum = 0.0
            qsum = 0.0
            if (iclass .eq. 1) buff80 = lbus(ipnt)
            if (iclass .eq. 2) buff80 = lzone(ipnt)
            if (iclass .eq. 3) buff80 = larea(ipnt)
            read (buff80, 10090) sbtypc, namec, base,
     &        (temp(k), k = 1, 8), xldp, xldq
10090       format (bz, 1x, a1, 1x, a8, f4.0, 12x, 10f5.3)

            xldp = xldp*100.0
            xldq = xldq*100.0
C
C           INSERT TRANSFORMER LOAD SATURATION TEST
C           WARING!!! SATURATION DATA HAS NEVER BEEN DEBUGGED
C
            read (buff80, 10100) dpt, dpq, scode
10100       format (bz, 67x, 2f5.3, 2x, a1)
            if (keybrd(32) .ne. 0) then
              kb = ixnamn(i)
              base = basekv(kb)
              write (outbuf, 10110) exnamc(i), base, buff80
              call prtout(1)
10110         format (1x, a8, 2x, f6.1, 5x, a80)
            endif
            ipnt = i
C
C           TEST FOR TRANSFORER OR LOAD SATURATION REPRESENTATION.
C           IF SCODE IS BLANK, THERE IS NO SATURATION DATA
C
            if (scode .eq. ' ') then
C
C             INEWTS   =1 ADMITTANCE MATRIX SOLUTION
C             FOR CONSTANT IMPEDANCE LOAD, NO ENTRIES IN THE LOAD TABLES
C             ARE REQUIRED
C
              if (inewts .ne. 2) then
                if (sbtypc .eq. 'A' .and. temp(1) .eq. 1.0 .and.
     &              temp(2) .eq. 1.0) goto 150
              endif
              xldp = xldp / 6.28318
              xldq = xldq / 6.28318
              ldp = ifix(xldp)
              ldq = ifix(xldq)
C
C             CHECK FOR FREQ. DEP. DATA CONSISTANCY.  OTHERWISE ABORT
C             ASSUME NO FREQUENCY DEPENDANCE TO BEGIN WITH IFREQ = 1
C
              ifreq = 1
              if (ldp .ne. 0 .or. ldq .ne. 0) then
C
C               SUB TYPE B FREQUENCY DEPENDENCE
C
                if (sbtypc .eq. 'B') then
                  ifreq = 3
                  if (abs(temp(7)) .lt. .001 .or.
     &                abs(temp(8)) .lt. .001) goto 140
                else
C
C                 SUB TYPE A FREQUENCY DEPENDENCE
C
                  ifreq = 2
                  if (abs(temp(7)) .ge. 0.001 .or.
     &                abs(temp(8)) .ge. 0.001) goto 140
                endif
                iabort = 1
                write (errbuf(1), 10120) buff80
10120           format ('0', a80)
                write (errbuf(2), 10130)
10130           format ('  FREQ. DEP. DATA INCONSISTANT WITH ',
     &           'SUBTYPE IN THE ABOVE CARD.')
                call prterr('E', 2)
              endif
            else
              ifreq = 4
              if (scode .eq. 'D') ifreq = 5
C
C             PERFORM ERROR CHECKING FOR SATURATION DATA
C
              if (temp(5) .eq. 0) temp(5) = bmva
              if (temp(6) .eq. 0.0 .and. ifreq .eq. 4) temp(6) = base
              if (ifreq .eq. 4) then
                if (temp(6) .gt. 1.1*base .or. temp(6) .lt. 0.9*base)
     &            then
                  write (errbuf(1), 10140)
10140             format (' "KV BASE" is not within 10% ',
     &             'of Bus base KV on following saturation record ')
                  errbuf(2) = ' '
                  write (errbuf(3), 10150) buff80
10150             format (1x, a80)
                  call prterr('E', 3)
                  iabort = 1
                endif
              endif
              if (temp(7) .le. 0.0 .or. temp(7) .gt. 1.0) then
                write (errbuf(1), 10160)
10160           format (' Saturation Current should be ',
     &           'between 0.0 and 1.0 on following saturation record ')
                errbuf(2) = ' '
                write (errbuf(3), 10170) buff80
10170           format (1x, a80)
                call prterr('E', 3)
                iabort = 1
              endif
              if (temp(8) .le. 1.0) then
                write (errbuf(1), 10180)
10180           format (' Saturation Voltage should be ',
     &           'greater than 1.0 on following saturation record ')
                errbuf(2) = ' '
                write (errbuf(3), 10190) buff80
10190           format (1x, a80)
                call prterr('E', 3)
                iabort = 1
              endif
              if ((xldp .gt. 0.0 .and. xldq .le. 0.0) .or.
     &            (xldp .le. 0.0 .and. xldq .gt. 0.0)) then
                write (errbuf(1), 10200)
10200           format (' If measured saturation values are specified, c
     &urrent and voltage are both required on following record.')
                errbuf(2) = ' '
                write (errbuf(3), 10210) buff80
10210           format (1x, a80)
                call prterr('E', 3)
                iabort = 1
              endif
            endif
  140       if (ifreq .ne. 1) ifqsw = 1
C
C           I1 COUNTS THE NUMBER OF BUSES WITH LOAD REPRESENTATION
C
            i1 = i1 + 1
            ldndxp(1, i1) = ldp
            ldndxp(2, i1) = ldq
            ldndxp(3, i1) = ifreq
            ldndxp(4, i1) = ipnt
  150       ityp1 = 0
            ityp2 = 0
            ityp3 = 0
            ityp4 = 0
            nitem = 0
C
C           CHECK FOR CONSTANT IMPEDANCE LOADS
C
            if (abs(temp(1)) .ge. .001 .or.
     &          abs(temp(2)) .ge. .001) then
              psum = psum + temp(1)
              qsum = qsum + temp(2)
              m = m + 1
              ityp1 = 1
              nitem = nitem + 1
              busldn(1, m) = temp(1)
              busldn(2, m) = temp(2)
              istsw = 1
            endif
C
C           CHECK FOR FREQUENCY SENSITIVE LOAD
C
            if (abs(temp(7)) .ge. .001 .or.
     &          abs(temp(8)) .ge. .001) then
C
C             BYPASS PSUM TALLY FOR SATURATION DATA
C
              if (ifreq .le. 3) then
                psum = psum + temp(7)
                qsum = qsum + temp(8)
              endif
              m = m + 1
              ityp2 = 2
              nitem = nitem + 1
              busldn(1, m) = temp(7)
              busldn(2, m) = temp(8)
            endif
C
C           CHECK FOR CONSTANT MVA LOADS
C
            if (abs(temp(5)) .ge. .001 .or.
     &          abs(temp(6)) .ge. .001) then
C
C             BYPASS PSUM TALLY FOR SATURATION DATA
C
              if (ifreq .le. 3) then
                psum = psum + temp(5)
                qsum = qsum + temp(6)
              endif
              m = m + 1
              ityp3 = 3
              nitem = nitem + 1
              busldn(1, m) = temp(5)
              busldn(2, m) = temp(6)
            endif
C
C           CHECK FOR CONSTANT CURRENT LOAD
C
            if (abs(temp(3)) .ge. .001 .or.
     &          abs(temp(4)) .ge. .001) then
              psum = psum + temp(3)
              qsum = qsum + temp(4)
              m = m + 1
              ityp4 = 4
              nitem = nitem + 1
              busldn(1, m) = temp(3)
              busldn(2, m) = temp(4)
            endif
            if (psum .lt. .999 .or. psum .gt. 1.001) then
              write (errbuf(1), 10220) buff80
              write (errbuf(2), 10230) psum
              call prterr('E', 2)
10220         format ('0', a80)
10230         format (5x, 'THIS CARD IS INCORRECT. SUM OF PER UNIT ',
     &         'REAL POWERS NOT EQUAL 1. ', 5x, f6.3)
              iabort = 1
            else
              psum = 1.0
            endif
            if (qsum .lt. .999 .or. qsum .gt. 1.001) then
              write (errbuf(1), 10220) buff80
              write (errbuf(2), 10240) qsum
              call prterr('E', 2)
10240         format (5x, 'THIS CARD IS INCORRECT. SUM OF PER UNIT ',
     &         'REACTIVE POWERS NOT EQUAL 1. ', 5x, f6.3)
              iabort = 1
            else
              qsum = 1.0
            endif
C
C           ZERO OUT LDIDXF BEFORE STORING DATA
C
            do k = 1, 6
              ldidxn(k, i3) = 0
            enddo
            ldidxn(1, i3) = ityp4
            ldidxn(2, i3) = ityp3
            ldidxn(3, i3) = ityp2
            ldidxn(4, i3) = ityp1
            ldidxn(5, i3) = nitem
            ldidxn(6, i3) = iecsl
            iecsl = iecsl + nitem
          endif
          i3 = i3 + 1
          kldin = kldint
        enddo
        kldid = kldidt
        kbusld = kbuslt
        lrep = i1
        if (lrep .eq. 0) ipwr = 1
        kshed = kbusld + 1
        keclst = kshed
      endif
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


program nc_colors
! @(#) Print out some ncurses colors.
!(LICENSE:PD)
!---------------------------------------------
! John S. Urban, 20141214
!---------------------------------------------
! based on a C program example by
!  Yoran Heling <projects@yorhel.nl>
!  Date: 2011-06-11
!  License: MIT
!  Web: http://dev.yorhel.nl/dump/nccolour
!---------------------------------------------
use M_ncurses
integer :: i

type color
  integer(C_SHORT)  :: code
  character(len=20) :: desc
end type color

type(color)         :: colours(0:8)
   colours(0) = color(-1_C_SHORT,    "default" )
   colours(1) = color(COLOR_BLACK,   "BLACK"   )
   colours(2) = color(COLOR_RED,     "RED"     )
   colours(3) = color(COLOR_GREEN,   "GREEN"   )
   colours(4) = color(COLOR_YELLOW,  "YELLOW"  )
   colours(5) = color(COLOR_BLUE,    "BLUE"    )
   colours(6) = color(COLOR_MAGENTA, "MAGENTA" )
   colours(7) = color(COLOR_CYAN,    "CYAN"    )
   colours(8) = color(COLOR_WHITE,   "WHITE"   )

  stdscr=initscr()
  if (.not.c_associated(stdscr)) then
     ierr=endwin()
     write(*,*)'could not open standard screen. input not a tty? ERRNO=',ierr
     write(*,'("Error initializing screen.")')
     stop 1
  endif

  if(.not.has_colors()) then
    write(*,'(a)')"This terminal does not support colors"
    stop 1
  endif
  ierr=start_color()

  ! this prevents ncurses from forcing a white-on-black color scheme,
  ! regardless of what scheme is used by the terminal.
  ierr=use_default_colors()

  ierr=attron(A_BOLD)
  ierr=mvaddstr(0, 0, " Color      FD  FB  FW  BG  BD    BFD BFB BFW BBG BBD"//C_NULL_CHAR)
  ierr=attroff(A_BOLD)

  do   i=0, size(colours)-1
    ierr=init_pair(int(i*6+1,C_SHORT), colours(i)%code, -1_C_SHORT)
    ierr=init_pair(int(i*6+2,C_SHORT), colours(i)%code, COLOR_BLACK)
    ierr=init_pair(int(i*6+3,C_SHORT), colours(i)%code, COLOR_WHITE)
    ierr=init_pair(int(i*6+4,C_SHORT), colours(i)%code, colours(i)%code)
    ierr=init_pair(int(i*6+5,C_SHORT), -1_C_SHORT, colours(i)%code)
    ierr=mvaddstr(i+2, 1, colours(i)%desc//C_NULL_CHAR)
    ierr=move(i+2, 8)
    call print_colorpart(0_C_LONG, i)
    call print_colorpart(A_BOLD, i)
  enddo

  ierr=mvaddstr(i+3, 2, "FD  = Front color on default background"//C_NULL_CHAR)
  ierr=mvaddstr(i+4, 2, "FB  = Front color on black background"//C_NULL_CHAR)
  ierr=mvaddstr(i+5, 2, "FW  = Front color on white background"//C_NULL_CHAR)
  ierr=mvaddstr(i+6, 2, "BG  = Front and background color"//C_NULL_CHAR)
  ierr=mvaddstr(i+7, 2, "BD  = Background color with default front color"//C_NULL_CHAR)
  ierr=mvaddstr(i+8, 2, "B?? = As above, with A_BOLD enabled"//C_NULL_CHAR)
  ierr=mvaddstr(i+10, 0, "Hit any key to exit."//C_NULL_CHAR)
  ierr=refresh()
  ierr=getch()
  ierr=endwin()
contains
subroutine print_colorpart(attr, i)
integer(C_LONG) :: attr
integer(C_INT)  :: i
!---------- WORKS TOO
!   ierr=  attron(ior(attr,COLOR_PAIR(i*6+2)))
!   ierr=  addstr("FB"//C_NULL_CHAR)
!   ierr=  attroff(ior(attr,COLOR_PAIR(i*6+2)))
!   ierr=  addstr("  "//C_NULL_CHAR)
!----------
!---------- JUST GET COLOR, NOT BOLD ???
!   ierr=  attron(COLOR_PAIR(i*6+2))
!   ierr=  addstr("FB"//C_NULL_CHAR)
!   ierr=  attroff(COLOR_PAIR(i*6+2))
!   ierr=  addstr("  "//C_NULL_CHAR)
!----------
!----------
   ierr=addstr("    "//C_NULL_CHAR)
   ierr=  attrset(ior(attr,COLOR_PAIR(i*6+1)))
   ierr=  addstr("FD"//C_NULL_CHAR)
   ierr=  attrset(attr)
   ierr=  addstr("  "//C_NULL_CHAR)
!----------
   ierr=  attrset(ior(attr,COLOR_PAIR(i*6+2)))
   ierr=  addstr("FB"//C_NULL_CHAR)
   ierr=  attrset(attr)
   ierr=  addstr("  "//C_NULL_CHAR)

   ierr=  attrset(ior(attr,COLOR_PAIR(i*6+3)))
   ierr=  addstr("FW"//C_NULL_CHAR)
   ierr=  attrset(attr)
   ierr=  addstr("  "//C_NULL_CHAR)

   ierr=  attrset(ior(attr,COLOR_PAIR(i*6+4)))
   ierr=  addstr("BG"//C_NULL_CHAR)
   ierr=  attrset(attr)
   ierr=  addstr("  "//C_NULL_CHAR)

   ierr=  attrset(ior(attr,COLOR_PAIR(i*6+5)))
   ierr=  addstr("BD"//C_NULL_CHAR)
   ierr=  attrset(attr)

end subroutine print_colorpart
end program nc_colors


text:
      subroutine qlbox9(Y,musq,Ires)
      implicit none
c  $I^{\{D=4-2 \epsilon\}}_4(0,p_2^2,p_3^2,m^2;s_{12},s_{23};0,0,0,m^2)$}

c                    [                       s12                ]
c                    [   0        0        - ---         0      ]
c                    [                        2                 ]
c                    [                                          ]
c                    [                       p2sq    msq - s23  ]
c                    [   0        0        - ----    ---------  ]
c                    [                        2          2      ]
c               y9 = [                                          ]
c                    [   s12     p2sq                msq - p3sq ]
c                    [ - ---   - ----        0       ---------- ]
c                    [    2       2                      2      ]
c                    [                                          ]
c                    [        msq - s23  msq - p3sq             ]
c                    [   0    ---------  ----------     msq     ]
c                    [            2          2                  ]
      include 'qlconstants.f'
      integer iep
      double complex Ires(-2:0),qllnrat,wlogt,wlog2,
     . dilog1,dilog2,qlLi2omrat,qlLi2omx2
      double precision Y(4,4),msq,musq,tabar,si,fac,mean,m3sqbar,mp2sq

      msq=Y(4,4)
      mean=sqrt(musq*msq)
      tabar=2d0*Y(2,4)
      si=2d0*Y(1,3)
      m3sqbar=2d0*Y(3,4)
      mp2sq=2d0*Y(2,3)
      fac=si*tabar

      wlogt=qllnrat(tabar,mean)
      wlog2=qllnrat(si,mp2sq)

      dilog1=qlLi2omx2(m3sqbar,tabar,mp2sq,msq)
      dilog2=qlLi2omrat(si,mp2sq)
      Ires(-2)=dcmplx(0.5d0)
      Ires(-1)=-wlogt-wlog2
      Ires( 0)=dilog1+2d0*dilog2+dcmplx(pisq/12d0)+(wlogt+wlog2)**2
      do iep=-2,0
      Ires(iep)=Ires(iep)/dcmplx(fac)
      enddo
      return
      end

ext: f


text:
!* Copyright (c) 2000, NVIDIA CORPORATION.  All rights reserved.
!*
!* Licensed under the Apache License, Version 2.0 (the "License");
!* you may not use this file except in compliance with the License.
!* You may obtain a copy of the License at
!*
!*     http://www.apache.org/licenses/LICENSE-2.0
!*
!* Unless required by applicable law or agreed to in writing, software
!* distributed under the License is distributed on an "AS IS" BASIS,
!* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!* See the License for the specific language governing permissions and
!* limitations under the License.
!
!	Test that array constructors are allowed as 'put' argument to
!	random_seed

	program pp
	integer n
	integer, allocatable :: s(:)
	integer result,expect
	data expect/1/

	call random_seed(size=n)

	allocate (s(1:n))
	s = (/ ( i, i=1,n ) /)
	call random_seed(put=s)

	call random_seed(put=(/ ( i, i=1,n ) /) )

	result = 1
	call check(result,expect,1)
	end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
      double precision verysmallnumber
      parameter(verysmallnumber=1d-9)

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


! PR fortran/44036
! { dg-do compile }
! { dg-options "-fopenmp" }
subroutine foo(a, b)
  integer, external :: a
  integer, external, pointer :: b
  integer, external :: c
  integer, external, pointer :: d
  integer :: x
  x = 6
!$omp parallel default(none) private (x)
  x = a(4)
!$omp end parallel
!$omp parallel default(none) private (x)	! { dg-error "enclosing 'parallel'" }
  x = b(5)					! { dg-error "not specified in" }
!$omp end parallel
!$omp parallel default(none) private (x)
  x = c(6)
!$omp end parallel
  d => a
!$omp parallel default(none) private (x)	! { dg-error "enclosing 'parallel'" }
  x = d(7)					! { dg-error "not specified in" }
!$omp end parallel
end


text:
module io_mod

use parameters_mod

contains

    subroutine output_vtk(Qin,nout,iam)

        implicit none
        real Qin(nx,ny,nz,nQ,nbasis)
        integer nout
        integer(I4P) , parameter :: nb=1*ngu,sprd=0*1

        real(R4P), dimension(nnx+1) :: x_xml_rect
        real(R4P), dimension(nny+1) :: y_xml_rect
        real(R4P), dimension(nnz+1) :: z_xml_rect
        real(R4P), dimension(nnx*nny*nnz) :: var_xml_val_x
        real(R4P), dimension(nnx*nny*nnz) :: var_xml_val_y
        real(R4P), dimension(nnx*nny*nnz) :: var_xml_val_z
        real(R4P), dimension(nnx,nny,nnz,nQ) :: qvtk
        real(R4P), dimension(nnx,nny,nnz) :: qvtk_dxvy,qvtk_dyvx
        real U, P, vx, vy, vz, dni, dxrh,dyrh,dxmy,dymx
        integer(I4P):: E_IO,i,j,k,l,num,iam,igrid,ir,jr,kr,ib,jb,kb,ieq
        character (50) :: out_name
        character (4) :: tname
        character (5) :: tname1
        character (4) :: pname
        character (5) :: pname1

        num=nout+10000

        write(tname1,'(i5)')num
        tname=tname1(2:5)
        tname = trim(tname)
        tname = adjustr(tname)

        num=iam+10000

        write(pname1,'(i5)')num
        pname=pname1(2:5)
        pname = trim(pname)
        pname = adjustr(pname)

        ! "outdir" is a global variable specifying the output directory
        out_name=''//outdir//'/perseus_p'//pname//'_t'//tname//'.vtr'

        if (iam .eq. print_mpi) print *, '  Data written to:  ', out_name
        out_name = trim(out_name)
        out_name = adjustr(out_name)

        E_IO = VTK_INI_XML(output_format = 'BINARY',              &
                           filename      = out_name, &
                           mesh_topology = 'RectilinearGrid',     &
                           nx1=1,nx2=nnx+1,ny1=1,ny2=nny+1,nz1=1,nz2=nnz+1)

        do i=1+nb,nnx-nb+1
            x_xml_rect(i-nb)=(xvtk(i) - 0.5*dxvtk)*1e-3
        enddo
        do j=1+nb,nny-nb+1
            y_xml_rect(j-nb)=(yvtk(j) - 0.5*dyvtk)*1e-3
        enddo
        do k=1+nb,nnz-nb+1
            z_xml_rect(k-nb)=(zvtk(k) - 0.5*dzvtk)*1e-3
        enddo

        E_IO = VTK_GEO_XML(nx1=1,nx2=nnx+1,ny1=1,ny2=nny+1,nz1=1,nz2=nnz+1, &
                             X=x_xml_rect,Y=y_xml_rect,Z=z_xml_rect)

        E_IO = VTK_DAT_XML(var_location     = 'cell', &
                           var_block_action = 'OPEN')

        do i=1+nb,nnx-nb

            ir = int(i/nvtk) + 1
            ib = i - nvtk*int(i/nvtk)
            if (ib .eq. 0) then
                ib = nvtk
                ir = ir - 1
            end if

            do j=1+nb,nny-nb

                jr = int(j/nvtk) + 1
                jb = j - nvtk*int(j/nvtk)
                if (jb .eq. 0) then
                    jb = nvtk
                    jr = jr - 1
                end if

                do k=1+nb,nnz-nb

                    kr = int(k/nvtk) + 1
                    kb = k - nvtk*int(k/nvtk)
                    if (kb .eq. 0) then
                        kb = nvtk
                        kr = kr - 1
                    end if
                    igrid = nvtk2*(ib-1) + nvtk*(jb-1) + kb
                    do ieq=1,nQ
                        qvtk(i,j,k,ieq) = sum(bfvtk(igrid,1:nbasis)*Qin(ir,jr,kr,ieq,1:nbasis))
                    end do

                    dxrh = sum(bfvtk_dx(igrid,1:nbasis)*Qin(ir,jr,kr,rh,1:nbasis))
                    dyrh = sum(bfvtk_dy(igrid,1:nbasis)*Qin(ir,jr,kr,rh,1:nbasis))
                    dxmy = sum(bfvtk_dx(igrid,1:nbasis)*Qin(ir,jr,kr,my,1:nbasis))
                    dymx = sum(bfvtk_dy(igrid,1:nbasis)*Qin(ir,jr,kr,mx,1:nbasis))
                    qvtk_dxvy(i,j,k) = (qvtk(i,j,k,rh)*dxmy - qvtk(i,j,k,my)*dxrh)/qvtk(i,j,k,rh)**2
                    qvtk_dyvx(i,j,k) = (qvtk(i,j,k,rh)*dymx - qvtk(i,j,k,mx)*dyrh)/qvtk(i,j,k,rh)**2
                end do
            end do
        end do

        do i=1+nb,nnx-nb
            do j=1+nb,nny-nb
                do k=1+nb,nnz-nb
                    l=(i-nb)+(j-nb-1)*(nnx-2*nb)+(k-nb-1)*(nnx-2*nb)*(nny-2*nb)
                    var_xml_val_x(l)=log(qvtk(i,j,k,rh)*n0)/log(10.)
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Log Density',                    &
                           var     = var_xml_val_x)

        do i=1+nb,nnx-nb
            do j=1+nb,nny-nb
                do k=1+nb,nnz-nb
                    l=(i-nb)+(j-nb-1)*(nnx-2*nb)+(k-nb-1)*(nnx-2*nb)*(nny-2*nb)
                    var_xml_val_x(l)=qvtk(i,j,k,rh)*n0
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Density',                    &
                           var     = var_xml_val_x)

        do i=1+nb,nnx-nb
            do j=1+nb,nny-nb
                do k=1+nb,nnz-nb
                    l=(i-nb)+(j-nb-1)*(nnx-2*nb)+(k-nb-1)*(nnx-2*nb)*(nny-2*nb)
                    dni = 1./qvtk(i,j,k,rh)
                    vx = qvtk(i,j,k,mx)*dni
                    vy = qvtk(i,j,k,my)*dni
                    vz = qvtk(i,j,k,mz)*dni
                    P = (aindex - 1.)*(qvtk(i,j,k,en) - 0.5*qvtk(i,j,k,rh)*(vx**2 + vy**2 + vz**2))
                    var_xml_val_x(l)=P*te0*dni
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Temperature',                    &
                           var     = var_xml_val_x)

       do i=1+nb,nnx-nb
           do j=1+nb,nny-nb
               do k=1+nb,nnz-nb
                   l=(i-nb)+(j-nb-1)*(nnx-2*nb)+(k-nb-1)*(nnx-2*nb)*(nny-2*nb)
                   dni = 1./qvtk(i,j,k,rh)
                   vx = qvtk(i,j,k,mx)*dni
                   vy = qvtk(i,j,k,my)*dni
                   vz = qvtk(i,j,k,mz)*dni
                   U = qvtk(i,j,k,en) - 0.5*qvtk(i,j,k,rh)*(vx**2 + vy**2 + vz**2)
                   P = (aindex - 1.)*U
                   var_xml_val_x(l)=P*(dni**aindm1)  ! Polytropic gas
                   ! var_xml_val_x(l)=(1./aindm1)*log(P/dni**aindm1)  ! I Do Like CFD
               enddo
           enddo
       enddo
       E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                          varname = 'Entropy density',                    &
                          var     = var_xml_val_x)

        do i=1+nb,nnx-nb
            do j=1+nb,nny-nb
                do k=1+nb,nnz-nb
                    l = (i-nb)+(j-nb-1)*(nnx-2*nb)+(k-nb-1)*(nnx-2*nb)*(nny-2*nb)
                    dni = 1./qvtk(i,j,k,rh)
                    vx = qvtk(i,j,k,mx)*dni
                    vy = qvtk(i,j,k,my)*dni
                    vz = qvtk(i,j,k,mz)*dni
                    ! from "viscosity" version
                    P = (aindex - 1.)*(qvtk(i,j,k,en) - 0.5*qvtk(i,j,k,rh)*(vx**2 + vy**2 + vz**2))
                    var_xml_val_x(l)=P
                    ! if(ieos .eq. 1) P = (aindex - 1.)*(qvtk(i,j,k,en) - 0.5*qvtk(i,j,k,rh)*(vx**2 + vy**2 + vz**2))
                    ! if(ieos .eq. 2)then
                    !     P = P_1*(qvtk(i,j,k,rh)**7.2 - 1.) + P_base
                    !     if(P < P_floor) P= P_floor
                    ! end if
                    ! var_xml_val_x(l) = P*P0
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'pressure',                    &
                           var     = var_xml_val_x)

        do i=1+nb,nnx-nb
            do j=1+nb,nny-nb
                do k=1+nb,nnz-nb
                    l=(i-nb)+(j-nb-1)*(nnx-2*nb)+(k-nb-1)*(nnx-2*nb)*(nny-2*nb)
                    var_xml_val_x(l)=v0*qvtk(i,j,k,mx)/qvtk(i,j,k,rh)
                    var_xml_val_y(l)=v0*qvtk(i,j,k,my)/qvtk(i,j,k,rh)
                    var_xml_val_z(l)=v0*qvtk(i,j,k,mz)/qvtk(i,j,k,rh)
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Velocity',                    &
                           varX     = var_xml_val_x, varY     = var_xml_val_y, varZ     = var_xml_val_Z )

        do i=1+nb,nnx-nb
            do j=1+nb,nny-nb
                do k=1+nb,nnz-nb
                    l=(i-nb)+(j-nb-1)*(nnx-2*nb)+(k-nb-1)*(nnx-2*nb)*(nny-2*nb)
                        var_xml_val_x(l)= -2.*(qvtk_dxvy(i,j,k)-qvtk_dyvx(i,j,k))*dxi/t0
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Vorticity',                    &
                           var     = var_xml_val_x)

        ! from "viscosity" version
        do i=1+nb,nnx-nb
            do j=1+nb,nny-nb
                do k=1+nb,nnz-nb
                    l=(i-nb)+(j-nb-1)*(nx-2*nb)+(k-nb-1)*(nx-2*nb)*(ny-2*nb)
                        var_xml_val_x(l)= -2*(Q_r2(i,j,k,my,2)-Q_r2(i,j,k,mx,3))/Q_r2(i,j,k,rh,1)*dxi/t0
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                          varname = 'Vorticity2',                    &
                          var     = var_xml_val_x)

        E_IO = VTK_DAT_XML(var_location     = 'cell', &
                           var_block_action = 'Close')
        E_IO = VTK_GEO_XML()
        E_IO = VTK_END_XML()

    end subroutine output_vtk

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

    subroutine output_vtk0(Qin,nout,iam)

        implicit none
        real Qin(nx,ny,nz,nQ,nbasis)
        integer nout
        integer(I4P) , parameter :: nb=1*ngu,sprd=0*1
        !    integer(I4P), parameter:: nnx=nx-2*nb,nny=ny-2*nb,nnz=nz-2*nb
        real(R4P), dimension(nnx+1):: x_xml_rect
        real(R4P), dimension(nny+1):: y_xml_rect
        real(R4P), dimension(nnz+1):: z_xml_rect
        real(R4P), dimension(nnx*nny*nnz):: var_xml_val_x
        real(R4P), dimension(nnx*nny*nnz):: var_xml_val_y
        real(R4P), dimension(nnx*nny*nnz):: var_xml_val_z
        real P, vx, vy, vz,  dni
        integer(I4P):: E_IO,i,j,k,l,num,iam
        character (50) :: out_name
        character (4) :: tname
        character (5) :: tname1
        character (4) :: pname
        character (5) :: pname1

        num=nout+10000

        write(tname1,'(i5)')num
        tname=tname1(2:5)
        tname = trim(tname)
        tname = adjustr(tname)

        num=iam+10000

        write(pname1,'(i5)')num
        pname=pname1(2:5)
        pname = trim(pname)
        pname = adjustr(pname)
        out_name='data/data3/perseus_p'//pname//'_t'//tname//'.vtr'
        ! out_name='data/data4/perseus_p'//pname//'_t'//tname//'.vtr'
        print *, '  Data written to:  ', out_name
        out_name = trim(out_name)
        out_name = adjustr(out_name)

        E_IO = VTK_INI_XML(output_format = 'BINARY',              &
                           filename      = out_name, &
                           mesh_topology = 'RectilinearGrid',     &
                           nx1=1,nx2=nnx+1,ny1=1,ny2=nny+1,nz1=1,nz2=nnz+1)

        do i=1+nb,nx-nb+1
            x_xml_rect(i)=(xc(i) - 0.5/dxi)*1e-3
        enddo
        do i=1+nb,ny-nb+1
            y_xml_rect(i)=(yc(i) - 0.5/dyi)*1e-3
        enddo
        do i=1+nb,nz-nb+1
            z_xml_rect(i)=(zc(i) - 0.5/dzi)*1e-3
        enddo

        E_IO = VTK_GEO_XML(nx1=1,nx2=nnx+1,ny1=1,ny2=nny+1,nz1=1,nz2=nnz+1, &
                             X=x_xml_rect,Y=y_xml_rect,Z=z_xml_rect)

        E_IO = VTK_DAT_XML(var_location     = 'cell', &
                           var_block_action = 'OPEN')

        do i=1+nb,nx-nb
            do j=1+nb,ny-nb
                do k=1+nb,nz-nb
                    l=(i-nb)+(j-nb-1)*(nx-2*nb)+(k-nb-1)*(nx-2*nb)*(ny-2*nb)
                    var_xml_val_x(l)=log(Qin(i,j,k,rh,1)*n0)/log(10.)

                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Log Density',                    &
                           var     = var_xml_val_x)

        do i=1+nb,nx-nb
            do j=1+nb,ny-nb
                do k=1+nb,nz-nb
                    l=(i-nb)+(j-nb-1)*(nx-2*nb)+(k-nb-1)*(nx-2*nb)*(ny-2*nb)
                    var_xml_val_x(l)=Qin(i,j,k,rh,1)*n0
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Density',                    &
                           var     = var_xml_val_x)

        do i=1+nb,nx-nb
            do j=1+nb,ny-nb
                do k=1+nb,nz-nb
                    l=(i-nb)+(j-nb-1)*(nx-2*nb)+(k-nb-1)*(nx-2*nb)*(ny-2*nb)
                    dni = 1./Qin(i,j,k,rh,1)
                    vx = Qin(i,j,k,mx,1)*dni
                    vy = Qin(i,j,k,my,1)*dni
                    vz = Qin(i,j,k,mz,1)*dni
                    P = (aindex - 1)*(Qin(i,j,k,en,1) - 0.5*Qin(i,j,k,rh,1)*(vx**2 + vy**2 + vz**2))*dni
                    var_xml_val_x(l)=P*te0
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Temperature',                    &
                           var     = var_xml_val_x)

        do i=1+nb,nx-nb
            do j=1+nb,ny-nb
                do k=1+nb,nz-nb
                    l=(i-nb)+(j-nb-1)*(nx-2*nb)+(k-nb-1)*(nx-2*nb)*(ny-2*nb)
                    dni = 1./Qin(i,j,k,rh,1)
                    vx = Qin(i,j,k,mx,1)*dni
                    vy = Qin(i,j,k,my,1)*dni
                    vz = Qin(i,j,k,mz,1)*dni

                    ! from "viscosity" version
                    P = (aindex - 1)*(Qin(i,j,k,en,1) - 0.5*Qin(i,j,k,rh,1)*(vx**2 + vy**2 + vz**2))
                    ! if(ieos .eq. 1)P = (aindex - 1)*(Qin(i,j,k,en,1) - 0.5*Qin(i,j,k,rh,1)*(vx**2 + vy**2 + vz**2))
                    var_xml_val_x(l)=P
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'pressure',                    &
                           var     = var_xml_val_x)

        do i=1+nb,nx-nb
            do j=1+nb,ny-nb
                do k=1+nb,nz-nb
                    l=(i-nb)+(j-nb-1)*(nx-2*nb)+(k-nb-1)*(nx-2*nb)*(ny-2*nb)
                    var_xml_val_x(l)=v0*Qin(i,j,k,mx,1)/Qin(i,j,k,rh,1)
                    var_xml_val_y(l)=v0*Qin(i,j,k,my,1)/Qin(i,j,k,rh,1)
                    var_xml_val_z(l)=v0*Qin(i,j,k,mz,1)/Qin(i,j,k,rh,1)
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Velocity',                    &
                           varX     = var_xml_val_x, varY     = var_xml_val_y, varZ     = var_xml_val_Z )

        do i=1+nb,nx-nb
            do j=1+nb,ny-nb
                do k=1+nb,nz-nb
                    l=(i-nb)+(j-nb-1)*(nx-2*nb)+(k-nb-1)*(nx-2*nb)*(ny-2*nb)
                        var_xml_val_x(l)= -2*(Q_r2(i,j,k,my,2)-Q_r2(i,j,k,mx,3))/Q_r2(i,j,k,rh,1)*dxi/t0
                enddo
            enddo
        enddo
        E_IO = VTK_VAR_XML(NC_NN   = nnx*nny*nnz, &
                           varname = 'Vorticity',                    &
                           var     = var_xml_val_x)

        E_IO = VTK_DAT_XML(var_location     = 'cell', &
                           var_block_action = 'Close')
        E_IO = VTK_GEO_XML()
        E_IO = VTK_END_XML()

    end subroutine output_vtk0

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

    subroutine writeQ(fprefix,irank,iddump,Qin,tnow,dtnow,noutnow,mpi_nxnow,mpi_nynow,mpi_nznow)

        implicit none
        real :: Qin(nx,ny,nz,nQ,nbasis),tnow,dtnow
        integer :: irank,iddump,noutnow,mpi_nxnow,mpi_nynow,mpi_nznow,nump,numd,qq,k,j,i,ir
        character (4) :: fprefix,pname,dname
        character (5) :: pname1,dname1
        character (30) :: fname2

        nump = iam + 10000

        write(pname1,'(i5)')nump
        pname=pname1(2:5)
        pname = trim(pname)
        pname = adjustr(pname)

        numd = iddump + 10000

        write(dname1,'(i5)')numd
        dname=dname1(2:5)
        dname = trim(dname)
        dname = adjustr(dname)

        fname2 = 'data/'//fprefix//'_p'//pname//'_d'//dname//'.dat'
        ! print *,'fname2 ',fname2

        open(unit=3,file=fname2)

        ! open(unit = 10, file = 'data/perseus_t'//dname//'_p'//pname//'.bin',form = 'unformatted',access = 'stream')

        do ir=1,nbasis
            do qq=1,nQ
                do k=1,nz
                    do j=1,ny
                        write(3,*) (Qin(i,j,k,qq,ir),i=1,nx)
                    enddo
                enddo
            enddo
        enddo

        write(3,*) tnow,dtnow,noutnow,mpi_nxnow,mpi_nynow,mpi_nznow
        close(3)

    end subroutine writeQ

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

    subroutine readQ(fprefix,irank,iddump,Qin,tnow,dtnow,noutnow,mpi_nxnow,mpi_nynow,mpi_nznow)

        implicit none
        real :: Qin(nx,ny,nz,nQ,nbasis),tnow,dtnow
        integer :: irank,iddump,noutnow,mpi_nxnow,mpi_nynow,mpi_nznow,nump,numd,qq,k,j,i,ir
        character (4) :: fprefix,pname,dname
        character (5) :: pname1,dname1
        character (30) :: fname2

        nump = irank + 10000

        write(pname1,'(i5)')nump
        pname=pname1(2:5)
        pname = trim(pname)
        pname = adjustr(pname)

        numd = iddump + 10000

        write(dname1,'(i5)')numd
        dname=dname1(2:5)
        dname = trim(dname)
        dname = adjustr(dname)

        fname2 = 'data/'//fpre//'_p'//pname//'_d'//dname//'.dat'

        open(unit=3,file=fname2,action='read')

        do ir=1,nbasis
        do qq=1,nQ
            do k=1,nz
                do j=1,ny
                    read(3,*) (Qin(i,j,k,qq,ir),i=1,nx)
                enddo
            enddo
        enddo
        enddo

        read(3,*) tnow,dtnow,noutnow,mpi_nxnow,mpi_nynow,mpi_nznow
        close(3)

    end subroutine readQ

end module io_mod

max_stars_repo_name: sseyler/HERMESHD

max_issues_repo_name: sseyler/HERMESHD

max_forks_repo_name: sseyler/HERMESHD


text:
// Generated Files
${PROJ_DIR}/llink/script/premade_examples/axi_st_d256_gen1_gen2/axi_st_d256_gen1_gen2_slave_top.sv
${PROJ_DIR}/llink/script/premade_examples/axi_st_d256_gen1_gen2/axi_st_d256_gen1_gen2_slave_concat.sv
${PROJ_DIR}/llink/script/premade_examples/axi_st_d256_gen1_gen2/axi_st_d256_gen1_gen2_slave_name.sv

// Logic Link files
-f ${PROJ_DIR}/llink/rtl/llink.f

// Common Files
-f ${PROJ_DIR}/common/rtl/common.f

max_stars_repo_name: chipsalliance/aib-protocols

max_issues_repo_name: chipsalliance/aib-protocols

max_forks_repo_name: chipsalliance/aib-protocols


text:
! { dg-do run }
! Option passed to avoid excess errors from obsolete warning
! { dg-options "-w" }
! PR18827
      integer i,j
      common /foo/ i,j
      assign 1000 to j
      j = 5
      goto j
 1000 continue
      end

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


MODULE mymodule

 IMPLICIT NONE
 CONTAINS

  FUNCTION running_max(arg)
    INTEGER, INTENT(IN) :: arg
    INTEGER :: running_max
    INTEGER :: current_max = 0

    current_max = MAX(current_max, arg)
    running_max = current_max
  END FUNCTION running_max

END MODULE mymodule

PROGRAM driver
  USE mymodule
  IMPLICIT NONE

  PRINT *, running_max(1)
  PRINT *, running_max(10)
  PRINT *, running_max(2)
  PRINT *, running_max(20)

END PROGRAM driver


text:
submodule(state_interface) state_implementation
  use kind_parameters, only : rkind
  implicit none

contains

  module procedure define
    use assertions_interface, only : assert, max_errmsg_len
    real(rkind) :: temperature_, pressure_
    namelist/state_list/ temperature_, pressure_

    this%time_ = time
    this%burn_depth_ = burn_depth

    block
      integer, parameter :: success = 0
      character(len=max_errmsg_len) error_message
      integer :: io_status, file_unit

      open(newunit=file_unit, file=input_file, status="old", iostat=io_status, iomsg=error_message)
      call assert(io_status == success, "state_t%define: io_status == success", diagnostic_data = error_message)
      read(file_unit, nml=state_list)
      close(file_unit)
    end block

    associate(V => (volume), T => (temperature_), p => (pressure_))
      this%mass_ = p*V/(R_gas*T)
      this%energy_ = this%mass_*c_v*T
    end associate
  end procedure

  module procedure absolute
    this_abs = abs([this%mass_, this%energy_, this%burn_depth_])
  end procedure

  module procedure new_state
    new_state%mass_ = mass
    new_state%energy_ = energy
    new_state%burn_depth_= burn_depth
    new_state%time_ = time
  end procedure

  module procedure row_vector
    row_vector = reshape([this%time_, this%mass_, this%energy_, this%burn_depth_], [1,4])
  end procedure

  module procedure time
    real(rkind) time
    time = this%time_
  end procedure

  module procedure burn_depth
    burn_depth = this%burn_depth_
  end procedure

  module procedure energy
    energy = this%energy_
  end procedure

  module procedure mass
    mass = this%mass_
  end procedure

  module procedure add
    total%time_       = lhs%time_       + rhs%time_
    total%burn_depth_ = lhs%burn_depth_ + rhs%burn_depth_
    total%mass_       = lhs%mass_       + rhs%mass_
    total%energy_     = lhs%energy_    + rhs%energy_
  end procedure

  module procedure multiply
    product_%time_       = lhs%time_       * rhs
    product_%burn_depth_ = lhs%burn_depth_ * rhs
    product_%mass_       = lhs%mass_       * rhs
    product_%energy_     = lhs%energy_     * rhs
  end procedure

end submodule

max_stars_repo_name: sourceryinstitute/rocket-science

max_issues_repo_name: sourceryinstitute/rocket-science

max_forks_repo_name: sourceryinstitute/rocket-science


text:
!****h* root/fortran/test/tH5L_F03.f90
!
! NAME
!  tH5L_F03.f90
!
! FUNCTION
!  Test FORTRAN HDF5 H5L APIs which are dependent on FORTRAN 2003
!  features.
!
! COPYRIGHT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   Copyright by The HDF Group.                                               *
!   Copyright by the Board of Trustees of the University of Illinois.         *
!   All rights reserved.                                                      *
!                                                                             *
!   This file is part of HDF5.  The full HDF5 copyright notice, including     *
!   terms governing use, modification, and redistribution, is contained in    *
!   the files COPYING and Copyright.html.  COPYING can be found at the root   *
!   of the source code distribution tree; Copyright.html can be found at the  *
!   root level of an installed copy of the electronic HDF5 document set and   *
!   is linked from the top-level documents page.  It can also be found at     *
!   http://hdfgroup.org/HDF5/doc/Copyright.html.  If you do not have          *
!   access to either file, you may request a copy from help@hdfgroup.org.     *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
! USES
!  liter_cb_mod
!
! CONTAINS SUBROUTINES
!  test_iter_group
!
!*****

MODULE liter_cb_mod

  USE HDF5
  USE ISO_C_BINDING
  IMPLICIT NONE

  TYPE iter_enum
     INTEGER RET_ZERO
     INTEGER RET_TWO
     INTEGER RET_CHANGE
     INTEGER RET_CHANGE2
  END TYPE iter_enum

  ! Custom group iteration callback data
  TYPE, bind(c) ::  iter_info
     CHARACTER(LEN=1), DIMENSION(1:10) :: name !  The name of the object
     INTEGER(c_int) :: TYPE    !  The TYPE of the object
     INTEGER(c_int) :: command ! The TYPE of RETURN value
  END TYPE iter_info

CONTAINS

!***************************************************************
!**
!**  liter_cb(): Custom link iteration callback routine.
!**
!***************************************************************

  INTEGER FUNCTION liter_cb(group, name, link_info, op_data) bind(C)

    USE HDF5
    USE ISO_C_BINDING
    IMPLICIT NONE

    INTEGER(HID_T), VALUE :: group
    CHARACTER(LEN=1), DIMENSION(1:10) :: name

    TYPE (H5L_info_t) :: link_info

    TYPE(iter_info) :: op_data

    INTEGER, SAVE :: count
    INTEGER, SAVE :: count2

!!$
!!$    iter_info *info = (iter_info *)op_data;
!!$    static int count = 0;
!!$    static int count2 = 0;

    op_data%name(1:10) = name(1:10)

    SELECT CASE (op_data%command)

    CASE(0)
       liter_cb = 0
    CASE(2)
       liter_cb = 2
    CASE(3)
       count = count + 1
       IF(count.GT.10) THEN
          liter_cb = 1
       ELSE
          liter_cb = 0
       ENDIF
    CASE(4)
       count2 = count2 + 1
       IF(count2.GT.10) THEN
          liter_cb = 1
       ELSE
          liter_cb = 0
       ENDIF
    END SELECT

  END FUNCTION liter_cb
END MODULE liter_cb_mod

! *****************************************
! ***        H 5 L   T E S T S
! *****************************************

!***************************************************************
!**
!**  test_iter_group(): Test group iteration functionality
!**
!***************************************************************
SUBROUTINE test_iter_group(total_error)

  USE HDF5
  USE ISO_C_BINDING
  USE liter_cb_mod
  IMPLICIT NONE

  INTEGER, INTENT(INOUT) :: total_error
  INTEGER(HID_T) :: fapl
  INTEGER(HID_T) :: file !  File ID
  INTEGER(hid_t) :: dataset          ! Dataset ID
  INTEGER(hid_t) :: datatype         ! Common datatype ID
  INTEGER(hid_t) :: filespace        ! Common dataspace ID
  INTEGER(hid_t) :: root_group,grp   ! Root group ID
  INTEGER i,j                  ! counting variable
  INTEGER(hsize_t) idx            ! Index in the group
  CHARACTER(LEN=11) :: DATAFILE = "titerate.h5"
  INTEGER, PARAMETER :: ndatasets = 50
  CHARACTER(LEN=10) :: name !  temporary name buffer
  CHARACTER(LEN=10), DIMENSION(1:ndatasets+2) :: lnames !  Names of the links created
!!$    char dataset_name[NAMELEN];   dataset name

  TYPE(iter_info), TARGET :: info

!!$    iter_info info;          Custom iteration information
!!$    H5G_info_t ginfo;        Buffer for querying object's info
!!$    herr_t ret;		     Generic return value

  INTEGER :: error
  INTEGER :: ret_value
  TYPE(C_PTR) :: f_ptr
  TYPE(C_FUNPTR) :: f1
  TYPE(C_PTR) :: f2
  CHARACTER(LEN=2) :: ichr2
  CHARACTER(LEN=10) :: ichr10

  ! Get the default FAPL
  CALL H5Pcreate_f(H5P_FILE_ACCESS_F, fapl, error)
  CALL check("h5pcreate_f", error, total_error)

  ! Set the "use the latest version of the format" bounds for creating objects in the file
  CALL H5Pset_libver_bounds_f(fapl, H5F_LIBVER_LATEST_F, H5F_LIBVER_LATEST_F, error)
  CALL check("H5Pset_libver_bounds_f",error, total_error)

  ! Create the test file with the datasets
  CALL h5fcreate_f(DATAFILE, H5F_ACC_TRUNC_F, file, error, H5P_DEFAULT_F, fapl)
  CALL check("h5fcreate_f", error, total_error)

  ! Test iterating over empty group
  idx = 0
  info%command = 0
  f1 = C_FUNLOC(liter_cb)
  f2 = C_LOC(info)

  CALL H5Literate_f(file, H5_INDEX_NAME_F, H5_ITER_INC_F, idx, f1, f2, ret_value, error)
  CALL check("H5Literate_f", error, total_error)

  CALL H5Tcopy_f(H5T_NATIVE_INTEGER, datatype, error)
  CALL check("H5Tcopy_f", error, total_error)

  CALL H5Screate_f(H5S_SCALAR_F, filespace, error)
  CALL check("H5Screate_f", error, total_error)

  DO i = 1, ndatasets
     WRITE(ichr2, '(I2.2)') i

     name = 'Dataset '//ichr2

     CALL h5dcreate_f(file, name, datatype, filespace, dataset, error)
     CALL check("H5dcreate_f", error, total_error)

     lnames(i) = name

     CALL h5dclose_f(dataset,error)
     CALL check("H5dclose_f", error, total_error)

  ENDDO

  !  Create a group and named datatype under root group for testing

  CALL H5Gcreate_f(file, "grp0000000", grp, error)
  CALL check("H5Gcreate_f", error, total_error)

  lnames(ndatasets+2) = "grp0000000"

!!$
!!$    lnames[NDATASETS] = HDstrdup("grp");
!!$    CHECK(lnames[NDATASETS], NULL, "strdup");
!!$

  CALL H5Tcommit_f(file, "dtype00000", datatype, error)
  CALL check("H5Tcommit_f", error, total_error)

  lnames(ndatasets+1) = "dtype00000"

  !  Close everything up

  CALL H5Tclose_f(datatype, error)
  CALL check("H5Tclose_f", error, total_error)

  CALL H5Gclose_f(grp, error)
  CALL check("H5Gclose_f", error, total_error)

  CALL H5Sclose_f(filespace, error)
  CALL check("H5Sclose_f", error, total_error)

  CALL H5Fclose_f(file, error)
  CALL check("H5Fclose_f", error, total_error)

  !  Iterate through the datasets in the root group in various ways
  CALL H5Fopen_f(DATAFILE, H5F_ACC_RDONLY_F, file, error, access_prp=fapl)
  CALL check("h5fopen_f", error, total_error)

  ! Test all objects in group, when callback always returns 0
  info%command = 0
  idx = 0
  CALL H5Literate_f(file, H5_INDEX_NAME_F, H5_ITER_INC_F, idx, f1, f2, ret_value, error)
  IF(ret_value.GT.0)THEN
     PRINT*,"ERROR: Group iteration function didn't return zero correctly!"
     CALL verify("H5Literate_f", error, -1, total_error)
  ENDIF

  !   Test all objects in group, when callback always returns 1
  !   This also tests the "restarting" ability, because the index changes

  info%command = 2
  idx = 0
  i = 0
  f1 = C_FUNLOC(liter_cb)
  f2 = C_LOC(info)
  DO
     CALL H5Literate_f(file, H5_INDEX_NAME_F, H5_ITER_INC_F, idx, f1, f2, ret_value, error)
     IF(error.LT.0) EXIT
     !  Verify return value from iterator gets propagated correctly
     CALL VERIFY("H5Literate", ret_value, 2, total_error)
     !  Increment the number of times "2" is returned
     i = i + 1
     ! Verify that the index is the correct value
     CALL VERIFY("H5Literate", INT(idx), INT(i), total_error)
     IF(idx .GT.ndatasets+2)THEN
        PRINT*,"ERROR: Group iteration function walked too far!"
     ENDIF

     ! Verify the correct name is retrieved
     DO j = 1, 10
        ichr10(j:j) = info%name(j)(1:1)
     ENDDO
     CALL verifystring("H5Literate_f", ichr10, lnames(INT(idx)), total_error)
     IF(i.EQ.52)EXIT ! prints out error message otherwise (for gcc/gfortran/g95) not intel (why) -FIXME- scot
  END DO

  ! put check if did not walk far enough -scot FIXME

  IF(i .NE. (NDATASETS + 2)) THEN
     CALL VERIFY("H5Literate_f", i, INT(NDATASETS + 2), total_error)
     PRINT*,"ERROR: Group iteration function didn't perform multiple iterations correctly"
  ENDIF

  ! Test all objects in group, when callback changes return value
  ! This also tests the "restarting" ability, because the index changes

  info%command = 3
  idx = 0
  i = 0

  f1 = C_FUNLOC(liter_cb)
  f2 = C_LOC(info)
  DO

     CALL H5Literate_f(file, H5_INDEX_NAME_F, H5_ITER_INC_F, idx, f1, f2, ret_value, error)
     IF(error.LT.0) EXIT
     CALL VERIFY("H5Literate_f", ret_value, 1, total_error)

     ! Increment the number of times "1" is returned
     i = i + 1

     ! Verify that the index is the correct value
     CALL VERIFY("H5Literate_f", INT(idx), INT(i+10), total_error)

     IF(idx .GT.ndatasets+2)THEN
        PRINT*,"Group iteration function walked too far!"
     ENDIF

     DO j = 1, 10
        ichr10(j:j) = info%name(j)(1:1)
     ENDDO
     ! Verify that the correct name is retrieved
     CALL verifystring("H5Literate_f", ichr10, lnames(INT(idx)), total_error)
     IF(i.EQ.42)EXIT ! prints out error message otherwise (for gcc/gfortran/g95) not intel (why) -FIX- scot
  ENDDO

  IF(i .NE. 42 .OR. idx .NE. 52)THEN
     PRINT*,"ERROR: Group iteration function didn't perform multiple iterations correctly!"
     CALL check("H5Literate_f",-1,total_error)
  ENDIF

  CALL H5Fclose_f(file, error)
  CALL check("H5Fclose_f", error, total_error)

END SUBROUTINE test_iter_group

max_stars_repo_name: pbasting/cactus

max_issues_repo_name: pbasting/cactus

max_forks_repo_name: pbasting/cactus


text:
      MODULE vname1
      USE Vname0
      REAL(rprec), PARAMETER :: one = 1, zero = 0
      INTEGER :: mpol, mpol1, mrho, nfp, ntheta, nphi, nphi2, n2,
     1           mpnt, ntor
      INTEGER :: mm(mu), nn(nv+1), m1(mnd), n1(mnd)
      REAL(rprec), DIMENSION(nv) :: r0n, z0n, raxis, zaxis
      REAL(rprec), DIMENSION(0:mu) :: dm1
      REAL(rprec), DIMENSION(0:mu+2) :: t1m, t2m
      REAL(rprec), DIMENSION(0:mu,2) :: xmpq
      REAL(rprec), DIMENSION(:), ALLOCATABLE :: xvec, gvec,
     1   xdot, xstore
      REAL(rprec) :: gnorm, specw, delt, deltf, dnorm, elongate,
     1   twopi, r10, HB_Parameter
      END MODULE vname1


text:
      SUBROUTINE MB03WA( WANTQ, WANTZ, N1, N2, A, LDA, B, LDB, Q, LDQ,
     $                   Z, LDZ, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To swap adjacent diagonal blocks A11*B11 and A22*B22 of size
C     1-by-1 or 2-by-2 in an upper (quasi) triangular matrix product
C     A*B by an orthogonal equivalence transformation.
C
C     (A, B) must be in periodic real Schur canonical form (as returned
C     by SLICOT Library routine MB03XP), i.e., A is block upper
C     triangular with 1-by-1 and 2-by-2 diagonal blocks, and B is upper
C     triangular.
C
C     Optionally, the matrices Q and Z of generalized Schur vectors are
C     updated.
C
C         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)',
C         Z(in) * B(in) * Q(in)' = Z(out) * B(out) * Q(out)'.
C
C     This routine is largely based on the LAPACK routine DTGEX2
C     developed by Bo Kagstrom and Peter Poromaa.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     WANTQ   LOGICAL
C             Indicates whether or not the user wishes to accumulate
C             the matrix Q as follows:
C             = .TRUE. :  The matrix Q is updated;
C             = .FALSE.:  the matrix Q is not required.
C
C     WANTZ   LOGICAL
C             Indicates whether or not the user wishes to accumulate
C             the matrix Z as follows:
C             = .TRUE. :  The matrix Z is updated;
C             = .FALSE.:  the matrix Z is not required.
C
C     Input/Output Parameters
C
C     N1      (input) INTEGER
C             The order of the first block A11*B11. N1 = 0, 1 or 2.
C
C     N2      (input) INTEGER
C             The order of the second block A22*B22. N2 = 0, 1 or 2.
C
C     A       (input/output) DOUBLE PRECISION array, dimension
C             (LDA,N1+N2)
C             On entry, the leading (N1+N2)-by-(N1+N2) part of this
C             array must contain the matrix A.
C             On exit, the leading (N1+N2)-by-(N1+N2) part of this array
C             contains the matrix A of the reordered pair.
C
C     LDA     INTEGER
C             The leading dimension of the array A. LDA >= MAX(1,N1+N2).
C
C     B       (input/output) DOUBLE PRECISION array, dimension
C             (LDB,N1+N2)
C             On entry, the leading (N1+N2)-by-(N1+N2) part of this
C             array must contain the matrix B.
C             On exit, the leading (N1+N2)-by-(N1+N2) part of this array
C             contains the matrix B of the reordered pair.
C
C     LDB     INTEGER
C             The leading dimension of the array B. LDB >= MAX(1,N1+N2).
C
C     Q       (input/output) DOUBLE PRECISION array, dimension
C             (LDQ,N1+N2)
C             On entry, if WANTQ = .TRUE., the leading
C             (N1+N2)-by-(N1+N2) part of this array must contain the
C             orthogonal matrix Q.
C             On exit, the leading (N1+N2)-by-(N1+N2) part of this array
C             contains the updated matrix Q. Q will be a rotation
C             matrix for N1=N2=1.
C             This array is not referenced if WANTQ = .FALSE..
C
C     LDQ     INTEGER
C             The leading dimension of the array Q. LDQ >= 1.
C             If WANTQ = .TRUE., LDQ >= N1+N2.
C
C     Z       (input/output) DOUBLE PRECISION array, dimension
C             (LDZ,N1+N2)
C             On entry, if WANTZ = .TRUE., the leading
C             (N1+N2)-by-(N1+N2) part of this array must contain the
C             orthogonal matrix Z.
C             On exit, the leading (N1+N2)-by-(N1+N2) part of this array
C             contains the updated matrix Z. Z will be a rotation
C             matrix for N1=N2=1.
C             This array is not referenced if WANTZ = .FALSE..
C
C     LDZ     INTEGER
C             The leading dimension of the array Z. LDZ >= 1.
C             If WANTZ = .TRUE., LDZ >= N1+N2.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             = 1:  the transformed matrix (A, B) would be
C                   too far from periodic Schur form; the blocks are
C                   not swapped and (A,B) and (Q,Z) are unchanged.
C
C     METHOD
C
C     In the current code both weak and strong stability tests are
C     performed. The user can omit the strong stability test by changing
C     the internal logical parameter WANDS to .FALSE.. See ref. [2] for
C     details.
C
C     REFERENCES
C
C     [1] Kagstrom, B.
C         A direct method for reordering eigenvalues in the generalized
C         real Schur form of a regular matrix pair (A,B), in M.S. Moonen
C         et al (eds.), Linear Algebra for Large Scale and Real-Time
C         Applications, Kluwer Academic Publ., 1993, pp. 195-218.
C
C     [2] Kagstrom, B., and Poromaa, P.
C         Computing eigenspaces with specified eigenvalues of a regular
C         matrix pair (A, B) and condition estimation: Theory,
C         algorithms and software, Numer. Algorithms, 1996, vol. 12,
C         pp. 369-407.
C
C     CONTRIBUTORS
C
C     D. Kressner, Technical Univ. Berlin, Germany, and
C     P. Benner, Technical Univ. Chemnitz, Germany, December 2003.
C
C     REVISIONS
C
C     V. Sima, May 2008 (SLICOT version of the HAPACK routine DTGPX2).
C
C     KEYWORDS
C
C     Eigenvalue, periodic Schur form, reordering
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
      DOUBLE PRECISION   TEN
      PARAMETER          ( TEN = 1.0D+01 )
      INTEGER            LDST
      PARAMETER          ( LDST = 4 )
      LOGICAL            WANDS
      PARAMETER          ( WANDS = .TRUE. )
C     .. Scalar Arguments ..
      LOGICAL            WANTQ, WANTZ
      INTEGER            INFO, LDA, LDB, LDQ, LDZ, N1, N2
C     .. Array Arguments ..
      DOUBLE PRECISION   A(LDA,*), B(LDB,*), Q(LDQ,*), Z(LDZ,*)
C     .. Local Scalars ..
      LOGICAL            DTRONG, WEAK
      INTEGER            I, LINFO, M
      DOUBLE PRECISION   BQRA21, BRQA21, DDUM, DNORM, DSCALE, DSUM, EPS,
     $                   F, G, SA, SB, SCALE, SMLNUM, SS, THRESH, WS
C     .. Local Arrays ..
      INTEGER            IWORK( LDST )
      DOUBLE PRECISION   AI(2), AR(2), BE(2), DWORK(32), IR(LDST,LDST),
     $                   IRCOP(LDST,LDST), LI(LDST,LDST),
     $                   LICOP(LDST,LDST), S(LDST,LDST),
     $                   SCPY(LDST,LDST), T(LDST,LDST), TAUL(LDST),
     $                   TAUR(LDST), TCPY(LDST,LDST)
C     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH
C     .. External Subroutines ..
      EXTERNAL           DGEMM, DGEQR2, DGERQ2, DLACPY, DLARTG, DLASET,
     $                   DLASSQ, DORG2R, DORGR2, DORM2R, DORMR2, DROT,
     $                   DSCAL, MB03YT, SB04OW
C     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, SQRT
C
C     .. Executable Statements ..
C
      INFO = 0
C
C     Quick return if possible.
C     For efficiency, the arguments are not checked.
C
      IF ( N1.LE.0 .OR. N2.LE.0 )
     $   RETURN
      M = N1 + N2
C
      WEAK = .FALSE.
      DTRONG = .FALSE.
C
C     Make a local copy of selected block.
C
      CALL DLASET( 'All', LDST, LDST, ZERO, ZERO, LI, LDST )
      CALL DLASET( 'All', LDST, LDST, ZERO, ZERO, IR, LDST )
      CALL DLACPY( 'Full', M, M, A, LDA, S, LDST )
      CALL DLACPY( 'Full', M, M, B, LDB, T, LDST )
C
C     Compute threshold for testing acceptance of swapping.
C
      EPS = DLAMCH( 'P' )
      SMLNUM = DLAMCH( 'S' ) / EPS
      DSCALE = ZERO
      DSUM = ONE
      CALL DLACPY( 'Full', M, M, S, LDST, DWORK, M )
      CALL DLASSQ( M*M, DWORK, 1, DSCALE, DSUM )
      CALL DLACPY( 'Full', M, M, T, LDST, DWORK, M )
      CALL DLASSQ( M*M, DWORK, 1, DSCALE, DSUM )
      DNORM = DSCALE*SQRT( DSUM )
      THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
C
      IF ( M.EQ.2 ) THEN
C
C        CASE 1: Swap 1-by-1 and 1-by-1 blocks.
C
C        Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
C        using Givens rotations and perform the swap tentatively.
C
         F =  S(2,2)*T(2,2) - T(1,1)*S(1,1)
         G = -S(2,2)*T(1,2) - T(1,1)*S(1,2)
         SB = ABS( T(1,1) )
         SA = ABS( S(2,2) )
         CALL DLARTG( F, G, IR(1,2), IR(1,1), DDUM )
         IR(2,1) = -IR(1,2)
         IR(2,2) =  IR(1,1)
         CALL DROT( 2, S(1,1), 1, S(1,2), 1, IR(1,1), IR(2,1) )
         CALL DROT( 2, T(1,1), LDST, T(2,1), LDST, IR(1,1), IR(2,1) )
         IF( SA.GE.SB ) THEN
            CALL DLARTG( S(1,1), S(2,1), LI(1,1), LI(2,1), DDUM )
         ELSE
            CALL DLARTG( T(2,2), T(2,1), LI(1,1), LI(2,1), DDUM )
            LI(2,1) = -LI(2,1)
         END IF
         CALL DROT( 2, S(1,1), LDST, S(2,1), LDST, LI(1,1), LI(2,1) )
         CALL DROT( 2, T(1,1), 1, T(1,2), 1, LI(1,1), LI(2,1) )
         LI(2,2) =  LI(1,1)
         LI(1,2) = -LI(2,1)
C
C        Weak stability test:
C           |S21| + |T21| <= O(EPS * F-norm((S, T))).
C
         WS = ABS( S(2,1) ) + ABS( T(2,1) )
         WEAK = WS.LE.THRESH
         IF ( .NOT.WEAK )
     $      GO TO 50
C
         IF ( WANDS ) THEN
C
C           Strong stability test:
C             F-norm((A-QL'*S*QR, B-QR'*T*QL)) <= O(EPS*F-norm((A,B))).
C
            CALL DLACPY( 'Full', M, M, A, LDA, DWORK(M*M+1), M )
            CALL DGEMM( 'No Transpose', 'No Transpose', M, M, M, ONE,
     $                  LI, LDST, S, LDST, ZERO, DWORK, M )
            CALL DGEMM( 'No Transpose', 'Transpose', M, M, M, -ONE,
     $                  DWORK, M, IR, LDST, ONE, DWORK(M*M+1), M )
            DSCALE = ZERO
            DSUM = ONE
            CALL DLASSQ( M*M, DWORK(M*M+1), 1, DSCALE, DSUM )
C
            CALL DLACPY( 'Full', M, M, B, LDB, DWORK(M*M+1), M )
            CALL DGEMM( 'No Transpose', 'No Transpose', M, M, M, ONE,
     $                  IR, LDST, T, LDST, ZERO, DWORK, M )
            CALL DGEMM( 'No Transpose', 'Transpose', M, M, M, -ONE,
     $                  DWORK, M, LI, LDST, ONE, DWORK(M*M+1), M )
            CALL DLASSQ( M*M, DWORK(M*M+1), 1, DSCALE, DSUM )
            SS = DSCALE*SQRT( DSUM )
            DTRONG = SS.LE.THRESH
            IF( .NOT.DTRONG )
     $         GO TO 50
         END IF
C
C        Update A and B.
C
         CALL DLACPY( 'All', M, M, S, LDST, A, LDA )
         CALL DLACPY( 'All', M, M, T, LDST, B, LDB )
C
C        Set  N1-by-N2 (2,1) - blocks to ZERO.
C
         A(2,1) = ZERO
         B(2,1) = ZERO
C
C        Accumulate transformations into Q and Z if requested.
C
         IF ( WANTQ )
     $      CALL DROT( 2, Q(1,1), 1, Q(1,2), 1, LI(1,1), LI(2,1) )
         IF ( WANTZ )
     $      CALL DROT( 2, Z(1,1), 1, Z(1,2), 1, IR(1,1), IR(2,1) )
C
C        Exit with INFO = 0 if swap was successfully performed.
C
         RETURN
C
      ELSE
C
C        CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
C                and 2-by-2 blocks.
C
C        Solve the periodic Sylvester equation
C                 S11 * R - L * S22 = SCALE * S12
C                 T11 * L - R * T22 = SCALE * T12
C        for R and L. Solutions in IR and LI.
C
         CALL DLACPY( 'Full', N1, N2, T(1,N1+1), LDST, LI, LDST )
         CALL DLACPY( 'Full', N1, N2, S(1,N1+1), LDST, IR(N2+1,N1+1),
     $                LDST )
         CALL SB04OW( N1, N2, S, LDST, S(N1+1,N1+1), LDST,
     $                IR(N2+1,N1+1), LDST, T, LDST, T(N1+1,N1+1), LDST,
     $                LI, LDST, SCALE, IWORK, LINFO )
         IF ( LINFO.NE.0 )
     $      GO TO 50
C
C        Compute orthogonal matrix QL:
C
C                    QL' * LI = [ TL ]
C                               [ 0  ]
C        where
C                    LI =  [      -L              ].
C                          [ SCALE * identity(N2) ]
C
         DO 10 I = 1, N2
            CALL DSCAL( N1, -ONE, LI(1,I), 1 )
            LI(N1+I,I) = SCALE
   10    CONTINUE
         CALL DGEQR2( M, N2, LI, LDST, TAUL, DWORK, LINFO )
         CALL DORG2R( M, M, N2, LI, LDST, TAUL, DWORK, LINFO )
C
C        Compute orthogonal matrix RQ:
C
C                    IR * RQ' =   [ 0  TR],
C
C         where IR = [ SCALE * identity(N1), R ].
C
         DO 20 I = 1, N1
            IR(N2+I,I) = SCALE
   20    CONTINUE
         CALL DGERQ2( N1, M, IR(N2+1,1), LDST, TAUR, DWORK, LINFO )
         CALL DORGR2( M, M, N1, IR, LDST, TAUR, DWORK, LINFO )
C
C        Perform the swapping tentatively:
C
         CALL DGEMM( 'Transpose', 'No Transpose', M, M, M, ONE, LI,
     $               LDST, S, LDST, ZERO, DWORK, M )
         CALL DGEMM( 'No Transpose', 'Transpose', M, M, M, ONE, DWORK,
     $               M, IR, LDST, ZERO, S, LDST )
         CALL DGEMM( 'No Transpose', 'No Transpose', M, M, M, ONE, IR,
     $               LDST, T, LDST, ZERO, DWORK, M )
         CALL DGEMM( 'No Transpose', 'No Transpose', M, M, M, ONE,
     $               DWORK, M, LI, LDST, ZERO, T, LDST )
         CALL DLACPY( 'All', M, M, S, LDST, SCPY, LDST )
         CALL DLACPY( 'All', M, M, T, LDST, TCPY, LDST )
         CALL DLACPY( 'All', M, M, IR, LDST, IRCOP, LDST )
         CALL DLACPY( 'All', M, M, LI, LDST, LICOP, LDST )
C
C        Triangularize the B-part by a QR factorization.
C        Apply transformation (from left) to A-part, giving S.
C
         CALL DGEQR2( M, M, T, LDST, TAUR, DWORK, LINFO )
         CALL DORM2R( 'Right', 'No Transpose', M, M, M, T, LDST, TAUR,
     $                S, LDST, DWORK, LINFO )
         CALL DORM2R( 'Left', 'Transpose', M, M, M, T, LDST, TAUR,
     $                IR, LDST, DWORK, LINFO )
C
C        Compute F-norm(S21) in BRQA21. (T21 is 0.)
C
         DSCALE = ZERO
         DSUM = ONE
         DO 30 I = 1, N2
            CALL DLASSQ( N1, S(N2+1,I), 1, DSCALE, DSUM )
   30    CONTINUE
         BRQA21 = DSCALE*SQRT( DSUM )
C
C        Triangularize the B-part by an RQ factorization.
C        Apply transformation (from right) to A-part, giving S.
C
         CALL DGERQ2( M, M, TCPY, LDST, TAUL, DWORK, LINFO )
         CALL DORMR2( 'Left', 'No Transpose', M, M, M, TCPY, LDST,
     $                TAUL, SCPY, LDST, DWORK, LINFO )
         CALL DORMR2( 'Right', 'Transpose', M, M, M, TCPY, LDST,
     $                TAUL, LICOP, LDST, DWORK, LINFO )
C
C        Compute F-norm(S21) in BQRA21. (T21 is 0.)
C
         DSCALE = ZERO
         DSUM = ONE
         DO 40 I = 1, N2
            CALL DLASSQ( N1, SCPY(N2+1,I), 1, DSCALE, DSUM )
   40    CONTINUE
         BQRA21 = DSCALE*SQRT( DSUM )
C
C        Decide which method to use.
C          Weak stability test:
C             F-norm(S21) <= O(EPS * F-norm((S, T)))
C
         IF ( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESH ) THEN
            CALL DLACPY( 'All', M, M, SCPY, LDST, S, LDST )
            CALL DLACPY( 'All', M, M, TCPY, LDST, T, LDST )
            CALL DLACPY( 'All', M, M, IRCOP, LDST, IR, LDST )
            CALL DLACPY( 'All', M, M, LICOP, LDST, LI, LDST )
         ELSE IF ( BRQA21.GE.THRESH ) THEN
            GO TO 50
         END IF
C
C        Set lower triangle of B-part to zero
C
         CALL DLASET( 'Lower', M-1, M-1, ZERO, ZERO, T(2,1), LDST )
C
         IF ( WANDS ) THEN
C
C           Strong stability test:
C              F-norm((A-QL*S*QR', B-QR*T*QL')) <= O(EPS*F-norm((A,B)))
C
            CALL DLACPY( 'All', M, M, A, LDA, DWORK(M*M+1), M )
            CALL DGEMM( 'No Transpose', 'No Transpose', M, M, M, ONE,
     $                  LI, LDST, S, LDST, ZERO, DWORK, M )
            CALL DGEMM( 'No Transpose', 'No Transpose', M, M, M, -ONE,
     $                  DWORK, M, IR, LDST, ONE, DWORK(M*M+1), M )
            DSCALE = ZERO
            DSUM = ONE
            CALL DLASSQ( M*M, DWORK(M*M+1), 1, DSCALE, DSUM )
C
            CALL DLACPY( 'All', M, M, B, LDB, DWORK(M*M+1), M )
            CALL DGEMM( 'Transpose', 'No Transpose', M, M, M, ONE,
     $                  IR, LDST, T, LDST, ZERO, DWORK, M )
            CALL DGEMM( 'No Transpose', 'Transpose', M, M, M, -ONE,
     $                  DWORK, M, LI, LDST, ONE, DWORK(M*M+1), M )
            CALL DLASSQ( M*M, DWORK(M*M+1), 1, DSCALE, DSUM )
            SS = DSCALE*SQRT( DSUM )
            DTRONG = ( SS.LE.THRESH )
            IF( .NOT.DTRONG )
     $         GO TO 50
C
         END IF
C
C        If the swap is accepted ("weakly" and "strongly"), apply the
C        transformations and set N1-by-N2 (2,1)-block to zero.
C
         CALL DLASET( 'All', N1, N2, ZERO, ZERO, S(N2+1,1), LDST )
C
C        Copy (S,T) to (A,B).
C
         CALL DLACPY( 'All', M, M, S, LDST, A, LDA )
         CALL DLACPY( 'All', M, M, T, LDST, B, LDB )
         CALL DLASET( 'All', LDST, LDST, ZERO, ZERO, T, LDST )
C
C        Standardize existing 2-by-2 blocks.
C
         CALL DLASET( 'All', M, M, ZERO, ZERO, DWORK, M )
         DWORK(1) = ONE
         T(1,1) = ONE
         IF ( N2.GT.1 ) THEN
            CALL MB03YT( A, LDA, B, LDB, AR, AI, BE, DWORK(1), DWORK(2),
     $                   T(1,1), T(2,1) )
            DWORK(M+1) = -DWORK(2)
            DWORK(M+2) = DWORK(1)
            T(N2,N2) = T(1,1)
            T(1,2) = -T(2,1)
         END IF
         DWORK(M*M) = ONE
         T(M,M) = ONE
C
         IF ( N1.GT.1 ) THEN
            CALL MB03YT( A(N2+1,N2+1), LDA, B(N2+1,N2+1), LDB, TAUR,
     $                   TAUL, DWORK(M*M+1), DWORK(N2*M+N2+1),
     $                   DWORK(N2*M+N2+2), T(N2+1,N2+1), T(M,M-1) )
            DWORK(M*M) = DWORK(N2*M+N2+1)
            DWORK(M*M-1 ) = -DWORK(N2*M+N2+2)
            T(M,M) = T(N2+1,N2+1)
            T(M-1,M) = -T(M,M-1)
         END IF
C
         CALL DGEMM( 'Transpose', 'No Transpose', N2, N1, N2, ONE,
     $               DWORK, M, A(1,N2+1), LDA, ZERO, DWORK(M*M+1), N2 )
         CALL DLACPY( 'All', N2, N1, DWORK(M*M+1), N2, A(1,N2+1), LDA )
         CALL DGEMM( 'Transpose', 'No Transpose', N2, N1, N2, ONE,
     $               T(1,1), LDST, B(1,N2+1), LDB, ZERO,
     $               DWORK(M*M+1), N2 )
         CALL DLACPY( 'All', N2, N1, DWORK(M*M+1), N2, B(1,N2+1), LDB )
         CALL DGEMM( 'No Transpose', 'No Transpose', M, M, M, ONE, LI,
     $               LDST, DWORK, M, ZERO, DWORK(M*M+1), M )
         CALL DLACPY( 'All', M, M, DWORK(M*M+1), M, LI, LDST )
         CALL DGEMM( 'No Transpose', 'No Transpose', N2, N1, N1, ONE,
     $               A(1,N2+1), LDA, T(N2+1,N2+1), LDST, ZERO,
     $               DWORK(M*M+1), M )
         CALL DLACPY( 'All', N2, N1, DWORK(M*M+1), M, A(1,N2+1), LDA )
         CALL DGEMM( 'No Transpose', 'No Transpose', N2, N1, N1, ONE,
     $               B(1,N2+1), LDB, DWORK(N2*M+N2+1), M, ZERO,
     $               DWORK(M*M+1), M )
         CALL DLACPY( 'All', N2, N1, DWORK(M*M+1), M, B(1,N2+1), LDB )
         CALL DGEMM( 'Transpose', 'No Transpose', M, M, M, ONE, T,
     $               LDST, IR, LDST, ZERO, DWORK, M )
         CALL DLACPY( 'All', M, M, DWORK, M, IR, LDST )
C
C        Accumulate transformations into Q and Z if requested.
C
         IF( WANTQ ) THEN
            CALL DGEMM( 'No Transpose', 'No Transpose', M, M, M, ONE, Q,
     $                  LDQ, LI, LDST, ZERO, DWORK, M )
            CALL DLACPY( 'All', M, M, DWORK, M, Q, LDQ )
         END IF
C
         IF( WANTZ ) THEN
            CALL DGEMM( 'No Transpose', 'Transpose', M, M, M, ONE, Z,
     $                  LDZ, IR, LDST, ZERO, DWORK, M )
            CALL DLACPY( 'Full', M, M, DWORK, M, Z, LDZ )
C
         END IF
C
C        Exit with INFO = 0 if swap was successfully performed.
C
         RETURN
C
      END IF
C
C     Exit with INFO = 1 if swap was rejected.
C
   50 CONTINUE
C
      INFO = 1
      RETURN
C *** Last line of MB03WA ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
module pyHazelMod
use iso_c_binding, only: c_int, c_double
use vars
use maths
use io
use SEE
use rt_coef
use synth
use allen
implicit none

contains
subroutine c_hazel(synModeInput, nSlabsInput, B1Input, B2Input, hInput, tau1Input, tau2Input, boundaryInput, &
    transInput, atomicPolInput, magopttInput, anglesInput, nLambdaInput, lambdaAxisInput, dopplerWidthInput, dopplerWidth2Input, dampingInput, &
    dopplerVelocityInput, dopplerVelocity2Input, ffInput, betaInput, beta2Input, nbarInput, omegaInput, normalizationInput, deltaCollInput, &
    wavelengthOutput, stokesOutput, epsOutput, etaOutput) bind(c)

    integer(c_int), intent(in) :: synModeInput, nSlabsInput, transInput, atomicPolInput, magopttInput
    integer(c_int), intent(in) :: nLambdaInput, normalizationInput
    real(c_double), intent(in), dimension(nLambdaInput) :: lambdaAxisInput
    real(c_double), intent(in), dimension(3) :: B1Input, anglesInput
    real(c_double), intent(in), dimension(3) :: B2Input
    real(c_double), intent(in), dimension(4,nLambdaInput) :: boundaryInput
    real(c_double), intent(in):: tau2Input, dopplerWidth2Input, dopplerVelocity2Input, ffInput
    real(c_double), intent(in), dimension(4) :: nbarInput, omegaInput
    real(c_double), intent(in) :: hInput, tau1Input, dopplerWidthInput, dampingInput, dopplerVelocityInput, betaInput, beta2Input, deltaCollInput
    real(c_double), intent(out), dimension(nLambdaInput) :: wavelengthOutput
    real(c_double), intent(out), dimension(4,nLambdaInput) :: stokesOutput
    real(c_double), intent(out), dimension(4,nLambdaInput) :: epsOutput
    real(c_double), intent(out), dimension(4,4,nLambdaInput) :: etaOutput

    integer :: n, nterml, ntermu, error

    real(c_double) :: ae, wavelength, reduction_factor, reduction_factor_omega, j10
    integer :: i, j

! Initialize the random number generator
!   call random_seed

! Initialize Allen's data
!   call read_allen_data

! Fill the factorial array
!   call factrl

! Read the configuration file
!   call read_config
    input_model_file = 'helium.mod'
    input_experiment = 'init_parameters.dat'
    verbose_mode = 0
    linear_solver = 0
    synthesis_mode = synModeInput
    working_mode = 0

    if (normalizationInput == 0) then
        observation%normalization = 'maximum'
    else
        observation%normalization = 'peak'
    endif

! Read the atomic model
!   call read_model_file(input_model_file)

! Set the variables for the experiment from the parameters of the subroutine
    isti = 1
    imag = 1
    idep = 0
    use_paschen_back = 1
    params%nslabs = nSlabsInput
    if (params%nslabs == 3 .or. params%nslabs == -2) then
        params%bgauss = B1Input(1)
        params%thetabd = B1Input(2)
        params%chibd = B1Input(3)
        params%bgauss2 = B2Input(1)
        params%thetabd2 = B2Input(2)
        params%chibd2 = B2Input(3)
    else
        params%bgauss = B1Input(1)
        params%thetabd = B1Input(2)
        params%chibd = B1Input(3)
    endif
    params%height = hInput
    if (params%nslabs == 1) then
        params%dtau = tau1Input
    endif

    if (params%nslabs == 2 .or. params%nslabs == 3) then
        params%dtau = tau1Input
        params%dtau2 = tau2Input
    endif

! Two components with filling factor
    if (params%nslabs == -2) then
        params%dtau = tau1Input
        params%dtau2 = tau2Input
        params%ff = ffInput
    endif

    params%beta = betaInput
    if (params%nslabs /= 1) then
        params%beta2 = beta2Input
    endif

    fixed%nemiss = transInput
    fixed%use_atomic_pol = atomicPolInput
    fixed%thetad = anglesInput(1)
    fixed%chid = anglesInput(2)
    fixed%gammad = anglesInput(3)

    if (params%nslabs == 1 .or. params%nslabs == 2) then
        params%vdopp = dopplerWidthInput
        params%damping = dampingInput
        fixed%damping_treatment = 0
    endif
    if (params%nslabs == 3 .or. params%nslabs == -2) then
        params%vdopp = dopplerWidthInput
        params%vdopp2 = dopplerWidth2Input
        params%damping = dampingInput
        fixed%damping_treatment = 0
    endif

! Read the wavelength of the transition that we want to synthesize
    fixed%wl = atom%wavelength(transInput)

! Set the values of nbar and omega in case they are given
    nbarExternal = nbarInput
    omegaExternal = omegaInput

! If nbar=0 or omega=0, use the numbers from Allen. If not, treat them as reduction factors
    do i = 1, atom%ntran
        if (nbarExternal(i) == 0) then
            nbarExternal(i) = 1.0
        endif
        if (omegaExternal(i) == 0) then
            omegaExternal(i) = 1.0
        endif
    enddo

    if (params%nslabs == 1) then
        params%vmacro = dopplerVelocityInput
    else
        params%vmacro = dopplerVelocityInput
        params%vmacro2 = dopplerVelocity2Input
    endif

    use_mag_opt_RT = magopttInput
    use_stim_emission_RT = 1

    if (deltaCollInput /= 0) then
        idep = 1
        params%delta_collision = deltaCollInput
    endif

!*********************************
!** SYNTHESIS MODE
!*********************************
    fixed%no = nLambdaInput
    observation%n = fixed%no

    if (associated(observation%wl)) deallocate(observation%wl)
    if (associated(inversion%stokes_unperturbed)) deallocate(inversion%stokes_unperturbed)
    if (associated(fixed%stokes_boundary)) deallocate(fixed%stokes_boundary)

    allocate(observation%wl(observation%n))
    allocate(inversion%stokes_unperturbed(0:3,fixed%no))
    allocate(fixed%stokes_boundary(0:3,observation%n))

    fixed%stokes_boundary(0:3,:) = boundaryInput

    observation%wl = lambdaAxisInput

    fixed%omax = minval(lambdaAxisInput)
    fixed%omin = maxval(lambdaAxisInput)

    call do_synthesis(params, fixed, observation, inversion%stokes_unperturbed, error)

    do i = 1, 4
        stokesOutput(i,:) = inversion%stokes_unperturbed(i-1,:)
    enddo

    wavelengthOutput = observation%wl + fixed%wl

! Fill the emission vector and absorption matrix
    do i = 1, 4
        epsOutput(i,:) = epsilon(i-1,:)
    enddo

    etaOutput(1,1,:) = eta(0,:) - eta_stim(0,:)
    etaOutput(2,2,:) = eta(0,:) - eta_stim(0,:)
    etaOutput(3,3,:) = eta(0,:) - eta_stim(0,:)
    etaOutput(4,4,:) = eta(0,:) - eta_stim(0,:)
    etaOutput(1,2,:) = eta(1,:) - eta_stim(1,:)
    etaOutput(2,1,:) = eta(1,:) - eta_stim(1,:)
    etaOutput(1,3,:) = eta(2,:) - eta_stim(2,:)
    etaOutput(3,1,:) = eta(2,:) - eta_stim(2,:)
    etaOutput(1,4,:) = eta(3,:) - eta_stim(3,:)
    etaOutput(4,1,:) = eta(3,:) - eta_stim(3,:)
    etaOutput(2,3,:) = mag_opt(3,:) - mag_opt_stim(3,:)
    etaOutput(3,2,:) = -mag_opt(3,:) - mag_opt_stim(3,:)
    etaOutput(2,4,:) = -mag_opt(2,:) - mag_opt_stim(2,:)
    etaOutput(4,2,:) = mag_opt(2,:) - mag_opt_stim(2,:)
    etaOutput(3,4,:) = mag_opt(1,:) - mag_opt_stim(1,:)
    etaOutput(4,3,:) = -mag_opt(1,:) - mag_opt_stim(1,:)

    if (allocated(epsilon)) deallocate(epsilon)
	if (allocated(epsilon_zeeman)) deallocate(epsilon_zeeman)
    if (allocated(eta)) deallocate(eta)
    if (allocated(eta_zeeman)) deallocate(eta_zeeman)
    if (allocated(eta_stim)) deallocate(eta_stim)
    if (allocated(eta_stim_zeeman)) deallocate(eta_stim_zeeman)
    if (allocated(mag_opt)) deallocate(mag_opt)
    if (allocated(mag_opt_zeeman)) deallocate(mag_opt_zeeman)
    if (allocated(mag_opt_stim)) deallocate(mag_opt_stim)
    if (allocated(mag_opt_stim_zeeman)) deallocate(mag_opt_stim_zeeman)

!   open(unit=31,file=input_model_file,action='write',status='replace')
!   close(31,status='delete')

end subroutine c_hazel

subroutine c_init() bind(c)

! Initialize the random number generator
    call random_seed

! Deallocate memory. Important if we change the number of wavelength points from Python
    if (fixed%no /= size(observation%wl)) then
        if (associated(observation%wl)) deallocate(observation%wl)
        if (associated(inversion%stokes_unperturbed)) deallocate(inversion%stokes_unperturbed)
        if (associated(fixed%stokes_boundary)) deallocate(fixed%stokes_boundary)
    endif

! Initialize Allen's data
    call read_allen_data

! Fill the factorial array
    call factrl

    input_model_file = 'helium.mod'

! Read the atomic model
    call read_model_file(input_model_file)

end subroutine c_init

end module pyHazelMod

max_stars_repo_name: aasensio/hazel

max_issues_repo_name: aasensio/hazel

max_forks_repo_name: aasensio/hazel


text:
      subroutine intbs6(jflag,jflgi,kmax,narn1a,narn2a,nbaspd,nbtd,
     $ nbti,nbtmax,ndrsrd,ndecsp,noutpt,nsta,nstmax,nttyo,uspeca,
     $ ubmtbi)
c
c     This subroutine interprets the data file basis species read from
c     the input file. "Data file" species to be created according to
c     directives read from the input file are ignored. This subroutine
c     also sets up the jflag array.
c
c     This subroutine is called by:
c
c       EQ6/eq6.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c
c     Principal output:
c
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer kmax,nbtmax,nstmax
c
      integer jflag(nstmax),jflgi(nbtmax),nbaspd(nbtmax),
     $ ndecsp(nbtmax),ndrsrd(2,nstmax)
c
      integer narn1a,narn2a,nbtd,nbti,noutpt,nsta,nttyo
c
      character*56 uspn56
      character*48 uspeca(nstmax),ubmtbi(nbtmax)
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer jlen,j2,nb,nbi,nerr,ns,nt
c
      integer ilnobl,nbasis
c
      character*8 ux8
c
c-----------------------------------------------------------------------
c
      nerr = 0
c
      do ns = 1,nsta
        jflag(ns) = 0
      enddo
c
      do ns = narn1a,narn2a
        jflag(ns) = 30
      enddo
c
      do nb = 1,nbtd
        ns = nbaspd(nb)
        nt = ndrsrd(2,ns) - ndrsrd(1,ns) + 1
        if (nt .lt. 2) jflag(ns) = -1
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      do nbi = 1,nbti
        do ns = 1,nsta
          if (uspeca(ns)(1:48) .eq. ubmtbi(nbi)(1:48)) then
            if (jflgi(nbi).ne.0 .and. jflgi(nbi).ne.30) then
c
c             Calling sequence substitutions:
c               ubmtbi(nbi) for unam48
c
              call fmspnx(jlen,ubmtbi(nbi),uspn56)
              write (ux8,'(i5)') jflgi(nbi)
              call lejust(ux8)
              j2 = ilnobl(ux8)
              write (noutpt,1000) uspn56(1:jlen),ux8(1:j2)
              write (nttyo,1000) uspn56(1:jlen),ux8(1:j2)
 1000         format(/' * Error - (EQ6/intbs6) The species ',a,
     $        /7x,'has an has illegal jflag value of ',a,' on the',
     $        ' input file',/7x,'(the only legal values are 0 and 30).')
              nerr = nerr + 1
            endif
            jflag(ns) = jflgi(nbi)
c
c           Calling sequence substitutions:
c             nbaspd for nbasp
c             nbtd for nbt
c
            nb = nbasis(nbaspd,nbtd,nbtmax,ns)
            if (nb .eq. 0) then
c
c             Add a species to the working basis set.
c
              nbtd = nbtd + 1
              if (nbtd .gt. nbtmax) then
c
c               Calling sequence substitutions:
c                 ubmtbi(nbi) for unam48
c
                call fmspnx(jlen,ubmtbi(nbi),uspn56)
                write (ux8,'(i5)') nbtmax
                call lejust(ux8)
                j2 = ilnobl(ux8)
                write (noutpt,1005) ux8(1:j2),uspn56(1:jlen)
                write (nttyo,1005) ux8(1:j2),uspn56(1:jlen)
 1005           format(/' * Error - (EQ6/intbs6) The maximum ',a,
     $          /7x,'basis species have been exceeded in',
     $          ' interpreting the',/7x,'input file while processing',
     $          ' data file basis species',/7x,a,'. Increase the',
     $          ' dimensioning parameter nbtpar.')
                nerr = nerr + 1
                go to 115
              endif
              nb = nbtd
              nbaspd(nb) = ns
            endif
            ndecsp(nbi) = nb
            go to 115
          endif
        enddo
c
c       Calling sequence substitutions:
c         ubmtbi(nbi) for unam48
c
        call fmspnx(jlen,ubmtbi(nbi),uspn56)
        write (noutpt,1010) uspn56(1:jlen)
        write (nttyo,1010) uspn56(1:jlen)
 1010   format(/' * Note - (EQ6/intbs6) The species "',a,'"',
     $  /7x,'appears as a basis species on the input file, but it',
     $  " wasn't",/7x,'read from the data file. If it is a species',
     $  ' to be created by the code,',/7x,'such as a generic ion',
     $  ' exchanger species, there is no problem.')
c
  115   continue
      enddo
c
      if (nerr .gt. 0) stop
c
      end