text:
      SUBROUTINE RNPERM (SIZE, PERM)

C  RNPERM returns a random permutation of integers 1-SIZE in PERM.
C  For example, if SIZE is specified as 4, PERM might = [2,3,1,4].

      INTEGER   SIZE, PERM(SIZE), IX, RANPOS, SWAP, RNDINT

      DO 100 IX = 1,SIZE
         PERM(IX) = IX
100   CONTINUE

      DO 200 IX = SIZE,2, -1
         RANPOS = RNDINT (1,IX)
         SWAP = PERM(IX)
         PERM(IX) = PERM(RANPOS)
         PERM(RANPOS) = SWAP
200   CONTINUE

      END

ext: f


text:
MODULE iarray
!
!  Contains all array var names.
!

    use constants, only : wp

    implicit none

    real(kind=wp), allocatable :: jmean(:,:,:), jmeanGLOBAL(:,:,:)
end MODULE iarray

max_stars_repo_name: lewisfish/signedMCRT

max_issues_repo_name: lewisfish/signedMCRT

max_forks_repo_name: lewisfish/signedMCRT


text:
      subroutine prreac(cdrs,ndrs,ndrsmx,ndrsr,nf,ns,nstmax,uspec)
c
c     This subroutine writes the n-th reaction in a set on the file
c     whose unit number isf nf.
c
c     This subroutine is called by:
c
c       EQLIB/alters.f
c       EQLIB/echolk.f
c       EQLIB/switch.f
c       EQLIB/swtchb.f
c       EQLIB/swtchk.f
c       EQ3NR/arrsim.f
c       EQ3NR/dawfix.f
c       EQ3NR/echox.f
c       EQ3NR/eq3nr.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       cdrs   = the array of reaction coefficients
c       ndrs   = the array of corresponding species indices
c       ndrsmx = dimension of the cdrs and ndrs arrays
c       ndrsr  = aray giving the range in the cdrs and ndrs arrays
c                  corresponding to the reaction for a given species
c       nf     = the unit number of the file to write on
c       ns     = the index of the species whose reaction is to
c                  be written
c       uspec  = array of species names
c
c     Principal output:
c
c       None
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ndrsmx,nstmax
c
      integer ndrs(ndrsmx),ndrsr(2,nstmax)
      integer nf,ns
c
      character*48 uspec(nstmax)
c
      real*8 cdrs(ndrsmx)
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer jlen,n,nr1,nr2,nss,nt
c
      logical qdtach,qfirst
c
      character*56 uspn56
c
      real*8 cx
c
c-----------------------------------------------------------------------
c
      write (nf,1000)
 1000 format(1x)
c
      nr1 = ndrsr(1,ns)
      nr2 = ndrsr(2,ns)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      nt = nr2 - nr1 + 1
      if (nt .lt. 2) then
c
c       Calling sequence substitutions:
c         uspec(ns) for unam48
c
        call fmspnx(jlen,uspec(ns),uspn56)
c
        write (nf,1010) uspn56(1:jlen)
 1010   format(3x,a,' is a strict basis species and has no reaction.',/)
        go to 999
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      qdtach = .false.
      do n = nr1 + 1,nr2
        nss = ndrs(n)
        if (nss .eq. 0) then
          qdtach = .true.
          go to 100
        endif
      enddo
  100 continue
c
      if (qdtach) then
c
c       Calling sequence substitutions:
c         uspec(ns) for unam48
c
        call fmspnx(jlen,uspec(ns),uspn56)
c
        write (nf,1020) uspn56(1:jlen)
 1020   format(3x,a,' is a detached auxiliary basis species and',
     $  ' in effect',/3x,'has no reaction.',/)
        go to 999
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Write the coefficients and names of the reactants.
c
      qfirst = .true.
      do n = nr1,nr2
        cx = -cdrs(n)
        nss = ndrs(n)
        if (cx .gt. 0.) then
c
c         Calling sequence substitutions:
c           uspec(nss) for unam48
c
          call fmspnx(jlen,uspec(nss),uspn56)
c
          if (qfirst) then
            write (nf,1030) cx,uspn56(1:jlen)
 1030       format(4x,f7.3,2x,a)
            qfirst = .false.
          else
            write (nf,1050) cx,uspn56(1:jlen)
 1050       format(2x,'+ ',f7.3,2x,a)
          endif
        endif
      enddo
c
      write (nf,1070)
 1070 format(10x,'==')
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Write the coefficients and names of the products.
c
      qfirst = .true.
      do n = nr1,nr2
        cx = cdrs(n)
        nss = ndrs(n)
        if (cx .gt. 0.) then
c
c         Calling sequence substitutions:
c           uspec(nss) for unam48
c
          call fmspnx(jlen,uspec(nss),uspn56)
c
          if (qfirst) then
            write (nf,1030) cx,uspn56(1:jlen)
            qfirst = .false.
          else
            write (nf,1050) cx,uspn56(1:jlen)
          endif
        endif
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
  999 continue
      end


SUBROUTINE FISHHABITAT(IOPENFISH)
USE GLOBAL;USE MAIN;USE SCREENC; USE KINETIC, ONLY:O2, CHLA, NO3,NH4,PO4, TP, GAMMA, SED, SATO; USE TVDC, ONLY: CONSTITUENTS; USE NAMESC, ONLY: CNAME2,CDNAME2; USE LOGICC
IMPLICIT NONE
INTEGER :: IFISH, N, NSEG,IOPENFISH,KKMAX,KSEG,JJW,JWFILE,JBFILE    !,JWFILE1,JBFILE1
REAL    :: VOLTOT,O2CORR,DOSAT
CHARACTER*80, ALLOCATABLE, DIMENSION(:) :: FISHNAME
CHARACTER*80 :: CONHAB,CONAVG,CONSURF,CONSOD
REAL, ALLOCATABLE, DIMENSION(:)   :: FISHTEMPL,FISHTEMPH,FISHDO,HABVOL,PHABVOL,CDO,CPO4,CNO3,CNH4,CCHLA,CTOTP,CDOS,CPO4S,CNO3S,CNH4S,CCHLAS,CTOTPS,CGAMMA,SSEDD,VOLTOTBR,VOLTOTWB
REAL, ALLOCATABLE, DIMENSION(:,:) :: HABVOLBR,HABVOLWB,PHABVOLBR,PHABVOLWB
INTEGER, ALLOCATABLE, DIMENSION (:) :: ISEGVOL
SAVE

        !JWFILE1=9549
        !JBFILE1=9749

IF(IOPENFISH .NE. 3)THEN
! read input file

if(nit == 1.or.iopenfish==0)then
open(FISHHABFN,file='w2_habitat.npt',status='old')
! skip 1st 2 lines
read(FISHHABFN,*)
read(FISHHABFN,*)
read(FISHHABFN,*)ifish,conhab
allocate (fishname(ifish),fishtempl(ifish),fishtemph(ifish),fishdo(ifish),habvol(ifish),phabvol(ifish),habvolbr(nbr,ifish),habvolwb(nwb,ifish),phabvolbr(nbr,ifish),phabvolwb(nwb,ifish),voltotbr(nbr),voltotwb(nwb))
read(FISHHABFN,*)
do i=1,ifish
  read(FISHHABFN,*)fishname(i),fishtempl(i),fishtemph(i),fishdo(i)
enddo
read(FISHHABFN,*)
read(FISHHABFN,*)nseg,conavg     ! volume weighted averages of critical WQ parameters
read(FISHHABFN,*)
allocate(isegvol(nseg),cdo(nseg),cpo4(nseg),cno3(nseg),cnh4(nseg),cchla(nseg),ctotp(nseg),cdos(nseg),cpo4s(nseg),cno3s(nseg),cnh4s(nseg),cchlas(nseg),ctotps(nseg),cgamma(nseg))
allocate(ssedd(imx))
read(FISHHABFN,*)(isegvol(i),i=1,nseg)
read(FISHHABFN,*)
read(FISHHABFN,*)kseg,consurf     ! # of layers for surface average
read(FISHHABFN,*)
read(FISHHABFN,*)consod
close(FISHHABFN)

