text:
!! R733 pointer-assignment-stmt
!   - F2008 addition
!   - a pointer function reference can denote a variable in any variable
!     definition context
!
subroutine pointer_function_ref(proc_ptr)
  INTERFACE
     SUBROUTINE int_arg_sub(i)
       integer :: i
     END SUBROUTINE
  END INTERFACE

  procedure(int_arg_sub), pointer :: proc_ptr

  proc_ptr => internal_proc

contains

  subroutine internal_proc(i)
    integer :: i
    i = 13
  end subroutine

end subroutine

max_stars_repo_name: ouankou/rose

max_issues_repo_name: WildeGeist/rose

max_forks_repo_name: WildeGeist/rose


text:
      SUBROUTINE SXLOOP (X,Y,N)
C
C     SINGLE PRECISION VERSION OF CXLOOP, CALLED ONLY BY CDCMPS
C
      REAL       X(1),   Y(1)
      REAL       XX(2) , YY(2), MPY(2)
      NN = N + N
      DO 10 I = 1,NN
   10 X(I) = Y(I)
      RETURN
C
      ENTRY SLOOP (XX,YY,MPY,M)
      MM = M + M
      DO 20 I = 1,MM,2
      XX(I  ) = XX(I  ) - MPY(1)*YY(I) + MPY(2)*YY(I+1)
   20 XX(I+1) = XX(I+1) - MPY(2)*YY(I) - MPY(1)*YY(I+1)
      RETURN
      END

ext: f


text: Anns Nursery School is a Schools Preschool for children aged three to five.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
submodule (io) mag

use sanity_check, only : check_finite_mag

implicit none (type, external)

interface ! mag_*.f90

module subroutine output_magfields_hdf5(outdir,ymd,UTsec,Br,Btheta,Bphi)
character(*), intent(in) :: outdir
integer, intent(in) :: ymd(3)
real(wp), intent(in) :: UTsec
real(wp), dimension(:), intent(in)  :: Br,Btheta,Bphi
end subroutine output_magfields_hdf5

module subroutine output_magfields_nc4(outdir,ymd,UTsec,Br,Btheta,Bphi)
character(*), intent(in) :: outdir
integer, intent(in) :: ymd(3)
real(wp), intent(in) :: UTsec
real(wp), dimension(:), intent(in)  :: Br,Btheta,Bphi
end subroutine output_magfields_nc4

module subroutine output_magfields_raw(outdir,ymd,UTsec,Br,Btheta,Bphi)
character(*), intent(in) :: outdir
integer, intent(in) :: ymd(3)
real(wp), intent(in) :: UTsec
real(wp), dimension(:), intent(in)  :: Br,Btheta,Bphi
end subroutine output_magfields_raw

end interface

contains

module procedure create_outdir_mag
! subroutine create_outdir_mag(outdir,fieldpointfile)
!! CREATES OUTPUT DIRECTORY FOR MAGNETIC FIELD CALCULATIONS

integer :: ierr

