text:
      subroutine compute_entropy(s)
! computes entropy at istep and pushes the stack down for previous
! steps needed to compute ds/dt via finite difference (for now).
! hardcoded for Burgers equation. More later when folded into CMT-nek
! for Burgers, s=energy=1/2 U^2
      include 'SIZE'
      include 'TOTAL'  ! tlag is lurking. be careful
      include 'CMTDATA'
! I've always seen lorder=3, but I still need three steps
!          s(:,               1,       1)  ! entropy at current step
!          s(:,               2,       1)  ! entropy at step n-1
!          s(:,               1,       2)  ! entropy at step n-2
      real s(lx1*ly1*lz1*lelt,lorder-1,*)
      real ntol
      integer e
      data icalld /0/
      save icalld

      n=lx1*ly1*lz1
      ntot=n*nelt
      ntol=1.0e-10

      if (icalld .eq. 0) then
         if (nio .eq. 0) write(6,*) 'zeroing out entropy stack',istep
         icalld=1
         call rzero(s,ntot)
         call rzero(s(1,1,2),ntot) ! s_{n-1}
         call rzero(s(1,2,1),ntot) ! s_n
      endif

! compute the current entropy. This actually needs to go back in the
! usr file because it's EOS-dependent
      rgam=rgasref/(gmaref-1.0)
      do i=1,ntot
         rho=max(vtrans(i,1,1,1,irho),ntol)
         s(i,1,1)=rgam*rho*log(pr(i,1,1,1)/(rho**gmaref))
      enddo

      if (stage .eq. 1) then
! push the stack
         call copy(s(1,1,2),s(1,2,1),ntot) ! s_{n-1}=s_n
         call copy(s(1,2,1),s(1,1,1),ntot) ! fill s_n
      endif

      return
      end

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

      subroutine entropy_viscosity
      include 'SIZE'
      include 'TOTAL'
      include 'CMTDATA'
      parameter (lxyz=lx1*ly1*lz1)
      common /scrns/ scrent(lxyz,lelt)
      integer e
      character*132 deathmessage

      pi=4.0*atan(1.0)

      n=lx1*ly1*lz1
      ntot=n*nelt

! entropy at this and lorder prior steps
      call compute_entropy(tlag)
! compute maxval(|S-<S>|)
      savg    =    glsc2(tlag,bm1,ntot)
      savg    = -savg/volvm1
      call cadd2(scrent,tlag,savg,ntot)
      maxdiff =     glamax(scrent,ntot)
      if (maxdiff.le.0.0) then
         write(6,*) 'zero maxdiff usually means NAN$'
!        write(deathmessage,*) 'zero maxdiff usually means NAN$'
!        call exittr(deathmessage,maxdiff,istep)
      endif
      call entropy_residual(tlag) ! fill res2
      call copy(res2(1,1,1,1,2),res2,ntot) ! raw residual in res2
      call wavevisc(t(1,1,1,1,3))
      call resvisc(res2) ! overwrite res2
      call evmsmooth(res2,t(1,1,1,1,3),.true.) ! endpoints=.false.
                                               ! is intended to
                                               ! preserve face states,
                                               ! but this is easier to
                                               ! test 1D
!     call evmsmooth(res2,t(1,1,1,1,3),.true.) ! And again.
      call dsavg(res2) ! you DEFINITELY don't want a min here

      return
      end

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

      subroutine piecewiseAV(shock_detector)
      include 'SIZE'
      include 'TOTAL'
      include 'CMTDATA'
      integer e
      character*132 deathmessage
      common /scrvh/ avmask(lelt)
      real avmask
      external shock_detector

      pi=4.0*atan(1.0)

      nxyz=lx1*ly1*lz1
      ntot=nxyz*nelt

! toggle shock detector with AV application Lv, See & Ihme (2016) JCP 322

      if (time4av) then
         call shock_detector(avmask)

! old nu_max
         call wavevisc(t(1,1,1,1,3))
! diagnostic
!        if (stage.gt.0)then
!        do e=1,nelt
!!           do i=1,nxyz
!              write(stage*100+nid,*)xm1(i,1,1,e),ym1(i,1,1,e),t(i,1,1,e,3)
!           enddo
!        enddo
!        endif
         do e=1,nelt
            call cmult(t(1,1,1,e,3),avmask(e),nxyz)
            if (avmask(e).ne.1.0) write(6,*) 'duh sir'
         enddo

      else
         call rzero(t(1,1,1,1,3),ntot)
!        call shock_detector(avmask)
      endif

! diagnostic
!      if (stage.gt.0)then
!      do e=1,nelt
!         do i=1,nxyz
!            write(stage*1000+nid,'(5e15.7)') xm1(i,1,1,e),ym1(i,1,1,e),
!!     >        t(i,1,1,e,3),avmask(e),epsebdg(e)
!         enddo
!      enddo
!      endif

!     call max_to_trilin(t(1,1,1,1,3))

      return
      end

c-----------------------------------------------------------------------

      subroutine entropy_residual(s)
! COMPUTE R=ds/dt + div F(s) for entropy pair s and (hardcoded) F
! and store its norm functional thing |R| in res2 where it will
! provide artificial viscosity according to the code in
! entropy_viscosity and the method of Guermond
      include 'SIZE'
      include 'TOTAL' ! tlag lurks
      include 'CMTDATA'
      integer e
      real dsdtcoef(3) ! put this in /TIMESTEPCOEF/ someday
      data dsdtcoef /1.0,1.0,0.5/
      real s(lx1*ly1*lz1*lelt,lorder-1,ldimt) ! because it's really tlag

      n=lx1*ly1*lz1
      ntot=n*nelt

      if (istep .eq. 1) return
      rdt=1.0/(dsdtcoef(stage)*dt)
      if (stage .eq. 1) then ! THE MOST SABOTAGEABLE PART OF THE CODE
         call sub3(res2,s(1,1,1),s(1,1,2),ntot) ! EVALUATE s_n-s_{n-1}
      else
         call sub3(res2,s(1,1,1),s(1,2,1),ntot) ! EVALUATE s^{(stage)}-s_n
      endif
      call cmult(res2,rdt,ntot)
! res2=ds/dt. now,

!-----------------------------------------------------------------------
! cons approach: strong-conservation form of flux divergence in entropy residual
!-----------------------------------------------------------------------
! get around to expanding totalh to store fluxes for whole fields and
! properly vectorize evaluate_*_h and flux_div_integral
      do e=1,nelt
         call evaluate_entropy_flux(e) ! diffh. zero it before diffusion
         call flux_div_mini(e) ! into res2 it goes.
      enddo

      return
      end

c-----------------------------------------------------------------------

      subroutine evaluate_entropy_flux(e)
! entropy flux function for entropy residual.
! just vel*s for now
      include 'SIZE'
      include 'SOLN'
      include 'INPUT'
      include 'CMTDATA'
      integer e

      call rzero(totalh,3*lxd*lyd*lzd)
      n=lx1*ly1*lz1

      call col3(totalh(1,1),vx(1,1,1,e),tlag(1,1,1,e,1,1),n)
      call col3(totalh(1,2),vy(1,1,1,e),tlag(1,1,1,e,1,1),n)
      if (if3d) call col3(totalh(1,3),vz(1,1,1,e),tlag(1,1,1,e,1,1),n)

      return
      end

c-----------------------------------------------------------------------

      subroutine flux_div_mini(e)
      include 'SIZE'
      include 'INPUT'
      include 'GEOM'
      include 'DXYZ'
      include 'SOLN'
      include 'CMTDATA'
      parameter (ldd=lx1*ly1*lz1)
      parameter (ldg=lx1**3,lwkd=2*ldg)
      common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ju(ldd),ud(ldd),tu(ldd)
      real ju

      integer e

      nrstd=ldd
      nxyz=lx1*ly1*lz1
      mdm1=lx1-1

      call rzero(ur,nrstd)
      call rzero(us,nrstd)
      call rzero(ut,nrstd)
      call rzero(ud,nrstd)
      call rzero(tu,nrstd)

      if (if3d) then
         call local_grad3(ur,us,ut,totalh(1,1),mdm1,1,dxm1,dxtm1)
         do i=1,nxyz
            ud(i) = jacmi(i,e)*(rxm1(i,1,1,e)*ur(i)
     $                        + sxm1(i,1,1,e)*us(i)
     $                        + txm1(i,1,1,e)*ut(i))
         enddo
         call local_grad3(ur,us,ut,totalh(1,2),mdm1,1,dxm1,dxtm1)
         do i=1,nxyz ! confirmed to have no effect in 1D
            ud(i)=ud(i)+jacmi(i,e)*(rym1(i,1,1,e)*ur(i)
     $                            + sym1(i,1,1,e)*us(i)
     $                            + tym1(i,1,1,e)*ut(i))
         enddo
         call local_grad3(ur,us,ut,totalh(1,3),mdm1,1,dxm1,dxtm1)
         do i=1,nxyz ! confirmed to have no effect in 1D
            ud(i)=ud(i)+jacmi(i,e)*(rzm1(i,1,1,e)*ur(i)
     $                            + szm1(i,1,1,e)*us(i)
     $                            + tzm1(i,1,1,e)*ut(i))
         enddo
      else
         call local_grad2(ur,us,totalh(1,1),mdm1,1,dxm1,dxtm1)
         do i=1,nxyz
            ud(i) = jacmi(i,e)*(rxm1(i,1,1,e)*ur(i)
     $                        + sxm1(i,1,1,e)*us(i))
         enddo
         call local_grad2(ur,us,totalh(1,2),mdm1,1,dxm1,dxtm1)
         do i=1,nxyz
            ud(i)=ud(i)+jacmi(i,e)*(rym1(i,1,1,e)*ur(i)
     $                            + sym1(i,1,1,e)*us(i))
         enddo
      endif
      call add2(res2(1,1,1,e,1),ud,nxyz)

      return
      end

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

      subroutine resvisc(residual)
! smooth residual-based entropy visc, defined by Guermond, Popov, whoever
! DSAVG assumes IFFLOW to work correctly. IFHEAT still doesn't work with CMT-nek
      include 'SIZE'
      include 'TOTAL'
      include 'CMTDATA'
      include 'DG'
      integer lfq,heresize,hdsize
      parameter (lxyz=lx1*ly1*lz1)
      real residual(lxyz,nelt)
      integer e,f

      nxyz =lx1*ly1*lz1
      nxz  =lx1*lz1
      nface=2*ldim
      nxzf =nxz*nface
      nfq  =nxzf*nelt

! ensure continuity at faces. doing this before |abs| causes some cancellation that,
! so far, appears to beneficially reduce spikiness at faces.

      do e=1,nelt
         do i=1,nxyz
            residual(i,e)=residual(i,e)*meshh(e)**2
         enddo
      enddo

      call dsavg(residual) ! signed, can cancel at faces. Hope it does

      do e=1,nelt
         do i=1,nxyz
            residual(i,e)=abs(residual(i,e))
         enddo
      enddo

      call cmult(residual,c_sub_e,nxyz*nelt)

      if (maxdiff .ne. 0) then
         const=1.0/maxdiff
         call cmult(residual,const,nxyz*nelt)
      endif

      return
      end

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

      subroutine evmsmooth(resvisc,wavevisc,endpoints)
      include 'SIZE'
      include 'INPUT'
      real resvisc(lx1,ly1,lz1,nelt),wavevisc(lx1,ly1,lz1,nelt)
      real rtmp
      common /ctmp1/ rtmp(lx1,ly1,lz1)