if(restart_in)then
        open(FISHHABFN,file=conhab,POSITION='APPEND')
        JDAY1=0.0
        REWIND (FISHHABFN)
        READ   (FISHHABFN,'(///)')
        DO I=1,IFISH
        READ(FISHHABFN,*)
        ENDDO
        READ(FISHHABFN,'(//)',END=106)

        DO WHILE (JDAY1 < JDAY)
          READ (FISHHABFN,'(F10.0)',END=106) JDAY1
        END DO
        BACKSPACE (FISHHABFN)
        106     JDAY1=0.0
        jbfile=jbfile1;jwfile=jwfile1
        do jw=1,nwb
        jwfile=jwfile+1
         WRITE (SEGNUM,'(I0)') JW
         SEGNUM = ADJUSTL(SEGNUM)
         L      = LEN_TRIM(SEGNUM)
         OPEN(JWFILE,FILE='fish_habitat_wb'//SEGNUM(1:L)//'.opt',POSITION='APPEND')
         JDAY1=0.0
         REWIND (JWFILE)
        READ   (JWFILE,'(///)')
        DO I=1,IFISH
        READ(JWFILE,*)
        ENDDO
        READ(JWFILE,'(//)',END=110)
        DO WHILE (JDAY1 < JDAY)
          READ (JWFILE,'(F10.0)',END=110) JDAY1
        END DO
        BACKSPACE (JWFILE)
        110     JDAY1=0.0

            do jb=bs(jw),be(jw)
            jbfile=jbfile+1
            WRITE (SEGNUM,'(I0)') JB
            SEGNUM = ADJUSTL(SEGNUM)
            L      = LEN_TRIM(SEGNUM)
            OPEN(JBFILE,FILE='fish_habitat_br'//SEGNUM(1:L)//'.opt',POSITION='APPEND')
            JDAY1=0.0
            REWIND (JBFILE)
            READ   (JBFILE,'(///)')
            DO I=1,IFISH
            READ(JBFILE,*)
            ENDDO
            READ(JBFILE,'(//)',END=111)
            DO WHILE (JDAY1 < JDAY)
            READ (JBFILE,'(F10.0)',END=111) JDAY1
            END DO
            BACKSPACE (JBFILE)
            111     JDAY1=0.0

            enddo
        enddo

        if(oxygen_demand)then
            open(FISHHABFN+1,file=conavg,POSITION='APPEND')
            REWIND (FISHHABFN+1)
            READ   (FISHHABFN+1,'(//)')
            DO WHILE (JDAY1 < JDAY)
            READ (FISHHABFN+1,'(F10.0)',END=107) JDAY1
            END DO
            BACKSPACE (FISHHABFN+1)
            107     JDAY1=0.0
            open(FISHHABFN+2,file=consurf,POSITION='APPEND')
            REWIND (FISHHABFN+2)
            READ   (FISHHABFN+2,'(//)')
            DO WHILE (JDAY1 < JDAY)
            READ (FISHHABFN+2,'(F10.0)',END=108) JDAY1
            END DO
            BACKSPACE (FISHHABFN+2)
            108     JDAY1=0.0

                do jjw=1,nwb
                IF (SEDIMENT_CALC(JJW))then
                open(FISHHABFN+3,file=consod,POSITION='APPEND')
                REWIND (FISHHABFN+3)
                READ   (FISHHABFN+3,'(/)')
                DO WHILE (JDAY1 < JDAY)
                READ (FISHHABFN+3,'(F10.0)',END=109) JDAY1
                END DO
                BACKSPACE (FISHHABFN+3)
                109     JDAY1=0.0
                exit
                ENDIF
                enddo
         endif
else

        open(FISHHABFN,file=conhab,status='unknown')
        write(FISHHABFN,*)'Fish habitat analysis: CE-QUAL-W2 model results'
        write(FISHHABFN,*)
        write(FISHHABFN,*)'Species, Temperature minimum, Temperature maximum, Dissolved oxygen minimum'
        do i=1,ifish
        write(FISHHABFN,"(a,',',t25,f8.2,',',f8.2,',',f8.2)")trim(fishname(i)),fishtempl(i),fishtemph(i),fishdo(i)
        enddo
        write(FISHHABFN,*)
        write(FISHHABFN,100)(trim(fishname(i)),trim(fishname(i)),i=1,ifish)
        100 format('JDAY,',<ifish>('%VOL-',A,',','HAB-VOL(m3)-',A,','))

        jbfile=jbfile1;jwfile=jwfile1
        do jw=1,nwb
        jwfile=jwfile+1
         WRITE (SEGNUM,'(I0)') JW
         SEGNUM = ADJUSTL(SEGNUM)
         L      = LEN_TRIM(SEGNUM)
         OPEN(JWFILE,FILE='fish_habitat_wb'//SEGNUM(1:L)//'.opt',STATUS='UNKNOWN')

                write(JWFILE,*)'Fish habitat analysis: CE-QUAL-W2 model results'
                write(JWFILE,*)'FOR WATERBODY:',JW
                write(JWFILE,*)'Species, Temperature minimum, Temperature maximum, Dissolved oxygen minimum'
                do i=1,ifish
                write(JWFILE,"(a,',',t25,f8.2,',',f8.2,',',f8.2)")trim(fishname(i)),fishtempl(i),fishtemph(i),fishdo(i)
                enddo
                write(JWFILE,*)
                write(JWFILE,100)(trim(fishname(i)),trim(fishname(i)),i=1,ifish)

            do jb=bs(jw),be(jw)
            jbfile=jbfile+1
            WRITE (SEGNUM,'(I0)') JB
            SEGNUM = ADJUSTL(SEGNUM)
            L      = LEN_TRIM(SEGNUM)
            OPEN(JBFILE,FILE='fish_habitat_br'//SEGNUM(1:L)//'.opt',STATUS='UNKNOWN')
                write(JBFILE,*)'Fish habitat analysis: CE-QUAL-W2 model results'
                write(JBFILE,*)'FOR BRANCH:',JB
                write(JBFILE,*)'Species, Temperature minimum, Temperature maximum, Dissolved oxygen minimum'
                do i=1,ifish
                write(JBFILE,"(a,',',t25,f8.2,',',f8.2,',',f8.2)")trim(fishname(i)),fishtempl(i),fishtemph(i),fishdo(i)
                enddo
                write(JBFILE,*)
                write(JBFILE,100)(trim(fishname(i)),trim(fishname(i)),i=1,ifish)
            enddo
        enddo

  if(oxygen_demand)then
    open(FISHHABFN+1,file=conavg,status='unknown')
    write(FISHHABFN+1,'(a,80(1x,i4))')'Volume weighted WQ parameters at segments:',(isegvol(i),i=1,nseg)
    write(FISHHABFN+1,101)(trim(cname2(NPO4)),isegvol(i),trim(cname2(NNH4)),isegvol(i),trim(cname2(NNO3)),isegvol(i),trim(cname2(NDO)),isegvol(i),trim(cdname2(12)),isegvol(i),trim(cdname2(14)),isegvol(i),i=1,nseg)   ! Chlor a and TP
    101 format('JDAY,',<nseg>(6((A,'-',i3,','))))
    102 format('JDAY,',<nseg>(7((A,'-',i3,','))))

    open(FISHHABFN+2,file=consurf,status='unknown')
    write(FISHHABFN+2,'(a,i4,a,80(1x,i4))')'Surface (upper',kseg,' model layers) Volume weighted WQ parameters at segments:',(isegvol(i),i=1,nseg)
    write(FISHHABFN+2,102)(trim(cname2(NPO4)),isegvol(i),trim(cname2(NNH4)),isegvol(i),trim(cname2(NNO3)),isegvol(i),trim(cname2(NDO)),isegvol(i),trim(cdname2(12)),isegvol(i),trim(cdname2(14)),isegvol(i),'Gamma(m-1)',isegvol(i),i=1,nseg)

    do jjw=1,nwb
    IF (SEDIMENT_CALC(JJW))then
    open(FISHHABFN+3,file=consod,status='unknown')
    write(FISHHABFN+3,"('JDAY,',1000(i3,','))")(((i,i=us(jb),ds(jb)),jb=bs(jw),be(jw)),jw=1,nwb)
    exit
    ENDIF
    enddo
  endif

endif
endif

! compute total volume and habitat volume
habvol=0.0
voltot=0.0
HAB=100.0
voltotbr=0.0
voltotwb=0.0
habvolbr=0.0
habvolwb=0.0
do jw=1,nwb
    do jb=bs(jw),be(jw)
        do i=cus(jb),ds(jb)
            do k=ktwb(jw),kb(i)
                voltot=voltot+vol(k,i)
                voltotbr(jb)=voltotbr(jb)+vol(k,i)
                voltotwb(jw)=voltotwb(jw)+vol(k,i)
                    do ii=IFISH,1,-1
                    if(oxygen_demand)then
                        if(t2(k,i)<=fishtemph(ii).and.t2(k,i)>fishtempl(ii).and.o2(k,i)>=fishdo(ii))then
                            habvol(ii)=habvol(ii)+vol(k,i)
                            habvolbr(jb,ii)=habvolbr(jb,ii)+vol(k,i)
                            habvolwb(jw,ii)=habvolwb(jw,ii)+vol(k,i)
                            hab(k,i)=ii
                        endif
                    else
                        if(t2(k,i)<=fishtemph(ii).and.t2(k,i)>fishtempl(ii))then
                            habvol(ii)=habvol(ii)+vol(k,i)
                            habvolbr(jb,ii)=habvolbr(jb,ii)+vol(k,i)
                            habvolwb(jw,ii)=habvolwb(jw,ii)+vol(k,i)
                            hab(k,i)=ii
                        endif
                    endif
                    enddo
            enddo
        enddo
    enddo
enddo

do ii=1,ifish
phabvol(ii)=habvol(ii)/voltot
    do jw=1,nwb
        phabvolwb(jw,ii)=habvolwb(jw,ii)/voltotwb(jw)
        do jb=bs(jw),be(jw)
        phabvolbr(jb,ii)=habvolbr(jb,ii)/voltotbr(jb)
        end do
    enddo
enddo

! write out results

write(FISHHABFN,210)jday,(100.*phabvol(i),habvol(i),i=1,ifish)
210 format(f10.3,',',<ifish>(f8.2,',',e12.4,','))
jbfile=jbfile1;jwfile=jwfile1
do jw=1,nwb
    jwfile=jwfile+1
    write(jwfile,210)jday,(100.*phabvolwb(jw,i),habvolwb(jw,i),i=1,ifish)
    do jb=bs(jw),be(jw)
        jbfile=jbfile+1
        write(jbfile,210)jday,(100.*phabvolbr(jb,i),habvolbr(jb,i),i=1,ifish)
    enddo
enddo

if(oxygen_demand)then

    cno3=0.0
    cdo=0.0
    cpo4=0.0
    cchla=0.0
    cnh4=0.0
    cgamma=0.0
    ctotp=0.0

    do n=1,nseg

    i=isegvol(n)

    ! Find waterbody associated with this segment
    do jjw=1,nwb
        if(i >= us(bs(jjw)) .and. i <= ds(be(jjw)))exit
    enddo

    voltot=0.0
    kkmax=min(kseg,kb(i)-ktwb(jjw))      ! kseg is the # of layers
    if(kkmax < 0)cycle
        do k=ktwb(jjw),kb(i)
        voltot=voltot+vol(k,i)
        cpo4(n)=cpo4(n)+po4(k,i)*vol(k,i)
        if(k <= ktwb(jjw)+kkmax)cgamma(n)=cgamma(n)+gamma(k,i)*vol(k,i)
        ! NOTE*********** No credit for superstauration - if DO > saturation, then set DO=100% saturation
        DOSAT=SATO(t2(k,i),0.d0,palt(i),SALT_WATER(jjw))
            if(o2(k,i) > DOSAT )then
            o2corr=DOSAT
            else
            o2corr=o2(k,i)
            endif

        cdo(n)=cdo(n)+o2corr*vol(k,i)
        cno3(n)=cno3(n)+no3(k,i)*vol(k,i)
        cchla(n)=cchla(n)+chla(k,i)*vol(k,i)
        ctotp(n)=ctotp(n)+tp(k,i)*vol(k,i)
        cnh4(n)=cnh4(n)+nh4(k,i)*vol(k,i)

        if(k == ktwb(jjw)+kkmax)then
        cdos(n)=cdo(n)/voltot
        cpo4s(n)=cpo4(n)/voltot
        cno3s(n)=cno3(n)/voltot
        cnh4s(n)=cnh4(n)/voltot
        cchlas(n)=cchla(n)/voltot
        ctotps(n)=ctotp(n)/voltot
        cgamma(n)=cgamma(n)/voltot
        endif

        enddo
    cpo4(n)=cpo4(n)/voltot
    cdo(n)=cdo(n)/voltot
    cno3(n)=cno3(n)/voltot
    cnh4(n)=cnh4(n)/voltot
    cchla(n)=cchla(n)/voltot
    ctotp(n)=ctotp(n)/voltot
    enddo

    write(FISHHABFN+2,211)jday,(cpo4s(n),cnh4s(n),cno3s(n),cdos(n),ctotps(n),cchlas(n),cgamma(n),n=1,nseg)
    write(FISHHABFN+1,211)jday,(cpo4(n),cnh4(n),cno3(n),cdo(n),ctotp(n),cchla(n),n=1,nseg)
    211 format(f10.3,',',<nseg>(7(f10.4,',')))

! write out sed for each segment
ssedd=0.0
do jw=1,nwb
   IF (SEDIMENT_CALC(JW))then
    do jb=bs(jw),be(jw)
        do i=us(jb),ds(jb)
           if(ktwb(jw)<=kb(i))then
            do k=ktwb(jw),kb(i)
                ssedd(i)=ssedd(i)+sed(k,i)*vol(k,i)
            end do
           else
            ssedd(i)=-99.
           endif
        enddo
    enddo
  ENDIF
enddo

do jjw=1,nwb
IF (SEDIMENT_CALC(JJW))then
write(FISHHABFN+3,"(f10.3,',',1000(e13.4,','))")jday,(((ssedd(i),i=us(jb),ds(jb)),jb=bs(jw),be(jw)),jw=1,nwb)    ! OUTPUT IS IN GRAMS
exit
ENDIF
enddo

endif

else

deallocate(isegvol,cdo,cpo4,cno3,cnh4,cchla,ctotp,cdos,cpo4s,cno3s,cnh4s,cchlas,ctotps,cgamma,ssedd,fishname,fishtempl,fishtemph,fishdo,habvol,phabvol,habvolbr,habvolwb,phabvolbr,phabvolwb,voltotbr,voltotwb)
close(FISHHABFN)
close(FISHHABFN+1)
CLOSE(FISHHABFN+2)
CLOSE(FISHHABFN+3)
jbfile=jbfile1;jwfile=jwfile1
    do jw=1,nwb
     jwfile=jwfile+1
    close(jwfile)
        do jb=bs(jw),be(jw)
        jbfile=jbfile+1
        close(jbfile)
        enddo
    enddo

endif

return
end subroutine fishhabitat


      INTEGER FUNCTION FCB_SKIP_WHITESPACE(TAPIN,LAST_CHAR,             &
		       FCB_BYTES_IN_REC,BYTE_IN_FILE,REC_IN_FILE,BUFFER,&
		       LINE,N,LINELEN,ICUR,FRESH_LINE)
!-----------------------------------------------------------------------
!     Skips forward on the current LINE of size N with data in
!     LINE(1:LINELEN) from the current position ICUR moving over
!     whitespace and comments, reading new lines into LINE if
!     needed. The flag FRESH_LINE indicates that a fresh line
!     should be read on entry.
!-----------------------------------------------------------------------
      IMPLICIT                   NONE
      INTEGER,      INTENT(IN):: TAPIN,FCB_BYTES_IN_REC,N
      INTEGER,   INTENT(INOUT):: BYTE_IN_FILE,REC_IN_FILE,LINELEN,ICUR, &
				 FRESH_LINE
      INTEGER(1),INTENT(INOUT):: BUFFER(FCB_BYTES_IN_REC),LINE(N),      &
				 LAST_CHAR
! Local variables
      INTEGER                    IC,COMMENT_LEVEL
! External functions called
      INTEGER                    FCB_READ_LINE,FCB_NBLEN_ARRAY
!-----------------------------------------------------------------------
      FCB_SKIP_WHITESPACE = 0
      IF (FRESH_LINE.NE.0) THEN
	 ICUR = LINELEN+1
	 RETURN
      END IF

      IC = ICUR
      COMMENT_LEVEL = 0
      DO
	 IF ((IC.LE.LINELEN).AND.(LINE(IC).NE.IACHAR(' ')).AND.         &
	     (LINE(IC).NE.Z'09').AND.(LINE(IC).NE.IACHAR('(')) )EXIT

	 IF (IC.GT.LINELEN) THEN
	    FCB_SKIP_WHITESPACE = FCB_READ_LINE (TAPIN,LAST_CHAR,       &
		       FCB_BYTES_IN_REC,BYTE_IN_FILE,REC_IN_FILE,BUFFER,&
			LINE,N,LINELEN)
	    IF (FCB_SKIP_WHITESPACE.NE.0) RETURN
	    IC = 1
	    IF ((LINELEN.EQ.0)                                          &
	      .OR.(FCB_NBLEN_ARRAY(LINE(1:LINELEN),LINELEN).EQ.0)       &
	      .OR.(FCB_NBLEN_ARRAY(LINE(1:1),1).NE.0)) THEN
	       FRESH_LINE = 1
	       ICUR = LINELEN+1
	       RETURN
	    END IF
	 ELSE
	    IF (LINE(IC) .EQ. IACHAR('(')) THEN
	       IC = IC+1
	       COMMENT_LEVEL = COMMENT_LEVEL+1
	       DO
		  IF (COMMENT_LEVEL.EQ.0)EXIT
		  IF (IC.GT.LINELEN) THEN
		     FCB_SKIP_WHITESPACE=FCB_READ_LINE(TAPIN,LAST_CHAR, &
			FCB_BYTES_IN_REC,BYTE_IN_FILE,REC_IN_FILE,      &
			BUFFER,LINE,N,LINELEN)
		     IF (FCB_SKIP_WHITESPACE.NE.0) RETURN
		     IC = 1
		     IF ((LINELEN.EQ.0).OR.                             &
		      (FCB_NBLEN_ARRAY(LINE(1:LINELEN),LINELEN).EQ.0)   &
			.OR. FCB_NBLEN_ARRAY(LINE(1:1),1).NE.0) THEN
			FRESH_LINE = 1
			ICUR = LINELEN+1
			RETURN
		     END IF
		  ELSE
		     SELECT CASE (LINE(IC))
		     CASE (Z'5C')      ! backslash
			IC = IC+1      ! force skip of next character
		     CASE (Z'28')      ! open paren
			COMMENT_LEVEL = COMMENT_LEVEL+1;
		     CASE (Z'29')      ! close paren
			COMMENT_LEVEL = COMMENT_LEVEL-1;
		     END SELECT
		     IC = IC+1;
		  END IF
	       END DO
	    ELSE
	       IC = IC+1;
	    END IF
	 END IF
      END DO
      FRESH_LINE = 0
      FCB_SKIP_WHITESPACE = 0
      ICUR = IC
      RETURN
      END FUNCTION FCB_SKIP_WHITESPACE


text:
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                       Fluides geophysiques  v14         (11.04.96)
!                       -------------------------
!   . Red Black
!   . Calcul des vitesses
!   . Sauvegarde tous les X pas de temps
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! ------------------------------------------------------------------

module GEOPHY

    integer, parameter :: DP  = KIND(1.0D0)            !< double precision

    real(DP), parameter :: U = 1D-1
    real(DP), parameter :: g = 1D-2
    real(DP), parameter :: f = 1D-4

    integer, parameter :: pasx  = 31
    integer, parameter :: pasy  = 31

end module GEOPHY

! ------------------------------------------------------------------
subroutine init_xy(x, y, dx, dy, Ly)
    use GEOPHY
    implicit none

    real(DP) :: x(pasx+1), y(pasx+1)   !! [RB] taille exacte necessaire sinon ca plante au runtime
    real(DP) :: dx, dy, Ly

    integer :: i, j

    do i = 2, pasx+1
        print *,"i=",i
        x(i) = (i-2)*dx
    end do

    x(1) = x(pasx)
    x(2) = x(pasx+1)

    do j = 1, pasy
        y(j) = -Ly + (j-1)*dy
    end do
end subroutine init_xy

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

program NIHOUL

    use GEOPHY
    implicit none

    ! Declarations

    integer :: Nmax
    logical :: fin
    integer :: i, j, n, compt, startj, start, nsave, fsave, saveu
    real(DP) :: pi, h0, L, T, Lx, Ly, Bu, Ro, Et, w
    real(DP) :: dx, dy, dt, tol, kmoy

    real(DP), allocatable :: x(:), y(:)

    real(DP) :: temp1, temp2, temp3, temp4, temp5, erreur
    real(DP) :: d1, d2, d3
    real(DP) :: aeta1, aeta2

    real(DP), allocatable :: eta(:,:), b(:,:), q(:,:)
    real(DP), allocatable :: jacob(:,:), deta(:,:)
    real(DP), allocatable :: leta(:,:), v1(:,:), v2(:,:)

    ! Donnees

    print *, 'L (5000)?'
    read(*,*) L
    print *, 'Lx/L (5)?'
    read(*,*) Lx

    pi    = acos(-1.0D0)
    h0    = 1D2
    T     = L/U
    Ly    = min((3.0D0/4.0D0 + g*h0/(2.0D0*f*U*L)), 4.0D0)
    Bu    = g*h0/(f*L)**2
    Ro    = U/(f*L)
    Et    = 1.0D0/(f*T)
    dx    = 2.0D0*Lx/(pasx-1)
    dt    = dx*dx/2.0D0
    dy    = 2.0D0*Ly/(pasy-1)
    kmoy  = 0.0D0
    nsave  = 1

    print *, 'Ly/L =', Ly, '(', (3.0D0/4.0D0 + g*h0/(2.0D0*f*U*L)), ')'
    print *, 'Bu   =', Bu, ' ; Ro   =', Ro
    print *, 'dx   =', dx, ' ; dy   =', dy
    print *, 'pasx =', pasx,' ; pasy =', pasy
    print *, 'T    =', T, ' ; dH   =', f*L*U/g
    print *, 'dt (max) =', dt, 'sec.'
    print *, 'Le canal fait ', 2.0D0*Lx*L, 'm de long'
    print *, '           et ', 2.0D0*Ly*L, 'm de large'
    print *, 'tol (1D-6)?'
    read(*,*) tol
    print *, 'w (1.825)?'
    read(*,*) w
    print *, 'dt (2D-2)?'
    read(*,*) dt
    print *, 'Nmax (100)?'
    read(*,*) Nmax
    print *, 'fsave (Nmax)?'
    read(*,*) fsave
    print *, 'save U (1 = oui)?'
    read(*,*) saveu

    if(fsave > Nmax) then
        fsave = Nmax
    endif

    print *
    print *, 'Un pas de temps =', T*dt, 'sec.'

    ! dynamic arrays
    allocate( x(pasx+1) )
    allocate( y(pasx) )

    allocate( eta(pasx+1, pasy) )
    allocate( b(pasx+1, pasy) )
    allocate( q(pasx+1, pasy) )

    allocate( jacob(pasx+1, pasy) )
    allocate( deta(pasx+1, 2) )

    allocate( leta(pasx+1, pasy) )
    allocate( v1(pasx+1, pasy) )
    allocate( v2(pasx+1, pasy) )

    ! Calcul du profil initial
    call init_xy(x, y, dx, dy, Ly)

    temp1 = g*h0/f/L/U
    do i = 1, pasx+1
        temp2 = 0.01*sin(pi*x(i)/Lx)
        do j = 1, pasy
            if (abs(y(j)) < 1.0D0) then
                eta(i,j) = temp1*temp2 - 0.5*y(j)**2 * (1 + temp2)
            else
                eta(i,j) = temp1*temp2 + (0.5-abs(y(j))) * (1 + temp2)
            endif
        end do
    end do

    ! Ouverture du fichier RES.M et sauvegarde du profil initial

    open (UNIT = 1, FILE = 'res.m', STATUS ='unknown')

    do i = 2, pasx+1
        do j = 1, pasy
            write(1,*) 'a(', (nsave-1)*pasx+i-1, ',', j, ')=', eta(i,j), ';'
        end do
    end do

    ! Boucle principale (progression dans le temps)

    do n = 1, Nmax
        print *, 'Pas de temps n :', n

        ! Calcul de Beta (b) et du laplacien (q)

        ! Interieur

        do i = 2, pasx
            do j = 2, pasy-1
                temp1 = (eta(i-1,j) + eta(i+1,j) - 2*eta(i,j))/(dx*dx)      &
                      + (eta(i,j-1) + eta(i,j+1) - 2*eta(i,j))/(dy*dy)
                temp2 = ((eta(i+1,j) - eta(i-1,j))/(2*dx))**2               &
                      + ((eta(i,j+1) - eta(i,j-1))/(2*dy))**2

                b(i,j) = (Bu+Ro*eta(i,j))*temp1 + Ro/2*temp2
                q(i,j) = temp1
            end do
        end do

        ! Bords y

        do i = 2, pasx
            j = 1
            temp1 = (eta(i-1,j) + eta(i+1,j) - 2*eta(i,j))   /(dx*dx)       &
                  + (eta(i,j)   + eta(i,j+2) - 2*eta(i,j+1)) /(dy*dy)
            temp2 = ((  eta(i+1,j) - eta(i-1,j)) / (2*dx))**2               &
                  + ((4*eta(i,j+1) - eta(i,j+2) - 3*eta(i,j)) / (2*dy))**2

            b(i,j) = (Bu+Ro*eta(i,j))*temp1 + Ro/2*temp2
            q(i,j) = temp1

            j = pasy
            temp1 = (eta(i-1,j) + eta(i+1,j) - 2*eta(i,j))   /(dx*dx)       &
                  + (eta(i,j)   + eta(i,j-2) - 2*eta(i,j-1)) /(dy*dy)
            temp2 = ((eta(i+1,j) - eta(i-1,j))/(2*dx))**2                   &
                  + ((-4*eta(i,j-1) + eta(i,j-2) + 3*eta(i,j)) / (2*dy))**2

            b(i,j) = (Bu+Ro*eta(i,j))*temp1 + Ro/2*temp2
            q(i,j) = temp1
        end do

        ! Ajout des lignes i = 1 et i = pasx+1

        do j = 1, pasy
            b(1,j) = b(pasx,j)
            b(pasx+1,j) = b(2,j)

            q(1,j) = q(pasx,j)
            q(pasx+1,j) = q(2,j)
        end do

        ! Calcul du Jacobien

        do i = 2, pasx
            do j = 2, pasy-1
                temp1 = ((eta(i+1,j) - eta(i-1,j))*(b(i,j+1) - b(i,j-1)) -         &
                (eta(i,j+1) - eta(i,j-1))*(b(i+1,j) - b(i-1,j)))/4/dx/dy
                temp2 = (eta(i+1,j)*(b(i+1,j+1) - b(i+1,j-1)) - eta(i-1,j)*        &
                    (b(i-1,j+1) - b(i-1,j-1)) - eta(i,j+1)*(b(i+1,j+1)             &
                    -b(i-1,j+1)) + eta(i,j-1)*(b(i+1,j-1) - b(i-1,j-1)))           &
                    /4/dx/dy
                temp3 = (b(i,j+1)*(eta(i+1,j+1) - eta(i-1,j+1)) - b(i,j-1)*        &
                    (eta(i+1,j-1) - eta(i-1,j-1)) - b(i+1,j)*(eta(i+1,j+1) -       &
                    eta(i+1,j-1)) + b(i-1,j)*(eta(i-1,j+1) - eta(i-1,j-1)))        &
                    /4/dx/dy
                jacob(i,j) = (temp1 + temp2 + temp3)/3.0D0
            end do
        end do

        ! Ajout des lignes i = 1 et i = pasx+1

        do j = 1,pasy
            jacob(1,j) = jacob(pasx,j)
            jacob(pasx+1,j) = jacob(2,j)
        end do

        ! Calcul de q

        do i = 1, pasx+1
            do j = 1, pasy
                q(i,j) = Bu*q(i,j) - eta(i,j)
                q(i,j) = q(i,j) - Ro/Et*dt*jacob(i,j)
            end do
        end do

        ! Calcul de d(eta)/dy pour les bords

        do i = 1, pasx+1
            deta(i,1) = (4*eta(i,2) - eta(i,3) - 3*eta(i,1)) / (2*dy)
            deta(i,2) = (-4*eta(i,pasy-1) + eta(i,pasy-2) + 3*eta(i,pasy)) / (2*dy)
        end do

        ! Red Black

        fin = .false.
        compt = 0

        do while (fin == 0)
            compt = compt + 1
            fin = .true.

            ! Deux passages alternes

            do startj = 1, 2
                start = startj

                ! Red Black

                do i = 2, pasx
                    do j = start+1, pasy-1, 2
                        temp1 = -2/dx**2 - 2/dy**2 - 1/Bu
                        temp1 = (q(i,j)/Bu - (eta(i+1,j) + eta(i-1,j))/dx**2   &
                                           - (eta(i,j+1) + eta(i,j-1))/dy**2  ) / temp1
                        temp1 = eta(i,j) + w*(temp1-eta(i,j))
                        if(abs(eta(i,j)) > 0.0D0) then
                            erreur = abs( (temp1-eta(i,j))/eta(i,j) )
                            if(erreur > tol) then
                                fin = .false.
                            endif
                        endif
                        eta(i,j) = temp1
                    end do

                    if(start == 2) then
                        start = 1
                    else
                        start = 2
                    endif
                end do

                do j = 1, pasy
                    eta(1,j) = eta(pasx,j)
                    eta(pasx+1,j) = eta(2,j)
                end do

            end do

            ! Conditions aux limites

            aeta1 = eta(1,1)
            aeta2 = eta(1,pasy)

            do i = 2, pasx
                j = 1
                temp1 = (eta(i+1,j) - aeta1) / (2*dx)
                temp2 = (eta(i,j+2) - 2*eta(i,j+1)) / dy**2
                d1 = (eta(i+1,j) - aeta1) / (2*dx)
                d2 = (eta(i+1,j+1) - eta(i-1,j+1))/(2*dx)
                d3 = (eta(i+1,j+2) - eta(i-1,j+2))/(2*dx)
                temp3 = (4*d2 - d3 - 3*d1)/(2*dy)
                temp4 = -3*Et/(2*dy*dt) + 3*Ro*temp3/(2*dy) + Ro*temp1/dy**2
                temp5 = (temp1  -Et*((4*eta(i,j+1) - eta(i,j+2)) / (2*dy*dt)    &
                        -deta(i,1)/dt) - Ro*temp1*temp2 + Ro*temp3*             &
                        (4*eta(i,j+1) - eta(i,j+2))/(2*dy) ) / temp4
                temp5 = eta(i,j) + w*(temp5-eta(i,j))
                if(abs(eta(i,j)) > 0.0D0) then
                    erreur = abs( (temp5-eta(i,j))/eta(i,j) )
                    if(erreur > tol) then
                        fin = .false.
                    endif
                endif
                aeta1 = eta(i,j)
                eta(i,j) = temp5

                j = pasy
                temp1 = (eta(i+1,j) - aeta2) / (2*dx)
                temp2 = (eta(i,j-2) - 2*eta(i,j-1)) / dy**2
                d1 = (eta(i+1,j) - aeta2) / (2*dx)
                d2 = (eta(i+1,j-1) - eta(i-1,j-1)) / (2*dx)
                d3 = (eta(i+1,j-2) - eta(i-1,j-2)) / (2*dx)
                temp3 = (-4*d2 + d3 + 3*d1) / (2*dy)
                temp4 = 3*Et/(2*dy*dt) - 3*Ro*temp3/(2*dy) + Ro*temp1/dy**2
                temp5 = (temp1 - Et*((-4*eta(i,j-1) + eta(i,j-2)) / (2*dy*dt)    &
                        - deta(i,2)/dt) - Ro*temp1*temp2 + Ro*temp3*             &
                        (-4*eta(i,j-1) + eta(i,j-2))/(2*dy) ) / temp4
                temp5 = eta(i,j) + w*(temp5-eta(i,j))
                if (abs(eta(i,j)) > 0.0D0) then
                    erreur = abs( (temp5-eta(i,j)) / eta(i,j) )
                    if (erreur > tol) then
                        fin = .false.
                    endif
                endif
                aeta2 = eta(i,j)
                eta(i,j) = temp5
            end do

            ! Periodicite
            do j = 1, pasy
                eta(1,j) = eta(pasx,j)
                eta(pasx+1,j) = eta(2,j)
            end do

        end do

        print *, '    no it. =', compt
        kmoy = kmoy + compt

        ! Copie de eta (futur calcul des vitesses)

        if(n == nsave*fsave-1) then
            do i = 1, pasx+1
                do j = 1, pasy
                    leta(i,j) = eta(i,j)
                end do
            end do
        endif

        ! Calcul des vitesses

        if(saveu == 1) then
            if(n == nsave*fsave) then
                do i = 2, pasx
                    do j = 2, pasy-1

                        ! u (interieur)

                        temp1 = (eta(i,j+1) - eta(i,j-1))/2/dy
                        temp2 = (eta(i+1,j) - eta(i-1,j))/2/dx
                        temp3 = (leta(i+1,j) - leta(i-1,j))/2/dx
                        temp4 = (eta(i+1,j+1) - eta(i-1,j+1) - eta(i+1,j-1)     &
                                + eta(i-1,j-1)) / (4*dy*dy)
                        temp5 = (eta(i+1,j) + eta(i-1,j) - 2*eta(i,j)) / (dx*dx)
                        v1(i,j) = -temp1 - Et*(temp2 - temp3)/dt           &
                                - Ro*(temp2*temp4 - temp1*temp5)

                        ! v (interieur)

                        temp3 = (leta(i,j+1) - leta(i,j-1))/2/dy
                        temp5 = (eta(i,j+1) + eta(i,j-1) - 2*eta(i,j))/(dy*dy)
                        v2(i,j) = temp2 - Et*(temp1 - temp3)/dt            &
                                - Ro*(temp2*temp5 - temp1*temp4)
                    end do

                    ! u (j=1)
                    j = 1
                    temp1 = (4*eta(i,j+1)  -eta(i,j+2) - 3*eta(i,j))/2/dy
                    temp2 = (eta(i+1,j) - eta(i-1,j))/2/dx
                    temp3 = (leta(i+1,j) - leta(i-1,j))/2/dx
                    d1 = temp2
                    d2 = (eta(i+1,j+1) - eta(i-1,j+1))/(2*dx)
                    d3 = (eta(i+1,j+2) - eta(i-1,j+2))/(2*dx)
                    temp4 = (4*d2 - d3 - 3*d1)/(2*dy)
                    temp5 = (eta(i+1,j) + eta(i-1,j) - 2*eta(i,j))/(dx*dx)
                    v1(i,j) = -temp1 - Et*(temp2-temp3)/dt            &
                              -Ro*(temp2*temp4 - temp1*temp5)

                    ! v (j=1)
                    temp3 = (4*leta(i,j+1) - leta(i,j+2) - 3*leta(i,j))/2/dy
                    temp5 = (eta(i,j+2) + eta(i,j) - 2*eta(i,j+1))/(dy*dy)
                    v2(i,j) = temp2-Et*(temp1 - temp3)/dt             &
                             - Ro*(temp2*temp5 - temp1*temp4)

                    ! u (j=pasy)
                    j = pasy
                    temp1 = (-4*eta(i,j-1) + eta(i,j-2) + 3*eta(i,j))/2/dy
                    temp2 = (eta(i+1,j) - eta(i-1,j))/2/dx
                    temp3 = (leta(i+1,j) - leta(i-1,j))/2/dx
                    d1 = temp2
                    d2 = (eta(i+1,j-1) - eta(i-1,j-1))/(2*dx)
                    d3 = (eta(i+1,j-2) - eta(i-1,j-2))/(2*dx)
                    temp4 = (-4*d2 + d3 + 3*d1)/(2*dy)
                    temp5 = (eta(i+1,j) + eta(i-1,j) - 2*eta(i,j))/(dx*dx)
                    v1(i,j) = -temp1 - Et*(temp2 - temp3)/dt            &
                               - Ro*(temp2*temp4 - temp1*temp5)

                    ! v (j=pasy)
                    temp3 = (-4*leta(i,j-1) + leta(i,j-2) + 3*leta(i,j))/2/dy
                    temp5 = (eta(i,j-2) + eta(i,j) - 2*eta(i,j-1))/(dy*dy)
                    v2(i,j) = temp2 - Et*(temp1 - temp3)/dt         &
                             - Ro*(temp2*temp5 - temp1*temp4)
                end do

                ! periodicite

                do j = 1, pasy
                    v1(pasx+1,j) = v1(2,j)
                    v2(pasx+1,j) = v2(2,j)
                end do
            endif
        endif

        ! Resultats (save -> Matlab)

        if(n == nsave*fsave) then
            do i = 2, pasx+1
                do j = 1, pasy
                    write(1,*) 'a(', nsave*pasx+i-1, ',', j, ')=', eta(i,j), ';'
                    if(saveu == 1) then
                        write(1,*) 'u(', nsave*pasx+i-1, ',', j, ')=', v1(i,j), ';'
                        write(1,*) 'v(', nsave*pasx+i-1, ',', j, ')=', v2(i,j), ';'
                    endif
                end do
            end do
            nsave = nsave + 1
        endif

    end do

    close (UNIT=1)

    !  Commentaires finaux
    print "('Calcul effectue sur :', F10.1,'s')", Nmax*T*dt
    print "('Nb iterations (moy) :', F5.1)", kmoy/Nmax
    print "('Resultats dans RES.M')"

    ! free memory
    deallocate( x )
    deallocate( y )

    deallocate( eta )
    deallocate( b )
    deallocate( q )

    deallocate( jacob )
    deallocate( deta )

    deallocate( leta )
    deallocate( v1 )
    deallocate( v2 )

end program NIHOUL

max_stars_repo_name: rboman/progs

max_issues_repo_name: rboman/progs

max_forks_repo_name: rboman/progs


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE LNSTOR (A, ISTEP, TYP, NWRDS, NVAR, NPT, NPTS,
     &   XLN, YLN, ZLN, DATA)
C=======================================================================

C   --*** LNSTOR *** (PATHLN) Read and store pathline data from database
C   --   Written by Amy Gilkey - revised 05/27/88
C   --
C   --LNSTOR reads variables from the database and stores any that are
C   --pathline data in the appropriate location.  It reads only one
C   --variable type (history, global, nodal, element) per call.
C   --
C   --Parameters:
C   --   A - IN - the dynamic memory base array
C   --   ISTEP - IN - the time step number
C   --   TYP - IN - the type of variable:
C   --      'H'istory, 'G'lobal, 'N'odal, 'E'lement
C   --   NWRDS - IN - the number of words in a data record
C   --   NVAR - IN - the number of records to read
C   --   NPT - IN - the XLN, YLN, ZLN time index to fill
C   --   NPTS - IN - the maximum XLN, YLN, ZLN time index
C   --   XLN, YLN, ZLN - IN/OUT - the pathline data array
C   --   DATA - SCRATCH - size = NWRDS
C   --
C   --Common Variables:
C   --   Uses NLNCRV, ILVNE, ILVID of /LNVARS/

      include 'lnvars.blk'
      include 'dbnums.blk'

      DIMENSION A(*)
      CHARACTER TYP
      REAL XLN(NPTS,NLNCRV), YLN(NPTS,NLNCRV), ZLN(NPTS,NLNCRV)
      REAL DATA(NWRDS)

      LOGICAL NEED
      CHARACTER T

      CALL DBVIX_BL (TYP, 1, ISID)
      CALL DBVIX_BL (TYP, NVAR, IEID)

      IF ((TYP .EQ. 'H') .OR. (TYP .EQ. 'G')) THEN

C      --Read history or global variables

         CALL GETVAR (A, ISID, -999, ISTEP, NWRDS, DATA)

C      --Scan plot variable information and store if a history or global
C      --point from this data is to be plotted

         DO 100 NP = 1, NLNCRV
            CALL DBVTYP_BL (ILVID(1,NP), T, IDUM)
            IF ((T .EQ. 'H') .OR. (T .EQ. 'G')) THEN
               CALL DBVTYP_BL (ILVID(1,NP), T, NE)
               XLN(NPT,NP) = DATA(NE)
               CALL DBVTYP_BL (ILVID(2,NP), T, NE)
               YLN(NPT,NP) = DATA(NE)
               IF (NDIM .GE. 3) THEN
                  CALL DBVTYP_BL (ILVID(3,NP), T, NE)
                  ZLN(NPT,NP) = DATA(NE)
               END IF
            END IF
  100    CONTINUE

      ELSE
         DO 130 ID = ISID, IEID

C         --Determine if variable is needed

            NEED = .FALSE.
            DO 110 NP = 1, NLNCRV
               NEED = NEED .OR. (LOCINT (ID, NDIM, ILVID(1,NP)) .GT. 0)
  110       CONTINUE

            IF (NEED) THEN

C            --Read nodal/element variable

               CALL GETVAR (A, ID, -1, ISTEP, NWRDS, DATA)

C            --Scan plot variable information and store if a node/element
C            --point from this data is to be plotted

               DO 120 NP = 1, NLNCRV
                  NE = ILVNE(NP)
                  IF (ID .EQ. ILVID(1,NP)) XLN(NPT,NP) = DATA(NE)
                  IF (ID .EQ. ILVID(2,NP)) YLN(NPT,NP) = DATA(NE)
                  IF (NDIM .GE. 3) THEN
                     IF (ID .EQ. ILVID(3,NP)) ZLN(NPT,NP) = DATA(NE)
                  END IF
  120          CONTINUE
            END IF

  130    CONTINUE
      END IF

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
program literals
  ENUM, BIND(C)
    ENUMERATOR :: RED = 4, BLUE = 9
    ENUMERATOR YELLOW
  END ENUM

  a = 1
  b = 1_8
  c = 1.1
  d = 1.1_4
  e = 1e4_4+2e+4_8
  f = 1_"Deidranna &
  &rules!"
  str = 'Very ''&\\ bad & &
  & &
  & string\'
  g = B"101"
  h = O"1765"
  i = Z'9A7f'
  j = .TRUE._1
  k = ( 10.0, +4_8)
  l = (/ 10, 9, 8, 7 /)
  m = [ 1, 2, 3 ]
endprogram

max_stars_repo_name: LChernigovskaya/fortran-plugin

max_issues_repo_name: LChernigovskaya/fortran-plugin

max_forks_repo_name: LChernigovskaya/fortran-plugin


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:58:48 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE ISEND_CR__genmod
          INTERFACE
            SUBROUTINE ISEND_CR(SENDBUF,SENDCOUNT,DEST,SENDTAG,REQ)
              INTEGER(KIND=4) :: SENDCOUNT
              REAL(KIND=4) :: SENDBUF(SENDCOUNT)
              INTEGER(KIND=4) :: DEST
              INTEGER(KIND=4) :: SENDTAG
              INTEGER(KIND=4) :: REQ
            END SUBROUTINE ISEND_CR
          END INTERFACE
        END MODULE ISEND_CR__genmod


text:
	program qm
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccc Coupled Morse oscillators ccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	implicit real*8 (a-h,o-z)
	integer*4 ni,nj
	parameter (ni=256, nj=256)
	real*8 ak2(ni,nj),mx,my
	complex*16 v(ni,nj)
	complex*16 im,psi(ni,nj),hpsi(ni,nj)
	complex*16 psi0(ni,nj),c1,c2
	common /params/ d,x0,y0,di,wi
	common /grid/ xmin,ymin,xmax,ymax,dx,dy
C
C
	im=(0.d0,1.d0)
	pi=4.d0*atan(1.d0)
	N2=ni*nj
	d=7d0
	X0=0d0
	y0=0d0
c	z0=1.04435d0
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	open(15,file='cor1')
	open(16,file='cor2')
	open(17,file='wf0')
	open(18,file='wft')
	open(19,file='xav')
	open(5,file='IN')
c time propagation parameters
	read(5,*) dt,kmax,nout
c grid spacing

	read(5,*) dx,dy

	read(5,*) mx,my
c initial wave packet parameters
	read(5,*) ax,ay
	read(5,*) qx0,qy0
	read(5,*) px0,py0
c potential surface parameters
        read(5,*)DI,WI
	close(5)
c inverse mass prefactor 1/2/m
c hydrogen
c deuterium m=2*mh
c	am=0.5d0
c	alfa=alfa*sqrt(2.d0)
cc tritium=3*mh
c	am=1.d0/3.d0
cc	alfa=alfa*sqrt(3.d0)
cc
	xmin=0.1d0
	xmax=xmin+ni*dx
	ymin=0.1d0
	ymax=ymin+nj*dy
	ds=dx*dy
        qx0=qx0+x0
	qy0=qy0+y0
        write(*,*)'initial state',ax,ay,qx0,qy0,px0,py0
c define vector ak2
	consI=2.d0*pi/dx/ni
	consJ=2.d0*pi/dy/nj
	ni2=ni/2
	nj2=nj/2
	do 10 j=1,nj
	 akj=(j-1)*consj
	 if (j.gt.nj2) akj=-(nj+1-j)*consj
	 do 20 i=1,ni
	  aki=(i-1)*consi
	  if (i.gt.ni2) aki=-(ni+1-i)*consi
20	ak2(i,j)=aki*aki/(2d0*mx)+akj*akj/(2d0*my)
10	continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c------------ set up the initial wavefunction -----------------------
        an0=0d0
        ch_n=0d0
	do j=1,nj
	 y=ymin+dy*(j-1)
	 do i=1,ni
	  x=xmin+dx*(i-1)
          psi(i,j)=exp(-ax*(x-qx0)**2-ay*(y-qy0)**2
     &          +im*px0*(x-qx0)+im*py0*(y-qy0))
	  an0=an0+abs(psi(i,j))**2*ds
	enddo
	enddo
	an0=1d0/dsqrt(an0)

	 do j=1,nj
         y=ymin+dy*(j-1)
          do i=1,ni
           x=xmin+dx*(i-1)
          psi(i,j)=an0*psi(i,j)
	  psi0(i,j)=psi(i,j)
	   !if(abs(psi(i,j)).gt.1d-3)
	 write(17,*) x,y,abs(psi(i,j))
         ch_n=ch_n+conjg(psi(i,j))*psi(i,j)*ds
	enddo
	write(17,*) ' '
	enddo
        !continue
        write(*,*) 'NORM=',ch_n
cccccccccccccccccc propagate psi ccccccccccccccccccccccccccccccccc
	call ham(ni,nj,N2,ak2,psi,v,hpsi,en0)
	call correl(ni,nj,psi0,psi,c1,c2)
	t=0.d0
	write(15,1000) t,c1,abs(c1)
	write(16,1000) 2d0*t,c2,abs(c2)
	call split(kmax,nout,ni,nj,N2,dt,t,v,ak2,psi0,psi)
c------------ print the final  wavefunction -----------------------
	do j=1,nj
         y=ymin+(j-1)*dy
	 do i=1,ni
	  x=xmin+dx*(i-1)
	  !if(abs(psi(i,j)).gt.1d-3)
	write(18,*) x,y,abs(psi(i,j))
	enddo
	write (18,*) ' '
	enddo

c--------------------------------------------------------------------
1000    format(20(e14.7,1x))
	stop
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	subroutine ham(ni,nj,N2,ak2,psi,v,hpsi,en0)
	implicit real*8(a-h,o-z)
	parameter(nij=256)
	complex*16 psi(nij,nij),hpsi(nij,nij),tpsi(nij,nij),vpsi(nij,nij)
	dimension ak2(nij,nij)
	complex*16 v(nij,nij),c
	common /grid/ xmin,ymin,xmax,ymax,dx,dy
	call potpsi(ni,nj,v,psi,vpsi)
	call aver(ni,nj,psi,vpsi,C)
	write(*,*) 'potential energy',C
	en0=dreal(c)
	call kinpsi(ni,nj,N2,ak2,psi,tpsi)
	call aver(ni,nj,psi,tpsi,C)
	write(*,*) 'kinetic energy', C
	en0=en0+dreal(c)
	write(*,*) 'TOTAL E', en0
	do 11 j=1,nj
	 do 11 i=1,ni
11	  hpsi(i,j)=tpsi(i,j)+vpsi(i,j)
	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	subroutine potpsi(ni,nj,v,psi,vpsi)
	implicit real*8(a-h,o-z)
	parameter(nij=256)
	complex*16 psi(nij,nij),vpsi(nij,nij),v(nij,nij),im
	common /params/ d,x0,y0,di,wi
	common /grid/ xmin,ymin,xmax,ymax,dx,dy
c------------- set up potetial -------------------------------------------
        zx=1.0435d0
        zy=1.0435d0
        im=(0d0,1d0)

        do 10 j=1,nj
         y=ymin+(j-1)*dy
         ry=exp(-zy*(y-y0))
         do 10 i=1,ni
          x=xmin+(i-1)*dx
          rx=exp(-zx*(x-x0))
          vr=d*(rx-1d0)**2+d*(ry-1d0)**2
          if(vr.gt.5d0*d) vr=5d0*d
          vi=0d0
          r=dsqrt(rx*rx+ry*ry)
          if(r.gt.wi) vi=vi+di*(r-wi)**2
          v(i,j)=vr-im*vi
          vpsi(i,j)=psi(i,j)*v(i,j)
10      continue
        return
        end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        subroutine kinpsi(ni,nj,N2,ak2,cx1,aux)
        implicit real*8(a-h,o-z)
	parameter(nij=256)
        dimension ak2(nij,nij)
        complex*16 cx1(nij,nij),aux(nij,nij)
        dimension nn(2)
        nn(1)=ni
        nn(2)=nj
        do 10 j=1,nj
         do 10 i=1,ni
10      aux(i,j)=cx1(i,j)
        call fourn(aux,nn,2,1)
        do 11 j=1,nj
         do 11 i=1,ni
11      aux(i,j)=aux(i,j)*ak2(i,j)
        isign=-1
        call fourn(aux,nn,2,-1)
        do 12 j=1,nj
         do 12 i=1,ni
12      aux(i,j)=aux(i,j)/N2
        return
        end

c-------------------------------------------------------
c-------------------------------------------------------
        subroutine split(kmax,nout,ni,nj,N2,dt,t,v,ak2,psi0,cwf)
        implicit real*8(a-h,o-z)
	parameter(ij=256)
        dimension ak2(ij,ij)
        dimension nn(2)
        complex*16 cwf(ij,ij),v(ij,ij)
        complex*16 psi0(ij,ij),hpsi(ij,ij),c1,c2
	common /grid/ xmin,ymin,xmax,ymax,dx,dy
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        dt2=dt*0.5
	n2=ni*nj
        nn(1)=ni
        nn(2)=nj
        do 11 k=1,kmax
	 t=t+dt
         call fourn(cwf,nn,2,1)
         call diff(cwf,ak2,dt2,ni,nj)
         call fourn(cwf,nn,2,-1)
         do 12 j=1,nj
          do 12 i=1,ni
12       cwf(i,j)=cwf(i,j)/N2

         call phase(cwf,v,dt2,ni,nj)

         call fourn(cwf,nn,2,1)
         call diff(cwf,ak2,dt2,ni,nj)
         call fourn(cwf,nn,2,-1)
         do 13 j=1,nj
          do 13 i=1,ni
13       cwf(i,j)=cwf(i,j)/N2

        call phase(cwf,v,dt2,ni,nj)
	call correl(ni,nj,psi0,cwf,c1,c2)
	call xex(ni,nj,cwf,xav,yav)
	write(15,1000) t,c1,abs(c1)
	write(16,1000) 2d0*t,c2,abs(c2)
	write(19,1000) t,xav,yav
11	enddo

1000    format(20(e14.7,1x))
15	return
	end

c-------------------------------------------------------
c-------------------------------------------------------
	subroutine diff(cwf,ak2,ts,ni,nj)
	implicit real*8(a-h,o-z)
	parameter(nij=256)
	complex*16 nim,cwf(nij,nij)
	dimension ak2(nij,nij)
	nim=(0.d0,-1.d0)
	do 11 j=1,nj
	 do 11 i=1,ni
11	cwf(i,j)=cwf(i,j)*exp(nim*ts*ak2(i,j))
	return
	end

	subroutine phase(cwf,v,ts,ni,nj)
	implicit real*8(a-h,o-z)
	complex*16 nim,cwf(ni,nj)
	complex*16 v(ni,nj)
	nim=(0.d0,-1.d0)
	do 11 j=1,nj
	 do 11 i=1,ni
	if (ts.ge.0) cwf(i,j)=cwf(i,j)*exp(nim*ts*v(i,j))
11	if (ts.lt.0) cwf(i,j)=cwf(i,j)*exp(nim*ts*conjg(v(i,j)))
	return
	end

c-------------------------------------------------------
c-------------------------------------------------------
	subroutine correl(ni,nj,psi0,psi,c1,c2)
	implicit real*8(a-h,o-z)
	parameter(ij=256)
	complex*16 psi0(ij,ij),psi(ij,ij),c1,c2
	common /grid/ xmin,ymin,xmax,ymax,dx,dy
ccc---------- compute correlation functions ------------------------
	c1=(0d0,0d0)
	c2=(0d0,0d0)
	do 10 j=1,nj
	 y=ymin+(j-1)*dy
	 do 10 i=1,ni
	  x=xmin+(i-1)*dx
          c1=c1+conjg(psi0(i,j))*psi(i,j)
          c2=c2+psi(i,j)**2
10	continue
        c1=c1*dx*dy
        c2=c2*dx*dy
	return
	end

c-------------------------------------------------------
c-------------------------------------------------------
        subroutine aver(ni,nj,psi0,psi,c)
        implicit real*8(a-h,o-z)
        parameter(ij=256)
        complex*16 c,psi0(ij,ij),psi(ij,ij)
	common /grid/ xmin,ymin,xmax,ymax,dx,dy
	c=(0.d0,0.d0)
	ct=0.d0
        do 11 j=1,nj
         do 11 i=1,ni
	ct=ct+abs(psi0(i,j))**2
11      c=c+psi(i,j)*conjg(psi0(i,j))
        c=c/dx/dy
        return
        end

      SUBROUTINE fourn(data,nn,ndim,isign)
      INTEGER isign,ndim,nn(ndim)
      REAL*8 data(*)
      INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1,
     *k2,n,nprev,nrem,ntot
      REAL*8 tempi,tempr
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      ntot=1
      do 11 idim=1,ndim
        ntot=ntot*nn(idim)
11    continue
      nprev=1
      do 18 idim=1,ndim
        n=nn(idim)
        nrem=ntot/(n*nprev)
        ip1=2*nprev
        ip2=ip1*n
        ip3=ip2*nrem
        i2rev=1
        do 14 i2=1,ip2,ip1
          if(i2.lt.i2rev)then
            do 13 i1=i2,i2+ip1-2,2
              do 12 i3=i1,ip3,ip2
                i3rev=i2rev+i3-i2
                tempr=data(i3)
                tempi=data(i3+1)
                data(i3)=data(i3rev)
                data(i3+1)=data(i3rev+1)
                data(i3rev)=tempr
                data(i3rev+1)=tempi
12            continue
13          continue
          endif
          ibit=ip2/2
1         if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then
            i2rev=i2rev-ibit
            ibit=ibit/2
          goto 1
          endif
          i2rev=i2rev+ibit
14      continue
        ifp1=ip1
2       if(ifp1.lt.ip2)then
          ifp2=2*ifp1
          theta=isign*6.28318530717959d0/(ifp2/ip1)
          wpr=-2.d0*sin(0.5d0*theta)**2
          wpi=sin(theta)
          wr=1.d0
          wi=0.d0
          do 17 i3=1,ifp1,ip1
            do 16 i1=i3,i3+ip1-2,2
              do 15 i2=i1,ip3,ifp2
                k1=i2
                k2=k1+ifp1
                tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
                tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
                data(k2)=data(k1)-tempr
                data(k2+1)=data(k1+1)-tempi
                data(k1)=data(k1)+tempr
                data(k1+1)=data(k1+1)+tempi
15            continue
16          continue
            wtemp=wr
            wr=wr*wpr-wi*wpi+wr
            wi=wi*wpr+wtemp*wpi+wi
17        continue
          ifp1=ifp2
        goto 2
        endif
        nprev=n*nprev
18    continue
      return
      END


text:
C
C***********************************************************************
C                                                                      *
C            NACHOS II - A Finite Element Computer Program             *
C                        for Incompressible Flow Problems              *
C                                                                      *
C     Copyright (c) 1986,2019   National Technology & Engineering      *
C                               Solutions of Sandia, LLC (NTESS)       *
C                                                                      *
C                            All rights reserved.                      *
C                                                                      *
C     This software is distributed under the BSD 3-Clause License.     *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE CONVBC (KIND,IN,VALUE,NGAUSS,GS,GSWT,X,Y,KTQ,FTQ)
C
C     ******************************************************************
C
C     SUBROUTINE TO CALCULATE FLUX TERMS FOR A CONVECTIVE BOUNDARY
C     CONDITION (CONSTANT H AND TREF)
C
C     ******************************************************************
C
      REAL KTQ
C
      COMMON /PRBDAT/ IAXSYM,ITMDEP,IFORCE,IFREE,IVAR1,IVAR2,IPFUNC,
     1                IPNLTY,PNLTY
      COMMON /ELMDAT/ NNELM(6),NNCOR(6),NNMID(6),NNCTR(6),NSIDET(5),
     1                NSIDEQ(7),NNSIDE(6,4,3)
      COMMON /BCDAT/  BCSET(2,20,2),IPELM,IPNODE,PBC
C
      DIMENSION X(*), Y(*)
      DIMENSION GS(*), GSWT(*)
      DIMENSION XX(3), YY(3), INODE(3)
      DIMENSION KTQ(9,9), FTQ(9)
C
C     STATEMENT FUNCTION - ONE-DIMENSIONAL SHAPE FUNCTION
C
      FS(I,S)=(S**2-S)*.5*(I-2)*(I-3)*.5-(1.-S**2)*(I-1)*(I-3)+(S**2+S)*
     1.5*(I-1)*(I-2)*.5
C
C     ******************************************************************
C
C     LOCATE DATA FOR ELEMENT SIDE
C
      NN=NNELM(KIND)
      DO 10 J=1,3
      INODE(J)=NNSIDE(KIND,IN,J)
   10 CONTINUE
      DO 20 J=1,3
      I=INODE(J)
      XX(J)=X(I)
      YY(J)=Y(I)
   20 CONTINUE
C
C     COMPUTATION OF FLUX VECTOR
C
      IVALUE=INT(VALUE)
      VALUE=BCSET(1,IVALUE,1)*BCSET(1,IVALUE,2)
      DO 70 I=1,3
      FD=0.
      DO 60 K=1,NGAUSS
      S=GS(K)
      WW=GSWT(K)
      A=FS(I,S)
      GO TO (30, 40, 30, 40, 30, 40), KIND
   30 CONTINUE
      RAD=(1.-S)*XX(1)*.5+(1.+S)*XX(3)*.5
      C=(XX(3)-XX(1))*.5
      D=(YY(3)-YY(1))*.5
      GO TO 50
   40 CONTINUE
      RAD=(S**2-S)*XX(1)*.5+(1.-S**2)*XX(2)+(S**2+S)*XX(3)*.5
      C=(2.*S-1.)*XX(1)*.5+(-2.*S)*XX(2)+(2.*S+1.)*XX(3)*.5
      D=(2.*S-1.)*YY(1)*.5+(-2.*S)*YY(2)+(2.*S+1.)*YY(3)*.5
   50 CONTINUE
      IF (IAXSYM.EQ.0) RAD=1.0
      DELTA=SQRT(C**2+D**2)
      FD=FD+WW*A*RAD*DELTA
   60 CONTINUE
      II=INODE(I)
      FTQ(II)=FTQ(II)+FD*VALUE
   70 CONTINUE
C
C     COMPUTATION OF DIFFUSION CORRECTION MATRIX
C
      DO 120 J=1,3
      DO 120 I=1,3
      FD=0.
      DO 110 K=1,NGAUSS
      S=GS(K)
      WW=GSWT(K)
      A=FS(I,S)
      E=FS(J,S)
      GO TO (80, 90, 80, 90, 80, 90), KIND
   80 CONTINUE
      RAD=(1.-S)*XX(1)*.5+(1.+S)*XX(3)*.5
      C=(XX(3)-XX(1))*.5
      D=(YY(3)-YY(1))*.5
      GO TO 100
   90 CONTINUE
      RAD=(S**2-S)*XX(1)*.5+(1.-S**2)*XX(2)+(S**2+S)*XX(3)*.5
      C=(2.*S-1.)*XX(1)*.5+(-2.*S)*XX(2)+(2.*S+1.)*XX(3)*.5
      D=(2.*S-1.)*YY(1)*.5+(-2.*S)*YY(2)+(2.*S+1.)*YY(3)*.5
  100 CONTINUE
      IF (IAXSYM.EQ.0) RAD=1.0
      DELTA=SQRT(C**2+D**2)
      FD=FD+WW*A*RAD*DELTA*E
  110 CONTINUE
      II=INODE(I)
      JJ=INODE(J)
      KTQ(II,JJ)=KTQ(II,JJ)+FD*BCSET(1,IVALUE,1)
  120 CONTINUE
C
      RETURN
      END

max_stars_repo_name: nachos-ii/nachos-ii

max_issues_repo_name: nachos-ii/nachos-ii

max_forks_repo_name: nachos-ii/nachos-ii


module numa_f

  use, intrinsic :: iso_c_binding
  implicit none

  interface
     subroutine move_data_pages(addr, sz, to_node, info) bind(c,name='move_data_pages_')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: addr
       integer(c_size_t), intent(in) :: sz
       integer(c_int), intent(in) :: to_node
       integer(c_long), intent(out) :: info
     end subroutine move_data_pages
  end interface

contains

  subroutine move_data_pages_i(addr, sz, to_node, info)
    use, intrinsic :: iso_c_binding
    implicit none
    integer, intent(in), target :: addr
    integer, intent(in) :: sz
    integer, intent(in) :: to_node
    integer, intent(out) :: info

    type(c_ptr) :: addr_
    integer(c_size_t) :: sz_
    integer(c_int) :: to_node_
    integer(c_long) :: info_

    addr_ = c_loc(addr)
    sz_ = int(sz, c_size_t) * c_sizeof(addr)
    to_node_ = int(to_node, c_int)
    call move_data_pages(addr_, sz_, to_node_, info_)
    info = int(info_)
  end subroutine move_data_pages_i

  subroutine move_data_pages_l(addr, sz, to_node, info)
    use, intrinsic :: iso_c_binding
    implicit none
    logical, intent(in), target :: addr
    integer, intent(in) :: sz
    integer, intent(in) :: to_node
    integer, intent(out) :: info

    type(c_ptr) :: addr_
    integer(c_size_t) :: sz_
    integer(c_int) :: to_node_
    integer(c_long) :: info_

    addr_ = c_loc(addr)
    sz_ = int(sz, c_size_t) * c_sizeof(addr)
    to_node_ = int(to_node, c_int)
    call move_data_pages(addr_, sz_, to_node_, info_)
    info = int(info_)
  end subroutine move_data_pages_l

  subroutine move_data_pages_s(addr, sz, to_node, info)
    use, intrinsic :: iso_c_binding
    implicit none
    real(c_float), intent(in), target :: addr
    integer, intent(in) :: sz
    integer, intent(in) :: to_node
    integer, intent(out) :: info

    type(c_ptr) :: addr_
    integer(c_size_t) :: sz_
    integer(c_int) :: to_node_
    integer(c_long) :: info_

    addr_ = c_loc(addr)
    sz_ = int(sz, c_size_t) * c_sizeof(addr)
    to_node_ = int(to_node, c_int)
    call move_data_pages(addr_, sz_, to_node_, info_)
    info = int(info_)
  end subroutine move_data_pages_s

  subroutine move_data_pages_d(addr, sz, to_node, info)
    use, intrinsic :: iso_c_binding
    implicit none
    real(c_double), intent(in), target :: addr
    integer, intent(in) :: sz
    integer, intent(in) :: to_node
    integer, intent(out) :: info

    type(c_ptr) :: addr_
    integer(c_size_t) :: sz_
    integer(c_int) :: to_node_
    integer(c_long) :: info_

    addr_ = c_loc(addr)
    sz_ = int(sz, c_size_t) * c_sizeof(addr)
    to_node_ = int(to_node, c_int)
    call move_data_pages(addr_, sz_, to_node_, info_)
    info = int(info_)
  end subroutine move_data_pages_d

  subroutine move_data_pages_c(addr, sz, to_node, info)
    use, intrinsic :: iso_c_binding
    implicit none
    complex, intent(in), target :: addr
    integer, intent(in) :: sz
    integer, intent(in) :: to_node
    integer, intent(out) :: info

    type(c_ptr) :: addr_
    integer(c_size_t) :: sz_
    integer(c_int) :: to_node_
    integer(c_long) :: info_

    addr_ = c_loc(addr)
    sz_ = int(sz, c_size_t) * c_sizeof(addr)
    to_node_ = int(to_node, c_int)
    call move_data_pages(addr_, sz_, to_node_, info_)
    info = int(info_)
  end subroutine move_data_pages_c

  subroutine move_data_pages_z(addr, sz, to_node, info)
    use, intrinsic :: iso_c_binding
    implicit none
    complex(c_double), intent(in), target :: addr
    integer, intent(in) :: sz
    integer, intent(in) :: to_node
    integer, intent(out) :: info

    type(c_ptr) :: addr_
    integer(c_size_t) :: sz_
    integer(c_int) :: to_node_
    integer(c_long) :: info_

    addr_ = c_loc(addr)
    sz_ = int(sz, c_size_t) * c_sizeof(addr)
    to_node_ = int(to_node, c_int)
    call move_data_pages(addr_, sz_, to_node_, info_)
    info = int(info_)
  end subroutine move_data_pages_z

end module numa_f


text:
      SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
     $                   LDZ, IFST, ILST, WORK, LWORK, INFO )
*
*  -- LAPACK routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*     .. Scalar Arguments ..
      LOGICAL            WANTQ, WANTZ
      INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
     $                   WORK( * ), Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  DTGEXC reorders the generalized real Schur decomposition of a real
*  matrix pair (A,B) using an orthogonal equivalence transformation
*
*                 (A, B) = Q * (A, B) * Z',
*
*  so that the diagonal block of (A, B) with row index IFST is moved
*  to row ILST.
*
*  (A, B) must be in generalized real Schur canonical form (as returned
*  by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
*  diagonal blocks. B is upper triangular.
*
*  Optionally, the matrices Q and Z of generalized Schur vectors are
*  updated.
*
*         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
*         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
*
*
*  Arguments
*  =========
*
*  WANTQ   (input) LOGICAL
*          .TRUE. : update the left transformation matrix Q;
*          .FALSE.: do not update Q.
*
*  WANTZ   (input) LOGICAL
*          .TRUE. : update the right transformation matrix Z;
*          .FALSE.: do not update Z.
*
*  N       (input) INTEGER
*          The order of the matrices A and B. N >= 0.
*
*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
*          On entry, the matrix A in generalized real Schur canonical
*          form.
*          On exit, the updated matrix A, again in generalized
*          real Schur canonical form.
*
*  LDA     (input)  INTEGER
*          The leading dimension of the array A. LDA >= max(1,N).
*
*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
*          On entry, the matrix B in generalized real Schur canonical
*          form (A,B).
*          On exit, the updated matrix B, again in generalized
*          real Schur canonical form (A,B).
*
*  LDB     (input)  INTEGER
*          The leading dimension of the array B. LDB >= max(1,N).
*
*  Q       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
*          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
*          On exit, the updated matrix Q.
*          If WANTQ = .FALSE., Q is not referenced.
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q. LDQ >= 1.
*          If WANTQ = .TRUE., LDQ >= N.
*
*  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
*          On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
*          On exit, the updated matrix Z.
*          If WANTZ = .FALSE., Z is not referenced.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z. LDZ >= 1.
*          If WANTZ = .TRUE., LDZ >= N.
*
*  IFST    (input/output) INTEGER
*  ILST    (input/output) INTEGER
*          Specify the reordering of the diagonal blocks of (A, B).
*          The block with row index IFST is moved to row ILST, by a
*          sequence of swapping between adjacent blocks.
*          On exit, if IFST pointed on entry to the second row of
*          a 2-by-2 block, it is changed to point to the first row;
*          ILST always points to the first row of the block in its
*          final position (which may differ from its input value by
*          +1 or -1). 1 <= IFST, ILST <= N.
*
*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK. LWORK >= 4*N + 16.
*
*          If LWORK = -1, then a workspace query is assumed; the routine
*          only calculates the optimal size of the WORK array, returns
*          this value as the first entry of the WORK array, and no error
*          message related to LWORK is issued by XERBLA.
*
*  INFO    (output) INTEGER
*           =0:  successful exit.
*           <0:  if INFO = -i, the i-th argument had an illegal value.
*           =1:  The transformed matrix pair (A, B) would be too far
*                from generalized Schur form; the problem is ill-
*                conditioned. (A, B) may have been partially reordered,
*                and ILST points to the first row of the current
*                position of the block being moved.
*
*  Further Details
*  ===============
*
*  Based on contributions by
*     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*     Umea University, S-901 87 Umea, Sweden.
*
*  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
*      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
*      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
*      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY
      INTEGER            HERE, LWMIN, NBF, NBL, NBNEXT
*     ..
*     .. External Subroutines ..
      EXTERNAL           DTGEX2, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Executable Statements ..
*
*     Decode and test input arguments.
*
      INFO = 0
      LWMIN = MAX( 1, 4*N+16 )
      LQUERY = ( LWORK.EQ.-1 )
      IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
         INFO = -9
      ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
         INFO = -11
      ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
         INFO = -12
      ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
         INFO = -13
      ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -15
      END IF
*
      IF( INFO.EQ.0 ) THEN
         WORK( 1 ) = LWMIN
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DTGEXC', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( N.LE.1 )
     $   RETURN
*
*     Determine the first row of the specified block and find out
*     if it is 1-by-1 or 2-by-2.
*
      IF( IFST.GT.1 ) THEN
         IF( A( IFST, IFST-1 ).NE.ZERO )
     $      IFST = IFST - 1
      END IF
      NBF = 1
      IF( IFST.LT.N ) THEN
         IF( A( IFST+1, IFST ).NE.ZERO )
     $      NBF = 2
      END IF
*
*     Determine the first row of the final block
*     and find out if it is 1-by-1 or 2-by-2.
*
      IF( ILST.GT.1 ) THEN
         IF( A( ILST, ILST-1 ).NE.ZERO )
     $      ILST = ILST - 1
      END IF
      NBL = 1
      IF( ILST.LT.N ) THEN
         IF( A( ILST+1, ILST ).NE.ZERO )
     $      NBL = 2
      END IF
      IF( IFST.EQ.ILST )
     $   RETURN
*
      IF( IFST.LT.ILST ) THEN
*
*        Update ILST.
*
         IF( NBF.EQ.2 .AND. NBL.EQ.1 )
     $      ILST = ILST - 1
         IF( NBF.EQ.1 .AND. NBL.EQ.2 )
     $      ILST = ILST + 1
*
         HERE = IFST
*
   10    CONTINUE
*
*        Swap with next one below.
*
         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
*
*           Current block either 1-by-1 or 2-by-2.
*
            NBNEXT = 1
            IF( HERE+NBF+1.LE.N ) THEN
               IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
     $                   LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               RETURN
            END IF
            HERE = HERE + NBNEXT
*
*           Test if 2-by-2 block breaks into two 1-by-1 blocks.
*
            IF( NBF.EQ.2 ) THEN
               IF( A( HERE+1, HERE ).EQ.ZERO )
     $            NBF = 3
            END IF
*
         ELSE
*
*           Current block consists of two 1-by-1 blocks, each of which
*           must be swapped individually.
*
            NBNEXT = 1
            IF( HERE+3.LE.N ) THEN
               IF( A( HERE+3, HERE+2 ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
     $                   LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               RETURN
            END IF
            IF( NBNEXT.EQ.1 ) THEN
*
*              Swap two 1-by-1 blocks.
*
               CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
     $                      LDZ, HERE, 1, 1, WORK, LWORK, INFO )
               IF( INFO.NE.0 ) THEN
                  ILST = HERE
                  RETURN
               END IF
               HERE = HERE + 1
*
            ELSE
*
*              Recompute NBNEXT in case of 2-by-2 split.
*
               IF( A( HERE+2, HERE+1 ).EQ.ZERO )
     $            NBNEXT = 1
               IF( NBNEXT.EQ.2 ) THEN
*
*                 2-by-2 block did not split.
*
                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
     $                         Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
     $                         INFO )
                  IF( INFO.NE.0 ) THEN
                     ILST = HERE
                     RETURN
                  END IF
                  HERE = HERE + 2
               ELSE
*
*                 2-by-2 block did split.
*
                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
     $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
                  IF( INFO.NE.0 ) THEN
                     ILST = HERE
                     RETURN
                  END IF
                  HERE = HERE + 1
                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
     $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
                  IF( INFO.NE.0 ) THEN
                     ILST = HERE
                     RETURN
                  END IF
                  HERE = HERE + 1
               END IF
*
            END IF
         END IF
         IF( HERE.LT.ILST )
     $      GO TO 10
      ELSE
         HERE = IFST
*
   20    CONTINUE
*
*        Swap with next one below.
*
         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
*
*           Current block either 1-by-1 or 2-by-2.
*
            NBNEXT = 1
            IF( HERE.GE.3 ) THEN
               IF( A( HERE-1, HERE-2 ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
     $                   LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
     $                   INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               RETURN
            END IF
            HERE = HERE - NBNEXT
*
*           Test if 2-by-2 block breaks into two 1-by-1 blocks.
*
            IF( NBF.EQ.2 ) THEN
               IF( A( HERE+1, HERE ).EQ.ZERO )
     $            NBF = 3
            END IF
*
         ELSE
*
*           Current block consists of two 1-by-1 blocks, each of which
*           must be swapped individually.
*
            NBNEXT = 1
            IF( HERE.GE.3 ) THEN
               IF( A( HERE-1, HERE-2 ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
     $                   LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
     $                   INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               RETURN
            END IF
            IF( NBNEXT.EQ.1 ) THEN
*
*              Swap two 1-by-1 blocks.
*
               CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
     $                      LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
               IF( INFO.NE.0 ) THEN
                  ILST = HERE
                  RETURN
               END IF
               HERE = HERE - 1
            ELSE
*
*             Recompute NBNEXT in case of 2-by-2 split.
*
               IF( A( HERE, HERE-1 ).EQ.ZERO )
     $            NBNEXT = 1
               IF( NBNEXT.EQ.2 ) THEN
*
*                 2-by-2 block did not split.
*
                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
     $                         Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
                  IF( INFO.NE.0 ) THEN
                     ILST = HERE
                     RETURN
                  END IF
                  HERE = HERE - 2
               ELSE
*
*                 2-by-2 block did split.
*
                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
     $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
                  IF( INFO.NE.0 ) THEN
                     ILST = HERE
                     RETURN
                  END IF
                  HERE = HERE - 1
                  CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
     $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
                  IF( INFO.NE.0 ) THEN
                     ILST = HERE
                     RETURN
                  END IF
                  HERE = HERE - 1
               END IF
            END IF
         END IF
         IF( HERE.GT.ILST )
     $      GO TO 20
      END IF
      ILST = HERE
      WORK( 1 ) = LWMIN
      RETURN
*
*     End of DTGEXC
*
      END


text:
      SUBROUTINE FRLG
C
C     FREQUENCE RESPONSE LOAD GENERATOR
C
C     INPUTS  - CASEXX,USETD,DLT,FRL,GMD,GOD,DIT,PHIDH
C
C     OUTPUTS - PPF,PSF,PDF,FOL,PHF
C
C     4 SCRATHCES
C
C
      EXTERNAL        ANDF
      INTEGER         CASEXX,USETD,DLT,FRL,GMD,GOD,DIT,PHIDH,PPF,PSF,
     1                PDF,FOL,PHF,SCR1,SCR2,SCR3,SCR4,MCB(7),ANDF,
     2                SINGLE,OMIT,IFREQ(2),ITRAN(2)
      COMMON /TWO   / ITWO(32)
      COMMON /BLANK / MODAL(2),NOTRD,IAPP(2)
      COMMON /BITPOS/ IUM,IUO,SKP(6),IUS
      DATA    CASEXX, USETD,DLT,FRL,GMD,GOD,DIT,PHIDH      /
     1        101   , 102  ,103,104,105,106,107,108        /
      DATA    PPF   , PSF,PDF,FOL,PHF, SCR1,SCR2,SCR3,SCR4 /
     1        201   , 202,203,204,205, 301 ,302 ,303 ,304  /
      DATA    MODA  / 4HMODA  /
      DATA    IFREQ , ITRAN   /4HFREQ,1H ,4HTRAN,1H /
C
C     DETERMINE USET DATA
C
      MCB(1) = USETD
      CALL RDTRL (MCB)
      LUSETD = MCB(2)
      MULTI  =-1
      IF (ANDF(MCB(5),ITWO(IUM)) .NE. 0) MULTI =  1
      SINGLE =-1
      IF (ANDF(MCB(5),ITWO(IUS)) .NE. 0) SINGLE =  1
      OMIT =-1
      IF (ANDF(MCB(5),ITWO(IUO)) .NE. 0) OMIT =  1
C
      IAPP(1) = IFREQ(1)
      IAPP(2) = IFREQ(2)
C
C     BUILD LOADS ON P SET
C
C     ORDER IS ALL FREQUENCIES FOR GIVEN LOAD TOGETHER
C
      CALL FRLGA (DLT,FRL,CASEXX,DIT,PPF,LUSETD,NFREQ,NLOAD,FRQSET,FOL,
     1            NOTRD)
      IF (NOTRD .EQ. -1) GO TO 10
      IAPP(1) = ITRAN(1)
      IAPP(2) = ITRAN(2)
  10  CONTINUE
C
C     REDUCE LOADS TO D OR H SET
C
      IF (MULTI.LT.0 .AND. SINGLE.LT.0 .AND. OMIT.LT.0 .AND.
     1    MODAL(1).NE.MODA) RETURN
      CALL FRLGB (PPF,USETD,GMD,GOD,MULTI,SINGLE,OMIT,MODAL,PHIDH,PDF,
     1            PSF,PHF,SCR1,SCR2,SCR3,SCR4)
      RETURN
      END


#INCLUDE 'MR_H_ALIGN_PADDING.H'
!***********************************************************************************************************************************
! UNIT:
!
!  (MODULE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  MODULE MR_MOD_SCALE

    USE MR_KINDS

    USE MR_DEF_FIELD_VARS_DSET_NAMES

    IMPLICIT NONE

    PRIVATE

    PUBLIC :: MR_SCALE_UV
    PUBLIC :: MR_SCALE_SS

    PUBLIC :: MR_SCALE_XY

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

  CONTAINS

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_DETER_SCALE_DSET( PATH_DSET , SCALE , SCALE_DSET )

    CHARACTER(   *   ) , INTENT(IN   ) :: PATH_DSET

    REAL   (PARD_KIND) , INTENT(IN   ) :: SCALE
    REAL   (PARD_KIND) , INTENT(OUT  ) :: SCALE_DSET

    IF( PATH_DSET == "" ) THEN
      RETURN
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_TBFUV) ) /= 0 ) THEN
      SCALE_DSET =  SCALE
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_TBUV)  ) /= 0 ) THEN
      SCALE_DSET =  SCALE
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_CSS)   ) /= 0 ) THEN
      SCALE_DSET =  1.0
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_QSBUV) ) /= 0 ) THEN
      SCALE_DSET =  SCALE**1.5
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_KI)    ) /= 0 ) THEN
      SCALE_DSET =  SCALE
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_DI)    ) /= 0 ) THEN
      SCALE_DSET =  SQRT( SCALE )
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_ZB)    ) /= 0 ) THEN
      SCALE_DSET =  SCALE
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_ZS)    ) /= 0 ) THEN
      SCALE_DSET =  SCALE
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_UVA)   ) /= 0 ) THEN
      SCALE_DSET =  SQRT( SCALE )
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_UV)    ) /= 0 ) THEN
      SCALE_DSET =  SQRT( SCALE )
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_WW)    ) /= 0 ) THEN
      SCALE_DSET =  SQRT( SCALE )
  ! THE DETECTION OF H MUST BE AFTER THE DETECTION OF UVA
  ! FOR THE NAME OF H IS CONTAINED IN THE NAME OF UVA
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_H)     ) /= 0 ) THEN
      SCALE_DSET =  SCALE
  ! SEE MR_DEF_FIELD_VARS_DSET_NAMES
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_R)     ) /= 0 ) THEN
      SCALE_DSET =  1.0
  ! THE DETECTION OF VZWW MUST BE AFTER THE DETECTION OF KI AND/OR DI
  ! FOR THE NAME OF VZWW IS CONTAINED IN THE PATH OF KI AND/OR DI
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_VZWW)  ) /= 0 ) THEN
      SCALE_DSET =  SCALE**1.5
  ! SEE MR_MOD_INIT_OUTPUT
    ELSE IF( INDEX( TRIM(PATH_DSET) , "/"//TRIM(DSET_NAME_DZWW)  ) /= 0 ) THEN
      SCALE_DSET =  SCALE**1.5
    END IF

  END SUBROUTINE MR_DETER_SCALE_DSET

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_UV( PATH_UV , SCALE , NI , NJ , UV , UU , VV , UVO )

    IMPLICIT NONE

    CHARACTER(   *   ) , INTENT(IN   ) :: PATH_UV

    REAL   (PARD_KIND) , INTENT(IN   ) :: SCALE

    INTEGER(IJID_KIND) , INTENT(IN   ) :: NI , NJ

    REAL   (FDRD_KIND) , INTENT(INOUT) , DIMENSION(1:NI1(NI,FDRD_KIND),1:NJ,1:2) :: UV
    REAL   (FDRD_KIND) , INTENT(INOUT) , DIMENSION(0:NI0(NI,FDRD_KIND),1:NJ,1:2) :: UU
    REAL   (FDRD_KIND) , INTENT(INOUT) , DIMENSION(1:NI1(NI,FDRD_KIND),0:NJ,1:2) :: VV
    REAL   (FDRD_KIND) , INTENT(INOUT) , DIMENSION(0:NI0(NI,FDRD_KIND),0:NJ,1:2) :: UVO

    REAL   (PARD_KIND) :: SCALE_UV

    CALL MR_DETER_SCALE_DSET( PATH_UV , SCALE , SCALE_UV )

   !DIR$ FORCEINLINE
    CALL MR_SCALE_UV_FOR_W_NODES
   !DIR$ FORCEINLINE
    CALL MR_SCALE_UV_FOR_U_NODES
   !DIR$ FORCEINLINE
    CALL MR_SCALE_UV_FOR_V_NODES
   !DIR$ FORCEINLINE
    CALL MR_SCALE_UV_FOR_O_NODES
   !END$ FORCEINLINE

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

  CONTAINS

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_UV_FOR_W_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J
    INTEGER            :: DIM

    DO DIM = 1 , 2

      DO J = 1 , NJ
       !DIR$ VECTOR ALIGNED
        DO I = 1 , NI
          UV( I , J ,DIM) = UV( I , J ,DIM) * SCALE_UV
        END DO
      END DO

    END DO

  END SUBROUTINE MR_SCALE_UV_FOR_W_NODES

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_UV_FOR_U_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J
    INTEGER            :: DIM

    DO DIM = 1 , 2

      DO J = 1 , NJ
       !DIR$ VECTOR ALIGNED
        DO I = 0 , NI
          UU( I , J ,DIM) = UU( I , J ,DIM) * SCALE_UV
        END DO
      END DO

    END DO

  END SUBROUTINE MR_SCALE_UV_FOR_U_NODES

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_UV_FOR_V_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J
    INTEGER            :: DIM

    DO DIM = 1 , 2

      DO J = 0 , NJ
       !DIR$ VECTOR ALIGNED
        DO I = 1 , NI
          VV( I , J ,DIM) = VV( I , J ,DIM) * SCALE_UV
        END DO
      END DO

    END DO

  END SUBROUTINE MR_SCALE_UV_FOR_V_NODES

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_UV_FOR_O_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J
    INTEGER            :: DIM

    DO DIM = 1 , 2

      DO J = 0 , NJ
       !DIR$ VECTOR ALIGNED
        DO I = 0 , NI
          UVO( I , J ,DIM) = UVO( I , J ,DIM) * SCALE_UV
        END DO
      END DO

    END DO

  END SUBROUTINE MR_SCALE_UV_FOR_O_NODES

  END SUBROUTINE MR_SCALE_UV

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_SS( PATH_SS , SCALE , NI , NJ , SS , SU , SV , SO )

    IMPLICIT NONE

    CHARACTER(   *   ) , INTENT(IN   ) :: PATH_SS

    REAL   (PARD_KIND) , INTENT(IN   ) :: SCALE

    INTEGER(IJID_KIND) , INTENT(IN   ) :: NI , NJ

    REAL   (FDRD_KIND) , INTENT(INOUT) , DIMENSION(1:NI1(NI,FDRD_KIND),1:NJ) :: SS
    REAL   (FDRD_KIND) , INTENT(INOUT) , DIMENSION(0:NI0(NI,FDRD_KIND),1:NJ) :: SU
    REAL   (FDRD_KIND) , INTENT(INOUT) , DIMENSION(1:NI1(NI,FDRD_KIND),0:NJ) :: SV
    REAL   (FDRD_KIND) , INTENT(INOUT) , DIMENSION(0:NI0(NI,FDRD_KIND),0:NJ) :: SO

    REAL   (PARD_KIND) :: SCALE_SS

    CALL MR_DETER_SCALE_DSET( PATH_SS , SCALE , SCALE_SS )

   !DIR$ FORCEINLINE
    CALL MR_SCALE_SS_FOR_W_NODES
   !DIR$ FORCEINLINE
    CALL MR_SCALE_SS_FOR_U_NODES
   !DIR$ FORCEINLINE
    CALL MR_SCALE_SS_FOR_V_NODES
   !DIR$ FORCEINLINE
    CALL MR_SCALE_SS_FOR_O_NODES
   !END$ FORCEINLINE

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

  CONTAINS

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_SS_FOR_W_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J

    DO J = 1 , NJ
     !DIR$ VECTOR ALIGNED
      DO I = 1 , NI
        SS( I , J ) = SS( I , J ) * SCALE_SS
      END DO
    END DO

  END SUBROUTINE MR_SCALE_SS_FOR_W_NODES

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_SS_FOR_U_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J

    DO J = 1 , NJ
     !DIR$ VECTOR ALIGNED
      DO I = 0 , NI
        SU( I , J ) = SU( I , J ) * SCALE_SS
      END DO
    END DO

  END SUBROUTINE MR_SCALE_SS_FOR_U_NODES

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_SS_FOR_V_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J

    DO J = 0 , NJ
     !DIR$ VECTOR ALIGNED
      DO I = 1 , NI
        SV( I , J ) = SV( I , J ) * SCALE_SS
      END DO
    END DO

  END SUBROUTINE MR_SCALE_SS_FOR_V_NODES

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_SS_FOR_O_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J

    DO J = 0 , NJ
     !DIR$ VECTOR ALIGNED
      DO I = 0 , NI
        SO( I , J ) = SO( I , J ) * SCALE_SS
      END DO
    END DO

  END SUBROUTINE MR_SCALE_SS_FOR_O_NODES

  END SUBROUTINE MR_SCALE_SS

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_XY( SCALE , NI , NJ , XYUV , XYUU , XYVV , XYOO )

    IMPLICIT NONE

    REAL   (PARD_KIND) , INTENT(IN   ) :: SCALE

    INTEGER(IJID_KIND) , INTENT(IN   ) :: NI , NJ

    REAL   (XYRD_KIND) , INTENT(INOUT) , DIMENSION(1:NI1(NI,XYRD_KIND),1:NJ,1:2) :: XYUV
    REAL   (XYRD_KIND) , INTENT(INOUT) , DIMENSION(0:NI0(NI,XYRD_KIND),1:NJ,1:2) :: XYUU
    REAL   (XYRD_KIND) , INTENT(INOUT) , DIMENSION(1:NI1(NI,XYRD_KIND),0:NJ,1:2) :: XYVV
    REAL   (XYRD_KIND) , INTENT(INOUT) , DIMENSION(0:NI0(NI,XYRD_KIND),0:NJ,1:2) :: XYOO

   !DIR$ FORCEINLINE
    CALL MR_SCALE_XY_FOR_W_NODES
   !DIR$ FORCEINLINE
    CALL MR_SCALE_XY_FOR_U_NODES
   !DIR$ FORCEINLINE
    CALL MR_SCALE_XY_FOR_V_NODES
   !DIR$ FORCEINLINE
    CALL MR_SCALE_XY_FOR_O_NODES
   !END$ FORCEINLINE

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

  CONTAINS

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_XY_FOR_W_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J
    INTEGER            :: DIM

    DO DIM = 1 , 2

      DO J = 1 , NJ
       !DIR$ VECTOR ALIGNED
        DO I = 1 , NI
          XYUV( I , J ,DIM) = XYUV( I , J ,DIM) * SCALE
        END DO
      END DO

    END DO

  END SUBROUTINE MR_SCALE_XY_FOR_W_NODES

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_XY_FOR_U_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J
    INTEGER            :: DIM

    DO DIM = 1 , 2

      DO J = 1 , NJ
       !DIR$ VECTOR ALIGNED
        DO I = 0 , NI
          XYUU( I , J ,DIM) = XYUU( I , J ,DIM) * SCALE
        END DO
      END DO

    END DO

  END SUBROUTINE MR_SCALE_XY_FOR_U_NODES

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_XY_FOR_V_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J
    INTEGER            :: DIM

    DO DIM = 1 , 2

      DO J = 0 , NJ
       !DIR$ VECTOR ALIGNED
        DO I = 1 , NI
          XYVV( I , J ,DIM) = XYVV( I , J ,DIM) * SCALE
        END DO
      END DO

    END DO

  END SUBROUTINE MR_SCALE_XY_FOR_V_NODES

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_SCALE_XY_FOR_O_NODES

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J
    INTEGER            :: DIM

    DO DIM = 1 , 2

      DO J = 0 , NJ
       !DIR$ VECTOR ALIGNED
        DO I = 0 , NI
          XYOO( I , J ,DIM) = XYOO( I , J ,DIM) * SCALE
        END DO
      END DO

    END DO

  END SUBROUTINE MR_SCALE_XY_FOR_O_NODES

  END SUBROUTINE MR_SCALE_XY

  END MODULE MR_MOD_SCALE


text:
subroutine gaussb(a, b, deter)
  implicit none

  real, dimension(:,:), intent(inout) :: a
  real, dimension(:), intent(inout) :: b
  real, intent(inout) :: deter

  real :: factor, piv

  integer :: n, k, i, j

  n = size(b)
  deter = 1

  do k=1,n-1
    piv=a(k,k)
    if(abs(piv) < 1.e-12) then
      print*, 'Pivote nulo na etapa: ', k
      stop
    end if

    deter = deter*piv

    a(k+1:n,k) = a(k+1:n,k)/piv
    do j=k+1,n
      a(k+1:n,j) = a(k+1:n,j)-a(k+1:n,k)*a(k,j)
    end do
    b(k+1:n) = b(k+1:n)-a(k+1:n,k)*b(k)

  end do

  if(abs(a(n,n)) < 1.e-12) then
    print*, 'Pivote nulo na etapa: ', n
    stop
  end if

  deter = deter*a(n,n)
end subroutine

max_stars_repo_name: UxioAndrade/Analisis-Numerico-Matricial

max_issues_repo_name: UxioAndrade/Analisis-Numerico-Matricial

max_forks_repo_name: UxioAndrade/Analisis-Numerico-Matricial


text:
Describe Users/ThrowAwayAccount100 here.
Welcome to the Wiki!  Why the anonymous account?  Folks here tend to be reasonably friendly, and theres a lot of value in establishing an Identity.  At a minimum, it gives more weight to your words, since the first thing youre telling readers with an anonymous account is Im not willing to stand behind what Im saying. Users/TomGarberson

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search