!NOTE HERE THAT WE INTERPRET OUTDIR AS THE BASE DIRECTORY CONTAINING SIMULATION OUTPUT
if ( mkdir(outdir//'/magfields/') /= 0 ) error stop 'could not create magfields output directory'
if ( mkdir(outdir//'/magfields/input/') /= 0 ) error stop 'could not create magfields/input directory'
if ( copyfile(fieldpointfile, outdir//'/magfields/input/magfieldpoints.dat') /=0 ) error stop 'could not create magfields dat file'

end procedure create_outdir_mag

module procedure output_magfields

select case (out_format)
case ('raw')
  call output_magfields_raw(outdir,ymd,UTsec,Br,Btheta,Bphi)
case ('h5', 'hdf5')
  call output_magfields_hdf5(outdir,ymd,UTsec,Br,Btheta,Bphi)
case ('nc', 'nc4')
  call output_magfields_nc4(outdir,ymd,UTsec,Br,Btheta,Bphi)
case default
  write(stderr,*) 'mag:output_magfields: unknown file format' // out_format
  error stop 6
end select

call check_finite_mag(Br, Btheta, Bphi)

end procedure output_magfields

end submodule mag

max_stars_repo_name: jklenzing/gemini

max_issues_repo_name: jklenzing/gemini

max_forks_repo_name: jklenzing/gemini


text:
      SUBROUTINE LIFRE(A,FA,AHEAD,ATAIL)
C     PREFIX THE FREE LIST WITH A LIST
C
      INTEGER*2 I
     = ,A(1)     ! THE ARRAY FOR THE POINTERS OF THE LIST
     = ,AHEAD    ! -> FIRST ELEMENT
     = ,AHIB     ! HBOUND(A,1)
     = ,ATAIL    ! -> LAST ELEMENT
     = ,FA       ! -> THE 1ST FREE ELEMENT
     = ,I1       ! -> THE ELEMENT BEFORE 'I2'
     = ,I2       ! -> CURRENT ELEMENT
C
      IF (AHEAD .EQ. ATAIL) GOTO 3
        A(1) = AHEAD
        I1 = 1
        I2 = A(I1)
11      IF(I2 .EQ. ATAIL) GOTO 2
            I1 = I2
            I2 = A(I1)
        GOTO 11
2       CONTINUE ! ALL ELEMENTS
        ! NOW PREFIX THE FREE LIST WITH 'A(AHEAD : I1)'
        A(I1) = FA
        FA = AHEAD
        AHEAD = ATAIL
3     CONTINUE ! AHEAD .EQ. ATAIL
      RETURN ! LIFRE
      END

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


subroutine new_vuf(vx,ux,fx,dd,xlat,constituent)
  implicit none
      ! Parameters
      integer, parameter   ::   maxfreq=170
      real,parameter       ::   pi=3.141592653589793
      complex,parameter    ::   i=(0.0,1.0)

      ! Outout variables
      real*8      ::   fx,vx,ux

      ! Input variables
      real*8      ::   xlat, dd
      character(len=5)     ::   constituent

      ! Variables for main constituents
      integer,dimension(maxfreq)            ::   ii=0,jj=0,kk=0,ll=0,mm=0,nn=0
      integer,dimension(maxfreq)            ::   ikmpr=0,nsat=0
      integer,dimension(maxfreq)            ::   nshallow=0,ishallow=0
      real,dimension(maxfreq)               ::   freq=0,df=0,semi=0
      character(len=5),dimension(maxfreq)   ::   name,kmpr,kon

      ! Variables for satellite data
      integer,dimension(180)                ::   ldel=0,mdel=0,ndel=0
      integer,dimension(180)                ::   ilatfac=0,iconst=0
      real,dimension(180)                   ::   amprat=0,phcorr=0

      ! Variables for shallow water constituents
      integer,dimension(320)                ::   shallow_iconst=0,shallow_iname=0
      real,dimension(320)                   ::   shallow_coef=0

      ! Temporary variables for file input
      integer                          ::   iii,jjj,kkk,lll,mmm,nnn,nnj
      real                             ::   nsemi
      real,dimension(4)                ::   ncoef
      character(len=5)                 ::   nkon
      character(len=5),dimension(4)    ::   nkonco

      ! Variables for calculation of output arguments
      real,dimension(146)              ::    f=0,u=0,v=0
      real,dimension(180)              ::    rr=0,uu=0
      complex,dimension(146)           ::    fsum=0
      ! Temporary variables used in the calculations
      real                             ::    dum,dumm,tmp2
      complex                          ::    tmp
      ! Loop variables
      integer   ::   t,j,k
      integer   ::   ntotal,ntidal,nsatellite,m
      integer   ::   jbase,j1,jend,jstart

      ! The astronomical arguments
      real,dimension(6)   ::   astro,ader

      ! Miscellaneous
      integer             ::   konx,kr,ierr,ier
      real                ::   slat
      character(len=5)    ::   kblank

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

      xlat=55
      kr=82
      DATA kblank/'     '/
      if (kr.gt.0) then
         open(unit=kr,file='constits.dat'                                &
              ,status='old'                                               &
              ,form='formatted')
      end if

      do k=1,maxfreq
         READ(kr,1010) name(k),freq(k),kmpr(k)
1010     FORMAT(4X,A5,3X,F13.10,4X,A5)

         if (name(k).eq.constituent) then
            konx=k
         elseif(name(k).eq.kblank) then
            m=k-1
            exit
         end if

      end do

      do k=1,m
         do t=1,m

            if (kmpr(k)==name(t))then
               ikmpr(k)=t
               df(k)=abs(freq(k)-freq(t))
               exit
            end if

         end do
      end do
      df(1)=0

      ierr = 0

!***********************************************************************
!  THE ASTRONOMICAL ARGUMENTS AND THEIR RATES OF CHANGE,
!  S0,H0,P0,ENP0,PP0,DS,DH,DP,DNP,DPP,  ARE READ FROM TWO RECORDS IN
!  THE FORMAT(5F13.10):
!
!	These values are no longer used though they are still
!	read in. More accurate polynomial approximations are
!	now employed.
      read(kr,50)dum,dum,dum,dum,dum,dum,dum,dum,dum,dum
50    FORMAT(5F13.10)

!***********************************************************************
!  HERE THE MAIN CONSTITUENTS AND THEIR DOODSON NUMBERS ARE READ IN
!  FORMAT (6X,A5,1X,6I3,F5.2,I4). THE VALUES ARE RESPECTIVELY
!     kon    = CONSTITUENT NAME
!  ii,jj,kk,ll,mm,nn = THE SIX DOODSON NUMBERS
!     semi   = PHASE CORRECTION
!     nj     = THE NUMBER OF SATELLITES FOR THIS CONSTITUENT.
!  THE END OF ALL MAIN CONSTITUENTS IS DENOTED BY A BLANK CARD.
!  There are 45 main constituents in the original package.   - bh
!
      nsatellite=0

      jbase=1
      jend=0

      do k=1,1000
         read(kr,60)nkon,iii,jjj,kkk,lll,mmm,nnn,nsemi,nnj
60       FORMAT(6X,A5,1X,6I3,F5.2,I4)

         if(nkon.eq.kblank)then
            ntidal=k-1
            exit
         end if

         do t=1,m
            if (nkon.eq.name(t)) then
               j1=t
               exit
            end if
         end do

         kon(j1)=nkon
         ii(j1)=iii
         jj(j1)=jjj
         kk(j1)=kkk
         ll(j1)=lll
         mm(j1)=mmm
         nn(j1)=nnn
         semi(j1)=nsemi
         nsat(j1)=nnj

         if(nsat(j1).ge.1)then
            nsatellite=nsatellite+nsat(j1)
            jend=jend+nsat(j1)
            !******************************************************************
            !  IF NJ>0, INFORMATION ON THE SATELLITE CONSTITUENTS IS READ ,
            !  THREE SATELLITES PER CARD, IN THE FORMAT
            ! (11X,3(3I3,F4.2,F7.4,1X,I1,1X)).
            !  FOR EACH SATELLITE THE VALUES READ ARE
            !  LDEL,MDEL,NDEL = THE CHANGES IN THE LAST THREE DOODSON NUMBERS
            !                   FROM THOSE OF THE MAIN CONSTITUENT.
            !     ph     = THE PHASE CORRECTION
            !     ee     = THE AMPLITUDE RATIO OF THE SATELLITE TIDAL POTENTIAL
            !              TO THAT OF THE MAIN CONSTITUENT.
            !     ir     = 1 IF THE AMPLITUDE RATIO HAS TO BE MULTIPLIED BY THE
            !              LATITUDE CORRECTION FACTOR FOR DIURNAL CONSTITUENTS
            !              2 IF THE AMPLITUDE RATIO HAS TO BE MULTIPLIED BY THE
            !              LATITUDE CORRECTION FACTOR FOR SEMI-DIURNAL CONSTI-
            !              TUENTS.
            !              OTHERWISE IF NO CORRECTION IS REQUIRED TO THE
            !              AMPLITUDE RATIO
            read(kr,80)(ldel(j),mdel(j),ndel(j),phcorr(j),amprat(j),ilatfac(j),j=jbase,jend)
80          FORMAT((11X,3(3I3,F4.2,F7.4,1X,I1,1X)))
            do t=jbase,jend
               iconst(t)=j1
            end do
            jbase=jend+1
         end if
      end do

!
!***********************************************************************
!  THE SHALLOW WATER CONSTITUENTS AND THE MAIN CONSTITUENTS FROM WHICH
!  THEY ARE DERIVED ARE READ IN HERE WITH THE FORMAT
!  (6X,A5,I1,2X,4(F5.2,A5,5X)). THE VALUES ARE RESPECTIVELY
!     KON    = NAME  OF THE SHALLOW WATER CONSTITUENT
!     NJ     = NUMBER OF MAIN CONSTITUENTS FROM WHICH IT IS DERIVED.
!     COEF,KONCO = COMBINATION NUMBER AND NAME OF THESE MAIN
!                  CONSTITUENTS.
!  THE END OF THESE CONSTITUENTS IS DENOTED BY A BLANK CARD.
!   There are 101 shallow water constituents in the original package.
!   Together with the main constituents we get 146 altogether.
!
      jbase=0
      do k=1,1000
         read(kr,130)nkon,nnj,(ncoef(j),nkonco(j),j=1,4)
130      FORMAT(6X,A5,I1,2X,4(F5.2,A5,5X))

         if(nkon.eq.kblank) then
            ntotal=k-1
           exit

        end if

        do t=1,m
           if (nkon.eq.name(t)) then
              j1=t
              exit
           end if
        end do

        nshallow(j1)=nnj
        ishallow(j1)=jbase+1

        do t=1,nnj
           jbase=jbase+1
           shallow_coef(jbase)=ncoef(t)
           shallow_iconst(jbase)=j1

           do kkk=1,m

              if (nkonco(t).eq.name(kkk)) then
                 shallow_iname(jbase)=kkk
                 exit
              end if

           end do

        end do

     end do
     close(unit=KR)

! Calculate astronomical arguments at mid-point of data time series

     call new_astro(dd,astro,ader)

! Phase relative to Greenwich (in units of cycles, presumeably).
! (This only returns values when we have doodson#s, i.e., not for the
! shallow water components, but these will be computed later.)
!v=rem( const.doodson*astro+const.semi, 1);

     do k=1,m
        dum=ii(k)*astro(1) + jj(k)*astro(2) + kk(k)*astro(3) +         &
             ll(k)*astro(4) + mm(k)*astro(5) + nn(k)*astro(6) +         &
             semi(k)
        v(k)=dum-floor(dum)
     end do

!     call new_astro(dd,astro,ader)
     ! Apparently the second-order terms in the tidal potential go to zero
     ! at the equator, but the third-order terms do not. Hence when trying
     ! to infer the third-order terms from the second-order terms, the
     ! nodal correction factors blow up. In order to prevent this, it is
     ! assumed that the equatorial forcing is due to second-order forcing
     ! OFF the equator, from about the 5 degree location. Latitudes are
     ! hence (somewhat arbitrarily) forced to be no closer than 5 deg to
     ! the equator.

     if (abs(xlat).lt.5.0) then
        if (xlat<0) then
           xlat=-5.0
        else
           xlat=5.0
        end if
     end if

     slat=sin(pi*xlat/180.0)

     ! satellite amplitude ratio adjustment for latitude

     rr=amprat  ! no amplitude correction

     do k=1,nsatellite
        if (ilatfac(k).eq.1) then
           ! latitude correction for diurnal constituents
           rr(k)=rr(k)*0.36309*(1.0-5.0*slat*slat)/slat
        elseif (ilatfac(k).eq.2) then
           ! latitude correction for semi-diurnal constituents
           rr(k)=rr(k)*2.59808*slat
        end if
     end do

     ! calculate nodal amplitude and phase corrections
     do k=1,nsatellite
        dum=ldel(k)*astro(4)+mdel(k)*astro(5)+ndel(k)*astro(6)+phcorr(k)
        uu(k)=dum-floor(dum)
     end do

     ! Sum up all of the satellite factors for all satellites
     fsum = 0  ! emsures that all elements of fsum equals zero before we start
     do k=1,nsatellite
        fsum(iconst(k))=fsum(iconst(k))+rr(k)*exp(i*2.0*pi*uu(k))
     end do

     fsum=fsum+(1,0)

     f=abs(fsum)

     do k=1,m
        tmp2=aimag(fsum(k))
        tmp=tmp2*i
        tmp=fsum(k)-tmp
        u(k)=atan2(tmp2,abs(tmp))/(2*pi)
     end do

     ! Compute amplitude and phase corrections for shallow water constituents
     do k=1,m

        if(nshallow(k).gt.0) then
           jbase=ishallow(k)
           jend=jbase+nshallow(k)-1
           tmp2=1
           dum=0
           dumm=0
           do t=jbase,jend
              tmp2=f(shallow_iname(t))**abs(shallow_coef(t))*tmp2
              dum=u(shallow_iname(t))*shallow_coef(t)+dum
              dumm=v(shallow_iname(t))*shallow_coef(t)+dumm
           end do
           f(k)=tmp2
           u(k)=dum
           v(k)=dumm
        end if

     end do
     fx=f(konx)
     ux=u(konx)
     vx=v(konx)

   end subroutine new_vuf


text:
      SUBROUTINE realft(data,n,isign)
      INTEGER isign,n
      REAL data(n)
CU    USES four1
      INTEGER i,i1,i2,i3,i4,n2p3
      REAL c1,c2,h1i,h1r,h2i,h2r,wis,wrs
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      theta=3.141592653589793d0/dble(n/2)
      c1=0.5
      if (isign.eq.1) then
        c2=-0.5
        call four1(data,n/2,+1)
      else
        c2=0.5
        theta=-theta
      endif
      wpr=-2.0d0*sin(0.5d0*theta)**2
      wpi=sin(theta)
      wr=1.0d0+wpr
      wi=wpi
      n2p3=n+3
      do 11 i=2,n/4
        i1=2*i-1
        i2=i1+1
        i3=n2p3-i2
        i4=i3+1
        wrs=sngl(wr)
        wis=sngl(wi)
        h1r=c1*(data(i1)+data(i3))
        h1i=c1*(data(i2)-data(i4))
        h2r=-c2*(data(i2)+data(i4))
        h2i=c2*(data(i1)-data(i3))
        data(i1)=h1r+wrs*h2r-wis*h2i
        data(i2)=h1i+wrs*h2i+wis*h2r
        data(i3)=h1r-wrs*h2r+wis*h2i
        data(i4)=-h1i+wrs*h2i+wis*h2r
        wtemp=wr
        wr=wr*wpr-wi*wpi+wr
        wi=wi*wpr+wtemp*wpi+wi
11    continue
      if (isign.eq.1) then
        h1r=data(1)
        data(1)=h1r+data(2)
        data(2)=h1r-data(2)
      else
        h1r=data(1)
        data(1)=c1*(h1r+data(2))
        data(2)=c1*(h1r-data(2))
        call four1(data,n/2,-1)
      endif
      return
      END
!
      SUBROUTINE four1(data,nn,isign)
      INTEGER isign,nn
      REAL data(2*nn)
      INTEGER i,istep,j,m,mmax,n
      REAL tempi,tempr
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      n=2*nn
      j=1
      do 11 i=1,n,2
        if(j.gt.i)then
          tempr=data(j)
          tempi=data(j+1)
          data(j)=data(i)
          data(j+1)=data(i+1)
          data(i)=tempr
          data(i+1)=tempi
        endif
        m=n/2
1       if ((m.ge.2).and.(j.gt.m)) then
          j=j-m
          m=m/2
        goto 1
        endif
        j=j+m
11    continue
      mmax=2
2     if (n.gt.mmax) then
        istep=2*mmax
        theta=6.28318530717959d0/(isign*mmax)
        wpr=-2.d0*sin(0.5d0*theta)**2
        wpi=sin(theta)
        wr=1.d0
        wi=0.d0
        do 13 m=1,mmax,2
          do 12 i=m,n,istep
            j=i+mmax
            tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
            tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
            data(j)=data(i)-tempr
            data(j+1)=data(i+1)-tempi
            data(i)=data(i)+tempr
            data(i+1)=data(i+1)+tempi
12        continue
          wtemp=wr
          wr=wr*wpr-wi*wpi+wr
          wi=wi*wpr+wtemp*wpi+wi
13      continue
        mmax=istep
      goto 2
      endif
      return
      END

max_stars_repo_name: BrunoSanchez/fisicaComp

max_issues_repo_name: BrunoSanchez/fisicaComp

max_forks_repo_name: BrunoSanchez/fisicaComp


text:
! SolveMatrixEquation   Routines for solution of R(X)=0
!
! Main subroutines
!   sme_setparam        Set parameters for solving of matrix equations. Required to call, but can be called without parameters
!   numjac              Calculate the numerical jacobian
!   solve_step          Solve a step, i.e. solve for dx so that A(x0+dx)=0
!   invert_jacobian     Calculate the inverted jacobian from the lu-factorized coming from solve_step
!   newton_raphson_num  Solve the equation system R(X)=0 using numerical jacobian
!   newton_raphson_ana  Solve the equation system R(X)=0 using analytical jacobian
!
! --------------------------------------------------------------------------------------------------------------
! Test subroutines (Use for testing and example)
!   test_residual       An example residual function, 3 unknowns
!   test_anajac         Analytical jacobian for test_residual
! End of SolveMatrixEquation description

module SolveMatrixEquation
implicit none
! Default parameters
integer, parameter          :: DEF_maxiter = 10         !Maximum number of iterations for newton-raphson iterations
double precision, parameter :: DEF_itertol = 1.d-6      !Tolerance on norm(residual) for newton-raphson iterations
double precision, parameter :: DEF_numdeps = 1.d-8      !Pertubation for numerical jacobian
logical, parameter          :: DEF_symnumjac = .false.  !Set true to use symmetrical calculation of jacobian (slower)
integer, parameter          :: DEF_line_search_start = 0!Default to start imediately
double precision, parameter :: DEF_alpha_min = 1.d-3     ! Default is to take minimum 1/100 of step size
double precision, parameter :: DEF_alpha_mult = 5.d-1     ! Default is to take half step at every new attempt

double precision, parameter :: infty = huge(1.d0)       !Set value to compare to check for infinity

! Define actual parameters
integer          :: SME_maxiter     !Maximum number of iterations for newton-raphson iterations
double precision :: SME_itertol     !Tolerance on norm(residual) for newton-raphson iterations
double precision :: SME_numdeps     !Pertubation for numerical jacobian
logical          :: SME_symnumjac   !Set true to use symmetrical calculation of jacobian (slower)
integer          :: SME_line_search_start   ! From what iteration to initiate line search if residual increases
double precision :: SME_alpha_min           ! Minimum value of alpha before failing
double precision :: SME_alpha_mult          ! Multiplication factor for reducing alpha

    contains
!Setup parameters for SolveMatrixEquation
subroutine sme_setparam(maxiter, itertol, pertubation, symnumjac, line_search_start, alpha_min, alpha_mult)
    !Set parameters for SolveMatrixEquation module (see specification statements of module)
    !Note: All options are optional, if not specified default values will be used
    !To call routine, use call sme_setparam(maxiter=NN, itertol=nn.d-n, pertubation=nn.d-n, symnumjac=.true.),
    !noting that you don't need to specify all
    implicit none
    integer, optional           :: maxiter
    double precision, optional  :: itertol
    double precision, optional  :: pertubation
    logical, optional           :: symnumjac
    integer, optional           :: line_search_start
    double precision, optional  :: alpha_min
    double precision, optional  :: alpha_mult

    SME_maxiter             = DEF_maxiter
    SME_itertol             = DEF_itertol
    SME_numdeps             = DEF_numdeps
    SME_symnumjac           = DEF_symnumjac
    SME_line_search_start   = DEF_line_search_start
    SME_alpha_min           = DEF_alpha_min
    SME_alpha_mult          = DEF_alpha_mult

    !Set maxiter
    if (present(maxiter)) then
        SME_maxiter = maxiter
    endif

    !Set itertol
    if (present(itertol)) then
        SME_itertol = itertol
    endif

    !Set pertubation (numdeps)
    if (present(pertubation)) then
        SME_numdeps = pertubation
    endif

    !Set condition for symmetric jacobian
    if (present(symnumjac)) then
        SME_symnumjac = symnumjac
    endif

    ! Line search start
    if (present(line_search_start)) then
        SME_line_search_start = line_search_start
    endif

    ! Line search minimum alpha value
    if (present(alpha_min)) then
        if (alpha_min<=0.d0) then  ! Check needed to prevent infinite loop
            write(*,*) 'alpha_min must be between 0 and 1 (strictly larger than 0!'
            stop
        endif

        SME_alpha_min = alpha_min
    endif

    ! Line search alpha reduction factor
    if (present(alpha_mult)) then
        if (alpha_mult>=1.d0) then  ! Check needed to prevent infinite loop
            write(*,*) 'alpha_mult must be strictly less than 1.0'
            stop
        endif
        SME_alpha_mult = alpha_mult
    endif

end subroutine sme_setparam

!Calculate numerical jacobian
subroutine numjac(residual_function, x0, jacobian, props, statev, param)
    !residual: Residual subroutine
    !x0      : Input variable to residual subroutine where the jacobian should be calculated
    !param   : Set of parameters used in the residual subroutine
    !jacobian: Output from subroutine, giving the numerical jacobian of the residual function at x0

    implicit none
    double precision, intent(in)                :: x0(:)
    double precision, optional, intent(in)      :: props(:), statev(:), param(:)!Input parameters
    double precision, allocatable, intent(out)  :: jacobian(:,:)                !Output arguments
    double precision, allocatable               :: dx(:), r0(:), rc(:)          !Internal allocatable float variables
    integer                                     :: nvar, k1, ninput             !Internal integers

    interface   !Interface for residual function
        subroutine residual_function(r, x, props, statev, param)
            double precision, intent(in)                :: x(:)
            double precision, optional, intent(in)      :: props(:), statev(:), param(:)      !Input arguments
            double precision, allocatable, intent(out)  :: r(:)
        end subroutine residual_function
    end interface

    nvar  = size(x0)  ! Size of problem

    !Allocate internal variables
    allocate(dx(nvar), r0(nvar), rc(nvar), jacobian(nvar,nvar))

    if (SME_symnumjac) then
        do k1=1,nvar
            dx = 0.d0
            dx(k1) = SME_numdeps/2.d0
            call residual_function(rc, x0+dx, props, statev, param)    !Get residual rc at x0+dx
            call residual_function(r0, x0-dx, props, statev, param)    !Get residual r0 at x0-dx
            jacobian(:,k1) = (rc-r0)/SME_numdeps   !Calculate k1 column in jacobian
        end do
    else
        call residual_function(r0, x0, props, statev, param)    !Get initial residual r0
        do k1=1,nvar
            dx = 0.d0
            dx(k1) = SME_numdeps
            call residual_function(rc, x0+dx, props, statev, param) !Get current residual rc
            jacobian(:,k1) = (rc-r0)/SME_numdeps   !Calculate k1 column in jacobian
        end do
    endif

end subroutine numjac

!Solve one step in newton iteration
subroutine solve_step(jacobian, r0, dx, lconv, ipiv, info)
! Solve jacobian*dx+r0=0 => jacobian*dx = -r0
! Note: jacobian is returned as lu-factorized, as well as pivot indicies ipiv
! usable in dgetri directly instead of calling dgetrf first
! dgesv documentation: https://software.intel.com/en-us/node/468876
    implicit none
    double precision, intent(inout) :: jacobian(:,:)    !Outputted ad LU-factorized
    double precision, intent(in)    :: r0(:)            !Residual at current x0
    double precision, intent(out)   :: dx(:)            !Suggested step
    logical, intent(out)            :: lconv            !Tells if solution succeeded
    integer, intent(out)            :: info             !Information about solution
    integer                         :: nvar             !Number of variables
    integer,intent(out)             :: ipiv(:)          !pivot indices for LU decomposition

    dx = -r0    !Input b in dgesv is overwritten to answer x
    nvar = size(r0)

    !    dgesv( n,   nrhs,        a,  lda, ipiv,  b,  ldb, info )
    call dgesv(nvar,    1, jacobian, nvar, ipiv, dx, nvar, info)
    lconv = (info==0)
end subroutine solve_step

!Invert matrix (typically jacobian)
subroutine invert_jacobian(lu_jacobian, ipiv)
    !lu_jacobian is the jacobian lu-factorized using either solve_step, or dgesv/dgetrf directly.
    !ipiv (indicies of pivot elements) are also output from solve_step, dgesv and dgetrf required
    !lu_jacobian is modified in this subroutine to become the inverse of the jacobian matrix
    implicit none
    double precision, intent(in) :: lu_jacobian(:,:)
    integer                      :: nvar, lwork, info
    integer, intent(out)         :: ipiv(:)
    double precision, allocatable:: work(:)
    nvar = size(lu_jacobian, dim=1)
    lwork = 2*nvar  !2 to be safe, to optimize see dgetri documentation
    !Probably not critical, as this routine is normally only called once per material execution

    allocate(work(lwork))
    !    dgetri(   n,           a,  lda, ipiv, work, lwork, info )
    call dgetri(nvar, lu_jacobian, nvar, ipiv, work, lwork, info)

end subroutine invert_jacobian

!Solve equation system using numerical jacobian
subroutine newton_raphson_num(residual_function, x0, jac_inv, lconv, info, props, statev, param, rvalues, xvalues, printres)
!Main routine for use to solve equation system with numerical jacobian:
!Specify the residual_function (by giving the actual name of the subroutine defining the residual function for the local problem)
!The inverse jacobian will be given as output for further use. lconv is true if convergence achieved, otherwise it is false
implicit none
    double precision,allocatable,intent(out):: jac_inv(:,:)         !Inverted jacobian at stationary point returned (different uses during subroutine)
    double precision, intent(inout)         :: x0(:)                !Start guess
    double precision, intent(in), optional  :: props(:), statev(:), param(:)    !Residual function parameters
    integer,          intent(in), optional  :: printres             !Logical if print residual or not
    logical,          intent(out)           :: lconv                !Logical if convergence or not
    integer                                 :: k1, k2               !Looping iterator
    integer                                 :: ipiv(size(x0))       !Indicies of pivot elements
    double precision                        :: dx(size(x0))         !Increment
    double precision, allocatable           :: r0(:)                !Residual
    integer                                 :: info                 !Information about solution of step (dgesv)
    logical                                 :: diagnostic           !Logical to describe if diagnostic information should be given or not
    integer                                 :: printinfo            !Integer to decide what info to print
    double precision, optional, allocatable :: rvalues(:,:)         !Optional diagnostics information - residual values
    double precision, optional, allocatable :: xvalues(:,:)         !Optional diagnostics information - x-values
    double precision, allocatable           :: rvalues_full(:,:)    !Used to fill calculate rvalues
    double precision, allocatable           :: xvalues_full(:,:)    !Used to fill calculate xvalues
    character(len=20)                       :: jacobian_format      !Used to print jacobian nicer

    interface   !interface for residual function
        subroutine residual_function(r, x, props, statev, param)
            double precision, intent(in)                :: x(:)
            double precision, intent(in), optional      :: props(:), statev(:), param(:)      !Residual function parameters
            double precision, allocatable, intent(out)  :: r(:)
        end subroutine residual_function
    end interface

    diagnostic = present(rvalues).or.present(xvalues)   !Check if diagnostic info should be saved
    ! If diagnostic, then allocate the full 2d arrays for saving information
    if (diagnostic) then
        allocate(rvalues_full(size(x0), SME_maxiter+1), xvalues_full(size(x0), SME_maxiter+1))
    endif

    ! Should residual norm be printed out for each iterations?
    if (present(printres)) then
        printinfo = printres       !Same if present
    else
        printinfo = 0              !Don't print if not specified
    endif

    call residual_function(r0, x0, props, statev, param)

    ! If diagnostic, then save the current state
    if (diagnostic) then
        rvalues_full(:, 1) = r0
        xvalues_full(:, 1) = x0
    endif

    ! If print residual, print
    if (printinfo==1) then
        write(*,*) 'Printing iteration nr and residual values:'
        write(*,"(I2,E15.5)") 0, sqrt(sum(r0*r0))
    elseif (printinfo==2) then
        write(*,*) 'Printing full variable and residual'
        write(*,*) 'x'
        write(*,*) x0
        write(*,"(A,E15.5)") 'r(x): ', sqrt(sum(r0*r0))
        write(*,*) r0
        write(*,*) ' '
    elseif (printinfo==3) then
        write(*,*) 'Printing full variable, residual and jacobian (transpose)'
        write(jacobian_format,"(I0)") size(r0)
        jacobian_format = '('//trim(adjustl(jacobian_format))//'E15.5)'
    endif

    do k1 = 1,SME_maxiter
        call numjac(residual_function, x0, jac_inv, props, statev, param)  !Get jacobian
        if (badcontent2(jac_inv)) then
            write(*,*) 'Jacobian NaN or Inf. Printing x'
            write(*,*) x0
            lconv = .false.
            return
        endif
        if (printinfo==3) then
            write(*,"(A,I2)") 'x_', k1
            write(*,*) x0
            write(*,"(A,E15.5)") 'r(x): ', sqrt(sum(r0*r0))
            write(*,*) r0
            write(*,"(A)") 'Jacobian: '
            do k2=1,size(r0)
                write(*,jacobian_format) jac_inv(:,k2)
            enddo
            write(*,*) ' '
        endif
        call solve_step(jac_inv, r0, dx, lconv, ipiv, info) !Get step increment dx
        x0 = x0 + dx                                        !Update step
        call residual_function(r0, x0, props, statev, param)               !Get new residual
        if (badcontent1(r0)) then
            write(*,*) 'Residual NaN or Inf. Printing x and r'
            write(*,*) x0
            write(*,*) r0
            lconv = .false.
            return
        endif
        ! If diagnostic, then save the current state
        if (diagnostic) then
            rvalues_full(:, 1+k1) = r0
            xvalues_full(:, 1+k1) = x0
        endif

        ! If print residual, print
        if (printinfo==1) then
            write(*,"(I2,E15.5)") k1, sqrt(sum(r0*r0))
        elseif (printinfo==2) then
            write(*,"(A,I2)") 'x_', k1
            write(*,*) x0
            write(*,"(A,E15.5)") 'r(x): ', sqrt(sum(r0*r0))
            write(*,*) r0
            write(*,*) ' '
        endif

        if (sqrt(sum(r0*r0))<SME_itertol) then              !Check for convergence
            lconv = .true.                                  !Solution converged
            call invert_jacobian(jac_inv, ipiv)             !Calculate inverse jacobian

            ! If diagnostic, allocate and save the diagnostic information to rvalues and xvalues as requested
            if (diagnostic) then
                if (present(rvalues)) then
                    allocate(rvalues(size(x0), 1+k1))
                    rvalues = rvalues_full(:, 1:(1+k1))
                endif
                if (present(xvalues)) then
                    allocate(xvalues(size(x0), 1+k1))
                    xvalues = xvalues_full(:, 1:(1+k1))
                endif
            endif

            return                                          !End subroutine when convergence achieved
        endif
    enddo
    lconv = .false. !If subroutine didn't end from convergence, r0>SME_itertol

end subroutine newton_raphson_num

!Solve equation system using analytical jacobian
subroutine newton_raphson_ana(residual_function, jacobian_function, x0, jac_inv, lconv, info, props, statev, param, rvalues, xvalues, printres)
!Main routine for use to solve equation system with analytical jacobian:
!Specify the residual_function (by giving the actual name of the subroutine defining the residual function for the local problem)
!Specify the jacobian_function (by giving the actual name of the subroutine defining the analytical jacobian for the local problem)
!The inverse jacobian will be given as output for further use. lconv is true if convergence achieved, otherwise it is false
implicit none
    double precision,allocatable,intent(out):: jac_inv(:,:)  !Inverted jacobian at stationary point returned (different uses during subroutine)
    double precision, intent(inout)         :: x0(:)         !Start guess
    double precision, intent(in), optional  :: props(:), statev(:), param(:)      !Residual function parameters
    integer,          intent(in), optional  :: printres             !Logical if print residual or not
    logical,          intent(out)           :: lconv         !Logical if convergence or not
    integer                                 :: k1, k2            !Looping iterator
    integer                                 :: ipiv(size(x0))!Indicies of pivot elements
    double precision                        :: dx(size(x0))  !Increment
    double precision, allocatable           :: r0(:)         !Residual
    integer                                 :: info          !Information about solution of step (dgesv)
    logical                                 :: diagnostic    !Logical to describe if diagnostic information should be given or not
    integer                                 :: printinfo
    double precision, optional, allocatable :: rvalues(:,:)  !Optional diagnostics information - residual values
    double precision, optional, allocatable :: xvalues(:,:)  !Optional diagnostics information - x-values
    double precision, allocatable           :: rvalues_full(:,:)       !Used to fill calculate rvalues
    double precision, allocatable           :: xvalues_full(:,:)       !Used to fill calculate xvalues
    double precision                        :: residual, residual_old   ! Current and previous residual (scalar = sqrt(sum(r0**2)) )
    double precision                        :: alpha
    character(len=20)                       :: jacobian_format      !Used to print jacobian nicer

    interface   !Interface for residual function
        subroutine residual_function(r, x, props, statev, param)
            double precision, intent(in)                :: x(:)
            double precision, intent(in), optional      :: props(:), statev(:), param(:)      !Residual function parameters
            double precision, allocatable, intent(out)  :: r(:)
        end subroutine residual_function
    end interface

    interface   !Interface for jacobian function
        subroutine jacobian_function(j, x, props, statev, param)
            double precision, intent(in)                :: x(:)
            double precision, intent(in), optional      :: props(:), statev(:), param(:)      !Residual function parameters
            double precision, allocatable, intent(out)  :: j(:,:)
        end subroutine jacobian_function
    end interface
    diagnostic = present(rvalues).or.present(xvalues)   !Check if diagnostic info should be saved

    ! Should residual norm be printed out for each iterations?
    if (present(printres)) then
        printinfo = printres       !Same if present
    else
        printinfo = 0              !Don't print if not specified
    endif

    ! If diagnostic, then allocate the full 2d arrays for saving information
    if (diagnostic) then
        allocate(rvalues_full(size(x0), SME_maxiter+1), xvalues_full(size(x0), SME_maxiter+1))
    endif

    !Get initial residual
    call residual_function(r0, x0, props, statev, param)
    residual_old = sqrt(sum(r0**2))

    ! If diagnostic, then save the current state
    if (diagnostic) then
        rvalues_full(:, 1) = r0
        xvalues_full(:, 1) = x0
    endif

    ! If print residual, print
    if (printinfo==1) then
        write(*,*) 'Printing iteration nr and residual values:'
        write(*,"(I2,E15.5)") 0, residual_old
    elseif (printinfo==2) then
        write(*,*) 'Printing full variable and residual'
        write(*,*) 'x'
        write(*,*) x0
        write(*,"(A,E15.5)") 'r(x): ', residual_old
        write(*,*) r0
        write(*,*) ' '
    elseif (printinfo==3) then
        write(*,*) 'Printing full variable, residual and jacobian (transpose)'
        write(jacobian_format,"(I0)") size(r0)
        jacobian_format = '('//trim(adjustl(jacobian_format))//'E15.5)'
    endif

    do k1 = 1,SME_maxiter
        call jacobian_function(jac_inv, x0, props, statev, param)      !Get jacobian
        if (diagnostic) then
            if (badcontent2(jac_inv)) then  ! Can take some time, so only check when in diagnostic mode
                write(*,*) 'Jacobian NaN or Inf. Printing x'
                write(*,*) x0
                lconv = .false.
                return
            endif
        endif
        if (printinfo==3) then
            write(*,"(A,I2)") 'x_', k1
            write(*,*) x0
            write(*,"(A,E15.5)") 'r(x): ', sqrt(sum(r0*r0))
            write(*,*) r0
            write(*,"(A)") 'Jacobian: '
            do k2=1,size(r0)
                write(*,jacobian_format) jac_inv(:,k2)
            enddo
            write(*,*) ' '
        endif
        call solve_step(jac_inv, r0, dx, lconv, ipiv, info) !Get step increment dx
        x0 = x0 + dx                                        !Update step

        call residual_function(r0, x0, props, statev, param)               !Get new residual

        residual = sqrt(sum(r0**2))
        if (badcontent0(residual)) then
            write(*,*) 'Residual NaN or Inf. Printing x and r'
            write(*,*) x0
            write(*,*) r0
            lconv = .false.
            return
        endif

        alpha = 1.d0
        do while ( (k1>SME_line_search_start).and.(residual>residual_old).and.(alpha>=(SME_alpha_min+1.d-14)) )
            x0 = x0 - alpha*dx                                      ! Reset step
            alpha = alpha*SME_alpha_mult                            ! Calculate new alpha
            if (alpha<SME_alpha_min) then
                alpha = SME_alpha_min
            endif
            x0 = x0 + alpha*dx                                      !Update step
            call residual_function(r0, x0, props, statev, param)    !Get new residual

            residual = sqrt(sum(r0**2))
            if (badcontent0(residual)) then
                write(*,*) 'Residual NaN or Inf'
                lconv = .false.
                return
            endif
        enddo
        residual_old = residual

        ! If diagnostic, then save the current state
        if (diagnostic) then
            rvalues_full(:, 1+k1) = r0
            xvalues_full(:, 1+k1) = x0
        endif

        ! If print residual, print
        if (printinfo==1) then
            write(*,"(I2,E15.5)") k1, residual
        elseif (printinfo==2) then
            write(*,"(A,I2)") 'x_', k1
            write(*,*) x0
            write(*,"(A,E15.5)") 'r(x): ', residual
            write(*,*) r0
            write(*,*) ' '
        endif

        if (residual<SME_itertol) then                      !Check for convergence
            lconv = .true.                                  !Solution converged
            call invert_jacobian(jac_inv, ipiv)             !Calculate inverse jacobian

            ! If diagnostic, allocate and save the diagnostic information to rvalues and xvalues as requested
            if (diagnostic) then
                if (present(rvalues)) then
                    allocate(rvalues(size(x0), 1+k1))
                    rvalues = rvalues_full(:, 1:(1+k1))
                endif
                if (present(xvalues)) then
                    allocate(xvalues(size(x0), 1+k1))
                    xvalues = xvalues_full(:, 1:(1+k1))
                endif
            endif

            return  !End subroutine when convergence achieved
        endif !End check of convergence
    enddo   !End loop
    lconv = .false. !If subroutine didn't end from convergence, r0>SME_itertol

end subroutine newton_raphson_ana

function badcontent0(scalar) result(isbad)
implicit none
    double precision scalar
    logical          isbad
    isbad = (scalar/=scalar).or.(abs(scalar)>=infty)
    !           NaN /= NaN      infty largest representable double precision (def as parameter)
end function

function badcontent1(vector) result(isbad)
implicit none
    double precision vector(:)
    logical          isbad
    isbad = (any(vector/=vector)).or.(any(abs(vector)>=infty))
    !           NaN /= NaN                          infty largest representable double precision (def as parameter)
end function

function badcontent2(matrix) result(isbad)
implicit none
    double precision matrix(:,:)
    logical          isbad
    isbad = (any(matrix/=matrix)).or.(any(abs(matrix)>=infty))  ! See badcontent1
end function

! ================================================================================================
! Test subroutines (use for reference and testing of modifications)
subroutine test_residual(r, x, props, statev, param)
    implicit none
    double precision, intent(in)                :: x(:)
    double precision, intent(in), optional      :: props(:), statev(:), param(:)      !Residual function parameters
    double precision, allocatable, intent(out)  :: r(:)
    allocate(r(size(x)))
    if (param(1)==1.d0) then
        r(1) = x(1)+2*x(2)
        r(2) = 3*x(1)+4*x(2)
        r(3) = x(3)
    else
        r(1) = x(1)*x(2)+1.0
        r(2) = x(1)**2-x(2)**2
        r(3) = x(3)
    endif
end subroutine test_residual

subroutine test_anajac(ajacobian, x, props, statev, param)
    implicit none
    double precision, intent(in)                :: x(:)
    double precision, intent(in), optional      :: props(:), statev(:), param(:)      !parameters
    double precision, allocatable, intent(out)  :: ajacobian(:,:)

    allocate(ajacobian(size(x),size(x)))
    if (param(1)==1.d0) then
        ajacobian(1,:) = (/1.d0, 2.d0, 0.d0/)
        ajacobian(2,:) = (/3.d0, 4.d0, 0.d0/)
        ajacobian(3,:) = (/0.d0, 0.d0, 1.d0/)
    else
        ajacobian(1,:) = (/x(2), x(1), 0.d0/)
        ajacobian(2,:) = (/2.d0*x(1), -2.d0*x(2), 0.d0/)
        ajacobian(3,:) = (/0.d0, 0.d0, 1.d0/)
    endif
end subroutine test_anajac

end module

max_stars_repo_name: kengwit/MaterialModels

max_issues_repo_name: kengwit/MaterialModels

max_forks_repo_name: kengwit/MaterialModels


text:
Court N Cedar is a collection of shops and religious and spiritual organizations religious centers at 231 G Street. It is an openair, twostory complex with a water fountains water fountain in the central courtyard. The second story spaces are mostly deserted. It is comprised of:

  Beautiful Massage
  Christian Science Reading Room
  Resler Brothers Barber Shop
  Smashwords
  State Farm Insurance (Heather Copeland)
  INDIGO/Hammond & Playle Architects, LLP


! =================================================
!   Main program.
! =================================================

program main

    use coordinates,    only: Nx, Ny
    use numeric_solver, only: calc_numeric => iterate
    use analytic_solver
    implicit none

    real*8 :: Tn(Nx,Ny), Tn_(Nx,Ny)
    real*8 :: Ta(Nx,Ny), Ta_(Nx,Ny)
    real*8 :: avg, std

    call init_global_vars()

    ! Numerical solution
    call calc_numeric(Tn, Tn_, method=3)
    print *, "Numerical solution:"
    call print_matrix(Tn_)

    ! Analytical solution
    call calc_analytic(Ta, Ta_)
    print *, "Analytical solution:"
    call print_matrix(Ta_)

    ! Compare numerical solution with analytical solution
    call comp_analytic(Tn_, Ta_, avg, std)
    print *, "Average error:     ", avg
    print *, "Standard deviation:", std

    pause

end program


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

C -*- Mode: fortran -*-
C=======================================================================
      SUBROUTINE COMAND (IDNPS, IDESS, IDNSET, IDESET,
     &   IDELB, NUMELB, NUMLNK, NAMELB, ELATTR,
     &  XN, YN, A, *)
C=======================================================================

C   $Id: comand.f,v 1.13 2005/01/17 23:09:58 gdsjaar Exp $

C   --*** COMAND *** (GENSHELL) Input and process commands
C   --
C   --COMAND inputs and executes an user command.
C   --
C   --Parameters:
C   --   IDNPS - IN - the IDs of existing node sets
C   --   IDESS - IN - the IDs of existing side sets
C   --   IDNSET - OUT - the IDs of the front and back surface node sets;
C   --      (0) = length
C   --   IDESET - OUT - the IDs of the front and back surface side sets;
C   --      (0) = length
C   --   IDELB - IN - the ids for each block
C   --   NUMELB - IN - the number of elements for each block
C   --   NUMLNK - IN - the number of nodes per element for each block
C   --   XN, YN - IN - the 2D nodal coordinates
C   --   * - return statement for QUIT
C   --
C   --Common Variables:
C   --   Sets NDIM, NUMNP, NUMEL, NELBLK,
C   --      NUMNPS, LNPSNL, NUMESS, LESSEL, LESSNL of /DBNUMS/
CCC   --   Sets ITRANT, NNREPL, NEREPL, DIM3, NRTRAN, D3TRAN, ZGRAD,
CCC   --      CENTER, NUMCOL of /PARAMS/
C   --   Sets XOFFS, YOFFS, ZOFFS of /XYZOFF/
C   --   Sets ROT3D, ROTMAT, ROTCEN of /XYZROT/

      PARAMETER (MAXFLD = 32)
      PARAMETER (MAXSET = 10)

      INCLUDE 'gs_dbnums.blk'
      INCLUDE 'gs_params.blk'
      INCLUDE 'gs_xyzoff.blk'
      INCLUDE 'gs_xyzrot.blk'
      INCLUDE 'gs_xyzmir.blk'
      INCLUDE 'gs_xyzero.blk'
      INCLUDE 'gs_xyzscl.blk'
      INCLUDE 'gs_splxyz.blk'

      INTEGER IDNPS(*), IDESS(*)
      INTEGER IDNSET(0:MAXSET,2), IDESET(0:MAXSET,2)
      INTEGER IDELB(*)
      CHARACTER*32 NAMELB(*)
      REAL ELATTR(7, NELBLK)
      INTEGER NUMELB(*)
      INTEGER NUMLNK(*)
      REAL XN(*), YN(*)
      REAL A(*)

      LOGICAL FFEXST, MATSTR, HELP

      CHARACTER*8 WORD, VERB, TMPWRD
      INTEGER     INTYP(MAXFLD+1)
      CHARACTER*8 CFIELD(MAXFLD)
      INTEGER     IFIELD(MAXFLD)
      REAL        RFIELD(MAXFLD)
      LOGICAL ISFRON
      LOGICAL ISHELP

c$$$      CHARACTER*8 CMDTBL(34)
c$$$      SAVE CMDTBL
c$$$C      --CMDTBL - the valid commands table
c$$$
c$$$C   --Command table follows.  Remember to change the dimensioned size when
c$$$C   --changing the table.
c$$$      DATA CMDTBL /
c$$$     1   'TRANSLAT', 'ROTATE  ', 'OFFSET  ', 'REVOLVE ', 'REVCEN  ',
c$$$     2   'WARP    ', 'BLOCK   ', 'CENTER  ', 'TUNNEL  ', 'SPECIAL ',
c$$$     3   'NSETS   ', 'NODESETS', 'SSETS   ', 'SIDESETS', 'MIRROR  ',
c$$$     4   'LIST    ', 'SHOW    ', 'HELP    ', 'ZERO    ', 'SCALE   ',
c$$$     5   'END     ', 'EXIT    ', 'QUIT    ', 'TWIST   ', 'PROJECT ',
c$$$     6   'SUMMARY ', 'INTERVAL', 'ROTCEN  ', 'SPAWN   ', 'SPLINE  ',
c$$$     7   'CHANGE  ', 'TRANSPLI', 'SHIFT   ', '        ' /
c$$$

      CHARACTER*8 CMDTBL(27)
      SAVE CMDTBL
C      --CMDTBL - the valid commands table

C   --Command table follows.  Remember to change the dimensioned size when
C   --changing the table.
      DATA CMDTBL /
     1   'TRANSLAT', 'OFFSET  ', 'REVOLVE ', 'REVCEN  ',
     2   'WARP    ',
     3   'NSETS   ', 'NODESETS', 'SSETS   ', 'SIDESETS', 'MIRROR  ',
     4   'LIST    ', 'SHOW    ', 'HELP    ', 'ZERO    ', 'SCALE   ',
     5   'END     ', 'EXIT    ', 'QUIT    ', 'RANDOMIZ',
     6   'SUMMARY ', 'SPAWN   ', 'SPLINE  ', 'ATTRIBUT',
     7   'CHANGE  ', 'SHIFT   ', 'THICKNES', '        ' /

C   --Initialize

      ITRANT = 0
      DIM3  = 1.0
      XOFFS = 0.0
      YOFFS = 0.0
      ZOFFS = 0.0
      XMIRR = 1.0
      YMIRR = 1.0
      ZMIRR = 1.0
      XZERO = 0.0
      YZERO = 0.0
      ZZERO = 0.0
      XSCAL = 1.0
      YSCAL = 1.0
      ZSCAL = 1.0
      XRAND = 0.0
      YRAND = 0.0
      ZRAND = 0.0

      ROT3D = .FALSE.
      CALL INIREA (3*3, 0.0, ROTMAT)
      DO 10 I = 1, 3
         ROTMAT(I,I) = 1.0
   10 CONTINUE

C     -- Initialize the elattr array here.  These will be used for the
C     -- attributes of any 3d BEAM blocks. User can change using
C     -- ATTRIBUTE command.
      do 15 i = 1, nelblk
        elattr(1,i) = 1.0
        elattr(2,i) = 1.0
        elattr(3,i) = 1.0
        elattr(4,i) = 1.0
        elattr(5,i) = 0.0
        elattr(6,i) = 0.0
        elattr(7,i) = 1.0
 15   continue

      CALL MINMAX (NUMNP, XN, XMIN, XMAX)
      ROTCEN(1) = XMIN
      CALL MINMAX (NUMNP, YN, YMIN, YMAX)
      ROTCEN(2) = YMIN
      ROTCEN(3) = 0.0

      IDNSET(0,1) = 0
      IDNSET(0,2) = 0
      IDESET(0,1) = 0
      IDESET(0,2) = 0

   20 CONTINUE

C   --Read command line

      WRITE (*, *)
      CALL FREFLD (0, 0, 'GenShell> ', MAXFLD,
     &   IOSTAT, NUMFLD, INTYP, CFIELD, IFIELD, RFIELD)
      IF (IOSTAT .LT. 0) GOTO 220
      IF (NUMFLD .EQ. 0) GOTO 20
      INTYP(MIN(MAXFLD,NUMFLD)+1) = -999

      IFLD = 1
      CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
      CALL ABRSTR (VERB, WORD, CMDTBL)
      IF (VERB .EQ. ' ') VERB = WORD

C   --Perform command

      IF (VERB .EQ. 'TRANSLAT') THEN
         ITRANT = 1
         CALL FFREAL (IFLD, INTYP, RFIELD,
     *      'translation thickness', 1.0, DIM3, *170)

      ELSE IF (VERB .EQ. 'THICKNES') THEN
         if (itrant .eq. 0) then
           call prterr('CMDERR',
     $       'A Mesh Transformation command must be specified')
            goto 170
         else
            CALL FFREAL (IFLD, INTYP, RFIELD,
     *           'thickness', 1.0, DIM3, *170)
         end if

      ELSE IF (VERB .EQ. 'SPLINE') THEN
         ITRANT = 64
         CALL FFREAL (IFLD, INTYP, RFIELD,
     *      'spline thickness', 1.0, DIM3, *170)

         CALL GETSPL (A)

C --- BEGINNING OF WARP
      ELSE IF (VERB .EQ. 'WARP') THEN
         ITRANT = 4
         CALL FFCHAR (IFLD, INTYP, CFIELD,
     *      'POINT', WORD)
         IF((.NOT. MATSTR(WORD, 'LINE',    1)) .AND.
     *      (.NOT. MATSTR(WORD, 'XAXIS',   1)) .AND.
     *      (.NOT. MATSTR(WORD, 'YAXIS',   1)) .AND.
     *      (.NOT. MATSTR(WORD, 'ELLIPSE', 1)) .AND.
     *      (.NOT. MATSTR(WORD, 'POINT',   1))) THEN
            CALL PRTERR ('CMDERR', 'Invalid WARP Option')
            GOTO 170
         END IF

         IF (INTYP(IFLD) .EQ. 0) THEN
            CALL FFCHAR (IFLD, INTYP, CFIELD, 'MAP', TMPWRD)
         ELSE
            TMPWRD = 'MAP'
         END IF

         CALL FFREAL (IFLD, INTYP, RFIELD,
     *      'warping thickness', 1.0, DIM3, *170)

         IF (MATSTR(WORD, 'LINE', 1)) THEN
            IWARP = 2

C ... For XAXIS and YAXIS warps, we can either map vertically onto
C     the cylinder, keeping both X and Y coordinates the same as in
C     the 2D mesh.  Or, we can map onto the surface keeping the length
C     of the mesh the same in the 2D and 3D meshes.

         ELSE IF (MATSTR(WORD, 'XAXIS', 1)) THEN
            IF (MATSTR(TMPWRD, 'VERTICAL', 1)) THEN
               IWARP = -3
            ELSE
               IWARP = -1
            END IF

         ELSE IF (MATSTR(WORD, 'YAXIS', 1)) THEN
            IF (MATSTR(TMPWRD, 'VERTICAL', 1)) THEN
               IWARP = -4
            ELSE
               IWARP = -2
            END IF

         ELSE IF (MATSTR(WORD, 'POINT', 1)) THEN
            IWARP = 1
         ELSE IF (MATSTR(WORD, 'ELLIPSE', 1)) THEN
            IWARP = 2
            CALL FFREAL (IFLD, INTYP, RFIELD,
     *         'ellipse Y axis', 0.0, XRAD, *170)
            CALL FFREAL (IFLD, INTYP, RFIELD,
     *         'ellipse Z axis', 0.0, YRAD, *170)
         ELSE
            IWARP = 0
            GO TO 170
         END IF
         CALL FFREAL (IFLD, INTYP, RFIELD,
     *      'warping distance', 0.0, DWARP, *170)

C --- END OF WARP

      ELSE IF (VERB .EQ. 'OFFSET' .OR. VERB .EQ. 'SHIFT') THEN

C     ... Originally offset just asked for NDIM values.  It was changed
C     to go by axis type (OFFSET Y 1.0). We do maintain compatibility
C     and check for both types of input.

         RMULT = 0.0
         IF (INTYP(IFLD) .EQ. 0) THEN
  930       CONTINUE
            IF (FFEXST (IFLD, INTYP)) THEN
               CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
               IF (MATSTR (WORD, 'RESET', 1)) THEN
                  XOFFS = 0.0
                  YOFFS = 0.0
                  ZOFFS = 0.0
                  RMULT = 0.0
               ELSE IF (MATSTR (WORD, 'ADD', 2)) THEN
C     ... Set for cumulative offsets/shifts
                  RMULT = 1.0
               ELSE IF (MATSTR (WORD, 'ALL', 2)) THEN
                  CALL FFREAL (IFLD, INTYP, RFIELD,
     &                 'model offset', 0.0, TOFFS, *170)
                  XOFFS = RMULT * XOFFS + TOFFS
                  YOFFS = RMULT * YOFFS + TOFFS
                  ZOFFS = RMULT * ZOFFS + TOFFS
               ELSE IF (WORD .EQ. 'X') THEN
                  CALL FFREAL (IFLD, INTYP, RFIELD,
     &                 'X coordinate offset', 0.0, TOFFS, *170)
                  XOFFS = RMULT * XOFFS + TOFFS
               ELSE IF (WORD .EQ. 'Y') THEN
                  CALL FFREAL (IFLD, INTYP, RFIELD,
     &                 'Y coordinate offset', 0.0, TOFFS, *170)
                  YOFFS = RMULT * YOFFS + TOFFS
               ELSE IF (WORD .EQ. 'Z') THEN
                  CALL FFREAL (IFLD, INTYP, RFIELD,
     &                 'Z coordinate offset', 0.0, TOFFS, *170)
                  ZOFFS = RMULT * ZOFFS + TOFFS
               ELSE
                  CALL PRTERR ('CMDERR',
     &              'Expected "X", "Y", "Z", "ALL", "ADD", or "RESET"')
                  GOTO 170
               END IF
               GOTO 930
            END IF
         ELSE
            CALL FFREAL (IFLD, INTYP, RFIELD,
     &           'X coordinate offset', XOFFS, XOFFS, *170)
            CALL FFREAL (IFLD, INTYP, RFIELD,
     &           'Y coordinate offset', YOFFS, YOFFS, *170)
            CALL FFREAL (IFLD, INTYP, RFIELD,
     &           'Z coordinate offset', ZOFFS, ZOFFS, *170)
         END IF

      ELSE IF (VERB .EQ. 'SCALE') THEN
   70    CONTINUE
         IF (FFEXST (IFLD, INTYP)) THEN
            CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
            IF (MATSTR (WORD, 'RESET', 1)) THEN
               XSCAL = 1.0
               YSCAL = 1.0
               ZSCAL = 1.0
            ELSE IF (MATSTR (WORD, 'ALL', 1)) THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &            'model scale factor', 1.0, XSCAL, *80)
               XSCAL = ABS(XSCAL)
               YSCAL = XSCAL
               ZSCAL = XSCAL
            ELSE IF (WORD .EQ. 'X') THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &            'X scale factor', 1.0, XSCAL, *80)
               XSCAL = ABS(XSCAL)
            ELSE IF (WORD .EQ. 'Y') THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &            'Y scale factor', 1.0, YSCAL, *80)
               YSCAL = ABS(YSCAL)
            ELSE IF (WORD .EQ. 'Z') THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &            'Z scale factor', 1.0, ZSCAL, *80)
               ZSCAL = ABS(ZSCAL)
            ELSE
               CALL PRTERR ('CMDERR',
     &            'Expected "X", "Y", "Z", "ALL", or "RESET"')
               GOTO 170
            END IF
            GOTO 70
         END IF

   80    CONTINUE

      ELSE IF (VERB .EQ. 'ZERO') THEN
   90    CONTINUE
         IF (FFEXST (IFLD, INTYP)) THEN
            CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
            IF (MATSTR (WORD, 'RESET', 1)) THEN
               XZERO = 0.0
               YZERO = 0.0
               ZZERO = 0.0
            ELSE IF (MATSTR (WORD, 'ALL', 2)) THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &              'minimum coordinate', 0.0, XZERO, *100)
               YZERO = XZERO
               ZZERO = XZERO
            ELSE IF (WORD .EQ. 'X') THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &            'minimum X coordinate', 0.0, XZERO, *100)
            ELSE IF (WORD .EQ. 'Y') THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &            'minimum Y coordinate', 0.0, YZERO, *100)
            ELSE IF (WORD .EQ. 'Z') THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &            'minimum Z coordinate', 0.0, ZZERO, *100)
            ELSE
               CALL PRTERR ('CMDERR',
     &            'Expected "X", "Y", "Z", "ALL" or "RESET"')
               GOTO 100
            END IF
            GOTO 90
         END IF

 100     CONTINUE

      ELSE IF (VERB .EQ. 'RANDOMIZ') THEN
 101     CONTINUE
         IF (FFEXST (IFLD, INTYP)) THEN
            CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
            IF (MATSTR (WORD, 'RESET', 1)) THEN
               XRAND = 0.0
               YRAND = 0.0
               ZRAND = 0.0
            ELSE IF (MATSTR (WORD, 'ALL', 1)) THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &              'model random factor', 0.0, TRAND, *141)
               XRAND = ABS(TRAND)
               YRAND = ABS(TRAND)
               ZRAND = ABS(TRAND)
            ELSE IF (WORD .EQ. 'X') THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &              'X random factor', 0.0, TRAND, *141)
               XRAND = ABS(TRAND)
            ELSE IF (WORD .EQ. 'Y') THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &              'Y random factor', 0.0, TRAND, *141)
               YRAND = ABS(TRAND)
            ELSE IF (WORD .EQ. 'Z') THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &              'Z random factor', 0.0, TRAND, *141)
               ZRAND = ABS(TRAND)
            ELSE
               CALL PRTERR ('CMDERR',
     &              'Expected "X", "Y", "Z", "ALL", or "RESET"')
               GOTO 141
            END IF
            GOTO 101
         END IF

 141     CONTINUE

      ELSE IF (VERB .EQ. 'MIRROR') THEN

  110    CONTINUE
         IF (FFEXST (IFLD, INTYP)) THEN
            CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
            IF (MATSTR (WORD, 'RESET', 1)) THEN
               XMIRR = 1.
               YMIRR = 1.
               ZMIRR = 1.
            ELSE IF (WORD .EQ. 'X') THEN
               XMIRR = -1.
            ELSE IF (WORD .EQ. 'Y') THEN
               YMIRR = -1.
            ELSE IF (WORD .EQ. 'Z') THEN
               ZMIRR = -1.
            ELSE
               CALL PRTERR ('CMDERR',
     &            'Expected "X", "Y", "Z" or "RESET"')
               GOTO 120
            END IF
            GOTO 110
         END IF

  120    CONTINUE

      ELSE IF (VERB .EQ. 'REVOLVE') THEN
         DEGANG = ATAN2(0.0, -1.0) / 180.0

  130    CONTINUE
         IF (FFEXST (IFLD, INTYP)) THEN
            CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
            IF (MATSTR (WORD, 'RESET', 1)) THEN
               ROT3D = .FALSE.
               CALL INIREA (3*3, 0.0, ROTMAT)
               DO 140 I = 1, 3
                  ROTMAT(I,I) = 1.0
  140          CONTINUE
            ELSE IF ((WORD .EQ. 'X') .OR. (WORD .EQ. 'Y')
     &         .OR. (WORD .EQ. 'Z')) THEN
               CALL FFREAL (IFLD, INTYP, RFIELD,
     &            'angle of rotation', 0.0, DEG, *150)
               ROT3D = .TRUE.
               CALL ROTXYZ (WORD(:1), DEG * DEGANG, ROTMAT)
            ELSE
               CALL PRTERR ('CMDERR',
     &            'Expected "X", "Y", "Z" or "RESET"')
               GOTO 150
            END IF
            GOTO 130
         END IF

  150    CONTINUE

      ELSE IF (VERB .EQ. 'REVCEN') THEN
         CALL MINMAX (NUMNP, XN, XMIN, XMAX)
         CALL MINMAX (NUMNP, YN, YMIN, YMAX)

         CALL FFREAL (IFLD, INTYP, RFIELD,
     &      'X revolution center', XMIN, ROTCEN(1), *170)
         CALL FFREAL (IFLD, INTYP, RFIELD,
     &      'Y revolution center', YMIN, ROTCEN(2), *170)
         CALL FFREAL (IFLD, INTYP, RFIELD,
     &      'Z revolution center', 0.0, ROTCEN(3), *170)

      ELSE IF ((VERB .EQ. 'NSETS') .OR. (VERB .EQ. 'NODESETS')
     &   .OR. (VERB .EQ. 'SSETS') .OR. (VERB .EQ. 'SIDESETS')) THEN
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         IF ((.NOT. MATSTR (WORD, 'FRONT', 1)) .AND.
     &      (.NOT. MATSTR (WORD, 'BACK', 1))) THEN
            CALL PRTERR ('CMDERR', 'Expected FRONT or BACK')
            GOTO 170
         END IF
         ISFRON = MATSTR (WORD, 'FRONT', 1)

         IF ((VERB .EQ. 'NSETS') .OR. (VERB .EQ. 'NODESETS')) THEN
            IF (ISFRON) THEN
               CALL USIDS (IFLD, INTYP, CFIELD, IFIELD,
     &            NUMNPS, IDNPS, IDNSET(0,2), IDNSET(1,2),
     &            IDNSET(0,1), IDNSET(1,1), *170)
            ELSE
               CALL USIDS (IFLD, INTYP, CFIELD, IFIELD,
     &            NUMNPS, IDNPS, IDNSET(0,1), IDNSET(1,1),
     &            IDNSET(0,2), IDNSET(1,2), *170)
            END IF
         ELSE
            IF (ISFRON) THEN
               CALL USIDS (IFLD, INTYP, CFIELD, IFIELD,
     &            NUMESS, IDESS, IDESET(0,2), IDESET(1,2),
     &            IDESET(0,1), IDESET(1,1), *170)
            ELSE
               CALL USIDS (IFLD, INTYP, CFIELD, IFIELD,
     &            NUMESS, IDESS, IDESET(0,1), IDESET(1,1),
     &            IDESET(0,2), IDESET(1,2), *170)
            END IF
         END IF

      ELSE IF (VERB .EQ. 'SUMMARY') THEN
         PRINT 160
  160    FORMAT (
     &        /' TRANSLATE,    thickness'
     &        /' WARP, POINT,  thickness, radius'
     &        /' WARP, axis,   thickness, radius'
     $        /' WARP, ELLIPSE,thickness, major, minor'
     $        /' SPLINE,       thickness'
     $        /' THICKNESS,    thickness'
     &        /' MIRROR,       axis, ... '
     &        /' OFFSET,       xoffset, yoffset, zoffset'
     $        /' OFFSET,       axis, offset, ...'
     &        /' REVCEN,       xcenter, ycenter, zcenter'
     &        /' REVOLVE,      axis, number_of_degrees, ...'
     $        /' SCALE,        axis, scale, ...'
     &        /' SHIFT,        axis, shift, ...'
     $        /' ZERO,         axis, minimum, ...'
     &        /' NSETS,        FRONT_or_BACK, id1, id2, ...'
     &        /' SSETS,        FRONT_or_BACK, id1, id2, ...'
     $        /' CHANGE,       MATERIAL|NODESET|SIDESET old_id, new_id'
     &        /' ATTRIBUTE,    id which_attribute new_value'
     &        )
         VERB = ' '

      ELSE IF ((VERB .EQ. 'SHOW') .OR. (VERB .EQ. 'LIST')) THEN
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         IF (MATSTR (WORD, 'COMMANDS', 3)) THEN
            CALL SHOCMD ('COMMANDS', CMDTBL)
         ELSE IF (MATSTR (WORD, 'BLOCK', 3)) THEN
            CALL SHOW ('BLOCK', 'BLOCK', IDNPS, IDESS, IDNSET, IDESET,
     &         IDELB, NAMELB, NUMELB, NUMLNK, ELATTR)
         ELSE
            CALL ABRSTR (VERB, WORD, CMDTBL)
            CALL SHOW (VERB, WORD, IDNPS, IDESS, IDNSET, IDESET,
     &         IDELB, NAMELB, NUMELB, NUMLNK, ELATTR)
         END IF
         VERB = ' '