! are faces included in smoothing? if not (say they're fixed by dsavg) then
! endpoints should be .false.
      logical endpoints
      integer e
! clip residual viscosity and smooth it.
! actually just smooth first and then clip it. Just average-with-my-neighbors for now
! weight the neighbors according to GLL-spacing instead of uniform-grid
! at a later date.

      nxyz=lx1*ly1*lz1
      kstart=1
      kend=lz1
      rldim=1.0/ldim

      if (endpoints) then
         istart=1
         jstart=1
         iend=lx1
         jend=ly1
      else
         istart=2
         jstart=2
         iend=lx1-1
         jend=ly1-1
         if (if3d) then
            kstart=2
            kend=lz1-1
         endif
      endif

! yes I know this loop sucks and I need to write a matrix or something
      do e=1,nelt
         do i=1,nxyz
            if (wavevisc(i,1,1,e).le. resvisc(i,1,1,e))
     >          resvisc(i,1,1,e)= wavevisc(i,1,1,e)
         enddo
         call copy(rtmp,resvisc(1,1,1,e),nxyz) ! really only for .false.
         do iz=kstart,kend
            if (if3d) then
               km1=iz-1
               kp1=iz+1
               izm=km1
!              if (km1 .lt. 1) izm=iz ! bias towards {{face point}}
               if (km1 .lt. 1) izm=kp1 ! Guermond symmetry
               izp=kp1
!              if (kp1 .gt. lz1) izp=iz ! bias towards {{face point}}
               if (kp1 .gt. lz1) izp=km1 ! Guermond symmetry
            else
               izm=iz
               izp=iz
            endif
            do iy=jstart,jend
               jm1=iy-1
               jp1=iy+1
               iym=jm1
!              if (jm1 .lt. 1) iym=iy ! bias towards {{face point}}
               if (jm1 .lt. 1) iym=jp1 ! Guermond symmetry
               iyp=jp1
!              if (jp1 .gt. ly1) iyp=iy ! bias toward {{face point}}
               if (jp1 .gt. ly1) iyp=jm1 ! Guermond symmetry
               do ix=istart,iend
                  im1=ix-1
                  ip1=ix+1
                  ixm=im1
!                 if (im1 .lt. 1) ixm=ix ! bias towards {{face point}}
                  if (im1 .lt. 1) ixm=ip1 ! Guermond symmetry
                  ixp=ip1
!                 if (ip1 .gt. lx1) ixp=ix ! bias towards {{face point}}
                  if (ip1 .gt. lx1) ixp=im1 ! Guermond symmetry
                  x0 = resvisc(ix ,iy ,iz ,e)
                  x1 = resvisc(ixm,iy ,iz ,e)
                  x2 = resvisc(ixp,iy ,iz ,e)
                  x3 = resvisc(ix ,iym,iz ,e)
                  x4 = resvisc(ix ,iyp,iz ,e)
                  if (if3d) then
                     x5 = resvisc(ix ,iy ,izm,e)
                     x6 = resvisc(ix ,iy ,izp,e)
                  else
                     x5=0.0
                     x6=0.0
                  endif
                  rtmp(ix,iy,iz)=0.25*(2.0*ldim*x0+x1+x2+x3+x4+x5+x6)
     >                               *rldim
               enddo
            enddo
         enddo
         call copy(resvisc(1,1,1,e),rtmp,nxyz)
      enddo

      return
      end

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

      subroutine wavevisc(numax)
! max entropy visc, defined by Guermond, Popov, whoever (chapter-verse)
! as
! numax = c_max*h*max(dH/dU)
! in a given element
      include 'SIZE'
      include 'TOTAL'
      include 'CMTDATA'
      parameter (lxyz=lx1*ly1*lz1)
      common /scrns/ wavespeed(lxyz)
      real wavespeed
      real maxeig
      real numax(lxyz,nelt)
      integer e

      nxyz=lx1*ly1*lz1

      do e=1,nelt
         do i=1,nxyz
            wavespeed(i)=csound(i,1,1,e)+
     >      sqrt(vx(i,1,1,e)**2+vy(i,1,1,e)**2+vz(i,1,1,e)**2)
         enddo
         maxeig=vlamax(wavespeed,nxyz)
! Zingan (2015) only. not long for this world
!        rhomax(e)=vlamax(vtrans(1,1,1,e,irho),nxyz)
         do i=1,nxyz
            numax(i,e)=c_max*maxeig*meshh(e)
         enddo
      enddo

      call max_to_trilin(t(1,1,1,1,3))

      return
      end

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

      subroutine max_to_trilin(field)
! stupid subroutine to take a stupid uniform field and compute a trilinear
! tent between maximum shared values at the vertices.
      include 'SIZE'
      include 'TOTAL'
      real field(lx1,ly1,lz1,nelt)
      integer e

      nxyz=lx1*ly1*lz1

! get maxima on faces
      call dsop(field,'MAX',lx1,ly1,lz1)

! trilinear interpolation. you should adapt xyzlin to your needs instead
      do e=1,nelt
         p000=field(1,  1,  1,  e)
         p100=field(lx1,1,  1,  e)
         p010=field(1,  ly1,1,  e)
         p110=field(lx1,ly1,1,  e)
         p001=field(1,  1,  lz1,e)
         p101=field(lx1,1,  lz1,e)
         p011=field(1,  ly1,lz1,e)
         p111=field(lx1,ly1,lz1,e)
         c1=p100-p000
         c2=p010-p000
         c3=p001-p000
         c4=p110-p010-p100+p000
         c5=p011-p001-p010+p000
         c6=p101-p001-p100+p000
         c7=p111-p011-p101-p110+p100+p001+p010-p000
         rdx=1.0/(xm1(lx1,1,1,e)-xm1(1,1,1,e)) ! cubes only!!!
         rdy=1.0/(ym1(1,ly1,1,e)-ym1(1,1,1,e))
         rdz=0.0
         if(if3d) rdz=1.0/(zm1(1,1,lz1,e)-zm1(1,1,1,e))
         do i=1,nxyz
            deltax=rdx*(xm1(i,1,1,e)-xm1(1,1,1,e)) ! cubes only!!!
            deltay=rdy*(ym1(i,1,1,e)-ym1(1,1,1,e))
            deltaz=0.0
            if (if3d) deltaz=rdz*(zm1(i,1,1,e)-zm1(1,1,1,e))
            field(i,1,1,e)=p000+c1*deltax+c2*deltay+c3*deltaz+
     >                          c4*deltax*deltay+c5*deltay*deltaz+
     >                       c6*deltaz*deltax+c7*deltay*deltaz*deltax
         enddo
      enddo

      return
      end

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

      subroutine perssonperaire(shkdet,var,shtmp)
! Peraire & Persson (2006) Eq. 7
      include 'SIZE'
      include 'TOTAL'
      include 'CMTDATA'
      parameter (lxyz=lx1*ly1*lz1)
      common /scrns/ fvar (lxyz,lelt)
!    $ ,             ytm1 (lx1,ly1,lz1,lelv)
      real shkdet(lxyz,nelt),var(lxyz,*),shtmp(lxyz,*)
      integer e
! this whole common block is atrocious
      common /CMTFILTERS/  intv(lx1,lx1),intt(lx1,lx1)
     $                  ,  intvd(lxd,lxd),inttd(lxd,lxd)
     $                  ,  wk1(lx1,lx1,lz1),wk2(lx1,lx1,lz1)
     $                  ,  wkd1(lxd,lxd,lzd),wkd2(lxd,lxd,lzd)
      real                 intv, intt, intvd, inttd
      real kappa
      integer icalld
      save    icalld
      data    icalld/0/

      if (icalld.eq.0) then
         ncut=-110
         icalld=1
         call userfilt(intv,zgm1,lx1,nio)
      endif

      n=lx1*ly1*lz1
      ntot=n*nelt

! JH082418 Hardcoded threshold
      cpp=0.001

      call copy(fvar,var,ntot)
      call filterq(fvar,intv,lx1,lz1,wk1,wk2,intt,if3d,dmax)

      do e=1,nelt
         call sub3(shkdet(1,e),var(1,e),fvar(1,e),n) ! store u-u^ in shkdet
         shock_or_not=vlsc3(shkdet(1,e),shkdet(1,e),bm1(1,1,1,e),n)
         denom=vlsc3(var(1,e),var(1,e),bm1(1,1,1,e),n)
         if (denom .gt. 0.0) then ! mask this somehow. ifs are bad
            shock_or_not=(shock_or_not/denom)*(lx1-1)**4
!           shock_or_not=log10(shock_or_not)
            if (shock_or_not .gt. cpp) then
!              call cfill(shkdet(1,e),shock_or_not,n)
               call rone(shkdet(1,e),n)
            else
               call rzero(shkdet(1,e),n)
            endif
         else
            write(6,*) 'nid,iel=',nid,e
            write(6,*) 'variables that change sign are bad for shox'
            call exitt
         endif
      enddo

      return
      end

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
program ch1206
    ! Easter
    implicit none

    integer ::  year, metcyc, century, error1, &
        error2, day
    integer :: epact, luna, temp

    print *, ' input the year for which easter'
    print *, ' is to be calculated'
    print *, ' enter the whole year, e.g. 1978 '
    read *, year

    metcyc = mod(year,19)+1
    if (year<=1582) then
        day = (5*year)/4
        epact = mod(11*metcyc-4,30)+1
    else
        century = (year/100)+1
        error1 = (3*century/4)-12
        error2 = ((8*century+5)/25)-5
        day = (5*year/4)-error1-10
        temp = 11*metcyc+20+error2-error1
        epact = mod(temp,30)
        if (epact<=0) then
            epact = 30+epact
        end if
        if ((epact==25 .and. metcyc>11) .or. &
            epact==24) then
            epact = epact+1
        end if
    end if
    luna = 44-epact
    if (luna<21) then
        luna = luna+30
    end if
    luna = luna+7-(mod(day+luna,7))
    if (luna>31) then
        luna = luna-31
        print *, ' for the year ', year
        print *, ' easter falls on April ', luna
    else
        print *, ' for the year ', year
        print *, ' easter falls on march ', luna
    end if
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
c     program DRDRANG
c>> 2001-05-22 DRDRANG  Krogh Minor change for making .f90 version.
c>> 1996-05-28 DRDRANG  Krogh Added external statement.
c>> 1994-10-19 DRDRANG  Krogh  Changes to use M77CON
c>> 1987-12-09 DRDRANG  Lawson  Initial Code.
c--D replaces "?": DR?RANG, ?RANG, ?STAT1, ?STAT2
c     Demonstration driver for DRANG.  Generates random numbers
c     from the Gaussian distribution.
c     ------------------------------------------------------------------
      integer NCELLS
      parameter(NCELLS = 12+2)
      external DRANG
      double precision     DRANG, STATS(5), Y1, Y2, YTAB(1), ZERO
      integer  I, IHIST(NCELLS), N
c
      parameter(ZERO = 0.0D0)
      data       N / 10000/
      data  Y1, Y2 / -3.0D0, +3.0D0/
c     ------------------------------------------------------------------
      STATS(1) = ZERO
      do 20 I=1,N
c                          GET RANDOM NUMBER
         YTAB(1) = DRANG()
c
c                          Accumulate statistics and histogram.
c
        call DSTAT1(YTAB(1), 1, STATS, IHIST, NCELLS, Y1, Y2)
   20 continue
c                          Print the statistics and histogram.
c
      write(*,'(13x,a//)') 'Gaussian random numbers from DRANG'
      call DSTAT2(STATS, IHIST, NCELLS, Y1, Y2)
      stop
      end

ext: f


text:
real function gaussian()
  implicit none
  !	real rd

  real x1,x2,w,y1
  real y2
  real n1,n2
  integer use_last
  integer ii,jj

  use_last=0
  y2=0
  w=2.0
  if(use_last.ne.0) then
    y1=y2
    use_last=0
  else
    do while (w.ge.1.0)
      call random_number(n1)
      call random_number(n2)
      x1=2.0*n1-1.0
      x2=2.0*n2-1.0
      w = x1 * x1 + x2 * x2
    enddo
    w=((-2.0*log(w))/w)**0.5
    y1=x1*w
    y2=x2*w
    use_last=1
  endif
  gaussian=y1
end function

max_stars_repo_name: ChetGoerzen/DSurfTomo

max_issues_repo_name: ChetGoerzen/DSurfTomo

max_forks_repo_name: ChetGoerzen/DSurfTomo


text:
       SUBROUTINE MA_CGID ( cgrpt, mszrpt, ipt, iret )
C************************************************************************
C* MA_CGID                                                              *
C*                                                                      *
C* This subroutine decodes the Coast Guard station id in one report.    *
C* The value is stored in common /cintfv/.				*
C*                                                                      *
C* MA_CGID  ( CGRPT, MSZRPT, IPT, IRET )                                *
C*                                                                      *
C* Input parameters:                                                    *
C*      CGRPT           CHAR*           Report array                    *
C*      MSZRPT          INTEGER         Length of report in bytes       *
C*                                                                      *
C* Input and Output parameters:                                         *
C*      IPT             INTEGER         Pointer to start of field.      *
C*                                                                      *
C* Output parameters:                                                   *
C*      CIVALS(ICSTID)  CHAR*           Station id                      *
C*      IRET            INTEGER         Return code                     *
C*                                        0 = Normal return             *
C*                                       -1 = ID is missing             *
C*                                       -2 = bad format                *
C**                                                                     *
C* Log:                                                                 *
C* C. Caruso Magee/NCEP  4/01	Original Author                         *
C* F. J. Yen/NCEP	 4/01	Cleaned up, reformatted, and renamed	*
C*				from CG_STID. Changed prologue.		*
C************************************************************************
       	INCLUDE       'macmn.cmn'
C*
       	CHARACTER*(*) cgrpt
C------------------------------------------------------------------------
       	iret = 0
C
C*  	First 3 characters of cgrpt will be id.
C*  	If ID consists of ???, set STID to missing.
C*  	IF 1st, 2nd, or 3rd char is blank, leave stid as missing (blank).
C
       	IF ( cgrpt(1:3) .ne. '???' ) THEN
            IF ( INDEX ( cgrpt(1:3), ' ' ) .eq. 0 ) THEN
        	civals(icstid)(1:1) = 'K'
        	civals(icstid)(2:4) = cgrpt(1:3)
              ELSE
                iret = -2
            END IF
          ELSE
            iret = -1
        END IF

        ipt = ipt + 4
C*
        RETURN
        END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Subroutine pyinki(chfram, chbeam, chtarg, win)
  Common /lujets/n, k(9000, 5), p(9000, 5), v(9000, 5)
  Save /lujets/
  Common /ludat1/mstu(200), paru(200), mstj(200), parj(200)
  Save /ludat1/
  Common /pysubs/msel, msub(200), kfin(2, -40:40), ckin(200)
  Save /pysubs/
  Common /pypars/mstp(200), parp(200), msti(200), pari(200)
  Save /pypars/
  Common /pyint1/mint(400), vint(400)
  Save /pyint1/
  Character chfram*8, chbeam*8, chtarg*8, chcom(3)*8, chalp(2)*26, chidnt(3)*8, chtemp*8, chcde(18)*8, chinit*76
  Dimension len(3), kcde(18)
  Data chalp/'abcdefghijklmnopqrstuvwxyz', 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
  Data chcde/'e-      ', 'e+      ', 'nue     ', 'nue~    ', 'mu-     ', 'mu+     ', 'numu    ', 'numu~   ', 'tau-    ', 'tau+    ', 'nutau   ', 'nutau~  ', 'pi+     ', 'pi-     ', 'n       ', 'n~      ', 'p       ', 'p~      '/
  Data kcde/11, -11, 12, -12, 13, -13, 14, -14, 15, -15, 16, -16, 211, -211, 2112, -2112, 2212, -2212/
  chcom(1) = chfram
  chcom(2) = chbeam
  chcom(3) = chtarg
  Do i = 1, 3
    len(i) = 8
    Do ll = 8, 1, -1
      If (len(i)==ll .And. chcom(i)(ll:ll)==' ') len(i) = ll - 1
      Do la = 1, 26
        If (chcom(i)(ll:ll)==chalp(2)(la:la)) chcom(i)(ll:ll) = chalp(1)(la:la)
      End Do
    End Do
    chidnt(i) = chcom(i)
    Do ll = 1, 6
      If (chidnt(i)(ll:ll+2)=='bar') Then
        chtemp = chidnt(i)
        chidnt(i) = chtemp(1:ll-1) // '~' // chtemp(ll+3:8) // '  '
      End If
    End Do
    Do ll = 1, 8
      If (chidnt(i)(ll:ll)=='_') Then
        chtemp = chidnt(i)
        chidnt(i) = chtemp(1:ll-1) // chtemp(ll+1:8) // ' '
      End If
    End Do
  End Do
  n = 2
  Do i = 1, 2
    k(i, 2) = 0
    Do j = 1, 18
      If (chidnt(i+1)==chcde(j)) k(i, 2) = kcde(j)
    End Do
    p(i, 5) = ulmass(k(i,2))
    mint(40+i) = 1
    If (iabs(k(i,2))>100) mint(40+i) = 2
    Do j = 1, 5
      v(i, j) = 0.
    End Do
  End Do
  If (k(1,2)==0) Write (mstu(11), 1000) chbeam(1:len(2))
  If (k(2,2)==0) Write (mstu(11), 1100) chtarg(1:len(3))
  If (k(1,2)==0 .Or. k(2,2)==0) Stop
  Do j = 6, 10
    vint(j) = 0.
  End Do
  chinit = ' '
  If (chcom(1)(1:2)=='cm') Then
    If (chcom(2)(1:1)/='e') Then
      loffs = (34-(len(2)+len(3)))/2
      chinit(loffs+1:76) = 'PYTHIA will be initialized for a ' // chcom(2)(1:len(2)) // '-' // chcom(3)(1:len(3)) // ' collider' // ' '
    Else
      loffs = (33-(len(2)+len(3)))/2
      chinit(loffs+1:76) = 'PYTHIA will be initialized for an ' // chcom(2)(1:len(2)) // '-' // chcom(3)(1:len(3)) // ' collider' // ' '
    End If
    s = win**2
    p(1, 1) = 0.
    p(1, 2) = 0.
    p(2, 1) = 0.
    p(2, 2) = 0.
    p(1, 3) = sqrt(((s-p(1,5)**2-p(2,5)**2)**2-(2.*p(1,5)*p(2,5))**2)/(4.*s))
    p(2, 3) = -p(1, 3)
    p(1, 4) = sqrt(p(1,3)**2+p(1,5)**2)
    p(2, 4) = sqrt(p(2,3)**2+p(2,5)**2)
  Else If (chcom(1)(1:3)=='fix') Then
    loffs = (29-(len(2)+len(3)))/2
    chinit(loffs+1:76) = 'PYTHIA will be initialized for ' // chcom(2)(1:len(2)) // ' on ' // chcom(3)(1:len(3)) // ' fixed target' // ' '
    p(1, 1) = 0.
    p(1, 2) = 0.
    p(2, 1) = 0.
    p(2, 2) = 0.
    p(1, 3) = win
    p(1, 4) = sqrt(p(1,3)**2+p(1,5)**2)
    p(2, 3) = 0.
    p(2, 4) = p(2, 5)
    s = p(1, 5)**2 + p(2, 5)**2 + 2.*p(2, 4)*p(1, 4)
    vint(10) = p(1, 3)/(p(1,4)+p(2,4))
    Call lurobo(0., 0., 0., 0., -vint(10))
  Else If (chcom(1)(1:3)=='use') Then
    loffs = (13-(len(1)+len(2)))/2
    chinit(loffs+1:76) = 'PYTHIA will be initialized for ' // chcom(2)(1:len(2)) // ' on ' // chcom(3)(1:len(3)) // 'user-specified configuration' // ' '
    p(1, 4) = sqrt(p(1,1)**2+p(1,2)**2+p(1,3)**2+p(1,5)**2)
    p(2, 4) = sqrt(p(2,1)**2+p(2,2)**2+p(2,3)**2+p(2,5)**2)
    Do j = 1, 3
      vint(7+j) = sngl((dble(p(1,j))+dble(p(2,j)))/dble(p(1,4)+p(2,4)))
    End Do
    Call lurobo(0., 0., -vint(8), -vint(9), -vint(10))
    vint(7) = ulangl(p(1,1), p(1,2))
    Call lurobo(0., -vint(7), 0., 0., 0.)
    vint(6) = ulangl(p(1,3), p(1,1))
    Call lurobo(-vint(6), 0., 0., 0., 0.)
    s = p(1, 5)**2 + p(2, 5)**2 + 2.*(p(1,4)*p(2,4)-p(1,3)*p(2,3))
  Else
    Write (mstu(11), 1800) chfram(1:len(1))
    Stop
  End If
  If (s<parp(2)**2) Then
    Write (mstu(11), 1900) sqrt(s)
    Stop
  End If
  mint(11) = k(1, 2)
  mint(12) = k(2, 2)
  mint(43) = 2*mint(41) + mint(42) - 2
  vint(1) = sqrt(s)
  vint(2) = s
  vint(3) = p(1, 5)
  vint(4) = p(2, 5)
  vint(5) = p(1, 3)
  If (mstp(82)<=1) vint(149) = 4.*parp(81)**2/s
  If (mstp(82)>=2) vint(149) = 4.*parp(82)**2/s
  Return
  1000 Format (1X, 'Error: unrecognized beam particle ''', A, '''.'/1X, 'Execution stopped!')
  1100 Format (1X, 'Error: unrecognized target particle ''', A, '''.'/1X, 'Execution stopped!')
  1800 Format (1X, 'Error: unrecognized coordinate frame ''', A, '''.'/1X, 'Execution stopped!')
  1900 Format (1X, 'Error: too low CM energy,', F8.3, ' GeV for event ', 'generation.'/1X, 'Execution stopped!')
End Subroutine pyinki

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
program underflow

    use, intrinsic :: iso_fortran_env, only : sp => REAL32, dp => REAL64
    implicit none
    integer, parameter :: n = 5
    real(kind=sp), dimension(n) :: x = [ 1.0e-31_sp, 1.0e-31_sp, 5.0_sp, &
                                         1.0e31_sp, 1.0e31_sp ]
    integer :: i;
    real(kind=sp) :: prod
    real(kind=sp) :: log_prod
    real(kind=dp) :: dprod
    prod = 1.0_sp
    log_prod = 0.0_sp
    dprod = 1.0_dp
    do i = 1, n
        prod = prod*x(i)
    end do
    do i = 1, n
        dprod = dprod*x(i)
    end do
    do i = 1, n
        log_prod = log_prod + log(x(i))
    end do
    print '(A, F20.15)', 'REAL32 product = ', prod
    print '(A, F20.15)', 'REAL64 product = ', dprod
    print '(A, F20.15)', 'log REAL32 product = ', exp(log_prod)
    prod = 1.0_sp;
    do i = 1, n
        prod = prod*x(i)
    end do
    print '(A, F20.15)', 'REAL32 product (reverse) = ', prod

end program underflow

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


    Program dsbevd_example

!     DSBEVD Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_file_print_matrix_real_gen
      Use lapack_interfaces, Only: dsbevd
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Integer :: i, ifail, info, j, kd, ldab, ldz, liwork, lwork, n
      Character (1) :: job, uplo
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: ab(:, :), w(:), work(:), z(:, :)
      Integer, Allocatable :: iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: max, min
!     .. Executable Statements ..
      Write (nout, *) 'DSBEVD Example Program Results'
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n, kd
      ldab = kd + 1
      ldz = n
      liwork = 5*n + 3
      lwork = 2*n*n + 5*n + 1
      Allocate (ab(ldab,n), w(n), work(lwork), z(ldz,n), iwork(liwork))

!     Read A from data file

      Read (nin, *) uplo
      If (uplo=='U') Then
        Do i = 1, n
          Read (nin, *)(ab(kd+1+i-j,j), j=i, min(n,i+kd))
        End Do
      Else If (uplo=='L') Then
        Do i = 1, n
          Read (nin, *)(ab(1+i-j,j), j=max(1,i-kd), i)
        End Do
      End If

      Read (nin, *) job

!     Calculate all the eigenvalues and eigenvectors of A
      Call dsbevd(job, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, iwork, &
        liwork, info)

      Write (nout, *)
      If (info>0) Then
        Write (nout, *) 'Failure to converge.'
      Else

!       Print eigenvalues and eigenvectors

        Write (nout, *) 'Eigenvalues'
        Write (nout, 100) w(1:n)
        Write (nout, *)
        Flush (nout)

!       Standardize the eigenvectors so that first elements are non-negative.
        Do i = 1, n
          If (z(1,i)<0.0_dp) Then
            z(1:n, i) = -z(1:n, i)
          End If
        End Do

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call nagf_file_print_matrix_real_gen('General', ' ', n, n, z, ldz, &
          'Eigenvectors', ifail)

      End If

100   Format (3X, (8F8.4))
    End Program


text:
program pointers
    INTEGER, POINTER :: val1
end program

max_stars_repo_name: phil-blain/fortls

max_issues_repo_name: phil-blain/fortls

max_forks_repo_name: phil-blain/fortls


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
C
C-----MATRIX-VECTOR MULTIPLY FOR REAL SPARSE NONSYMMETRIC MATRICES---------
C
      SUBROUTINE OPNONSYMt(n,p,W,U,ad,au,jq,irow)
!
      implicit none
!
C-----------------------------------------------------------------------
      integer  IROW(*),JQ(*),n,l,i,j,llast
      real*8   U(*),W(*),Au(*),AD(*),p(*)
C-----------------------------------------------------------------------
C>    SPARSE MATRIX-VECTOR MULTIPLY FOR LANCZS  U = A^T*W
C>    SEE USPEC SUBROUTINE FOR DESCRIPTION OF THE ARRAYS THAT DEFINE
C>    THE MATRIX
c>    the vector p is not needed but is kept for compatibility reasons
c>    with the calling program
C-----------------------------------------------------------------------
C
C     COMPUTE THE DIAGONAL TERMS
      DO 10 I = 1,N
         U(I) = AD(I)*W(I)+U(I)
 10   CONTINUE
C
C     COMPUTE BY COLUMN
      LLAST = 0
      DO 30 J = 1,N
C
         DO 20 L = JQ(J),JQ(J+1)-1
            I = IROW(L)
C
            U(j) = U(j) + Au(L)*W(i)
C
 20      CONTINUE
C
 30   CONTINUE
C
      RETURN
      END

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
C ************************************************************************
C  FILE: ser_pi_calc.f
C  DESCRIPTION:
C    Serial pi Calculation - Fortran Version
C    This program calculates pi using a "dartboard" algorithm.  See
C    Fox et al.(1988) Solving Problems on Concurrent Processors, vol.1
C    page 207.
C    Note: Requires Fortran90 compiler due to random_number() function
C  AUTHOR: unknown
C  LAST REVISED: 02/23/12 Blaise Barney
C ************************************************************************
C Explanation of constants and variables used in this program:
C   DARTS               = number of throws at dartboard
C   ROUNDS              = number of times "DARTS" is iterated
C   pi                  = average of pi for this iteration
C   avepi               = average pi value for all iterations
C ************************************************************************

	program pi_calc

	integer  DARTS, ROUNDS
	parameter(DARTS = 10000)
	parameter(ROUNDS = 100)
        real*4    seednum
        real*8    pi, avepi, dboard

        print *,'Starting serial version of pi calculation...'

 	avepi = 0
	do 40 i = 1, ROUNDS
	  pi = dboard(DARTS)
          avepi = ((avepi*(i-1)) + pi)/ i
          write(*,32) DARTS*i, avepi
 32       format('   After',i8,' throws, average value of pi = ',
     &           f10.8,$)
 40	continue
        print *, ' '
        print *,'Real value of PI: 3.1415926535897'
        print *, ' '
	end

C ************************************************************************
C  dboard.f
C ************************************************************************
C Explanation of constants and variables used in this function:
C   darts       = number of throws at dartboard
C   score       = number of darts that hit circle
C   n           = index variable
C   r           = random number between 0 and 1
C   x_coord     = x coordinate, between -1 and 1
C   x_sqr       = square of x coordinate
C   y_coord     = y coordinate, between -1 and 1
C   y_sqr       = square of y coordinate
C   pi          = computed value of pi
C ************************************************************************
      real*8    function dboard(darts)

      integer   darts, score, n
      real*4    r
      real*8    x_coord, x_sqr, y_coord, y_sqr, pi

      score = 0

C  Throw darts at board.  Done by generating random numbers
C  between 0 and 1 and converting them to values for x and y
C  coordinates and then testing to see if they "land" in
C  the circle."  If so, score is incremented.  After throwing the
C  specified number of darts, pi is calculated.  The computed value
C  of pi is returned as the value of this function, dboard.
C  Note:  the seed value for rand() is set in pi_calc.

C    Note: Requires Fortran90 compiler due to random_number() function
      print *, ' '
      do 10 n = 1, darts
        call random_number(r)
        x_coord = (2.0 * r) - 1.0
        x_sqr = x_coord * x_coord

        call random_number(r)
        y_coord = (2.0 * r) - 1.0
        y_sqr = y_coord * y_coord

        if ((x_sqr + y_sqr) .le. 1.0) then
          score = score + 1
        endif

 10   continue

      pi = 4.0 * score / darts
      dboard = pi
      end

max_stars_repo_name: hixio-mh/HPC-Tutorials

max_issues_repo_name: hixio-mh/HPC-Tutorials

max_forks_repo_name: hixio-mh/HPC-Tutorials


text:
!!# MODULE: SUBROUTINE <SUB_levQuadratureReflect>
MODULE SUB_levQuadratureReflect

!!## EXTERNAL KINDS
USE KND_IntrinsicTypes,ONLY: KIND_Rsp,KIND_Rdp !!((01-A-KND_IntrinsicTypes.f90))

!!## EXTERNAL PARAMETERS

!!## EXTERNAL PROCEDURES
USE FUN_Reverse                                !!((04-A-FUN_Reverse.f90))
USE FUN_SIZEa                                  !!((06-B-FUN_SIZEa.f90))
USE SUB_aQuadratureReflect                     !!((05-B-SUB_aQuadratureReflect.f90))

!!## DEFAULT IMPLICIT
IMPLICIT NONE

!!## DEFAULT ACCESS
PRIVATE

!!## PROCEDURE LISTING
INTERFACE levQuadratureReflect
 MODULE PROCEDURE levQuadratureReflect_Rsp
 MODULE PROCEDURE levQuadratureReflect_Rdp
END INTERFACE

!!## PUBLIC ACCESS LIST
PUBLIC :: levQuadratureReflect

!!## MODULE PROCEDURES
CONTAINS

!!### PURE SUBROUTINE: levQuadratureReflect_Rsp
PURE SUBROUTINE levQuadratureReflect_Rsp( xa,xp,wa,wp )
!!#### LOCAL KINDS
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rsp  !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "07-B-SUB_levQuadratureReflect.f90.hdr"
!!--begin--
INCLUDE "07-B-SUB_levQuadratureReflect.f90.bdy"
!!--end--
END SUBROUTINE

!!### PURE SUBROUTINE: levQuadratureReflect_Rdp
PURE SUBROUTINE levQuadratureReflect_Rdp( xa,xp,wa,wp )
!!#### LOCAL KINDS
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rdp  !!((01-A-KND_IntrinsicTypes.f90))
INCLUDE "07-B-SUB_levQuadratureReflect.f90.hdr"
!!--begin--
INCLUDE "07-B-SUB_levQuadratureReflect.f90.bdy"
!!--end--
END SUBROUTINE

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
*> \brief \b DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLADIV + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLADIV( A, B, C, D, P, Q )
*
*       .. Scalar Arguments ..
*       DOUBLE PRECISION   A, B, C, D, P, Q
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLADIV performs complex division in  real arithmetic
*>
*>                       a + i*b
*>            p + i*q = ---------
*>                       c + i*d
*>
*> The algorithm is due to Michael Baudin and Robert L. Smith
*> and can be found in the paper
*> "A Robust Complex Division in Scilab"
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] A
*> \verbatim
*>          A is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] C
*> \verbatim
*>          C is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is DOUBLE PRECISION
*>          The scalars a, b, c, and d in the above expression.
*> \endverbatim
*>
*> \param[out] P
*> \verbatim
*>          P is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[out] Q
*> \verbatim
*>          Q is DOUBLE PRECISION
*>          The scalars p and q in the above expression.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date January 2013
*
*> \ingroup auxOTHERauxiliary
*
*  =====================================================================
      SUBROUTINE DLADIV( A, B, C, D, P, Q )
*
*  -- LAPACK auxiliary routine (version 3.5.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     January 2013
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   A, B, C, D, P, Q
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   BS
      PARAMETER          ( BS = 2.0D0 )
      DOUBLE PRECISION   HALF
      PARAMETER          ( HALF = 0.5D0 )
      DOUBLE PRECISION   TWO
      PARAMETER          ( TWO = 2.0D0 )
*
*     .. Local Scalars ..
      DOUBLE PRECISION   AA, BB, CC, DD, AB, CD, S, OV, UN, BE, EPS
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           DLADIV1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX
*     ..
*     .. Executable Statements ..
*
      AA = A
      BB = B
      CC = C
      DD = D
      AB = MAX( ABS(A), ABS(B) )
      CD = MAX( ABS(C), ABS(D) )
      S = 1.0D0

      OV = DLAMCH( 'Overflow threshold' )
      UN = DLAMCH( 'Safe minimum' )
      EPS = DLAMCH( 'Epsilon' )
      BE = BS / (EPS*EPS)

      IF( AB >= HALF*OV ) THEN
         AA = HALF * AA
         BB = HALF * BB
         S  = TWO * S
      END IF
      IF( CD >= HALF*OV ) THEN
         CC = HALF * CC
         DD = HALF * DD
         S  = HALF * S
      END IF
      IF( AB <= UN*BS/EPS ) THEN
         AA = AA * BE
         BB = BB * BE
         S  = S / BE
      END IF
      IF( CD <= UN*BS/EPS ) THEN
         CC = CC * BE
         DD = DD * BE
         S  = S * BE
      END IF
      IF( ABS( D ).LE.ABS( C ) ) THEN
         CALL DLADIV1(AA, BB, CC, DD, P, Q)
      ELSE
         CALL DLADIV1(BB, AA, DD, CC, P, Q)
         Q = -Q
      END IF
      P = P * S
      Q = Q * S
*
      RETURN
*
*     End of DLADIV
*
      END

      SUBROUTINE DLADIV1( A, B, C, D, P, Q )
*
*  -- LAPACK auxiliary routine (version 3.5.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     January 2013
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   A, B, C, D, P, Q
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE
      PARAMETER          ( ONE = 1.0D0 )
*
*     .. Local Scalars ..
      DOUBLE PRECISION   R, T
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLADIV2
      EXTERNAL           DLADIV2
*     ..
*     .. Executable Statements ..
*
      R = D / C
      T = ONE / (C + D * R)
      P = DLADIV2(A, B, C, D, R, T)
      A = -A
      Q = DLADIV2(B, A, C, D, R, T)
*
      RETURN
*
*     End of DLADIV1
*
      END

      DOUBLE PRECISION FUNCTION DLADIV2( A, B, C, D, R, T )
*
*  -- LAPACK auxiliary routine (version 3.5.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     January 2013
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   A, B, C, D, R, T
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D0 )
*
*     .. Local Scalars ..
      DOUBLE PRECISION   BR
*     ..
*     .. Executable Statements ..
*
      IF( R.NE.ZERO ) THEN
         BR = B * R
         if( BR.NE.ZERO ) THEN
            DLADIV2 = (A + BR) * T
         ELSE
            DLADIV2 = A * T + (B * T) * R
         END IF
      ELSE
         DLADIV2 = (A + D * (B / C)) * T
      END IF
*
      RETURN
*
*     End of DLADIV12
*
      END

max_stars_repo_name: hpanderson/OpenBLAS

max_issues_repo_name: hpanderson/OpenBLAS

max_forks_repo_name: hpanderson/OpenBLAS


text: #redirect wiki:Sacramento:Kaveri Madras Cuisine

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
The UC Davis Campus Judicial Board (CJB) is a University body which has two primary responsibilities:  1) Educational outreach to faculty and students, 2) Hearing and deciding cases of contested student misconduct. The CJB consists of 12 student members and 8 faculty members. The student members are appointed under the authority of the Vice Chancellor for Student Affairs, but are selected through an interview and recruitment process conducted by the current/outgoing student members and the CJB Staff Advisor (who is also a SJA officer). The faculty members of the Board are appointed by the Academic Senate and serve on the Board for 2 year periods (IIRC). Because their mission and purpose is closely related to the duties of the Office of Student Judicial Affairs, the two groups often work sidebyside. The students of the CJB work out of the Office of Student Judicial Affairs in Dutton Hall

 Duties
CJB can also refer only to the student members, as they will meet and work together to conduct campus outreach programs in order to promote academic and social integrity amongst students. Other responsibilities of CJB students include serving as advisors or advocates for students or faculty who are going into formal hearings (for contested misconduct). CJB students also conduct followup meetings with students who have been through the informal disciplinary process to obtain valuable assessment information, as well as help facilitate a more comprehensive learning experience for students who have been referred for misconduct. Responsibilities and duties of CJB students have been known to vary, as well.

Typically, Campus Judicial Board hearings will be heard by a panel consisting of 2 student members and 1 faculty member. One of the student members chairs the formal hearing. Previously, more serious cases under SJA jurisdiction would go to the Student Conduct Committee, however, that setup has recently been eliminated. All cases that warrant a formal hearing now go to the Campus Judicial Board. Incidents that can result in a referral to SJA can occur oncampus and offcampus and can range from Academic Dishonesty to assault.

 Titles Held By Members of the CJB
  ALL Members are Campus Judicial Board Appointees.
  ALL Members are Vested Judicial Authorities of ASUCD, in accordance with the ASUCD Constitution.
  ALL Members are UC Davis Campus Officials as defined by the Administration of Student Discipline.
  TWO Members are appointed as Cochairs of the Board each academic school year.

 Relation to ASUCD

As of a few years ago, ASUCD amended its Government Codes to include the Campus Judicial Board in its procedures. Complaints against the Elections Committee, for example, are to be brought before the Campus Judicial Board. In recent years, this has happened on a few occasions: the socalled Election Statement Controversy, as well as an issue regarding the Move to Division I in Winter 2002. Under this authority, CJB heard and decided the Election Statement Controversy of Fall 2004. CJB eventually ruled against the Elections Committee, but declined to recommend disciplinary action (instead leaving this to the discretion of the Senate). Shortly afterwards, ASUCD removed CJB from its role in procedural review, returning such powers to the ASUCD Court. Although the CJB appears powerless now within ASUCD functions, it is still mentioned in the ASUCD Constitution as having the ASUCDs judicial power vested within it. Hence, purely in name, the CJB remains a figurehead of ASUCD judicial authority, alongside the powerenabled ASUCD Court. Since this relation still exists, campaigns within student government are prohibited from utilizing endorsements from CJB members, as a way of promoting the impartial image of the CJB.

It should be noted that if ASUCD cases are referred to Student Judicial Affairs, and those cases are contested by all parties involved, the matter will go to a formal hearing. Most likely, the Campus Judicial Board would preside over that hearing as opposed to a single Hearing Officer which sometimes might be done.

 History

The Campus Judicial Board was created in 1976 under a modified honor code which UC Davis still upholds today. It was created in conjunction with the Academic Code of Conduct along with Student Judicial Affairs. Although the board itself is a fairly recent creation, the involvement of students in the disciplinary process at UC Davis can be traced to the beginning of the 20th century when the first student government was formed.

More information can be found on the CJB website listed above.

 Current Campus Judicial Board Appointees

 All CJB members hold office hours in SJA at varying times.

 Joining the Board

Student Judicial Affairs recruits students for open positions on the following years board during Winter Quarter. Generally, applications are available in January and are due in February. Interested students should contact sja@ucdavis.edu to be entered in the SJA listserve for prospective applicants. Everyone from undergraduate students to UCD professional school students are encouraged to apply.


text:
c ---------------------------------------------------------------------
c
c Subroutine that computes and assign Gaia errors to observables
c Errors are computed following Gaia Science Performance web page:
c http://www.cosmos.esa.int/web/gaia/science-performance
c
c Universitat de Barcelona
c Contact: Merce Romero-Gomez, Josep Manel Carrasco
c
c ---------------------------------------------------------------------

c Updated: October 2014

c Inputs:
c ---------

c   * External tables:
c     ----------------
c	'gfactor-Oct2012.dat': geometrical factors (Table 2 of the web)
c	'gfactor-Jun2013.dat': geometrical factors (Table 6 of the web)
c 	'TableVr-Oct2014.dat': Radial Velocity coefficients

c	Input files for the definition of constants:
c   		const_math.h : mathematical constants
c   		const_ast.h  : astronomical constants

c   * Parameters from the star:
c     -------------------------

c     jflag:
c     	The code offers two options depending of the type of errors
c     	you want to apply (see website):
c     		jflag=1, errors computed from mean geometrical factors (see website)
c     		jflag=-1,errors computed considering the scanning law of the satellite
c             		(the geometrical factors are then computed from the ecliptic
c             		latitude and the number transits)

c     V: V apparent magnitude
c     VI: observed (V-I) colour index (affected by extinction), in magnitudes

c     a(i): Astrometric actual values for the star (not affected by errors)
c         a(1): Equatorial heliocentric right ascension (units: radians)
c         a(2): Equatorial heliocentric declination: delta (units: radians)
c         a(3): Parallax: pi (units: mas)
c         a(4): Equatorial proper motions in right ascension in true arcs on the sky: mua_{*}=mua*cos(delta) (units: mas/yr)
c         a(5): Equatorial proper motions: mud (units: mas/yr)
c         a(6): Radial velocity: Vr (units: km/s)

c     ap(i): Atmospheric parameters actual values for the star (no affected by errors)
c         ap(1): Effective temperature (units: kelvin)
c         ap(2): surface gravitye (logg) (units: dex)
c         ap(3): [Fe/H] (units: dex)
c         ap(4): Absorption A0 (Av~A0, see Bailer-Jones et al., 2011)

c Output values:
c -------------
c Parameters of the star as observed by Gaia (affected by errors)
c Gaia errors assigned to each parameter:

c The output parameters of the star as observed by Gaia are:

c   Astrometric (+Vr) data affected by errors
c     ao(i): Astrometric values for the star affected by Gaia errors
c         ao(1): Observed Equatorial heliocentric right ascension (units: radians)
c         ao(2): Observed Equatorial heliocentric declination: delta (units: radians)
c         ao(3): Observed Parallax: pi (units: mas)
c         ao(4): Observed Equatorial proper motions in right ascension in true arcs on the sky: mua_{*}=mua*cos(delta) (units: mas/yr)
c         ao(5): Observed Equatorial proper motions: mud (units: mas/yr)
c         ao(6): Observed Radial velocity: Vr (units: km/s)
c     ae(i): Gaia standard deviation (1-sigma)
c         ae(1): Standard deviation in Equatorial heliocentric right ascension in true arcs on the sky:: alpha_{*}=alpha*cos(delta)  (units: mas)
c         ae(2): Standard deviation in Equatorial heliocentric declination: delta (units: mas)
c         ae(3): Standard deviation in Parallax: pi (units: mas)
c         ae(4): Standard deviation in Equatorial proper motions in right ascension in true arcs on the sky: mua_{*}=mua*cos(delta) (units: mas/yr)
c         ae(5): Standard deviation in Equatorial proper motions: mud (units: mas/yr)
c         ae(6): Standard deviation in Radial velocity: Vr (units: km/s)

c   Photometric data affected by errors (G, G_BP, G_RP, G_RVS)
c     p(i): Photometric actual values for the star (no affected by errors)
c         p(1): G magnitude (units: magnitudes)
c         p(2): G_BP-G_RP   (units: magnitudes)
c     po(i): Photometric values of the star affected by Gaia errors
c         po(1): Observed G magnitude (units: magnitudes)
c         po(2): Observed G_BP-G_RP   (units: magnitudes)
c     pe(i): Gaia standard deviation (errors, 1-sigma)
c         pe(1): Standard deviation in G magnitude (units: magnitudes)
c         pe(2): Standard deviation in G_BP-G_RP   (units: magnitudes)

c   Atmospheric parameters data affected by errors:
c     apo(i): Atmospheric parameters of the star affected by Gaia errors
c         apo(1): Observed Effective temperature (units: kelvin)
c         apo(2): Observed surface gravity (logg) (units: dex)
c         apo(3): Observed [Fe/H] (units: dex)
c         apo(4): Observed Absorption A0 (Av~A0, see Bailer-Jones et al., 2011) (units: magnitudes)
c     ape(i): Gaia standard deviation (1-sigma)
c         ap(1): Standard deviation in Effective temperature (units: kelvin)
c         ap(2): Standard deviation in surface gravity (logg) (units: dex)
c         ap(3): Standard deviation in [Fe/H] (units: dex)
c         ap(4): Standard deviation in Absorption A0 (Av~A0, see Bailer-Jones et al., 2011) (units: magnitudes)

c **********************************************************************
        subroutine Gaia_errors(jflag,V,VI,a,ao,ae,p,po,pe,
     &             ap,apo,ape,GRVS,Llavor)

	    implicit real*8 (a-h,o-z)
        dimension a(6),ao(6),ae(6),p(2),po(2),pe(2),g(5)
        dimension ap(4),apo(4),ape(4)
        dimension xvi(11),xavr(11),xbvr(11)
        dimension xsb(20),xNOBS(20),xga(20),xgd(20),xgpi(20),xgmua(20),
     &            xgmud(20)
        real*8 gasdev
        character*10 xchar
        logical ifirst
        data ifirst/.true./
c        save xsb,xNOBS,xga,xgd,xgpi,xgmua,xgmud,xnmean,gmeana,gmeand,
c     &       gmeanpi,gmeanmua,gmeanmud
        save xsb,xga,xgd,xgpi,xgmua,xgmud
        INCLUDE 'const_math.h'
        INCLUDE 'const_ast.h'
        COMMON / tableVr / xvi, xavr, xbvr
        integer Llavor

        if(ifirst)then
c Table 2 from the Gaia Science Performance Webpage. Only geometrical factors
c NOT including the effect of the number of passages (including 6.2% dead time)
c$$$        open(3,file='gfactor-Oct2012.dat')
c$$$            read(3,*)
c$$$            do i=1,20
c$$$              read(3,*)xsb(i),xk,xk,xNOBS(i),xga(i),xgd(i),xgpi(i),
c$$$     &                 xgmua(i),xgmud(i)
c$$$            enddo
c$$$            read(3,*)xchar,xk,xk,xnmean,gmeana,gmeand,gmeanpi,gmeanmua,
c$$$     &               gmeanmud
c$$$            close(3)

c Table 6 from the Gaia Science Performance Webpage. It includes the
c geometrical factors and the effect of the number of passages
c (including 0% dead time, which affects the number of passages,N_obs, and
c also the mean number. Then it applies the same formula as above)
        open(3,file='gfactor-Jun2013.dat')
            read(3,*)
            do i=1,20
              read(3,*)xsb(i),xga(i),xgd(i),xgpi(i),xgmua(i),xgmud(i)
            enddo
            read(3,*)xchar,gmeana,gmeand,gmeanpi,gmeanmua,gmeanmud
            close(3)

       open(3,file='TableVr-Oct2014.dat')
            read(3,*)
            do i=1,11
              read(3,*)xvi(i),xavr(i),xbvr(i)
            enddo
            close(3)
        ifirst=.false.
        endif

c Get Gaia errors in parallax,magnitudes and Vr
c        write(*,*)'aki'
        call ErrorsPiMag(v,vi,p,pe,empi,a(3))
        call ErrorsVr(v,vi,ae(6),GRVS)

c       Output errors on parallax and proper motion will be in mas
c        empi = empi /1e3

c       Mean errors (from mean g-geometrical factor)
        if(jflag.gt.0)then
          ae(3)=empi
          ae(1)=gmeana*ae(3)
          ae(2)=gmeand*ae(3)
          ae(4)=gmeanmua*ae(3)
          ae(5)=gmeanmud*ae(3)
        endif

c       Errors depending on the scanning law
c       (number of passages and the geometrical factor)
        if(jflag.lt.0)then
c Version using Table 2
c$$$          sbeta=dabs(0.9175*dsin(a(2))-0.3978*dcos(a(2))*dsin(a(1)))
c$$$
c$$$          call lininter(20,xsb,xNobs,sbeta,gns)
c$$$          call lininter(20,xsb,xga,sbeta,g(1))
c$$$          call lininter(20,xsb,xgd,sbeta,g(2))
c$$$          call lininter(20,xsb,xgpi,sbeta,g(3))
c$$$          call lininter(20,xsb,xgmua,sbeta,g(4))
c$$$          call lininter(20,xsb,xgmud,sbeta,g(5)) !!!change here
c$$$
c$$$          xnob=dsqrt(xnmean/gns)
c$$$
c$$$c          write(*,*)sbeta,g(1),xnob*g(1)
c$$$
c$$$          do i=1,5
c$$$          ae(i)=xnob*g(i)*empi
c$$$          enddo

c Version using Table 6
          sbeta=dabs(0.9175*dsin(a(2))-0.3978*dcos(a(2))*dsin(a(1)))

          call lininter(20,xsb,xga,sbeta,g(1))
          call lininter(20,xsb,xgd,sbeta,g(2))
          call lininter(20,xsb,xgpi,sbeta,g(3))
          call lininter(20,xsb,xgmua,sbeta,g(4))
          call lininter(20,xsb,xgmud,sbeta,g(5))

          do i=1,5
          ae(i)=g(i)*empi
          enddo
        endif

c       Computation of the errors in the atmospheric parameters
c       See Liu et al., 2012
c       A second order polinomial on G has been fitted to the
c       |Teff_eaneas_pq - Teff_real|

        call errorsAp(p(1),ape)

c       Computation of the observed astrometric quantities

c       The error in right ascension ae(1) denotes true arc on the sky
c       so the right ascension shall be converted to that before the
c       random error is assigned
c       alpha_{*}=alpha*cos(delta)

        a(1)=a(1)*dcos(a(2))

c	The errors in alpha* and  delta are in mas whereas alpha and delta are
c       in radiants, so (alpha*,delta)  are converted from radiants to mas
c       before the random error is assigned

	a(1)=a(1)/mas
	a(2)=a(2)/mas

        do i=1,6
        ao(i)=gasdev(Llavor,a(i),ae(i))
        enddo

c       The true and observed (alpha*,delta) are converted from mas to radiants
	a(1)=a(1)*mas
	a(2)=a(2)*mas
	ao(1)=ao(1)*mas
	ao(2)=ao(2)*mas

c       and the alpha_{*} is converted to alpha
        ao(1)=ao(1)/dcos(a(2))
c       we also convert back the input value
        a(1)=a(1)/dcos(a(2))

c       Computation of the observed photometric quantities
        do i=1,2
        po(i)=gasdev(Llavor,p(i),pe(i))
        enddo

c       Computation of the observed atmospheric parameters
        do i=1,4
        apo(i)=gasdev(Llavor,ap(i),ape(i))
        enddo
        return
	end

c******************************************************************
c     Subroutine ErrorsAp
c
c     A second order polinomial on G has been fitted to the
c     |Teff_eaneas_pq - Teff_real|
c     Commented, there is the option to fitting a second order
c     polinomial to the sigmas of the Aeneas pq-model
c     (see tableA4 of Liu & Bailer-Jones 2012)
c
c     input data
c     p(1)=G           Gaia apparent magnitude (G)

c     output data
c     ape(1):          Uncertainty in Teff (Kelvin)
c     ape(2):          Uncertainty in logg (dex)
c     ape(3):          Uncertainty in [Fe/H] (dex)
c     ape(4):          Uncertainty in A0 (mag)

c------------------------------------------------------------------

      subroutine errorsAp(G,ape)

      implicit real*8 (a-h,o-z)
      dimension ape(4)

c     Fitting |Teff_eaneas_pq - Teff_real|
      ape(1)= 630.3d0   - 107.54d0  *G + 4.9682d0    *G**2
      ape(2)= 0.49865d0 - 0.04432d0 *G + 0.0017055d0 *G**2
      ape(3)= 0.81192d0 - 0.12226d0 *G + 0.0056669d0 *G**2
      ape(4)= 0.56387d0 - 0.093435d0*G + 0.0042024d0 *G**2

c     Fitting sigTeff (Eaneas pq-model)
c      ape(1)=67.289d0   - 7.8492d0  *G + 0.34687d0   *G**2
c      ape(2)= 0.31915d0 - 0.044319d0*G + 0.0020501d0 *G**2
c      ape(3)= 0.71709d0 - 0.115d0   *G + 0.0050649d0 *G**2
c      ape(4)= 0.17787d0 - 0.026371d0*G + 0.0011578d0 *G**2

      return
      end

c******************************************************************
c     Subroutine ErrorsPiMag
c
c     input data
c     v=V           Johnson V apparent magnitude
c     vi=(V-I)c     Johnson-Cousins V-I colour

c     output data
c     p(1)=G                true Gaia G apparent magnitude
c     p(2)=G_BP-G_RP        true Gaia colour
c     pe(1)=sigG            uncertainty in G (in mag units)
c     pe(2)=sigBPRP         uncertainty in G_BP-G_RP colour (in mag units)
c     ae(3)=empi            uncertainty in parallax (in microarcsec units)

c------------------------------------------------------------------
      subroutine ErrorsPiMag(v,vi,p,pe,empi,xpi)

      implicit real*8 (a-h,o-z)
      dimension p(2),pe(2)

c---- Gaia photometry,
      p(2)=-0.0660d0+1.2061d0*vi-0.0614d0*vi*vi+0.0041d0*vi**3
      p(1)=V-0.0257d0-0.0924d0*vi-0.1623d0*vi*vi+0.0090d0*vi**3
      if(p(1).gt.20.) then
         empi=999999.999
         pe(1)=99.999
         pe(2)=99.999
         return
      endif

c---- Gaia error in parallax (assuming 70 transits)
c---- A - before commissioning error model
      z=MAX(10**(0.4*(12.0-15.0)), 10**(0.4*(p(1)-15.0)))
c      empi=dsqrt(9.3+658.1*z+4.568*z**2)*(0.986+(1.0-0.986)*vi)

c-----B - after commissioning
c----- new fitting provided by Rygl,Antoja,DeBruijne et al (October 2014)
c----- it assumes no dependance on the colour
c      empi=dsqrt(-11.5+706.1*z+32.6*z**2)*(0.986+(1.0-0.986)*vi)
c      empi=0.1*xpi
      empi=0.07*xpi
c      empi=0.16*xpi
c      empi=0.20*xpi
c      write(*,*)empi,xpi
      pe(1)=0.001*dsqrt(0.02076*z**2+2.7224*z+0.004352)
      pe(1)=pe(1)/dsqrt(70.d0)

      z=MAX(10**(0.4*(11.0-15.0)), 10**(0.4*(p(1)-15.0)))

      a=-0.003201*vi**3+0.0589*vi**2+0.3353*vi+0.7927
      b=-0.001019*vi**3+0.0244*vi**2+0.1756*vi+1.4684
      c=-0.004093*vi**3+0.0740*vi**2+0.2834*vi-3.4772
      sigBP=0.001*dsqrt((10.**a)*(z**2)+(10.**b)*z+10.**c)
      sigBP=sigBP/dsqrt(70.d0)

      a=-0.006560*vi**3+0.1080*vi**2-0.6296*vi+1.4470
      b=-0.003280*vi**3+0.0540*vi**2-0.3148*vi+1.7856
      c=-0.007992*vi**3+0.1482*vi**2-0.7544*vi-3.7232
      sigRP=0.001*dsqrt((10.**a)*(z**2)+(10.**b)*z+10.**c)
      sigRP=sigRP/dsqrt(70.d0)

c     We assume there is no correlation among colours
      pe(2)=dsqrt(sigBP**2+sigRP**2)

      return
      end

c******************************************************************
c     Subroutine ErrorsVr
c
c     input data
c     v=V           Johnson V apparent magnitude
c     vi=(V-I)c     Johnson-Cousins V-I colour

c     output data
c     ae(6)=Sigr            uncertainty in Vr (in km/s units)

c------------------------------------------------------------------
      subroutine ErrorsVr(v,vi,sigVr,GRVS)

      implicit real*8 (a-h,o-z)
      dimension xavr(11),xbvr(11),xvi(11)
      COMMON / tableVr / xvi, xavr, xbvr
c
c     Vr errors
c
c---- Before commissioning errors
      GRVS=V-0.0119D0-1.2092D0*VI+0.0188D0*VI*VI+0.0005D0*VI*VI*VI
      call lininter(11,xvi,xavr,vi,avr)
      call lininter(11,xvi,xbvr,vi,bvr)
      sigVr=1.d0+bvr*exp(avr*(V-14))

      return
      end

c******************************************************************

c     subroutine get_fi
c input: (alfa,delta) in radians
c output: fi parallactic angle in radians

C  --------------------------------------------------------------------
      subroutine get_fi(alfa,delta,fi)
      implicit real*8(a-h,o-z)
      include 'const_math.h'
      include 'const_ast.h'

      x=Dcos(Dg)*Dsin(Delta)*Dsin(Alfa-Ag)+Dsin(Dg)*Dcos(Delta)
      y=Dcos(Dg)*Dcos(Alfa-Ag)
c Parallactic angle
      Fi=Datan2(y,x)
      return
      end

c **********************************************************************
c          subroutine lininter(n,xa,ya,x,y)
c input:
c   xa(n) and ya(n) such that ya(i)=f(xa(i)), and xa ordered in
c                   increasing sense
c   x, point where you want to know the corresponding y=f(x)
c output:
c   y
C  --------------------------------------------------------------------

           subroutine lininter(n,xa,ya,x,y)
           implicit real*8(a-h,o-z)
           dimension xa(n),ya(n)

           if(x.gt.xa(n))then
             y=ya(n)
             return
           else if(x.lt.xa(1))then
             y=ya(1)
             return
           else if(xa(1).le.x.and.x.le.xa(n))then
             i=1
             do while (xa(i).lt.x)
               i=i+1
             enddo
             i=i-1     !!!change here
             x0=xa(i)
             x1=xa(i+1)
             y0=ya(i)
             y1=ya(i+1)
             y=((x-x0)/(x1-x0))*y1+((x1-x)/(x1-x0))*y0
             return
           endif

           end

c **********************************************************************
c Coordinate transformation: from galactic to equatiorial
c subroutine galactic_to_equatorial
c input: xi(l,b,pi,mul*cosb,mub,vr)
c output: xo(alpha,delta,pi,mua*cos(delta),mud,vr)
c Units:
c angles l,b,alpha,delta in radians
c proper motions in arcsec/yr
c note that pi(xi(3)) and vr(xi(6)) do not change
c------------------------------------------------------------------------

      subroutine galactic_to_equatorial(xi,xo)
      implicit real*8(a-h,o-z)
      dimension xi(6),xo(6)
      include 'const_math.h'
      include 'const_ast.h'

c l,b to alpha,delta
      cb=Dcos(xi(2))
      sb=Dsin(xi(2))
      cl=Dcos(xi(1))
      sl=Dsin(xi(1))

      sd=cb*Dsin(xi(1)-a2)*Dsin(a1)+sb*Dcos(a1)
      xo(2)=Dasin(sd)

      y=cb*Dsin(xi(1)-a2)*Dcos(a1)-sb*Dsin(a1)
      x=cb*Dcos(xi(1)-a2)
      xo(1)=Datan2(y,x) +Ag

c mulcosb,mub to muas,mud
      call get_fi(xo(1),xo(2),fi)
      cfi=dcos(fi)
      sfi=dsin(fi)
      xo(4)=cfi*xi(4)-sfi*xi(5)
      xo(5)=sfi*xi(4)+cfi*xi(5)
c
      xo(3)=xi(3)
      xo(6)=xi(6)
      return
      end

c **********************************************************************
c Coordinate transformation: from equatiorial to galactic
c subroutine equatorial_to_galactic(xi,xo)
c input: xi(alpha,delta,pi,mua*cos(delta),mud,vr)
c output: xo(l,b,pi,mul*cosb,mub,vr)
c Units:
c angles l,b,alpha,delta in radians
c proper motions in arcsec/yr
c note that pi(xi(3)) and vr(xi(6)) do not change
c -------------------------------------------------------------

      subroutine equatorial_to_galactic(xi,xo)
      implicit real*8(a-h,o-z)
      dimension xi(6),xo(6)
      include 'const_math.h'
      include 'const_ast.h'

c alpha,delta to l,b
      sa=dsin(xi(1)-Ag)
      ca=dcos(xi(1)-Ag)
      sd=dsin(xi(2))
      cd=dcos(xi(2))
      sb=sd*dcos(a1)-cd*dsin(a1)*sa
      xo(2)=dasin(sb)
      xo(1)=datan((sd-dcos(a1)*sb)/(dsin(a1)*cd*ca))+a2

c muastar,mud to mulstar,mub
      call get_fi(xi(1),xi(2),fi)
      cfi=dcos(fi)
      sfi=dsin(fi)
      xo(4)= cfi*xi(4)+sfi*xi(5)
      xo(5)=-sfi*xi(4)+cfi*xi(5)
c
      xo(3)=xi(3)
      xo(6)=xi(6)
      return
      end
c -------------------------------------------------------------


text:
c
!> For new fine grids (grids on level **level**), allocate space for
!! saving fluxes at boundary (even if it's physical boundary) after
!! each integration step.
!! These fluxes are for future conservative fixing of coarse grids
!! (level **level**-1 grids).
!! The address of this space is stored in node(ffluxptr, mkid) for grid
!! mkid.
!!
!! Note that if the refinment ratio is r, fluxes from every r cells
!! on grid **mkid** are saved to one slot in this space, since every r
!! cells on grid **mkid** border one cell on level **level**-1 grids.
!! \param[in] level boudnary lists of grids on this level get updated
!! \param[in] nvar number of equations for the system
!! \param[in] naux number of auxiliary variables for the system
c ----------------------------------------------------------
c
      subroutine prepf(level,nvar,naux)
c
      use amr_module
      implicit double precision (a-h,o-z)

c
c ::::::::::::::::::::::::::: PREPF :::::::::::::::::::::::::::::
c
c prepare new fine grids to save fluxes after each integration step
c for future conservative fixing of coarse grids
c save all boundary fluxes of fine grid (even if phys. bndry.) -
c but only save space for every intrat. (remember - 4 fluxes).
c
c :::::::::::::::::::::::::::::::::::;:::::::::::::::::::::::::::
c
      mkid = lstart(level)
 10   if (mkid .eq. 0) go to 99
          nx    = node(ndihi,mkid)-node(ndilo,mkid) + 1
          ny    = node(ndjhi,mkid)-node(ndjlo,mkid) + 1
          ikeep = nx/intratx(level-1)
          jkeep = ny/intraty(level-1)
          lenbc = 2*(ikeep+jkeep)
c
c         get twice the storage, one for plus or minus fluxes, the
c         other for coarse solution for wave fixing. also for aux vars.
c
          node(ffluxptr,mkid) = igetsp(2*nvar*lenbc+naux*lenbc)
          ist   = node(ffluxptr,mkid)

          do 20 i = 1, 2*nvar*lenbc + naux*lenbc
 20          alloc(ist+i-1) = 0.d0
          mkid = node(levelptr,mkid)
          go to 10

 99    return
       end

max_stars_repo_name: navravi/amrclaw


  subroutine dfxlinop(j,n19)
! --- calculations of more complicated line profiles, i.e. for hydrogen lines, and Stark profile lines.

    use types
    use dfcomm
    use hprof

    implicit none

    include 'common.rhoxbl'
    include 'common.tempbl'
    include 'common.turbpr'
    include 'common.stateb'
    include 'common.depart'
    include 'common.opsblk'
    include 'common.ionsbl'

      integer, intent(in) :: j, n19
!---- Local variables -------------------------------------------!

      integer n, k,  iline,  itemp1
      integer(kind=4) nlast,  ifstrong, minred, maxred, ibuff, maxblue, i, ib
!these are read as integers, but processed later as real

      integer(kind=4)  nblo, nbup, nelion, type, ncon, nelionx, nbuff
      real(kind=4) gammar, gammas, gammaw
      real(kind=4) bshore,  g,  cgf, ashore
      real(kind=4) gf

      equivalence (gammas,ashore),(gammaw,bshore)
      equivalence (gf,g,cgf),(type,nlast)!,(gammar,gaunt)

      real(kind=4) nstark,ndopp,nmerge,inglis
      real(kind=4) kappa,kapmin,kappa0,kapcen
      real(kind=4) kappa0red,kappared,kappa0blue,kappablue
      real(kind=4) oldkappa,newkappa,dkappa
      real(kind=8) dopph(maxd)

      real(kind=dp)  wave,wcon,wmerge,wshift,contx(26,17),wtail
      real(kind=dp)  wl,emergeh(maxd)
      real(kind=dp)  bluecut,wlplus1,wlplus2,redcut,wlminus1,wlminus2,vacair
      real(kind=dp)  ehyd(100),conth(15),alphahyd(99)
!     default implicit rule
      real(kind=dp) frelin, freq, epsil,  edgeblue, adamp, dopwl, vvoigt
      integer(kind=4)  nbuff1, nbuff2, nbuff3, ixwl
      real(kind=dp) dnbuff, xsectg, tail
!---------- External ----------------------------------------------!
      real(kind=dp) voigt !

      logical whilep, anotherp

  conth = (/ 109678.764,  27419.659,  12186.462,  6854.871,  4387.113,  &  !1.00
     &         3046.604,   2238.320,   1713.711,  1354.044,  1096.776,  &
     &          906.426,    761.650,    648.980,   559.579,   487.456/)

  contx= 0.0d0
  contx(1:10,1) = (/109678.764,  27419.659,  12186.462,  6854.871,  4387.113,  &   !1.00
     &         3046.604,   2238.320,   1713.711,  1354.044,  1096.776/)
  contx(1:6,2) =(/198310.760,  38454.691,  32033.214,  29223.753,27175.760,15073.868/) !2.00
  contx(1:6,3) =(/438908.850, 109726.529,  48766.491,  27430.925,17555.715,12191.437/) !2.01
  contx(1:10,4) =(/90883.840,  90867.420,  90840.420,  90820.420,90804.000,  &   !6.00
     &        90777.000,  80691.180,  80627.760,  69235.820,69172.400/)
  contx(1:4,6)= (/61671.020,  39820.615,  39800.556,  39759.842/) !12.00
  contx(1:2,8) =(/48278.370,  48166.309/) !13.00
  contx(1:10,10) = (/66035.000,  65957.885,  65811.843,  65747.550,65670.435,  &   !14.00
     &        65524.393,59736.150,59448.700,50640.630,50553.180/)

   itemp1 = 0
!
!

  if(itemp.ne.itemp1) then

    ehyd(1)=0.
    ehyd(2)=82259.105
    ehyd(3)=97492.302
    ehyd(4)=102823.893
    ehyd(5)=105291.651
    ehyd(6)=106632.160
    ehyd(7)=107440.444
    ehyd(8)=107965.051

    do  n=9,100
     ehyd(n)=109678.764d0-109677.576d0/n**2
    end do

    do n=1,99
      alphahyd(n)=1.d7/(ehyd(n+1)-ehyd(n))
    end do

    do k=1,nrhox
!     empirical
      inglis=1600./xne(k)**(2./15.)
      nmerge=inglis-1.5
      emerge(k)=109737.312/nmerge**2
      emergeh(k)=109677.576/nmerge**2
    end do
      itemp1=itemp
  end if

      rewind 19
!
!--- open big loop over all ntle lines ---!
  do iline= 1, n19

    read(19)wl,gf,nblo,nbup,nelion,type,ncon,nelionx,gammar,gammas,gammaw,nbuff

!--- Note that He lines do not work yet! -----
!!  ---- the following types are considered:
!
!     IF(TYPE.EQ.2)    coronal lines not implemented
!
!     IF(TYPE.LT.-2) OR IF(TYPE.EQ.0) or IF(TYPE.EQ.3) -> normal line
!
!      IF(TYPE.EQ.-1)GO TO 600 Hydrogen line
!      IF(TYPE.EQ.-2)GO TO 600 Hydrogen line
!
!      IF(TYPE.EQ.1)GO TO 700 autoionizing line
!--------------------------------------------------------------------

    if ((type .eq. 0 ) .or. (type .lt.-2) .or. (type .eq. 3)) then

!------     normal line ------------------------

      kappa0=cgf*xnfdop(nelion)
      if(kappa0 .ge. continuum(nbuff)) then
      ifstrong=0
!     since cont=cont/1000, cont*500 is 1/2 the actual continuum.  lines are
!     strong if the line center is 1/2 the continuum.
      if( kappa0 .gt. 500.*continuum(nbuff)) ifstrong=1

      mlines=mlines+1
      wcon=0.
      wtail=0.

      if(ncon.gt.0)then
        wcon=1.e7/(contx(ncon,nelionx)-emerge(j))
        wtail=1.e7/(contx(ncon,nelionx)-emerge(j)-500.)
      endif

      adamp=(gammar+gammas*xne(j)+gammaw*txnxn(j))/dopple(nelion)
      kapcen=kappa0*voigt(0.0d0,adamp)
      dopwl=dopple(nelion)*wl

        if(wl .le. wlend) then

!------     red wing
          minred=max(1,nbuff)
          maxred=min(nbuff+1000,length)
          if(ifstrong .eq. 1) maxred=length
          wave=wbegin*ratio**(minred-1)

          ibuff = minred -1
          whilep = .true.

          do while ((ibuff .lt. maxred) .and. whilep )
            ibuff = ibuff +1
            if(wave.ge.wcon) then
              vvoigt=abs(wave-wl)/dopwl
              kappa=kappa0*voigt(vvoigt,adamp)
              if(wave.lt.wtail)kappa=kappa*(wave-wcon)/(wtail-wcon)
              buffer(ibuff)=buffer(ibuff)+kappa

              if(kappa.lt.continuum(ibuff)) then
                whilep = .false.
              else
                wave=wave*ratio
              end if
            end if

        end do

       end if
       whilep = .false.

       if (wl .le. wlend) then
         if(minred.ne.1) then
          if(wl.ge.wlbeg) then
            whilep = .true.
          end if
         end if
       else
          whilep = .true.
       end if

       if (whilep) then

!------     blue wing
         ibuff=min(length+1,nbuff)
         maxblue=ibuff-1
         if(ifstrong.eq.0) maxblue=min(maxblue,1000)

         wave=wbegin*ratio**(ibuff-1)
         i = 0

         do while ((i .lt. maxblue) .and. (whilep))
            i= i+1
            ibuff=ibuff-1
            wave=wave/ratio

            if(wave .ge. wcon) then
              vvoigt=abs(wave-wl)/dopwl
              kappa=kappa0*voigt(vvoigt,adamp)
              if(wave.lt.wtail)kappa=kappa*(wave-wcon)/(wtail-wcon)
              buffer(ibuff)=buffer(ibuff)+kappa

              if(kappa.lt.continuum(ibuff)) then
                whilep = .false.
              end if
            end if

         end do !end do while

       end if
      end if !top if

    else if ((type .eq. -1) .or. (type .eq. -2) ) then

!-----  hydrogen line

      kappa0=cgf*xnfdop(1)
      if(kappa0.ge. continuum(nbuff)) then

        dopph(j)=dopple(1)
!--     deuterium

        if(type.eq.-2)dopph(j)=dopph(j)/1.4142

        mlines=mlines+1

!      if(ncon.eq.0)go to 620
!      if(nbup.eq.nblo+1)go to 620
!      if(nbup.eq.nblo+2)go to 630
!------------------------------------------------------------------------------
!----- three conditions 1) non of them are true, 2) either 1st or 2nd is true, 3) 3rd is true:

      if((ncon .ne. 0) .and. (nbup .ne. nblo+1) .and. (nbup .ne. nblo+2)) then

        wshift=1.0d7/(conth(ncon)-109677.576d0/81.d0**2)
        wmerge=1.0d7/(conth(ncon)-emergeh(j))

        if(wmerge.lt.0.)wmerge=wshift+wshift

        wcon=max(wshift,wmerge)
        wtail=1.0d7/(1.0d7/wcon-500.0d0)
        wcon=min(wshift+wshift,wcon)

        if(wtail .lt. 0.) wtail=wcon+wcon
        wtail=min(wcon+wcon,wtail)

        whilep = .false.

        if(wl  .le. wlend) then
!----------------------------------------------------------------------
!     computes the first 20 points, then every 20th point and log interpolates
!     the 19 missing points
!
!     red wing
!----------------------------------------------------------------------

          if(wlbeg .gt. alphahyd(nblo)) then
            whilep = .true.
          else

            redcut=1.d7/(109678.764d0-109677.576d0/(nbup-0.8d0)**2-ehyd(nblo))
            wlminus1=1.d7/(ehyd(nbup-1)-ehyd(nblo))
            wlminus2=1.d7/(ehyd(nbup-2)-ehyd(nblo))
            kappa0red=kappa0*hfnm(nblo,nbup-2)/hfnm(nblo,nbup)/ &
     &                (ehyd(nbup-2)-ehyd(nblo))*(ehyd(nbup)-ehyd(nblo))

            minred=max(1,nbuff)
            wave=wbegin*ratio**(minred-1)
            ib=0

            ibuff = minred-1
            whilep = .true.

            do while ((ibuff .lt. length) .and. whilep )

              ibuff = ibuff + 1
              if(wave.lt.wcon) then    ! 1
              ! nothing
              else
                if(wave.gt.wlminus1) then ! 2
                  whilep = .false.
                else

                  ib=ib+1
                  if((ib.gt.21) .and. (mod(ib,20) .ne. 1)) then !3
                    wave = wave*ratio
                  else

                    kappa=kappa0*hprof4(nblo,nbup,j,wl,wave-wl,dopph)

                    if(wave.lt.wtail)kappa=kappa*(wave-wcon)/(wtail-wcon)

                    if(wave.gt.redcut)then
                      kappared=kappa0red*hprof4m(nblo,nbup-2,j,wl,wave-wlminus2,dopph)
                      if(wave.lt.wtail)kappared=kappared*(wave-wcon)/(wtail-wcon)
                      if(kappared.ge.kappa) whilep = .false.
                    endif

                    if (whilep) then !4
                      buffer(ibuff)=buffer(ibuff)+kappa
                      newkappa=log(max(kappa,1.e-20))

                      if(ib.gt.21)then
                         dkappa=(newkappa-oldkappa)/20.
                         do  k=1,19
                           buffer(ibuff-k)=buffer(ibuff-k)+exp(newkappa-dble(k)*dkappa)
                         end do
                      endif

                      oldkappa=newkappa

                      if(kappa.lt.continuum(ibuff)) then
                        whilep = .false.
                      end if
                    end if ! 4

                    if (whilep) wave = wave*ratio

                  end if !3
                end if ! 2
              end if! 1

            end do ! while do loop

            whilep = .false.
           end if

        end if
!------- this is to decide if the blue wing has to be done :
        if (wl .le. wlend ) then
          if (minred .ne. 1) then
            if(wl.ge.wlbeg) whilep = .true.
          end if
        else
          whilep = .true.

        end if

        if (whilep) then
!     blue wing

          bluecut=1.d7/(109678.764d0-109677.576d0/(nbup+0.8d0)**2- &
     &             ehyd(nblo))
          wlplus1=1.d7/(ehyd(nbup+1)-ehyd(nblo))
          wlplus2=1.d7/(ehyd(nbup+2)-ehyd(nblo))
          kappa0blue=kappa0*hfnm(nblo,nbup+2)/hfnm(nblo,nbup)/     &
     &            (ehyd(nbup+2)-ehyd(nblo))*(ehyd(nbup)-ehyd(nblo))
          ibuff=min(length+1,nbuff)
          maxblue=ibuff-1
          wave=wbegin*ratio**(ibuff-1)

          i = 0
          do while ((i .lt. maxblue) .and. (whilep)) !614
            i = i + 1
            ibuff=ibuff-1
            wave=wave/ratio

           if((wave.lt.wcon) .or. wave .lt. wlplus1) then ! 1
             whilep = .false.

           else

             if(( i .gt. 21) .and. (mod(i,20) .ne. 1))  then ! 2

             else

               kappa=kappa0*hprof4(nblo,nbup,j,wl,wave-wl,dopph)
               if(wave.lt.wtail)kappa=kappa*(wave-wcon)/(wtail-wcon)

               if(wave .lt. bluecut)then
                 kappablue=kappa0blue*hprof4p(nblo,nbup+2,j,wl,wave-wlplus2,dopph)
                 if(wave.lt.wtail)kappablue=kappablue*(wave-wcon)/(wtail-wcon)
                 if(kappablue.gt.kappa) whilep = .false.

               endif

               if (whilep ) then
                 buffer(ibuff)=buffer(ibuff)+kappa
                 newkappa=log(max(kappa,1.e-20))

                 if(i.gt.21)then
                   dkappa=(newkappa-oldkappa)/20.

                  do  k=1,19
                   buffer(ibuff+k)=buffer(ibuff+k)+exp(newkappa-dble(k)*dkappa)
                  end do
                 endif

               oldkappa=newkappa

               if(kappa.lt.continuum(ibuff)) whilep = .false.

               end if
             end if !2
            end if !  1
          end do !do while loop 614

        end if

!--- second contitions

      else if ((ncon .eq. 0) .or. (nbup .eq. nblo+1) ) then

!--------------------------------------------------------
!     alpha presumed isolated lines and beta blue wings
!---------------------------------------------------------
       if(wl .le. wlend) then
!----     red wing-------!
         minred=max(1,nbuff)
         ib=0
         wave=wbegin*ratio**(minred-1)

         whilep = .true.
         ibuff = minred-1

         do while ((ibuff .lt. length) .and. whilep)
          ibuff = ibuff +1
          ib=ib+1
          if((ib.gt.21 ) .and. (mod(ib,20) .ne.1)) then
            wave = wave*ratio
          else

          kappa=kappa0*hprof4(nblo,nbup,j,wl,wave-wl,dopph)
          buffer(ibuff)=buffer(ibuff)+kappa
          newkappa=log(max(kappa,1.e-20))

          if(ib.gt.21)then
            dkappa=(newkappa-oldkappa)/20.
            do  k=1,19
             buffer(ibuff-k)=buffer(ibuff-k)+exp(newkappa-dble(k)*dkappa)
            end do
          endif

          oldkappa=newkappa

            if(kappa.lt.continuum(ibuff)) then
               whilep = .false.
            else
               wave = wave*ratio
            end if
          end if
         end do
       end if

       whilep = .false.

       if (wl .le. wlend) then
         if(minred.ne.1) then
          if(wl.ge.wlbeg) then
            whilep = .true.
          end if
         end if
       else
          whilep = .true.
       end if

        if (whilep ) then
!-------- alpha or beta blue wing
          ibuff=min(length+1,nbuff)

          maxblue=ibuff-1
          wave=wbegin*ratio**(ibuff-1)

          i = 0

          do while ((i .lt. maxblue) .and. whilep)
            i = i +1
            ibuff=ibuff-1
            wave=wave/ratio

            if((i.gt.21) .and. (mod(i,20).ne. 1)) then

            else

              kappa=kappa0*hprof4(nblo,nbup,j,wl,wave-wl,dopph)
              buffer(ibuff)=buffer(ibuff)+kappa
              newkappa=log(max(kappa,1.e-20))
              if(i.gt.21)then
                dkappa=(newkappa-oldkappa)/20.
                do k=1,19
                  buffer(ibuff+k)=buffer(ibuff+k)+exp(newkappa-dble(k)*dkappa)
                end do
              endif

              oldkappa=newkappa

              if(kappa.lt.continuum(ibuff)) then
               whilep = .false.
              end if
            end if
          end do
        end if

!--- third condition:

      else if (nbup .eq. nblo+2) then
       anotherp = .false.
!-----------------------------------------------------------
!     beta lines red wing

      if( wl .le. wlend) then
!
!-----  red wing

        minred=max(1,nbuff)
        ib=0
        wave=wbegin*ratio**(minred-1)

        ibuff = minred -1
        whilep = .true.

        do while ((ibuff .lt. length) .and. whilep) ! 631
          ibuff = ibuff +1

         if(wave.gt.alphahyd(nblo)) then ! (1)
            anotherp = .true.

           if(nbuff.ge.1) then !(2)
             ibuff=nbuff
             maxblue=ibuff-1
             wave=wbegin*ratio**(ibuff-1)

             i = 0
             do while ((i.lt.maxblue) .and. whilep  )
                i = i +1
                ibuff=ibuff-1
                wave=wave/ratio

                if((i.gt.21) .and. (mod(i,20).ne.1)) then
                   whilep = .false.
                else

                 kappa=kappa0*hprof4(nblo,nbup,j,wl,wave-wl,dopph)
                 buffer(ibuff)=buffer(ibuff)+kappa
                 newkappa=log(max(kappa,1.e-20))

                 if(i.gt.21)then
                  dkappa=(newkappa-oldkappa)/20.
                  do  k=1,19
                      buffer(ibuff+k)=buffer(ibuff+k)+ &
     &                              exp(newkappa-dble(k)*dkappa)
                  end do
                 endif

                 oldkappa=newkappa
                 if(kappa.lt.continuum(ibuff)) then
                  whilep = .false.
                 end if
                end if
              end do
           else
             whilep = .false.
           end if ! (2)
!
           whilep = .false.
! else from (1)
         else
           ib=ib+1

           if((ib.gt.21) .and. mod(ib,20) .ne. 1) then !(3)
              wave = wave*ratio
           else

             kappa=kappa0*hprof4(nblo,nbup,j,wl,wave-wl,dopph)
             buffer(ibuff)=buffer(ibuff)+kappa
             newkappa=log(max(kappa,1.e-20))

             if(ib.gt.21)then
               dkappa=(newkappa-oldkappa)/20.
               do k=1,19
                  buffer(ibuff-k)=buffer(ibuff-k)+exp(newkappa-dble(k)*dkappa)
               end do
             endif

             oldkappa=newkappa

             if (kappa .lt. continuum(ibuff)) then
               whilep = .false.
             else
               wave=wave*ratio
             end if

           end if ! (close (3))
         end if ! close(1)
        end do
!---------------
       end if

       if (wl .gt. wlend  ) then
        whilep = .true.
       else
        if ((minred .ne. 1 ) .and. (wl .ge. wlbeg)) whilep = .true.
       end if

       if (anotherp) whilep = .false.

       if (whilep ) then
!-------- (alpha or) beta blue wing
          ibuff=min(length+1,nbuff)

          maxblue=ibuff-1
          wave=wbegin*ratio**dble(ibuff-1)

          i = 0

          do while ((i .lt. maxblue) .and. whilep)
            i = i + 1
            ibuff=ibuff-1
            wave=wave/ratio

            if((i.gt.21) .and. (mod(i,20).ne. 1)) then
            !  i = i +1
            else

              kappa=kappa0*hprof4(nblo,nbup,j,wl,wave-wl,dopph)
              buffer(ibuff)=buffer(ibuff)+kappa
              newkappa=log(max(kappa,1.e-20))

              if(i.gt.21)then
                dkappa=(newkappa-oldkappa)/20.
                do k=1,19
                  buffer(ibuff+k)=buffer(ibuff+k)+exp(newkappa-dble(k)*dkappa)
                end do
              endif

              oldkappa=newkappa

              if(kappa.lt.continuum(ibuff)) then
               whilep = .false.
              end if
            end if
          end do

      end if

      end if ! three conditions

      end if ! -- top if

!--------------------------------------------------------------!
     else if (type .eq. 1) then
!-------------------------------
!     autoionizing line
!-------------------------------

       kappa0=bshore*g*xnfpel(nelion)
       if(kappa0.ge.continuum(nbuff)) then

         mlines=mlines+1
         frelin=2.99792458e17/wl

         if(wl .le. wlend) then
!          red wing

           minred=max(1,nbuff)
           freq=2.99792458e17/(wbegin*ratio**(minred-1))

           ibuff = minred-1
           whilep = .true.
           do while (whilep .and. (ibuff .lt. length) )
              ibuff = ibuff +1

              epsil=2.*(freq-frelin)/gammar
              kappa=kappa0*(ashore*epsil+bshore)/(epsil**2+1.)/bshore
              buffer(ibuff)=buffer(ibuff)+kappa
              freq=freq/ratio

              if (kappa .lt. continuum(ibuff)) then
                whilep = .false.
              end if

           end do
         end if

         whilep= .false.

         if (wl .le. wlend) then
           if (nbuff .ne. 1 ) then
             whilep = .true.
           end if
         else
           whilep = .true.
         end if

         if (whilep ) then
!          blue wing
           ibuff=min(length+1,nbuff)
           maxblue=ibuff-1
           freq=2.99792458e17/(wbegin*ratio**(ibuff-1))

           i = 0
           whilep = .true.
           do while ((i .lt. maxblue) .and. whilep)
             i = i +1
             ibuff=ibuff-1
             freq=freq*ratio
             epsil=2.*(freq-frelin)/gammar
             kappa=kappa0*(ashore*epsil+bshore)/(epsil**2+1.)/bshore
             buffer(ibuff)=buffer(ibuff)+kappa
            if(kappa.lt.continuum(ibuff)) then
             whilep = .false.
            end if
           end do

         end if
       end if

!----------------------------------------------------------------!
    else if (type .eq. 2 ) then
!------------------------------
!   coronal line
!   500 go to 900
!----------------------------------------------------------------!
    else
!----------------------------------
!     merged continuum
!     edge wavelengths are in vacuum
!-----------------------------------
      wshift=1.d7/(1.d7/wl-109737.312d0/nlast**2)
      wmerge=1.d7/(1.d7/wl-emerge(j))

      if(nelion.eq.1)then
        wshift=1.d7/(1.d7/wl-109677.576d0/nlast**2)
        wmerge=1.d7/(1.d7/wl-emergeh(j))
      endif

      if(wmerge .lt. 0.0d0) wmerge=wshift+wshift

      wmerge=max(wmerge,wshift)
      wmerge=min(wshift+wshift,wmerge)
      wtail=1.d7/(1.d7/wmerge-500.)

      if(wtail .lt. 0.0d0) wtail=wmerge+wmerge

      wtail=min(wmerge+wmerge,wtail)
      ixwl= log(wl)/ratiolg
      edgeblue=exp(ixwl*ratiolg)

      if(edgeblue .gt. wl) ixwl=ixwl-1

      nbuff1=ixwl+1-ixwlbeg+1
      ixwl= log(wmerge)/ratiolg+.5
      nbuff2=ixwl-ixwlbeg+1
      ixwl= log(wtail)/ratiolg+.5
      nbuff3=ixwl-ixwlbeg+1

      if(nbuff1.le.length) then

        if(nbuff3.ge.1) then

          dnbuff=dble(nbuff3-nbuff2)
          nbuff1=max(nbuff1,1)
          xsectg=dble(gf)
          kappa=xsectg*xnfpel(nelion)
          tail=1.

         do ibuff=nbuff1, min(nbuff3,length)
          if(ibuff.gt.nbuff2) tail=(nbuff3-ibuff)/dnbuff
          buffer(ibuff) = buffer(ibuff)+kappa*tail
         end do
        end if
      end if
    end if

!
   end do  ! 900 loop

  end

  double precision function voigt(v,a)
!
!
      use types
      use dfcomm

      implicit none
      real(kind=dp) v, a
      real(kind=dp) aa, vv, u, aau, vvu
      integer iv, ia
!     fast voigt

      if(a .ge. 0.20d0) then

        if((a .gt. 1.4d0) .or. ((a+v) .gt. 3.2d0)) then

        aa=a*a
        vv=v*v
        u=(aa+vv)*1.4142d0
        voigt=a*.79788d0/u
        if(a.gt. 100.0d0)return
        aau=aa/u
        vvu=vv/u
        voigt=voigt* &
     &  ((((aau-10.*vvu)*aau+5.*vvu*vvu)/u+vvu-aau/3.)*3./u+1.)

        return
        else

          iv=v*200.+1.5
          ia=a*200.+1.5
          voigt=atab(ia,iv)

        end if

      else

        if(v .le. 10.) then
          iv=v*200.+1.5
          voigt=(h2tab(iv)*a+h1tab(iv))*a+h0tab(iv)
        else
          voigt=0.5642d0*a/(v**2)
        end if

      end if

      return
  end