C ... ATTRIBUTE {id} {which } {value}
      ELSE IF (VERB .EQ. 'ATTRIBUT') THEN
         CALL FFINTG (IFLD, INTYP, IFIELD,
     *      'block ID', 0, ID, *170)
         CALL FFINTG (IFLD, INTYP, IFIELD,
     *      'which attribute', 0, IWHICH, *170)
         CALL FFREAL (IFLD, INTYP, RFIELD,
     *      'element attribute', 0.0, ATTRIB, *170)
         if (id .lt. 0) then
            call  prterr( 'ERROR', 'Invalid Block ID')
            go to 170
         end if
         if (iwhich .lt. 1 .or. iwhich .gt. 7) then
            call  prterr( 'ERROR',
     *       'Invalid Attribute Location. Range = 1..7')
            go to 170
         end if

         imat = locint (id, nelblk, idelb)
         IF (IMAT .EQ. 0) THEN
            CALL PRTERR ('ERROR', 'Invalid Block ID')
         else if (namelb(imat)(:3) .eq. 'BAR'
     &        .or. namelb(imat)(:4) .eq. 'BEAM'
     *        .or. namelb(imat)(:4) .eq. 'TRUS') then
            ELATTR(iwhich, imat) = ATTRIB
            write (*, 165) iwhich, id, attrib
 165        FORMAT(1x, 'Attribute number ',i5,' for block ',i5,
     *        ' set to ',1pe10.3)
         else
            CALL PRTERR ('ERROR', 'Block is not a beam, bar, or truss')
         end if
         verb = ' '

      ELSE IF (VERB .EQ. 'CHANGE') THEN
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         CALL FFINTG (IFLD, INTYP, IFIELD,
     *      'old ID', 0, IDOLD, *170)
         CALL FFINTG (IFLD, INTYP, IFIELD,
     *      'new ID', 0, IDNEW, *170)

         IF (MATSTR(WORD,'MATERIAL',1)) THEN
            CALL NEWID ('M', IDELB, NELBLK, IDNEW, IDOLD)
         ELSE IF (MATSTR(WORD,'NODESETS',1) .OR.
     $           MATSTR(WORD,'NSETS',1)) THEN
            CALL NEWID ('N', IDNPS, NUMNPS, IDNEW, IDOLD)
         ELSE IF (MATSTR(WORD,'SIDESETS',1) .OR.
     $           MATSTR(WORD,'SSETS',1)) THEN
            CALL NEWID ('S', IDESS, NUMESS, IDNEW, IDOLD)
         END IF
         VERB = ' '

      ELSE IF (VERB .EQ. 'HELP') THEN
         ISHELP = HELP ('GENSHELL', 'COMMANDS', CFIELD(IFLD))
         IF (.NOT. ISHELP) CALL SHOCMD ('COMMANDS', CMDTBL)
         VERB = ' '

      ELSE IF (VERB .EQ. 'SPAWN') THEN
          CALL PRTERR ('CMDSPEC', 'Spawn not implemented')

      ELSE IF ((VERB .EQ. 'END') .OR. (VERB .EQ. 'EXIT')) THEN
         CALL SCNEOF
         GOTO 220

      ELSE IF (VERB .EQ. 'QUIT') THEN
         CALL SCNEOF
         RETURN 1

      ELSE
         CALL PRTERR ('CMDERR', '"' // VERB(:LENSTR(VERB))
     &      // '" is an invalid command')
         VERB = ' '
      END IF

  170 CONTINUE

      IF (VERB .NE. ' ') THEN
         CALL SHOW (VERB, ' ', IDNPS, IDESS, IDNSET, IDESET,
     &      IDELB, NAMELB, NUMELB, NUMLNK, ELATTR)
      END IF

      GOTO 20

  220 CONTINUE
      IF (ITRANT .EQ. 0) ITRANT = 1

      RETURN
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20081211 08:57:08 nbsp Welcome to the Wiki Howdy Victoria, and Welcome to the Wiki!    My names Evan, pleased to meet you.  You dont need to worry about editing your own comments, so no need for the plea please allow change.  The only reason your last edit was partially reverted is because you removed the reviews of other people.  In general, reviews  both positive and negative  are retained.  You might want to read Welcome to the Wiki/Business Owner the informal, community written guide for business owners for some tips others have contributed.  Users/JabberWokky Evan JabberWokky Edwards

20081211 06:31:33 nbsp Thank you, Evan. I am new to all of this, so I appreciate the guidance. As long as I have your ear (or screen, as the case may be) I have a question about DavisWiki. I am the owner of Cartridge World Davis, but I did not make Cartridge Worlds listing on DavisWiki. My feeling is that if I chose to list my business on DavisWiki, then I am accepting that postings will be made about my business, negative or positive. If a business owner doesnt put their business out there, it seems that a page should be taken down at the owners request. I am new to this Wiki World, but it seems obvious that negative comments can be posted by any Tom, Dick or Harry without a business even deserving it. I understand that we can retaliate in kind, but the original comment is already out there. In this difficult economy, businesses are having a hard enough time making it. This other stuff seems to make it a hostile business climate...I thought I was avoiding that by choosing to open my business in the lovely town of Davis rather than Sacramento. Finally, the question...can I just have my Cartridge World page removed all together? Users/victoriarohl
   Well, its not mine to remove.  Its a creation of the community.  Its a bit like a newspaper article, only there are multiple contributors, all listed in the Info tab, and unlike newspaper article, everybody (including the owner and supporters) have an ability to respond.  Cartridge World is part of the community of Davis, so somebody wrote an entry here on the wiki about it (just like somebody wrote an entry about the Mayor, and another person wrote one about their church, and another person wrote up an entry about their kids softball team).  One key thing that might help is understanding that the wiki is not a business listing site  old businesses that closed a decade ago are written up in the same way as new ones (since they are part of the history of Davis), and there are entries about old murals, bike races and town events like local elections.  The Davis Wiki is a collection of things about Davis itself, and since Cartridge World is part of Davis, it has an entry that somebody in the community wrote and others added to.  That said, you are welcome to contribute: the wiki is open to the entire community of Davis, and that includes you.  The reason the negative reviews were restored was to give the same respect to those editors.  If they deleted your reply somebody would restore your words.  Once again, I encourage you to read Welcome to the Wiki/Business Owner the community written tips for business owners.  And of course, if you have any questions, there are a bunch of people who can help you  the entire wiki was created by volunteers all working together, and were pretty supportive of each other, especially new editors like you.  Users/JabberWokky Evan JabberWokky Edwards


text:
subroutine dimiwe(lunmd     ,lundia    ,error     ,nrrec     ,iweflg    , &
                & kmxdt     ,kmax      ,npiwe     ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: dimiwe.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/input/dimiwe.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: - Reads the dimension for Internal Wave Energy
!                dim. from an attribute file
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer , pointer :: itis
!
! Global variables
!
    integer               :: kmax   !  Description and declaration in esm_alloc_int.f90
    integer               :: kmxdt  !  Description and declaration in dimens.igs
    integer               :: lundia !  Description and declaration in inout.igs
    integer               :: lunmd  !  Description and declaration in inout.igs
    integer               :: npiwe  !  Description and declaration in dimens.igs
    integer               :: nrrec  !!  Record counter keeping the track of the last record read
    logical , intent(out) :: error  !!  Flag=TRUE if an error is encountered
    logical , intent(out) :: iweflg !  Description and declaration in procs.igs
!
! Local variables
!
    integer              :: ibeg
    integer              :: idef   ! Default value for NPIWE
    integer              :: iend
    integer              :: ier
    integer              :: iocond
    integer              :: lenc   ! Help var. (length of var. cvar to be looked for in the MD-file)
    integer              :: lfile  ! Number of non blank characters of file name
    integer              :: lkw
    integer              :: luntmp ! Unit number of (input par. for) FILIWE
    integer              :: newlun
    integer              :: nlook  ! Help var.: nr. of data to look for in the MD-file
    integer              :: ntrec  ! Help. var to keep track of NRREC
    logical              :: ex     ! Flag if file exists
    logical              :: found  ! Flag is true if KEYWRD is found
    logical              :: lerror ! Flag=TRUE if an error is encountered
    logical              :: newkw  ! Logical var. specifying whether a new recnam should be read from the MD-file or just new data in the continuation line
    character(12)        :: cdef   ! Default value when CHULP not found
    character(256)       :: filiwe ! File name for the IWE parameters
    character(300)       :: mdfrec ! Record read from the MD-file 300 = 256 + a bit (field, =, ##, etc.)
    character(6)         :: keyw   ! Name of record to look for in the MD-file (usually KEYWRD or RECNAM)
!
!! executable statements -------------------------------------------------------
!
    !
    itis  => gdp%gdrdpara%itis
    !
    lerror = .false.
    newkw  = .true.
    found  = .false.
    nlook  = 1
    mdfrec = ' '
    !
    ! locate 'Filiwe' containing parameter values for calculations of Internal Wave Energy;
    !
    keyw = 'Filiwe'
    ntrec = nrrec
    cdef = ' '
    lkw = 6
    call search(lunmd     ,lerror    ,newkw     ,nrrec     ,found     , &
              & ntrec     ,mdfrec    ,itis      ,keyw      ,lkw       , &
              & 'NO'      )
    lerror = .false.
    !
    ! found ?
    !
    if (found) then
       lenc = 12
       call read2c(lunmd     ,lerror    ,keyw      ,newkw     ,nlook     , &
                 & mdfrec    ,filiwe    ,cdef      ,lenc      ,nrrec     , &
                 & ntrec     ,lundia    ,gdp       )
       !
       ! reading error?
       !
       if (lerror) then
          lerror = .false.
          filiwe = cdef
       endif
       !
       ! File name defined (Y)
       !
       if (filiwe/=cdef) then
          !
          ! Define flag IWEFLG and recommended value for NPIWE
          !
          iweflg = .true.
          npiwe = 100
          !
          ! Test file existence and if open if found
          !
          lfile = len(filiwe)
          !
          inquire (file = filiwe(1:lfile), exist = ex)
          if (ex) then
             luntmp = newlun(gdp)
             open (luntmp, file = filiwe(:lfile), form = 'formatted',           &
                  & status = 'old')
             !
             ! Free formatted file
             ! Skip first record and then read dimensions from mdfrec
             ! for KMXDT default value not allowed incase IER < 0
             ! if NFREQS and/or NPIWE not found then defined for recommended
             ! value (100 according to Rob Uittenbogaard)
             !
             lenc = 300
             !
             read (luntmp, '(a)', iostat = iocond) mdfrec
             if (iocond/=0) then
                if (iocond<0) then
                   call prterr(lundia    ,'G006'    ,filiwe(1:lfile)      )
                else
                   call prterr(lundia    ,'G007'    ,filiwe(1:lfile)      )
                endif
                close (luntmp)
                error = .true.
                goto 9999
             endif
             !
             read (luntmp, '(a)', iostat = iocond) mdfrec
             if (iocond/=0) then
                if (iocond<0) then
                   call prterr(lundia    ,'G006'    ,filiwe(1:lfile)      )
                else
                   call prterr(lundia    ,'G007'    ,filiwe(1:lfile)      )
                endif
                close (luntmp)
                error = .true.
                goto 9999
             endif
             !
             ! Read KMAXTD from record
             !
             iend = 0
             ibeg = iend + 1
             call read1i(mdfrec    ,lenc      ,ibeg      ,iend      ,kmxdt     , &
                       & kmax      ,ier       )
             !
             ! End of record or no value defined (IER <= 0) not allowed
             !
             if (ier<=0) then
                error = .true.
                call prterr(lundia    ,'G006'    ,filiwe(1:lfile)      )
                !
                goto 9999
             endif
             !
             ! Read NPIWE from record
             !
             idef = npiwe
             ibeg = iend + 1
             call read1i(mdfrec    ,lenc      ,ibeg      ,iend      ,npiwe     , &
                       & idef      ,ier       )
             !
             ! End of record or no value defined (IER <= 0) allowed
             !
             if (ier<0) npiwe = idef
             !
             ! Close input file
             !
             close (luntmp)
          !
          ! File does not exist
          !
          else
             call prterr(lundia    ,'G004'    ,filiwe(1:lfile)      )
             !
             error = .true.
          endif
       endif
    endif
    !
 9999 continue
end subroutine dimiwe

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      MODULE snrc_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.3E  06:33:54  12/28/06
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
      SUBROUTINE snrc (IS, KAPS, KS, ND1, ND2, NE1, NE2, IBRD, IBRE)
      INTEGER, DIMENSION(4), INTENT(IN) :: IS
      INTEGER, DIMENSION(4), INTENT(IN) :: KAPS
      INTEGER, DIMENSION(4), INTENT(IN) :: KS
      INTEGER, INTENT(OUT) :: ND1
      INTEGER, INTENT(OUT) :: ND2
      INTEGER, INTENT(OUT) :: NE1
      INTEGER, INTENT(OUT) :: NE2
      INTEGER, INTENT(OUT) :: IBRD
      INTEGER, INTENT(OUT) :: IBRE
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
      program sigma2_to_sigma0
      implicit none
c
c     usage:  echo dens2 saln | sigma2_to_sigma0
c
c     input:  dens2 saln
c     output: saln pot.temp dens2 dens0
c
c             dens2  = Sigma-2 potential density
c             saln   = Salinity (psu)
c
c     alan j. wallcraft, naval research laboratory, may 2005.
c
      real*8  dens0,dens2,saln,temp
      integer ios
c
      write(6,'(a)') '#     saln  pot.temp     dens2     dens0'
      do
        read(5,*,iostat=ios) dens2,saln
        if     (ios.ne.0) then
          exit
        endif
        call tofsig2(temp,  dens2,saln)
        call  sigma0(dens0, temp, saln)
        write(6,'(4f10.3)') saln,temp,dens2,dens0
      enddo
      end
      SUBROUTINE TOFSIG2(TSEA, RSEA,SSEA)
      IMPLICIT NONE
C
      REAL*8  TSEA,SSEA,RSEA
C
C     CALCULATE POTENTIAL TEMPERATURE FROM SIGMA-2 DENSITY AND SALINITY.
C
C     STATEMENT FUNCTIONS.
C
      REAL*8     DZERO,DTHIRD
      PARAMETER (DZERO=0.D0,DTHIRD=1.D0/3.D0)
C
      REAL*8 C1,C2,C3,C4,C5,C6,C7
C --- coefficients for sigma-2 (based on Brydon & Sun fit)
      PARAMETER (C1= 9.77093e+00, C2=-2.26493e-02, C3= 7.89879E-01,
     .           C4=-6.43205E-03, C5=-2.62983E-03,
     .           C6= 2.75835E-05, C7= 3.15235E-05)
C
      REAL*8  R,S,T
      REAL*8  A0,A1,A2,CUBQ,CUBR,CUBAN,CUBRL,CUBIM,TOFSIG
C
C --- auxiliary statements for finding root of 3rd degree polynomial
      A0(S)=(C1+C3*S)/C6
      A1(S)=(C2+C5*S)/C6
      A2(S)=(C4+C7*S)/C6
      CUBQ(S)=DTHIRD*A1(S)-(DTHIRD*A2(S))**2
      CUBR(R,S)=DTHIRD*(0.5D0*A1(S)*A2(S)-1.5D0*(A0(S)-R/C6))
     .   -(DTHIRD*A2(S))**3
C --- if q**3+r**2>0, water is too dense to yield real root at given
C --- salinitiy. setting q**3+r**2=0 in that case is equivalent to
C --- lowering sigma until a double real root is obtained.
      CUBAN(R,S)=DTHIRD*ATAN2(SQRT(MAX(DZERO,
     .   -(CUBQ(S)**3+CUBR(R,S)**2))),CUBR(R,S))
      CUBRL(R,S)=SQRT(-CUBQ(S))*COS(CUBAN(R,S))
      CUBIM(R,S)=SQRT(-CUBQ(S))*SIN(CUBAN(R,S))
C
C --- temp (deg c) as a function of sigma and salinity (mil)
      TOFSIG(R,S)=-CUBRL(R,S)+SQRT(3.)*CUBIM(R,S)-DTHIRD*A2(S)
C
      TSEA = TOFSIG( RSEA,SSEA )
      RETURN
      END
      SUBROUTINE SIGMA0(RSEA,TSEA,SSEA)
      IMPLICIT NONE
C
      REAL*8  TSEA,SSEA,RSEA
C
C     CONVERT FROM SIGMA-2 TO SIGMA-0.
C
C     STATEMENT FUNCTIONS.
C
      REAL*8     DZERO,DTHIRD
      PARAMETER (DZERO=0.D0,DTHIRD=1.D0/3.D0)
C
      REAL*8     C1,C2,C3,C4,C5,C6,C7
C --- coefficients for sigma-0 (based on Brydon & Sun fit)
      DATA C1,C2,C3,C4,C5,C6,C7/
     . -1.36471E-01, 4.68181E-02, 8.07004E-01,-7.45353E-03,-2.94418E-03,
     .  3.43570E-05, 3.48658E-05/
C
      REAL*8  S,T
      REAL*8  SIG
C
C --- sigma-theta as a function of temp (deg c) and salinity (mil)
C --- (friedrich-levitus 3rd degree polynomial fit)
      SIG(T,S)=(C1+C3*S+T*(C2+C5*S+T*(C4+C7*S+C6*T)))
C
C     T AND S TO SIG-2,
C
      RSEA = SIG( TSEA, SSEA )
      RETURN
      END

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
C  Copyright (c) Mark J. Kilgard, 1994.

C  This program is freely distributable without licensing fees
C  and is provided without guarantee or warrantee expressed or
C  implied.  This program is -not- in the public domain.

C  GLUT Fortran program to draw red light sphere.

	subroutine display
#include "GL/fgl.h"
	call fglclear(GL_COLOR_BUFFER_BIT + GL_DEPTH_BUFFER_BIT)
	call fglcalllist(1)
	call glutswapbuffers
	end

	subroutine gfxinit
#include "GL/fgl.h"
	real diffuse(4),pos(4)
	data diffuse /1.0, 0.0, 0.0, 1.0/
	data pos /1.0, 1.0, 1.0, 0.0/
	call fglnewlist(1, GL_COMPILE)
	call glutsolidsphere(dble(1.0), 20, 20)
	call fglendlist
	call fgllightfv(GL_LIGHT0, GL_DIFFUSE, diffuse)
	call fgllightfv(GL_LIGHT0, GL_POSITION, pos)
	call fglenable(GL_LIGHTING)
	call fglenable(GL_LIGHT0)
	call fglenable(GL_DEPTH_TEST)
	call fglmatrixmode(GL_PROJECTION)
	call fgluperspective(dble(40.0), dble(1.0),
     2                       dble(1.0), dble(10.0))
	call fglmatrixmode(GL_MODELVIEW)
	call fglulookat(dble(0.0), dble(0.0), dble(5.0),
     2                  dble(0.0), dble(0.0), dble(0.0),
     3                  dble(0.0), dble(1.0), dble(1.0))
	call fgltranslatef(0.0, 0.0, -1.0)
	end

	program main
#include "GL/fglut.h"
	external display
	external reshape
	external submenu
	external mainmenu
	integer win
	call glutinit
	call glutinitdisplaymode(GLUT_DOUBLE+GLUT_RGB+GLUT_DEPTH)
	win =  glutcreatewindow('Fortran GLUT sphere')
	call gfxinit
	call glutdisplayfunc(display)
	call glutmainloop
	end


module propagate_cn2d

  include "./propagate_cn2d_src/moduse.src"

  implicit none

  private

  include "./propagate_cn2d_src/modvars.src"

  public :: propagate_cn2d_init
  public :: propagate_cn2d_cleanup
  public :: propagate_cn2d_splitop

contains
  subroutine propagate_cn2d_init()
    integer(ip) :: i_x, i_y

    dtau = dt

    include "./propagate_cn2d_src/init.src"

  end subroutine propagate_cn2d_init

  subroutine propagate_cn2d_cleanup()

    include "./propagate_cn2d_src/cleanup.src"

  end subroutine propagate_cn2d_cleanup

  subroutine propagate_cn2d_splitop(psi_arr, i_t)
    complex(dp), intent(inout) :: psi_arr(:,:)
    integer(ip), intent(in) :: i_t

    include "./propagate_cn2d_src/splitop.src"

  end subroutine propagate_cn2d_splitop

end module propagate_cn2d


text: Tim Coady was a former Chief Justice of the ASUCD Court. His term began in the Winter 2007 quarter at UC Davis, and ended Spring 2008.


text:
c******************************************************
      SUBROUTINE ESTRACEN(lidd,nlidd,step,e,nmm,denf,ff,xn,gis,ierr)
      implicit none
      include "tomo.h"
      include "line.h"
c TYPE definitions ++++++++++++++++++++++++++++++++++++++
      real*4    denf(NRAZ)
      real*4    ff(NRAZ+(NRAZA*NAZIPL-1)),gis(IAZIM,NRAZA)
      real*8    grid(4,3),w(4),e(3),step,xn(3),worazi,azms
      real*4    RBIMOD
      integer*4 lidd(NRAZ),nclo(4),ngr_out(4)
      integer*4 i,j,k,ierr,ncl,ncoin,nlidd,nmm,mo,ncl1
c END TYPE definitions ++++++++++++++++++++++++++++++++++
      ierr=0
      CALL E2PNTS(e,ngr_out,grid)
      CALL TRIAW(ngr_out,grid,e,w,ierr)
      do 20 j=1,4
      nclo(j)=ilocr(ngr_out(j))
      if(nclo(j).eq.0) then
      ierr=1
      return
      endif
   20 continue
      do 10 j=1,4
      if(j.eq.1) ncoin=nclo(1)*nmm+nclo(4)
      ncl=nclo(j)
      if(ncl.gt.nm) then
        write(*,*)'ERROR with local number',nm,ncl
        ierr=1
        return
      endif
      CALL SAZI(e,xn,azms)
      worazi=step*R/RBIMOD(e)
      do 30 i=1,nlidd
      if(lidd(i).eq.ncl) goto 40
   30 continue
      nlidd=nlidd+1
      if(nlidd.gt.NRAZ) then
      write(*,*) '(R0014) TOO SMALL SIZE OF lidd',nlidd-1
      STOP
      endif
      lidd(nlidd)=ncl
   40 do mo=1,1
      ncl1=nmm*(mo-1)+ncl
      if(mo.eq.1) ff(ncl1)=ff(ncl1)+w(j)*worazi
cxx   if(mo.eq.2) ff(ncl1)=ff(ncl1)-w(j)*worazi*DCOS(2.0d0*azms)
cxx   if(mo.eq.3) ff(ncl1)=ff(ncl1)-w(j)*worazi*DSIN(2.0d0*azms)
cxx   if(mo.eq.4) ff(ncl1)=ff(ncl1)-w(j)*worazi*DCOS(4.0d0*azms)
cxx   if(mo.eq.5) ff(ncl1)=ff(ncl1)-w(j)*worazi*DSIN(4.0d0*azms)
      enddo
      if(ncoin0.ne.ncoin) then
        denf(ncl)=1.0
          CALL AZI(j,grid,e,xn,azms)
          k=azms/const/(180.0/FLOAT(IAZIM))+1.0d0
          if(k.lt.1) k=1
          if(k.gt.IAZIM) k=IAZIM
        gis(k,ncl)=gis(k,ncl)+1.0
      endif
   10 continue
      ncoin0=ncoin
      return
      end
C******************************************************
      SUBROUTINE EINTEGRN(step,e,xn,nmm,sum,fsl,ierr)
      implicit none
      include "tomo.h"
      include "line.h"
c TYPE definitions ++++++++++++++++++++++++++++++++++++++
      integer*4 nmm,ierr
      real*4    fsl(nmm),sum
      real*8    grid(4,3),w(4),e(3),step,dsum,smult
      real*8    xn(3),scell
      real*4    RBIMOD
      integer*4 ngr_out(4),nclo,j,mo
c END TYPE definitions ++++++++++++++++++++++++++++++++++
      ierr=0
      CALL E2PNTS(e,ngr_out,grid)
      sum=0.0
      dsum=0.0d0
      smult=step*R/RBIMOD(e)
      do 20 mo=1,1
      do 10 j=1,4
      nclo=ilocr(ngr_out(j))
      if(nclo.lt.1) then
      ierr=1
      return
      endif
   10 w(j)=fsl(nmm*(mo-1)+nclo)
      CALL TRIAI(grid,e,w,scell)
      if(mo.eq.1) dsum=dsum+scell*smult
cxx   if(mo.eq.2) dsum=dsum-scell*smult*DCOS(2.0d0*azms)
cxx   if(mo.eq.3) dsum=dsum-scell*smult*DSIN(2.0d0*azms)
cxx   if(mo.eq.4) dsum=dsum-scell*smult*DCOS(4.0d0*azms)
cxx   if(mo.eq.5) dsum=dsum-scell*smult*DSIN(4.0d0*azms)
   20 continue
      sum=dsum
      return
      end

max_stars_repo_name: NoiseCIEI/RayTomo

max_issues_repo_name: NoiseCIEI/RayTomo

max_forks_repo_name: NoiseCIEI/RayTomo