text:
module root_module
   use precision_module, only : wp
   use error_module
   use constants_module

   implicit none
   public

contains

function root(cf,limitL,limitR)
   implicit none
   real(wp) :: cf(5), limitL, limitR
   real(wp) :: rcf1(2),rcf2(3),rcf3(4), rcf4(5)
   real(wp) :: root

   if     (cf(1).eq.0.0d0 .and. &
           cf(2).eq.0.0d0 .and. &
           cf(3).eq.0.0d0) then !linear  potentials
       rcf1=cf(4:5) !get rid of first three elements
       rcf1=rcf1/rcf1(1) !get reduced polynomial
       root=poly1_root(rcf1,limitL,limitR) !find root of cubic polynomial

   elseif (cf(1).eq.0.0d0 .and. &
           cf(2).eq.0.0d0) then !quadratic polynomial
       rcf2=cf(3:5) !get rid of first element
       rcf2=rcf2/rcf2(1) !get reduced polynomial
       root=poly2_root(rcf2,limitL,limitR) !find root of cubic polynomial

   elseif (cf(1).eq.0.0d0) then !cubic polynomial
       rcf3=cf(2:5) !get rid of the first element
       rcf3=rcf3/rcf3(1) !get reduced polynomial
       root=poly3_root(rcf3,limitL,limitR) !find root of cubic polynomial

   else !quiartic polynomial
       rcf4=cf/cf(1) ! get reduced polynomial
       root=poly4_root(rcf4,limitL,limitR) !find root of quartic polynomial
   endif
end function root

!---------------------------------------------------------------------
!---------------------------------------------------------------------
function poly1_root(rcf,limitL,limitR)
   implicit none
   real(wp) :: rcf(2), limitL, limitR
   real(wp) :: x, poly1_root

     x=-rcf(2)
     if (x.gt.limitL .and. x.lt.limitR) then
         poly1_root=x
     else
        call error_message(8)!error code in parentheses
     endif
end function poly1_root

!---------------------------------------------------------------------
!---------------------------------------------------------------------
function poly2_root(rcf,limitL,limitR)
   implicit none
   real(wp) :: rcf(3), limitL, limitR
   real(wp) :: b,c,dd,x1,x2, poly2_root
   integer          :: i

     b=rcf(2)
     c=rcf(3)
     dd=b**2.0d0-4.0d0*c !calculate determinant of reduced polynomial

     i=0 !number of roots
     if (dd.gt.0.0d0) then !two real root
         x1=(-b+sqrt(dd))/2
         x2=(-b-sqrt(dd))/2
         if (x1.gt.limitL .and. x1.lt.limitR) then
            poly2_root=x1
            i=i+1
         elseif (x2.gt.limitL .and. x2.lt.limitR) then
            poly2_root=x1
            i=i+1
         endif
     elseif (dd.eq.0.0d0) then
         x1=-b/2
         if (x1.gt.limitL .and. x1.lt.limitR) then
            poly2_root=x1
            i=i+1
         endif
     endif

   if (i.gt.1) then !multiple crossings in the region of interest
        call error_message(9)
   elseif (i.eq.0) then !no crossings in the region of interest
        call error_message(8)
   endif
end function poly2_root

!---------------------------------------------------------------------
!---------------ROOTS ARE COMPUTED USING CARDANO METHOD---------------
!----see details at http://mathworld.wolfram.com/CubicFormula.html-----
function poly3_root(rcf,limitL,limitR)
   implicit none
   real(wp) :: rcf(4), limitL, limitR
   real(wp) :: poly3_root
   real(wp) :: y1, y2, y3, factor, phi
   real(wp) :: q, r, det, s, t
   complex(wp)   :: y2i, y3i
   integer          :: i

!   rcf=cf/cf(1) ! coefficient in front of x^3 is set to one
   q=(3.0d0*rcf(3)-rcf(2)**2.0d0)/9.0d0
   r=(-2.0d0*rcf(2)**3.0d0+9.0d0*rcf(2)*rcf(3)-27.0d0*rcf(4))/54.0d0
   det=q**3.0d0+r**2.0d0

   if (det.gt.0.0d0) then !Only one real root
    s=r+sqrt(det)
     if (s.lt.0.0d0) then
        s=-((ABS(s))**(1.0d0/3.0d0))
     else
        s=s**(1.0d0/3.0d0)
     endif
    t=r-sqrt(det)
     if (t.lt.0.0d0) then
        t=-(ABS(t)**(1.0d0/3.0d0))
     else
        t=t**(1.0d0/3.0d0)
     endif
    y1=-rcf(2)/3.0d0+s+t

   !!!!!!!!!!!!!SKIP COMPEX ROOTS!!!!!!!!!!!!!!!!!!!!!!!!!
   !   re=-rcf(2)/3.0d0-(s+t)/2.0d0
   !   im=(sqrt(3.0d0)/2.0d0)*(s-t)
   !  y2i=cmplx(re,im)
   !  y3i=cmplx(re,-im)
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   elseif (det.eq.0.0d0) then ! three real roots (two are equal)
     if (s.lt.0.0d0) then
        s=-((abs(r))**(1.0d0/3.0d0))
     else
        s=r**(1.0d0/3.0d0)
     endif
    y1=-rcf(2)/3.0d0+s*2.0d0
    y2=-rcf(2)/3.0d0-s
    y3=y2

   elseif (det.lt.0.0d0) then ! three real unequal roots
      phi=r/sqrt(abs(q**3.0d0))
      phi=1.0d0/(cos(phi))
      factor=2.0d0*sqrt(-q)
    y1=factor*cos(phi/3.0d0)-rcf(2)/3.0d0
    y2=factor*cos((phi+2.0d0*pi)/3.0d0)-rcf(2)/3.0d0
    y3=factor*cos((phi+4.0d0*pi)/3.0d0)-rcf(2)/3.0d0
   endif

   i=0 !number of roots
   if (det.gt.0.0d0) then !
     if (y1.gt.limitL .and. y1.lt.limitR) then
        poly3_root=y1
        i=i+1
     else
        call error_message(7)!error code in parentheses
     endif
   elseif (det.eq.0.0d0) then !
     if (y1.gt.limitL .and. y1.lt.limitR) then
        poly3_root=y1
        i=i+1
     elseif (y2.gt.limitL .and. y2.lt.limitR) then
        poly3_root=y2 ! y3 is equal to y2, so there's no need to track y3
        i=i+1
     endif
   elseif (det.lt.0.0d0) then!
     if (y1.gt.limitL .and. y1.lt.limitR) then
        poly3_root=y1
        i=i+1
     elseif (y2.gt.limitL .and. y2.lt.limitR) then
        poly3_root=y2
        i=i+1
     elseif (y3.gt.limitL .and. y3.lt.limitR) then
        poly3_root=y3
        i=i+1
     else
        call error_message(8)!error code in parentheses
     endif
   endif

   if (i.gt.1) then !you have multiple crossings in region of interest
        call error_message(9)!error code in parentheses
   endif
end function poly3_root

!---------------------------------------------------------------------
!---------------ROOTS ARE COMPUTED USING FERRARI'S METHOD-------------
!----see details at https://en.wikipedia.org/wiki/Quartic_function----
function poly4_root(rcf,limitL,limitR)
   implicit none
   real(wp) :: rcf(5), limitL, limitR
   real(wp) :: poly4_root
   real(wp) :: a,b,c,d,e
   real(wp) :: p,q, dd,dd0,dd1,T,tmp
   complex(wp)   :: S,x(4)
   integer          :: i,j

   a=rcf(1) !simplify nottion for polynomial coefficients
   b=rcf(2)
   c=rcf(3)
   d=rcf(4)
   e=rcf(5)

   dd=256d0*e**3.0d0-192d0*b*d*e**2.0d0-128d0*c**2.0d0*e**2.0d0
   dd=dd+144*c*d**2.0d0*e-27d0*d**4.0d0+144*b**2.0d0*c*e**2.0d0
   dd=dd-6.0d0*b**2.0d0*d**2.d0*e-80.0d0*b*c**2.0d0*d*e
   dd=dd+18.0d0*b*c*d**3.0d0+16.0d0*c**4.0d0*e-4.0d0*c**3.0d0*d**2.0d0
   dd=dd-27.0d0*b**4.0d0*e**2.0d0+18.0d0*b**3.0d0*c*d*e
   dd=dd-4.0d0*b**3.0d0*d**3.0d0-4.0d0*b**2.0d0*c**3.0d0*e
   dd=dd+b**2.0d0*c**2.0d0*d**2.0d0

   dd0=c**2.0d0-3.0d0*b*d+12.0d0*e
   dd1=2.0d0*c**3.0d0-9.0d0*b*c*d+27.0d0*b**2.0d0*e
   dd1=dd1+27.0d0*d**2.0d0-72.0d0*c*e
   p=(8.0d0*c-3.0d0*b**2.0d0)/8
   q=(b**3.0d0+8.0d0*d-4.0d0*b*c)/8

!   if (dd.lt.0) then

     T=(sqrt(-27.0d0*dd)+dd1)/2
       if (T.lt.0.0d0) then
           T=-(abs(T)**(1.0d0/3.0d0))
       else
           T=T**(1.0d0/3.0d0)
       endif

     S=-2.0d0/3.0d0*p+(T+dd0/T)/3.0d0
     S=0.5d0*sqrt(S)

!print *,'a,b,c,d,e',a,b,c,d,e
!print *,'dd,dd1,dd0',dd,dd1,dd0
!print *,'p,q,T,S',p,q,T,S

   x(1)=-4.0d0*S**2.0d0-2.0d0*p+q/S
   x(1)=-b/4.0d0-S+0.5d0*sqrt(x(1))

   x(2)=-4.0d0*S**2.0d0-2.0d0*p+q/S
   x(2)=-b/4.0d0-S-0.5d0*sqrt(x(2))

   x(3)=-4.0d0*S**2.0d0-2.0d0*p-q/S
   x(3)=-b/4.0d0+S+0.5d0*sqrt(x(3))

   x(4)=-4.0d0*S**2.0d0-2.0d0*p-q/S
   x(4)=-b/4.0d0+S-0.5d0*sqrt(x(4))

!print *,'ROOOOOTS',x(1),x(2),x(3),x(4)
!print *,'IMAG PART',aimag(x(3)),aimag(x(4))
   i=0
   do j=1,4
      if (aimag(x(j)).eq.0.0d0) then !real root
          tmp=real(x(j))
!print *,'tmp',tmp
          if (tmp.gt.limitL .and. tmp.lt.limitR) then
              poly4_root=tmp
              i=i+1
          endif
      endif
   enddo
!print *,'iiiiiiii',i
!print *,'root',poly4_root

   if (i.gt.1) then !multiple crossings in the region of interest
        call error_message(9)
   elseif (i.eq.0) then !no crossings in the region of interest
        call error_message(8)
   endif
end function poly4_root
!---------------------------------------------------------------------
end module root_module

max_stars_repo_name: svarganov/NAST

max_issues_repo_name: svarganov/NAST

max_forks_repo_name: svarganov/NAST


text:
What are the top seven things you would tell an outoftown visitor to see in Davis?

 1. Pence Gallery
 1. John Natsoulas Art Gallery
 1. Mondavi Center
 1. Buehler Alumni & Visitors Center
 1. Davis Food Coop
 1. Village Homes
 1. Arboretum

Other candidates: The Domes, Stevenson Bridge, Delta of Venus, Memorial Union, Whymcycles?, companies?, churches?

20060417 10:51:32 nbsp Mondavi, The Domes, Toad Tunnel, The Eggheads,  Wed. Farmers Market in Summer, Whole Earth Fest (May Only though..), Dutton Fountain The saddest fountain   Those werent in any particular orderUsers/JulienBiewerElstob

20070227 15:16:48 nbsp Julien, those are all at ucd except for the Farmers Market and the Toad Tunnel... Stop being a poser. Users/NickSchmalenberger

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE PDLAUUM( UPLO, N, A, IA, JA, DESCA )
*
*  -- ScaLAPACK auxiliary routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            IA, JA, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * )
      DOUBLE PRECISION   A( * )
*     ..
*
*  Purpose
*  =======
*
*  PDLAUUM computes the product U * U' or L' * L, where the triangular
*  factor U or L is stored in the upper or lower triangular part of
*  the distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1).
*
*  If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
*  overwriting the factor U in sub( A ).
*  If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
*  overwriting the factor L in sub( A ).
*
*  This is the blocked form of the algorithm, calling Level 3 PBLAS.
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Arguments
*  =========
*
*  UPLO    (global input) CHARACTER*1
*          Specifies whether the triangular factor stored in the
*          distributed matrix sub( A ) is upper or lower triangular:
*          = 'U':  Upper triangular
*          = 'L':  Lower triangular
*
*  N       (global input) INTEGER
*          The number of rows and columns to be operated on, i.e. the
*          order of the triangular factor U or L. N >= 0.
*
*  A       (local input/local output) DOUBLE PRECISION pointer into the
*          local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
*          On entry, the local pieces of the triangular factor L or U.
*          On exit, if UPLO = 'U', the upper triangle of the distributed
*          matrix sub( A ) is overwritten with the upper triangle of the
*          product U * U'; if UPLO = 'L', the lower triangle of sub( A )
*          is overwritten with the lower triangle of the product L' * L.
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      DOUBLE PRECISION   ONE
      PARAMETER          ( ONE = 1.0D+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J, JB, JN
*     ..
*     .. External Subroutines ..
      EXTERNAL           PDGEMM, PDLAUU2, PDTRMM, PDSYRK
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICEIL
      EXTERNAL           ICEIL, LSAME
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( N.EQ.0 )
     $   RETURN
*
      JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 )
      IF(  LSAME( UPLO, 'U' ) ) THEN
*
*        Compute the product U * U'.
*
*        Handle first block separately
*
         JB = JN-JA+1
         CALL PDLAUU2( 'Upper', JB, A, IA, JA, DESCA )
         IF( JB.LE.N-1 ) THEN
            CALL PDSYRK( 'Upper', 'No transpose', JB, N-JB, ONE, A, IA,
     $                   JA+JB, DESCA, ONE, A, IA, JA, DESCA )
         END IF
*
*        Loop over remaining block of columns
*
         DO 10 J = JN+1, JA+N-1, DESCA( NB_ )
            JB = MIN( N-J+JA, DESCA( NB_ ) )
            I = IA + J - JA
            CALL PDTRMM( 'Right', 'Upper', 'Transpose',  'Non-unit',
     $                   J-JA, JB, ONE, A, I, J, DESCA, A, IA, J,
     $                   DESCA )
            CALL PDLAUU2( 'Upper', JB, A, I, J, DESCA )
            IF( J+JB.LE.JA+N-1 ) THEN
               CALL PDGEMM( 'No transpose', 'Transpose', J-JA, JB,
     $                      N-J-JB+JA, ONE, A, IA, J+JB, DESCA, A, I,
     $                      J+JB, DESCA, ONE, A, IA, J, DESCA )
               CALL PDSYRK( 'Upper', 'No transpose', JB, N-J-JB+JA, ONE,
     $                      A, I, J+JB, DESCA, ONE, A, I, J, DESCA )
            END IF
   10    CONTINUE
      ELSE
*
*        Compute the product L' * L.
*
*        Handle first block separately
*
         JB = JN-JA+1
         CALL PDLAUU2( 'Lower', JB, A, IA, JA, DESCA )
         IF( JB.LE.N-1 ) THEN
            CALL PDSYRK( 'Lower', 'Transpose', JB, N-JB, ONE, A, IA+JB,
     $                   JA, DESCA, ONE, A, IA, JA, DESCA )
         END IF
*
*        Loop over remaining block of columns
*
         DO 20 J = JN+1, JA+N-1, DESCA( NB_ )
            JB = MIN( N-J+JA, DESCA( NB_ ) )
            I = IA + J - JA
            CALL PDTRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', JB,
     $                   J-JA, ONE, A, I, J, DESCA, A, I, JA, DESCA )
            CALL PDLAUU2( 'Lower', JB, A, I, J, DESCA )
            IF( J+JB.LE.JA+N-1 ) THEN
               CALL PDGEMM( 'Transpose', 'No transpose', JB, J-JA,
     $                      N-J-JB+JA, ONE, A, I+JB, J, DESCA, A, I+JB,
     $                      JA, DESCA, ONE, A, I, JA, DESCA )
               CALL PDSYRK( 'Lower', 'Transpose', JB, N-J-JB+JA, ONE,
     $                      A, I+JB, J, DESCA, ONE, A, I, J, DESCA )
            END IF
   20    CONTINUE
      END IF
*
      RETURN
*
*     End of PDLAUUM
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text: Single Transferable Vote is the technical name for a Voting voting system known locally and more commonly as Choice Voting.  Preference Voting is yet another name.  Single Transferable Voting can be used in both singlewinner and multiwinner elections.  In multiwinner elections, it is a system of proportional representation.  For singlewinner elections, it reduces to the system more commonly known as Instant Runoff Voting.


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

* Vectorizer - stripmine, incorrect placement of temp stores of variant
*              expressions
*              iscbugs 911004a.f -- 911004f.f
*              same as kv02.f + use of pointers

	program p
	parameter (N=512)
	parameter (M=N/2)

	double precision x(N), y(N)
	double precision expect(N)

	pointer (px,idum)
	pointer (py,jdum)

	data x /N*1.0/
	data y /N*2.0/
	data (expect(i),i=1,N,2) / M*2.0 /
	data (expect(i),i=2,N,2) / M*1.0 /

	incx = 2
	incy = -2
	px = %loc(x)
	py = %loc(y)
	call dswap (N/2,px,incx,py,incy)
	call checkd(x, expect, N)
	end

      subroutine dswap (n,px,incx,py,incy)
      integer n, incx, incy
      double precision x, y
      pointer (px, x(n))
      pointer (py, y(n))
c
c     Swap vectors.
c     y <==> x
c
      integer i, ix, iy
      double precision t
c
   20 ix = 1
      iy = 1
      if (incx .lt. 0) ix = (-n+1)*incx + 1
      if (incy .lt. 0) iy = (-n+1)*incy + 1
c
c  need temps for computing the initial array addresses for x & y
c  these temps should be stored prior to the stripmine loop
c
      do 30 i = 1, n
         t = y(iy)
         y(iy) = x(ix)
         x(ix) = t
         ix = ix + incx
         iy = iy + incy
   30 continue
      return
      end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


module mo_profile
    USE omp_lib
    REAL(8) :: START, END

contains

integer function init()
    open  (130, file = "profile.log") ! , ACCESS='APPEND'
    init = 130
end function

subroutine close()
    close (130)
end subroutine

subroutine tic(prefix)
    character(len=*), optional, intent(in) :: prefix

    if (present(prefix)) then
        write (130, *) prefix
    endif
    START = omp_get_wtime ( )
end subroutine

subroutine toc()
    ! CHARACTER(256) outfile = 'profile.log'
    END = omp_get_wtime ()
    write (130, "(A10, 3F18.6)") "elapsed:", START, END, END - START
end subroutine

end module mo_profile


text:
       SUBROUTINE DECODE1
C*****THIS SUBROUTINE TRANSLATES AN ACCOS DATA FILE TO A CODEV LENS DECK
C
C** COMMON BLOCK ACCOSV DECLARATIONS
      CHARACTER*6 UNITS
      INTEGER*4 IORDER(66),IWT1,IWT2,IWT3,IWT4,IWT5
      REAL*4 WVW1,WVW2,WVW3,WVW4,WVW5,WT1,WT2,WT3,WT4,WT5,THETA(66)
      REAL*4 RN1
      REAL*8 N1(99),N2(99),N3(99),N4(99),N5(99),F2(225)
C** COMMON BLOCK BOTH DECLARATIONS
C BBC 9/20/84. ADDED THE VARIABLES NCONF(),ACVAL(),CPARM(), AND I11 TO HANDLE
C ALTERNATE CONFIGURATIONS.
      CHARACTER*4 SOLVE(40)
      CHARACTER*3 CPARM(30)
      CHARACTER*13 GLASS(0:225)
      CHARACTER*64 TITLE
      INTEGER*4 I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,M,ISTOP,IREF,ISURF
      INTEGER*4 NAPERT(225),NDECE(100),NSOLVE(40),NUMA,MESSAGE(20)
      REAL*4 DIF(66)
      REAL*8 ADE(100),BDE(100),CDE(100),XDE(100),YDE(100),CURVEY(0:225)
      REAL*8 THICKNS(0:225),RK(66),A(66),B(66),C(66),D(66),SOLVF1(40)
      REAL*8 F1(225),REF_OBJ_HT
      LOGICAL AFO,IC(66),CUF(66)
C** COMMON BLOCK ACONLY DECLARATIONS
      CHARACTER*4 TORIC(10),COMMAND(100),PIKUP(200),SHAPE(225),CASE(10)
      CHARACTER*4 SPECSF(225)
      CHARACTER*36 SUBMES1(20)
      CHARACTER*40 SUBMES2(20)
      CHARACTER*48 SUBMES3(20)
      INTEGER*4 INDEX(20),NTORIC(10),NCASE(10),NUMBER(10),IPIKTO(200)
      INTEGER*4 IPIKFR(200),NASP(66),J10,K1,NGRT(10),NSHAPE(225)
      REAL*4 PIKA(200),DDIST(10),RNUMBER(10),RIC1(225),RIC2(225)
      REAL*4 RCUF1(225),RCUF2(225)
      REAL*8 TCVR(10),PY,PUCY,RADIUS(0:225),PIKB(200),NAO,YAN
      LOGICAL OB(225)
C** LOCAL DECLARATIONS
      CHARACTER*4 CVAR		! JPL CHANGE.
      INTEGER*4 NMODELS,CFG1,REF_SURF
      CHARACTER*110 LINE
      INTEGER*4 DUMMY(100)
      LOGICAL SAME
C
      COMMON /ACCOSV/IORDER,WVW1,WVW2,WVW3,WVW4,WVW5,WT1,WT2,WT3,WT4,
     *        WT5,IWT1,IWT2,IWT3,IWT4,IWT5,THETA,F2,RN1,UNITS,N1,
     *        N2,N3,N4,N5
      COMMON /BOTH/I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,GLASS,MESSAGE,M,
     *        ISTOP,IREF,ISURF,AFO,NAPERT,NDECE,F1,DIF,ADE,BDE,CDE,XDE,
     *        YDE,CURVEY,THICKNS,RK,A,B,C,D,IC,CUF,NSOLVE,SOLVF1,SOLVE,
     *        TITLE,NUMA,REF_OBJ_HT
      COMMON /ACONLY/TCVR,PY,PUCY,INDEX,NTORIC,DDIST,SUBMES1,SUBMES2,
     *        SUBMES3,TORIC,NCASE,COMMAND,NUMBER,RADIUS,PIKA,PIKB,
     *        IPIKTO,IPIKFR,PIKUP,SHAPE,CASE,NASP,OB,SPECSF,NAO,YAN,
     *        J10,K1,NGRT,RIC1,RIC2,RCUF1,RCUF2,RNUMBER,NSHAPE
C
C******************************************************************************
C*****FORMAT STATEMENTS
C HAD TO CHANGE SOME OF THE FORMATS. VAX WAS ADDING A ZERO TO THE RIGHT END OF
C SOME OF THE INTEGERS BECAUSE AN EXTRA BLANK SPACE WAS READ. BBC 9/84.
  2   FORMAT(A110)
  3   FORMAT(5X,E15.6,1X,E15.6,3X,A13,3X,E8.6)
  4   FORMAT(5X,F15.6,1X,F15.6,3X,A13)
  5   FORMAT(1X,I3)
  6   FORMAT(1X,I3,3X,F9.6,2X,F9.6,2X,F9.6,2X,F9.6,2X,F9.6)
  7   FORMAT(19X,F13.6,1X,F13.6,3X,A13,3X,F8.6)
  8   FORMAT(1X,I3,1X,F13.5,F13.5,F13.5,F13.5,F13.5)
  9   FORMAT(19X,F13.6,1X,F13.6,3X,A13)
 10   FORMAT(2X,I3,4X,A4,17X,F6.3,16X,I4)
 12   FORMAT(1X,I3,4X,A4,F14.6)
 14   FORMAT(1X,I3,1X,A4,3X,A4,6X,F10.4,4X,F8.4)
 16   FORMAT(1X,I3,3X,A4,2X,F11.5,1X,F11.5,1X,F10.4,1X,F10.4,1X,F10.4)
 18   FORMAT(1X,I3,4X,A4,2X,I4,2X,F10.6,2X,F10.6)
 20   FORMAT(1X,I3,4X,A4,28X,F10.6)
 22   FORMAT(16X,F8.5,3X,F8.5,3X,F8.5,3X,F8.5,3X,F8.5)
 24   FORMAT(22X,I3)
 26   FORMAT(16X,A6)
 28   FORMAT(20X,A4)
 30   FORMAT(7X,F12.7,30X,F12.7)
 31   FORMAT(11X,A3,7X,I2,1X,E13.6)
 32   FORMAT(E14.6,40X,I3)
C******************************************************************************
C*****INITIALIZE
C
      L=0
      I1=0
      I2=0
      I3=0
      I4=0
      I6=0
      I7=0
      I8=0
      I9=0
      I10=0
      SAME=.FALSE.
      DO 100,I=1,225
        OB(I)=.FALSE.
100   CONTINUE
      M=0
      AFO=.FALSE.
      DO 110,I=1,20
        N1(I)=0.0
        N2(I)=0.0
        N3(I)=0.0
        N4(I)=0.0
        N5(I)=0.0
110   CONTINUE
      WT1=0.0
      WT2=0.0
      WT3=0.0
      WT4=0.0
      WT5=0.0
      WVW1=0.0
      WVW2=0.0
      WVW3=0.0
      WVW4=0.0
      WVW5=0.0
      DO 120,I=1,66
        RK(I)=0.0
        A(I)=0.0
        B(I)=0.0
        C(I)=0.0
        D(I)=0.0
120   CONTINUE
      DO 130,I=1,200
        PIKA(I)=1.0
130   CONTINUE
C******************************************************************************
C*****MAIN PROGRAM
C
C*****BA:`8NSBDATA
200   READ(6,2)LINE
      IF (LINE(2:6) .NE. '     ') THEN
        TITLE=LINE(2:64)
      ELSE
        GOTO 200
      END IF
205   READ(6,2,END=206)LINE
      IF (LINE(2:6) .NE. 'BASIC') GOTO 205
      READ(6,2)LINE
      IF (LINE(13:14) .EQ. 'CV') GOTO 220
      GOTO 210
C ACCOUNT FOR HITTING END-OF-FILE. BBC 9/84.
206   PRINT *,'COULD NOT FIND KEYWORD: BASIC IN COLUMN 2. PROGRAM HALT.'
      CALL EXIT
210   READ(6,2)LINE
      IF (LINE(4:5) .EQ. '  ') GOTO 210
      IF (LINE(5:5) .NE. ' ') GOTO 215
      IF (L .EQ. 0) THEN
        READ(LINE,3)RADIUS(L),THICKNS(L),GLASS(L),RN1
      ELSE
        READ(LINE,4)RADIUS(L),THICKNS(L),GLASS(L)
      END IF
      L=L+1
      GOTO 210
215   ISURF=L-1
      GOTO 240
220   READ(6,2)LINE
      IF (LINE(4:4) .EQ. ' ') GOTO 220
      IF (LINE(5:5) .NE. ' ') GOTO 230
      IF (L .EQ. 0) THEN
        READ(LINE,7)RADIUS(L),THICKNS(L),GLASS(L),RN1
      ELSE
        READ(LINE,9)RADIUS(L),THICKNS(L),GLASS(L)
      END IF
      L=L+1
      GOTO 220
230   ISURF=L-1
C******************************************************************************
C*****DETERMINE NEXT CATEGORY
240   IF (LINE(2:11) .EQ. 'REFRACTIVE') THEN
        GOTO 250
      ELSE IF (LINE(2:16) .EQ. 'CC AND ASPHERIC') THEN
        GOTO 300
      ELSE IF (LINE(2:11) .EQ. 'ASYMMETRIC') THEN
        GOTO 350
      ELSE IF (LINE(2:10) .EQ. 'SYMMETRIC') THEN
        GOTO 400
      ELSE IF (LINE(2:11) .EQ. 'ANAMORPHIC') THEN
        GOTO 425
      ELSE IF (LINE(2:7) .EQ. 'TORICS') THEN
        GOTO 450
      ELSE IF (LINE(3:9) .EQ. 'SPECIAL') THEN
        GOTO 500
      ELSE IF (LINE(2:5) .EQ. 'TILT') THEN
        GOTO 550
      ELSE IF (LINE(2:6) .EQ. 'CLEAR') THEN
        GOTO 600
      ELSE IF (LINE(2:8) .EQ. 'PICKUPS') THEN
        GOTO 650
      ELSE IF (LINE(2:7) .EQ. 'SOLVES') THEN
        GOTO 700
C READ 'REF SURF' JUST IN CASE WE GET A 'NO APERTURE STOP'. BBC 10/84.
      ELSE IF (LINE(53:60) .EQ. 'REF SURF') THEN
        READ(6,2)LINE
        READ(LINE,32) REF_OBJ_HT,REF_SURF
        READ(6,2)LINE	! SHOULD BE A BLANK LINE.
        GOTO 245
      ELSE IF (LINE(3:3) .EQ. ' ') THEN
245     READ(6,2)LINE
        GOTO 240
      ELSE
        GOTO 750
      END IF
C******************************************************************************
C*****REFRACTIVE INDICES
250   READ(6,2)LINE  ! SKIP 2 LINES
      READ(6,2)LINE
260   READ(6,2)LINE
      SAME=.FALSE.
      IF (LINE(2:5) .EQ. '    ') GOTO 245
      I7=I7+1
      READ(LINE,6)INDEX(I7),N1(I7),N2(I7),N3(I7),N4(I7),N5(I7)
      DO 265,K=1,I7-1
C IF N1 PARAM = 1.0 THEN KEEP UNIQUE GLASS NAME FOR THE PRINTOUT. BBC 10/84.
        IF ((N1(K) .EQ. N1(I7)) .AND. (N2(K) .EQ. N2(I7)) .AND.
     *      (N3(K) .EQ. N3(I7)) .AND. (N4(K) .EQ. N4(I7)) .AND.
     *      (N5(K) .EQ. N5(I7)) .AND. (N1(K) .NE. 1.0)) THEN
           GLASS(INDEX(I7))=GLASS(INDEX(K))
           SAME = .TRUE.
	ELSE IF ((GLASS(INDEX(I7))(1:5).EQ.'MODEL').AND.
     *          (GLASS(INDEX(K)) .EQ. GLASS(INDEX(I7)))) THEN
C**THE FOLLOWING STATEMENTS ARE NECESSARY BECAUSE THE 'CHAR' FUNCTION
C DOES NOT WORK. JPL CHANGE.
	   NMODELS=INDEX(I7)
	   IF(NMODELS.EQ.0) CVAR='0   '
	   IF(NMODELS.EQ.1) CVAR='1   '
	   IF(NMODELS.EQ.2) CVAR='2   '
	   IF(NMODELS.EQ.3) CVAR='3   '
	   IF(NMODELS.EQ.4) CVAR='4   '
	   IF(NMODELS.EQ.5) CVAR='5   '
	   IF(NMODELS.EQ.6) CVAR='6   '
	   IF(NMODELS.EQ.7) CVAR='7   '
	   IF(NMODELS.EQ.8) CVAR='8   '
	   IF(NMODELS.EQ.9) CVAR='9   '
	   IF(NMODELS.EQ.10) CVAR='10   '
	   IF(NMODELS.EQ.11) CVAR='11  '
	   IF(NMODELS.EQ.12) CVAR='12  '
	   IF(NMODELS.EQ.13) CVAR='13  '
	   IF(NMODELS.EQ.14) CVAR='14  '
	   IF(NMODELS.EQ.15) CVAR='15  '
	   IF(NMODELS.EQ.16) CVAR='16  '
	   IF(NMODELS.EQ.17) CVAR='17  '
	   IF(NMODELS.EQ.18) CVAR='18  '
	   IF(NMODELS.EQ.19) CVAR='19  '
	   IF(NMODELS.EQ.20) CVAR='20  '
	   IF(NMODELS.EQ.21) CVAR='21  '
	   IF(NMODELS.EQ.22) CVAR='22  '
	   IF(NMODELS.EQ.23) CVAR='23  '
	   IF(NMODELS.EQ.24) CVAR='24  '
	   IF(NMODELS.EQ.25) CVAR='25  '
	   IF(NMODELS.EQ.26) CVAR='26  '
	   IF(NMODELS.EQ.27) CVAR='27  '
	   IF(NMODELS.EQ.28) CVAR='28  '
	   IF(NMODELS.EQ.29) CVAR='29  '
	   IF(NMODELS.EQ.30) CVAR='30  '
C     CVAR='    '
C     CVAR=CHAR(INTS(INDEX(I7)))
	   GLASS(INDEX(I7))=GLASS(INDEX(I7))(1:5)//CVAR
	END IF
C***** DROP DUPLICATE PRIVATE CATALOG ENTRIES
265   CONTINUE
      IF (SAME) I7=I7-1
      IF ((GLASS(INDEX(I7))(1:6) .EQ. 'SCHOTT') .OR.
     *    (GLASS(INDEX(I7))(1:4) .EQ. 'HOYA') .OR.
     *    (GLASS(INDEX(I7))(1:5) .EQ. 'OHARA')) I7=I7-1
      GOTO 260
C******************************************************************************
C*****CC AND ASPHERIC DATA
300   READ(6,2)LINE  ! SKIP 2 LINES
      READ(6,2)LINE
320   READ(6,2)LINE
      IF (LINE(2:5) .EQ. '    ') GOTO 245
      I4=I4+1
      READ(LINE,8)NASP(I4),RK(I4),A(I4),B(I4),C(I4),D(I4)
      GOTO 320
C******************************************************************************
C*****ASYMMETRIC SPLINE
350   M=M+1
      MESSAGE(M)=1
      SUBMES1(M)=LINE(2:37)
      GOTO 245
C******************************************************************************
C*****SYMMETRIC SPLINE
400   M=M+1
      MESSAGE(M)=1
      SUBMES1(M)=LINE(2:37)
      GOTO 245
C******************************************************************************
C*****ANAMORPHIC ASPHERIC
425   M=M+1
      MESSAGE(M)=1
      SUBMES1(M)=LINE(2:37)
      READ(6,2)LINE !SKIP BLANK LINE
430   READ(6,2)LINE  ! SKIP ANAMORPHIC DATA
      IF (LINE(2:5) .NE. '    ') GOTO 430
      GOTO 245
C******************************************************************************
C*****TORICS
450   READ(6,2)LINE  ! SKIP A LINE
      READ(6,2)LINE  ! SKIP A LINE.	BBC. 9/18/84.
460   READ(6,2)LINE
      IF (LINE(2:5) .EQ. '   ') GOTO 245
      I8=I8+1
      READ(LINE,12)NTORIC(I8),TORIC(I8),TCVR(I8)
      IF (TORIC(I8) .EQ. '  X ') THEN
        M=M+1
        MESSAGE(M)=2
        I8=I8-1  ! DROP X-TORIC
      END IF
      GOTO 460
C******************************************************************************
C*****SPECIAL CONDITIONS
500   READ(6,2)LINE  ! SKIP 2 LINES
      READ(6,2)LINE
520   READ(6,2)LINE
      IF (LINE(2:5) .EQ. '    ') GOTO 245
      I3=I3+1
      READ(LINE,10)NCASE(I3),CASE(I3),DDIST(I3),NUMBER(I3)
      GOTO 520
C******************************************************************************
C*****TILT AND DECENTER DATA
550   READ(6,2)LINE  ! SKIP 2 LINES
      READ(6,2)LINE
560   READ(6,2)LINE
      IF (LINE(2:5) .EQ. '    ') GOTO 245
      I2=I2+1
      READ(LINE,16)NDECE(I2),COMMAND(I2),YDE(I2),XDE(I2),ADE(I2),
     *  BDE(I2),CDE(I2)
      GOTO 560
C******************************************************************************
C*****CLEAR APERTURES AND OBSTRUCTIONS
600   READ(6,2)LINE
      READ(6,2)LINE
620   READ(6,2)LINE
      IF (LINE(2:5) .EQ. '    ') GOTO 245
      I1=I1+1
      READ(LINE,14)NSHAPE(I1),FLAG1,SHAPE(I1),F2(I1),F1(I1)
      IF (FLAG1 .EQ. '(OB)') OB(I1)=.TRUE.
      GOTO 620
C******************************************************************************
C*****PICKUPS
650   READ(6,2)LINE  ! SKIP 2 LINES
      READ(6,2)LINE
660   READ(6,2)LINE
      IF (LINE(2:5) .EQ. '    ') GOTO 245
      I10=I10+1
      READ(LINE,18)IPIKFR(I10),PIKUP(I10),IPIKTO(I10),PIKA(I10),
     *  PIKB(I10)
      IF (PIKB(I10) .NE. 0.0) THEN
        M=M+1
        MESSAGE(M)=3
        SUBMES2(M)=LINE(3:42)
        I10=I10-1  ! DROP PICKUP WITH B VALUE
      END IF
      IF ((ABS(PIKA(I10)) .NE. 1.0) .AND. (PIKA(I10) .NE. 0.0)) THEN
        M=M+1
        MESSAGE(M)=6
        SUBMES2(M)=LINE(3:42)
        I10=I10-1
      END IF  ! DROP PICKUP WITH A VALUE NOT EQUAL TO +-1.0
      IF ((PIKUP(I10) .EQ. 'CLAP') .OR. (PIKUP(I10) .EQ. 'COBS')) THEN
        M=M+1
        MESSAGE(M)=4
        SUBMES2(M)=LINE(3:42)
        I10=I10-1  ! DROP CLAP OR COBS PICKUP
      END IF
      GOTO 660
C******************************************************************************
C*****SOLVES
700   READ(6,2)LINE ! SKIP 2 LINES
      READ(6,2)LINE
710   READ(6,2)LINE
      IF (LINE(2:5) .EQ. '    ') GOTO 245
      I6=I6+1
      READ(LINE,20)NSOLVE(I6),SOLVE(I6),SOLVF1(I6)
      IF ((SOLVE(I6) .EQ. 'CAY ') .OR. (SOLVE(I6) .EQ. 'PWRY') .OR.
     *    (SOLVE(I6) .EQ. 'PWRX') .OR. (SOLVE(I6) .EQ. 'CAX ')) THEN
        M=M+1
        MESSAGE(M)=5
        SUBMES3(M)=LINE(3:50)
        I6=I6-1  ! DROP DATA LINE
      END IF
      GOTO 710
C******************************************************************************
C*****ALTERNATE CONFIGURATIONS. BBC 9/20/84.
725   I11 = 0
      READ (6,2) LINE		! READ NEXT FOUR LINES.
      WRITE (5,244) LINE(1:78)
      READ (6,2) LINE
      WRITE (5,244) LINE(1:78)
      READ (6,2) LINE
      WRITE (5,244) LINE(1:78)
727   READ (6,2) LINE
      IF (LINE(2:14) .EQ. '             ') GOTO 750
      WRITE (5,244) LINE(1:78)
C      IF (LINE(2:4) .EQ. 'CFG') THEN
C      READ (LINE,32) CFG1
C      ELSE
C        I11 = I11 + 1
C        READ (LINE,31) CPARM(I11),NCONF(I11),ACVAL(I11)
C        CFG(I11) = CFG1
C      END IF
      GOTO 727		! DECODE THE NEXT LINE
C******************************************************************************
C*****GENERAL DATA
750   READ(6,2)LINE
      IF (LINE(2:11) .EQ. 'WAVELENGTH') THEN
        READ(LINE,22)WVW1,WVW2,WVW3,WVW4,WVW5
      ELSE IF (LINE(2:9) .EQ. 'SPECTRAL') THEN
        READ(LINE,22)WT1,WT2,WT3,WT4,WT5
      ELSE IF (LINE(2:9) .EQ. 'APERTURE') THEN
        READ(LINE,24)ISTOP
C IF 'NO APERTURE STOP' LINE, THEN 'S' DEFAULTS TO 'REF SURF'. BBC 10/84.
      ELSE IF (LINE(2:17) .EQ. 'NO APERTURE STOP') THEN
        ISTOP = REF_SURF
      ELSE IF (LINE(7:11) .EQ. 'UNITS') THEN
        READ(LINE,26)UNITS
      ELSE IF (LINE(13:16) .EQ. 'MODE') THEN
        READ(LINE,28)FLAG1
        IF (FLAG1 .EQ. 'AFOC') AFO=.TRUE.
      ELSE IF (LINE(2:8) .EQ. 'CONTROL') THEN
        READ(LINE,24)IREF
      ELSE IF (LINE(2:4) .EQ. 'PXT') THEN
        GOTO 800
C BBC 9/20/84. INSERT CODE TO HANDLE 'ALTERNATE CONFIGURATIONS'.
      ELSE IF (LINE(2:10) .EQ. 'ALTERNATE') THEN
        WRITE (5,244) LINE(1:78)
244     FORMAT (' :',A78)
        GOTO 725
      ELSE
        GOTO 750
      END IF
      GOTO 750
C******************************************************************************
C*****PXTY ALL: APERTURE AND FIELD SPECIFICATION
800   READ(6,2)LINE
      IF (LINE(4:4) .NE. '1') GOTO 800
      READ(LINE,30)PY,PUCY
C*****TEST PRINT
C      PRINT*, 'ISURF=',ISURF
C      DO 850,I=0,ISURF
C        PRINT*, I,' RADIUS',RADIUS(I),'THICKNS',THICKNS(I),
C     *     'GLASS',GLASS(I)
C850   CONTINUE
C      PRINT*, 'RN1',RN1
C      PRINT*, 'REFRACTIVE INDICES'
C      DO 860,I=1,I7
C      PRINT*, INDEX(I),N1(I),N2(I),N3(I),N4(I),N5(I)
C860   CONTINUE
C      PRINT*, 'ASPHERIC'
C      DO 870,I=1,I4
C      PRINT*, NASP(I),RK(I),A(I),B(I),C(I),D(I)
C870   CONTINUE
C      PRINT*, 'TORICS'
C      DO 880,I=1,I8
C        PRINT*, NTORIC(I),TORIC(I),TCVR(I)
C880   CONTINUE
C      PRINT*, 'SPECIAL CONDITIONS'
C      DO 890,I=1,I3
C        PRINT*, NCASE(I),CASE(I),DDIST(I),NUMBER(I)
C890   CONTINUE
C      PRINT*, 'TILT AND DECENTER'
C      DO 900,I=1,I2
C        PRINT*, NDECE(I),COMMAND(I),YDE(I),XDE(I),ADE(I),BDE(I),CDE(I)
C900   CONTINUE
C      PRINT*, 'APERTURES'
C      DO 910, I=1,I1
C        PRINT*, NAPERT(I),OB(I),SHAPE(I),F2(I),F1(I)
C910   CONTINUE
C      PRINT*, 'PICKUPS'
C      DO 920,I=1,I10
C        PRINT*, IPIKFR(I),PIKUP(I),IPIKTO(I),PIKA(I),PIKB(I)
C920   CONTINUE
C      PRINT*, 'SOLVES'
C      DO 930,I=1,I6
C        PRINT*, NSOLVE(I),SOLVE(I),SOLVF1(I)
C930   CONTINUE
C      PRINT*, 'WAVELENGTHS'
C      PRINT*, WVW1,WVW2,WVW3,WVW4,WVW5
C      PRINT*, 'WEIGHTS',WT1,WT2,WT3,WT4,WT5
C      PRINT*, 'ISTOP=',ISTOP
C      PRINT*, 'UNITS=',UNITS
C      PRINT*, 'AFOCAL=',AFO
C      PRINT*, 'IREF=',IREF
C      PRINT*, 'PY=',PY,'PUCY=',PUCY
C
      RETURN
      END

ext: for


text:
For the record I am 28 years old, college educated but lacking certification, I enjoy roller skating, bike rides, jogging, nature, science, music of various genres and walking aimlessly around town.

20080909 22:24:51 nbsp Re: Robin Williams:  You might want to tell Users/ScottMiller that.  I think he just thinks its an insult from somebody making fun of Episcopalians.  On the other hand, its probably not a big deal.  :)  Regardless, I enjoyed it. Users/JabberWokky Evan JabberWokky Edwards

20090125 00:17:34 nbsp Snobs unite. Users/CoffeeSnobDavis

20090205 20:04:13 nbsp hello CarrieB   I want to extend a meaningful  apology to you . If we have harmed you or not satisfied you , what could we do to make it up to you?  Users/jerryp

20090617 10:03:36 nbsp I keep a text file of neat quotes I find in books and historical essays.  The second paragraph from your Little Prague comment just was added to it, mostly due to the second and third sentences.  Well said! Users/JabberWokky

20091225 11:55:10 nbsp Hi Carrie Users/StevenDaubert


text:
Subroutine cren(px, py, pz, srt, i1, i2, iblock)
  Parameter (maxstr=150001, maxr=1, amn=0.939457, amp=0.93828, ap1=0.13496, ap2=0.13957, am0=1.232, pi=3.1415926, cutoff=1.8966, avmass=0.9383)
  Parameter (aka=0.498, ala=1.1157, asa=1.1974)
  Parameter (mx=4, my=4, mz=8, mpx=4, mpy=4, mpz=10, mpzp=10)
  Common /aa/r(3, maxstr)
  Common /bb/p(3, maxstr)
  Common /cc/e(maxstr)
  Common /ee/id(maxstr), lb(maxstr)
  Common /input1/masspr, massta, iseed, iavoid, dt
  Common /rndf77/nseed
  Save
  px0 = px
  py0 = py
  pz0 = pz
  ntag = 0
  iblock = 7
  ianti = 0
  If (lb(i1)<0 .Or. lb(i2)<0) Then
     ianti = 1
     iblock = -7
  End If
  kaonc = 0
  If (pnlka(srt)/(pnlka(srt)+pnska(srt))>ranart(nseed)) kaonc = 1
  If (e(i1)<=0.6) Then
     lb(i1) = 23
     e(i1) = aka
     If (kaonc==1) Then
        lb(i2) = 14
        e(i2) = ala
     Else
        lb(i2) = 15 + int(3*ranart(nseed))
        e(i2) = asa
     End If
     If (ianti==1) Then
        lb(i1) = 21
        lb(i2) = -lb(i2)
     End If
  Else
     lb(i2) = 23
     e(i2) = aka
     If (kaonc==1) Then
        lb(i1) = 14
        e(i1) = ala
     Else
        lb(i1) = 15 + int(3*ranart(nseed))
        e(i1) = asa
     End If
     If (ianti==1) Then
        lb(i2) = 21
        lb(i1) = -lb(i1)
     End If
  End If
  em1 = e(i1)
  em2 = e(i2)
  pr2 = (srt**2-em1**2-em2**2)**2 - 4.0*(em1*em2)**2
  If (pr2<=0.) pr2 = 1.E-09
  pr = sqrt(pr2)/(2.*srt)
  c1 = 1.0 - 2.0*ranart(nseed)
  t1 = 2.0*pi*ranart(nseed)
  s1 = sqrt(1.0-c1**2)
  ct1 = cos(t1)
  st1 = sin(t1)
  pz = pr*c1
  px = pr*s1*ct1
  py = pr*s1*st1
  Return
End Subroutine cren

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


! { dg-do compile }
! { dg-options "-fcoarray=lib -lcaf_single" }

! Check that pr70696 is really fixed.

  use iso_fortran_env
  type(event_type) :: x[*]

  ! exchange must not be called or the link problem before the patch
  ! does not occur.
contains
  subroutine exchange
    event post (x[1])
  end subroutine
end


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail,com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE EMG ( INT_ELEM_ID, OPT, WRITE_WARN, CALLING_SUBR, WRT_BUG_THIS_TIME )

! Main driver routine for calculation of matrices for all elements. This routine initializes appropriate arrays and calls
! other routines to calculate element matrices:

!  1) ME             = element mass matrix                             , if OPT(1) = 'Y'
!  2) PTE            = element thermal load vectors                    , if OPT(2) = 'Y'
!  3) SEi, STEi, BEi = element stress and strain data recovery matrices, if OPT(3) = 'Y'
!  4) KE             = element linea stiffness matrix                  , if OPT(4) = 'Y'
!  5) PPE            = element pressure load matrix                    , if OPT(5) = 'Y'
!  6) KED            = element differen stiff matrix calc              , if OPT(6) = 'Y'

! Also, calculate:

!  TE        = Coord transformation matrix (basic to elem)
!  ZS        = Stress recovery coeff's
!  FCONV     = Constants to convert stress to engineering force
!              NOTE: may need to calc KE to get this

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_BUG, WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR, MBUG, MEDAT0_CUSERIN, MELDOF, MEMATC, MOFFSET, NSUB, NTSUB
      USE SCONTR, ONLY                :  DEDAT_Q4_MATANG_KEY, DEDAT_Q4_POFFS_KEY, DEDAT_Q4_SHELL_KEY, DEDAT_Q4_THICK_KEY,          &
                                         DEDAT_T3_MATANG_KEY, DEDAT_T3_POFFS_KEY, DEDAT_T3_SHELL_KEY, DEDAT_T3_THICK_KEY,          &
                                         WARN_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE DEBUG_PARAMETERS, ONLY      :  DEBUG
      USE PARAMS, ONLY                :  SUPINFO, SUPWARN
      USE SUBR_BEGEND_LEVELS, ONLY    :  EMG_BEGEND
      USE CONSTANTS_1, ONLY           :  CONV_DEG_RAD, CONV_RAD_DEG, ZERO
      USE MODEL_STUF, ONLY            :  CAN_ELEM_TYPE_OFFSET, EDAT, EID, EPNT, ETYPE, ISOLID, MATANGLE, NUM_EMG_FATAL_ERRS,       &
                                         PCOMP_PROPS, PLY_NUM, TE_IDENT, THETAM, TYPE, XEL

      USE EMG_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'EMG'
      CHARACTER(LEN=*)  , INTENT(IN)  :: CALLING_SUBR       ! Name of subr that called this one
      CHARACTER( 1*BYTE), INTENT(IN)  :: OPT(6)             ! Array of EMG option indicators explained above
      CHARACTER(LEN=*)  , INTENT(IN)  :: WRITE_WARN         ! If 'Y" write warning messages, otherwise do not
      CHARACTER( 1*BYTE), INTENT(IN)  :: WRT_BUG_THIS_TIME  ! If 'Y' then write to BUG file if WRT_BUG array says to
      CHARACTER( 2*BYTE)              :: LOC                ! Location where THETAM is calculated (for DEBUG output purposes)
      CHARACTER( 1*BYTE)              :: FIX_EDAT      = 'N'! If 'Y', run code to change order of grids in EDAT for 3D elems
      CHARACTER( 1*BYTE)              :: RED_INT_SHEAR = 'N'! If 'Y', use Gaussian weighted average of B matrices for shear terms

      INTEGER(LONG), INTENT(IN)       :: INT_ELEM_ID        ! Internal element ID for which
      INTEGER(LONG)                   :: CASE_NUM    = 0    ! Can be subcase number (e.g. for UEL, PEL outout)
      INTEGER(LONG)                   :: DUM_BUG(0:MBUG-1)  ! Values from WRT_BUG sent to subr ELMOUT in a particular call
      INTEGER(LONG)                   :: EPNTK              ! Value from array EPNT at the row for this internal elem ID. It is the
!                                                             row number in array EDAT where data begins for this element.
      INTEGER(LONG)                   :: I                  ! DO loop index
      INTEGER(LONG)                   :: INT_ORDER   = 0    ! Gaussian integration order for element
      INTEGER(LONG)                   :: IORD_IJ            ! Integration order in the triangular plane for PENTA elements
      INTEGER(LONG)                   :: IORD_K             ! Integration order in Z direction for PENTA elements
      INTEGER(LONG)                   :: INT41,INT42        ! An integer used in getting MATANGLE
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = EMG_BEGEND

! **********************************************************************************************************************************
      EPNTK = EPNT(INT_ELEM_ID)
      EID   = EDAT(EPNTK)
      TYPE  = ETYPE(INT_ELEM_ID)
      CALL IS_ELEM_PCOMP_PROPS ( INT_ELEM_ID )

      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC,EID,TYPE, OPT
 9001    FORMAT(1X,A,' BEGN ',F10.3,5X,'Element No. ',I8,' ,Type ',A, ' OPT = ', 6A2)
      ENDIF
! **********************************************************************************************************************************
      NUM_EMG_FATAL_ERRS = 0

      CALL EMG_INIT

! **********************************************************************************************************************************
! Call ELMDAT1 subr to get some of the data needed for this elem.

      IF ((TYPE == 'ELAS1   ') .OR. (TYPE == 'ELAS2   ') .OR. (TYPE == 'ELAS3   ') .OR. (TYPE == 'ELAS4   ') .OR.                  &
          (TYPE == 'ROD     ') .OR. (TYPE == 'BAR     ') .OR. (TYPE == 'BEAM    ') .OR. (TYPE == 'BUSH    ') .OR.                  &
          (TYPE == 'HEXA8   ') .OR. (TYPE == 'HEXA20  ') .OR.                                                                      &
          (TYPE == 'PENTA6  ') .OR. (TYPE == 'PENTA15 ') .OR.                                                                      &
          (TYPE == 'TETRA4  ') .OR. (TYPE == 'TETRA10 ') .OR.                                                                      &
          (TYPE == 'USER1   ') .OR. (TYPE == 'USERIN  ') .OR. (TYPE == 'PLOTEL  ') .OR.                                            &
          (TYPE == 'SHEAR   ') .OR. (TYPE(1:5) == 'TRIA3') .OR. (TYPE(1:5) == 'QUAD4'   )) THEN
         CALL ELMDAT1 ( INT_ELEM_ID, WRITE_WARN )
      ELSE
         WRITE(ERR,1916) SUBR_NAME,EID,TYPE
         WRITE(F06,1916) SUBR_NAME,EID,TYPE
         FATAL_ERR = FATAL_ERR + 1
         CALL OUTA_HERE ( 'Y' )                            ! Coding error, so quit
      ENDIF

      IF (NUM_EMG_FATAL_ERRS > 0)   CALL EMG_QUIT

! **********************************************************************************************************************************
! For all elements but ELASi, USERIN, get routine to check element geometry, get element grid point coords in local elem
! coord system, and generate the coord transformation matrix to be used in transforming the stiffness matrices in
! elem system to basic sys. If elem is 'ELAS' the calcd stiff matrix IS in global coord's, so TE matrix is not needed.

!xx   IF      ((TYPE      == 'ROD     ') .OR. (TYPE == 'BAR     ') .OR. (TYPE == 'BEAM    ') .OR. (TYPE == 'BUSH    ') .OR.        &
      IF      ((TYPE      == 'ROD     ') .OR. (TYPE == 'BAR     ') .OR. (TYPE == 'BEAM    ') .OR.                                  &
               (TYPE(1:5) == 'TRIA3'   ) .OR.                                                                                &
               (TYPE      == 'PENTA6  ') .OR. (TYPE == 'PENTA15 ') .OR.                                                            &
               (TYPE      == 'TETRA4  ') .OR. (TYPE == 'TETRA10 ')) THEN
         CALL ELMGM1 ( INT_ELEM_ID, WRITE_WARN )
         FIX_EDAT = 'N'

         IF (TYPE(1:5) == 'PENTA') THEN
            IF (XEL(4,3) < ZERO) THEN
               FIX_EDAT = 'Y'
            ENDIF
         ELSE IF (TYPE(1:5) == 'TETRA') THEN
            IF (XEL(4,3) < ZERO) THEN
               FIX_EDAT = 'Y'
            ENDIF
         ENDIF

         IF (FIX_EDAT == 'Y') THEN
            CALL EDAT_FIXUP
            CALL ELMDAT1 ( INT_ELEM_ID, WRITE_WARN )
            CALL ELMGM1 ( INT_ELEM_ID, WRITE_WARN )
         ENDIF

      ELSE IF (TYPE == 'BUSH    ') THEN
         CALL ELMGM1_BUSH ( INT_ELEM_ID, WRITE_WARN )

      ELSE IF ((TYPE(1:5) == 'QUAD4') .OR. (TYPE == 'SHEAR   ')) THEN
         CALL ELMGM2 ( WRITE_WARN )

      ELSE IF ((TYPE == 'HEXA8   ') .OR. (TYPE == 'HEXA20  ')) THEN
         CALL ELMGM3 ( WRITE_WARN )
         FIX_EDAT = 'N'

         IF      (TYPE(1:4) == 'HEXA') THEN                ! See if we need to reorder grids to get local z in proper direction
            IF (XEL(5,3) < ZERO) THEN
               FIX_EDAT = 'Y'
            ENDIF
         ENDIF

         IF (FIX_EDAT == 'Y') THEN
            CALL EDAT_FIXUP
            CALL ELMDAT1 ( INT_ELEM_ID, WRITE_WARN )
            CALL ELMGM3 ( WRITE_WARN )
         ENDIF

      ELSE IF  (TYPE(1:4) == 'ELAS') THEN
         TE_IDENT = 'Y'

      ELSE IF ((TYPE == 'USER1   ') .OR. (TYPE == 'USERIN  ')) THEN
         CONTINUE

      ENDIF

      IF (NUM_EMG_FATAL_ERRS > 0)   CALL EMG_QUIT

! **********************************************************************************************************************************
! Generate element material matrices. Material matrices for shell elements with PCOMP props are generated elsewhere.
! Matrices of material props are not generated for 1-D elements
! --------

      IF ((TYPE(1:5) == 'TRIA3') .OR. (TYPE(1:5) == 'QUAD4') .OR. (TYPE == 'SHEAR   ')) THEN

         IF (PCOMP_PROPS == 'N') THEN                      ! SHEAR elem does not use PCOMP props

            THETAM = ZERO

            IF      (TYPE(1:5) == 'QUAD4') THEN
               INT41 = EDAT(EPNTK+DEDAT_Q4_MATANG_KEY)     ! Key to say whether matl angle is an actual angle or a coord ID
               INT42 = EDAT(EPNTK+DEDAT_Q4_MATANG_KEY+1)   ! Key to say (if INT41 is neg) whether coord ID is basic or otherwise
               IF      (INT41 >  0) THEN                   ! Angle is defined in array MATANGLE at row INT41
                  LOC = '#1'
                  THETAM = CONV_DEG_RAD*MATANGLE( INT41 )
               ELSE IF (INT41 <  0) THEN                   ! Angle is defined by a coord sys ID whose value is -INT41
                  LOC = '#2'
                  CALL GET_MATANGLE_FROM_CID ( -INT41 )
               ELSE IF (INT41 ==  0) THEN                  ! Angle is either specified as defined by basic coord sys or angle is 0.
                  IF      (INT42 == 1) THEN
                     LOC = '#3'
                     CALL GET_MATANGLE_FROM_CID ( 0 )
                  ELSE IF (INT42 == 0) THEN
                     LOC = '#4'
                     THETAM = ZERO
                  ENDIF
               ENDIF

            ELSE IF (TYPE(1:5) == 'TRIA3') THEN
               INT41 = EDAT(EPNTK+DEDAT_T3_MATANG_KEY)     ! Key to say whether matl angle is an actual angle or a coord ID
               INT42 = EDAT(EPNTK+DEDAT_T3_MATANG_KEY+1)   ! Key to say (if INT41 is neg) whether coord ID is basic or otherwise
               IF      (INT41 >  0) THEN                   ! Angle is defined in array MATANGLE at row INT41
                  LOC = '#1'
                  THETAM = CONV_DEG_RAD*MATANGLE( INT41 )
               ELSE IF (INT41 <  0) THEN                   ! Angle is defined by a coord sys ID whose value is -INT41
                  LOC = '#2'
                  CALL GET_MATANGLE_FROM_CID ( -INT41 )
               ELSE IF (INT41 ==  0) THEN                  ! Angle is either specified as defined by basic coord sys or angle is 0.
                  IF      (INT42 == 1) THEN
                     LOC = '#3'
                     CALL GET_MATANGLE_FROM_CID ( 0 )
                  ELSE IF (INT42 == 0) THEN
                     LOC = '#4'
                     THETAM = ZERO
                  ENDIF
               ENDIF

            ENDIF
                                                           ! Use WRITE_WARN even though the following is not a warning message
!                                                            this will allow THETAM to be printed out in only 1 call to EMG
            IF ((DEBUG(112) > 0) .AND. (WRITE_WARN == 'Y')) THEN
               IF (INT41 > 0) THEN
                  WRITE(F06,1001) TYPE, EID, CONV_RAD_DEG*THETAM, LOC, INT41, INT42
               ELSE
                  WRITE(F06,1002) TYPE, EID, CONV_RAD_DEG*THETAM, LOC, INT41, INT42
               ENDIF
            ENDIF

            CALL MATERIAL_PROPS_2D ( WRITE_WARN )
            CALL ROT_AXES_MATL_TO_LOC ( WRITE_WARN )

         ELSE

            IF (TYPE(1:5) == 'QUAD4') THEN
               INT41 = EDAT(EPNTK+DEDAT_Q4_MATANG_KEY)
               INT42 = EDAT(EPNTK+DEDAT_Q4_MATANG_KEY+1)
            ENDIF

            IF (TYPE(1:5) == 'TRIA3') THEN
               INT41 = EDAT(EPNTK+DEDAT_T3_MATANG_KEY)
               INT42 = EDAT(EPNTK+DEDAT_T3_MATANG_KEY+1)
            ENDIF

            IF (WRITE_WARN == 'Y') THEN
               IF ((INT41 /= 0) .OR. (INT42 /= 0)) THEN

                  WARN_ERR = WARN_ERR + 1
                  WRITE(ERR,1003) TYPE, EID
                  IF (SUPWARN == 'N') THEN
                     WRITE(F06,1003) TYPE, EID
                  ENDIF
               ENDIF
            ENDIF

         ENDIF

      ENDIF

      IF ((TYPE == 'HEXA8   ') .OR. (TYPE == 'HEXA20  ') .OR.                                                                      &
          (TYPE == 'PENTA6  ') .OR. (TYPE == 'PENTA15 ') .OR.                                                                      &
          (TYPE == 'TETRA4  ') .OR. (TYPE == 'TETRA10 ')) THEN
         CALL MATERIAL_PROPS_3D ( WRITE_WARN )
         CALL ROT_AXES_MATL_TO_LOC ( WRITE_WARN )
      ENDIF

! Call ELMOUT to output element data for item 0, 1

      IF (WRT_BUG(0) > 0) THEN
         CASE_NUM = 0
         DO I=0,MBUG-1
            DUM_BUG(I) = 0
         ENDDO
         DUM_BUG(0) = 1
         CALL ELMOUT ( INT_ELEM_ID, DUM_BUG, CASE_NUM, OPT )
      ENDIF

      IF ((TYPE(1:5) == 'TRIA3') .OR. (TYPE(1:5) == 'QUAD4') .OR. (TYPE(1:6) == 'SHEAR') .OR. (TYPE == 'USER1   ')) THEN
         CALL SHELL_ABD_MATRICES ( INT_ELEM_ID, WRITE_WARN )
      ENDIF

      IF (WRT_BUG(1) > 0) THEN
         CASE_NUM = 0
         DO I=0,MBUG-1
            DUM_BUG(I) = 0
         ENDDO
         DUM_BUG(1) = 1
         CALL ELMOUT ( INT_ELEM_ID, DUM_BUG, CASE_NUM, OPT )
      ENDIF

! Quick return if all OPT are 'N' (if we only need ELMDAT, ELMGMi called, and not all the other subr's which generate elem matrices

      IF ((OPT(1) == 'N') .AND. (OPT(2) == 'N') .AND. (OPT(3) == 'N') .AND. (OPT(4) == 'N') .AND. (OPT(5) == 'N') .AND.            &
          (OPT(6) == 'N')) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
         RETURN
      ENDIF

! **********************************************************************************************************************************
! For all but USERIN elem, call ELMDAT2 subr to get the rest of the data needed to calculate the matrices for this element.

      IF ((TYPE(1:4) == 'ELAS'    ) .OR. (TYPE      == 'ROD     ') .OR. (TYPE == 'BAR     ') .OR. (TYPE == 'BEAM    ') .OR.        &
          (TYPE(1:5) == 'TRIA3'   ) .OR. (TYPE(1:5) == 'QUAD4'   ) .OR. (TYPE == 'SHEAR   ') .OR. (TYPE == 'USER1   ') .OR.        &
          (TYPE      == 'HEXA8   ') .OR. (TYPE      == 'HEXA20  ') .OR.                                                            &
          (TYPE      == 'PENTA6  ') .OR. (TYPE      == 'PENTA15 ') .OR.                                                            &
          (TYPE      == 'TETRA4  ') .OR. (TYPE      == 'TETRA10 ')) THEN
         CALL ELMDAT2 ( INT_ELEM_ID, OPT, WRITE_WARN )
      ENDIF

      IF (NUM_EMG_FATAL_ERRS > 0)   CALL EMG_QUIT

! **********************************************************************************************************************************
! Now get the individual elem routines to calc the required elem matrices: ME and/or PTE and/or (SE1, SE2, STE1,STE2)
! and/or KE).

      IF (TYPE(1:4) == 'ELAS') THEN
         CALL ELAS1 ( OPT, WRITE_WARN )

      ELSE IF ((TYPE == 'ROD     ') .OR. (TYPE == 'BAR     ') .OR. (TYPE == 'BEAM    ')) THEN
         CALL BREL1 ( OPT, WRITE_WARN )
         IF (NUM_EMG_FATAL_ERRS > 0)   CALL EMG_QUIT

      ELSE IF (TYPE(1:4) == 'BUSH') THEN
         CALL BUSH ( INT_ELEM_ID, OPT, WRITE_WARN )
         IF (NUM_EMG_FATAL_ERRS > 0)   CALL EMG_QUIT

      ELSE IF (TYPE(1:5) == 'TRIA3') THEN
         CALL TREL1 ( OPT, WRITE_WARN )
         IF (NUM_EMG_FATAL_ERRS > 0)   CALL EMG_QUIT

      ELSE IF ((TYPE(1:5) == 'QUAD4') .OR. (TYPE == 'SHEAR   ')) THEN
         CALL QDEL1 ( OPT, WRITE_WARN )
         IF (NUM_EMG_FATAL_ERRS > 0)   CALL EMG_QUIT

      ELSE IF ((TYPE == 'HEXA8   ') .OR. (TYPE == 'HEXA20  ') .OR.                                                                 &
               (TYPE == 'PENTA6  ') .OR. (TYPE == 'PENTA15 ') .OR.                                                                 &
               (TYPE == 'TETRA4  ') .OR. (TYPE == 'TETRA10 ')) THEN

         IF (ISOLID(6) == 0) THEN                          ! Integration scheme
            RED_INT_SHEAR = 'Y'
         ELSE
            RED_INT_SHEAR = 'N'
         ENDIF

         IF ((TYPE == 'HEXA8   ') .OR. (TYPE == 'HEXA20  ')) THEN
            INT_ORDER = ISOLID(4)                          ! Integration order
            CALL HEXA ( OPT, INT_ELEM_ID, INT_ORDER, RED_INT_SHEAR, WRITE_WARN )
         ELSE IF ((TYPE == 'PENTA6  ') .OR. (TYPE == 'PENTA15 ')) THEN
            IF (ISOLID(4) == 2) THEN                       ! Integration order
               IORD_K  = 2
               IORD_IJ = 3
            ELSE
               IORD_K  = 3
               IORD_IJ = 7
            ENDIF
            CALL PENTA ( OPT, INT_ELEM_ID, IORD_IJ, IORD_K, RED_INT_SHEAR, WRITE_WARN )
         ELSE IF ((TYPE == 'TETRA4  ') .OR. (TYPE == 'TETRA10 ')) THEN
            IF (ISOLID(4) == 2) THEN                       ! Integration order
               INT_ORDER  = 1
            ELSE
               INT_ORDER  = 4
            ENDIF
            CALL TETRA ( OPT, INT_ELEM_ID, INT_ORDER, RED_INT_SHEAR, WRITE_WARN )
         ENDIF

      ELSE IF (TYPE == 'USER1   ') THEN
         CALL KUSER1 ( OPT, WRITE_WARN )

      ELSE IF (TYPE == 'USERIN  ') THEN
         CALL USERIN ( INT_ELEM_ID, OPT, CALLING_SUBR, WRITE_WARN )

      ENDIF

! **********************************************************************************************************************************
! For plate elements, process offsets (since they are specified in local element coordinates)

      IF ((TYPE(1:5) == 'TRIA3') .OR. (TYPE(1:5) == 'QUAD4')) THEN
         CALL ELMOFF ( OPT, WRITE_WARN )
      ENDIF

! **********************************************************************************************************************************
! Call ELMOUT to output for data items 2-5

      IF (WRT_BUG_THIS_TIME == 'Y') THEN

         IF ((WRT_BUG(2) > 0) .OR. (WRT_BUG(3) > 0) .OR. (WRT_BUG(4) > 0) .OR. (WRT_BUG(5) > 0)) THEN
            CASE_NUM = 0
            DO I=0,MBUG-1
               DUM_BUG(I) = 0
            ENDDO
            DUM_BUG(2) = WRT_BUG(2)
            DUM_BUG(3) = WRT_BUG(3)
            DUM_BUG(4) = WRT_BUG(4)
            DUM_BUG(5) = WRT_BUG(5)
            CALL ELMOUT ( INT_ELEM_ID, DUM_BUG, CASE_NUM, OPT )
         ENDIF

      ENDIF

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

      RETURN

! **********************************************************************************************************************************
 1916 FORMAT(' *ERROR  1916: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' ELEMENT',I8,' IS TYPE ',A,'. NO CODE WRITTEN FOR THIS TYPE.')

 1001 FORMAT(' THETAM material angle for ',a,i8,' = ',f8.3,' deg (at code location ',a,' in subr EMG: Row number in array MATANGLE'&
            ,' INT42 = ',I8,I3)

 1002 FORMAT(' THETAM material angle for ',a,i8,' = ',f8.3,' deg (at code location ',a,' in subr EMG: Material angle coord system' &
            ,' INT42 = ',I8,I3)

 1003 FORMAT(' WARNING     : MATERIAL ANGLE FOR ',A,I8,' IS NOT USED FOR PLATE ELEMENTS WITH PCOMP PROPERTIES')

! ##################################################################################################################################

      CONTAINS

! ##################################################################################################################################

      SUBROUTINE EMG_INIT

      USE SCONTR, ONLY                :  LSUB, MELDOF, MELGP, MDT, MPRESS, MAX_STRESS_POINTS
      USE MODEL_STUF, ONLY            :  AGRID, BE1, BE2, BE3, BGRID, DOFPIN, DT, ELAS_COMP, FCONV, KE, KED, ME,                   &
                                         OFFDIS, OFFSET, PEB, PEG, PEL, PPE, PRESS, PTE, SE1, SE2, SE3, STE1, STE2, STE3,          &
                                         UEB, UEG, UEL, UGG, XEB, XEL
      IMPLICIT NONE

      INTEGER(LONG)                   :: II,JJ,KK

! **********************************************************************************************************************************
                                                           !  1
      IF (ALLOCATED(AGRID)) THEN
         DO II=1,MELGP+1
            AGRID(II) = 0
         ENDDO
      ENDIF
                                                           !  2
      IF (ALLOCATED(BE1)) THEN
         DO II=1,3
            DO JJ=1,MELDOF
               DO KK=1,MAX_STRESS_POINTS+1
                  BE1(II,JJ,KK) = ZERO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
                                                           !  3
      IF (ALLOCATED(BE2)) THEN
         DO II=1,3
            DO JJ=1,MELDOF
               DO KK=1,MAX_STRESS_POINTS+1
                  BE2(II,JJ,KK) = ZERO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
                                                           !  4
      IF (ALLOCATED(BE3)) THEN
         DO II=1,3
            DO JJ=1,MELDOF
               DO KK=1,MAX_STRESS_POINTS+1
                  BE3(II,JJ,KK) = ZERO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
                                                           !  5
      IF (ALLOCATED(BGRID)) THEN
         DO II=1,MELGP+1
            BGRID(II) = 0
         ENDDO
      ENDIF
                                                           !  6
      IF (ALLOCATED(DOFPIN)) THEN
         DO II=1,MELDOF
            DOFPIN(II) = 0
         ENDDO
      ENDIF
                                                           !  7
      IF (ALLOCATED(DT)) THEN
         DO II=1,MDT
            DO JJ=1,LSUB
               DT(II,JJ) = 0
            ENDDO
         ENDDO
      ENDIF
                                                           !  8
      IF (ALLOCATED(KE)) THEN
         DO JJ=1,MELDOF
            DO II=1,MELDOF
               KE(II,JJ) = ZERO
            ENDDO
         ENDDO
      ENDIF
                                                           !  9
      IF (ALLOCATED(KED)) THEN
         DO JJ=1,MELDOF
            DO II=1,MELDOF
               KED(II,JJ) = ZERO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 10
      IF (ALLOCATED(ME)) THEN
         DO JJ=1,MELDOF
            DO II=1,MELDOF
               ME(II,JJ) = ZERO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 11
      IF (ALLOCATED(OFFDIS)) THEN
         DO II=1,MOFFSET
            DO JJ=1,3
               OFFDIS(II,JJ) = ZERO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 12
      IF (ALLOCATED(OFFSET)) THEN
         DO II=1,MOFFSET
            OFFSET(II) = ' '
         ENDDO
      ENDIF
                                                           ! 13
      IF (ALLOCATED(PEB)) THEN
         DO II=1,MELDOF
            PEB(II) = ZERO
         ENDDO
      ENDIF
                                                           ! 14
      IF (ALLOCATED(PEG)) THEN
         DO I=1,MELDOF
            PEG(I) = ZERO
         ENDDO
      ENDIF
                                                           ! 15
      IF (ALLOCATED(PEL)) THEN
         DO II=1,MELDOF
            PEL(II) = ZERO
         ENDDO
      ENDIF
                                                           ! 16
      IF (ALLOCATED(PPE)) THEN
         DO JJ=1,NSUB
            DO II=1,MELDOF
               PPE(II,JJ) = ZERO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 17
      IF (ALLOCATED(PRESS)) THEN
         DO II=1,MPRESS
            DO JJ=1,LSUB
               PRESS(II,JJ) = ZERO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 18
      IF (ALLOCATED(PTE)) THEN
         DO II=1,MELDOF
            DO JJ=1,LSUB
               PTE(II,JJ) = ZERO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 19
      IF (ALLOCATED(SE1)) THEN
         DO II=1,3
            DO JJ=1,MELDOF
               DO KK=1,MAX_STRESS_POINTS+1
                  SE1(II,JJ,KK) = ZERO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 20
      IF (ALLOCATED(SE2)) THEN
         DO II=1,3
            DO JJ=1,MELDOF
               DO KK=1,MAX_STRESS_POINTS+1
                  SE2(II,JJ,KK) = ZERO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 21
      IF (ALLOCATED(SE3)) THEN
         DO II=1,3
            DO JJ=1,MELDOF
               DO KK=1,MAX_STRESS_POINTS+1
                  SE3(II,JJ,KK) = ZERO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 22
      IF (ALLOCATED(STE1)) THEN
         DO II=1,3
            DO JJ=1,LSUB
               DO KK=1,MAX_STRESS_POINTS+1
                  STE1(II,JJ,KK) = ZERO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 23
      IF (ALLOCATED(STE2)) THEN
         DO II=1,3
            DO JJ=1,LSUB
               DO KK=1,MAX_STRESS_POINTS+1
                  STE2(II,JJ,KK) = ZERO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 24
      IF (ALLOCATED(STE3)) THEN
         DO II=1,3
            DO JJ=1,LSUB
               DO KK=1,MAX_STRESS_POINTS+1
                  STE3(II,JJ,KK) = ZERO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 25
      IF (ALLOCATED(UEB)) THEN
         DO II=1,MELDOF
            UEB(II) = ZERO
         ENDDO
      ENDIF
                                                           ! 26
      IF (ALLOCATED(UEG)) THEN
         DO II=1,MELDOF
            UEG(II) = ZERO
         ENDDO
      ENDIF
                                                           ! 27
      IF (ALLOCATED(UEL)) THEN
         DO II=1,MELDOF
            UEL(II) = ZERO
         ENDDO
      ENDIF
                                                           ! 28
      IF (ALLOCATED(UGG)) THEN
         DO II=1,MELDOF
            UGG(II) = ZERO
         ENDDO
      ENDIF
                                                           ! 29
      IF (ALLOCATED(XEB)) THEN
         DO II=1,MELGP+1
            DO JJ=1,3
               XEB(II,JJ) = ZERO
            ENDDO
         ENDDO
      ENDIF
                                                           ! 30
      IF (ALLOCATED(XEL)) THEN
         DO II=1,MELGP
            DO JJ=1,3
               XEL(II,JJ) = ZERO
            ENDDO
         ENDDO
      ENDIF

      DO I=1,3
         FCONV(I) = ZERO
      ENDDO

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

      END SUBROUTINE EMG_INIT

! ##################################################################################################################################

      SUBROUTINE EDAT_FIXUP

      USE MODEL_STUF, ONLY            :  ELGP

      IMPLICIT NONE

      INTEGER(LONG)                   :: DUM_AGRID(ELGP)   ! List og grids to which an elem is attached

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

      EPNTK = EPNT(INT_ELEM_ID)
      TYPE  = ETYPE(INT_ELEM_ID)
      EID   = EDAT(EPNTK)

      DO I=1,ELGP
         DUM_AGRID(I) = EDAT(EPNTK+I)
      ENDDO

      IF      (TYPE == 'HEXA8   ') THEN

         DUM_AGRID(1)  = EDAT(EPNTK+2)
         DUM_AGRID(2)  = EDAT(EPNTK+3)
         DUM_AGRID(3)  = EDAT(EPNTK+4)
         DUM_AGRID(4)  = EDAT(EPNTK+5)
         DUM_AGRID(5)  = EDAT(EPNTK+6)
         DUM_AGRID(6)  = EDAT(EPNTK+7)
         DUM_AGRID(7)  = EDAT(EPNTK+8)
         DUM_AGRID(8)  = EDAT(EPNTK+9)
         EDAT(EPNTK+2) = DUM_AGRID(5)
         EDAT(EPNTK+3) = DUM_AGRID(6)
         EDAT(EPNTK+4) = DUM_AGRID(7)
         EDAT(EPNTK+5) = DUM_AGRID(8)
         EDAT(EPNTK+6) = DUM_AGRID(1)
         EDAT(EPNTK+7) = DUM_AGRID(2)
         EDAT(EPNTK+8) = DUM_AGRID(3)
         EDAT(EPNTK+9) = DUM_AGRID(4)
         WRITE(ERR,9902) TYPE, EID, (EDAT(EPNTK+I+1),I=1,8)
         IF (SUPINFO == 'N') THEN
            WRITE(F06,9902) TYPE, EID, (EDAT(EPNTK+I+1),I=1,8)
         ENDIF

      ELSE IF (TYPE == 'PENTA6  ') THEN

         DUM_AGRID(1)  = EDAT(EPNTK+2)
         DUM_AGRID(2)  = EDAT(EPNTK+3)
         DUM_AGRID(3)  = EDAT(EPNTK+4)
         DUM_AGRID(4)  = EDAT(EPNTK+5)
         DUM_AGRID(5)  = EDAT(EPNTK+6)
         DUM_AGRID(6)  = EDAT(EPNTK+7)
         EDAT(EPNTK+2) = DUM_AGRID(4)
         EDAT(EPNTK+3) = DUM_AGRID(5)
         EDAT(EPNTK+4) = DUM_AGRID(6)
         EDAT(EPNTK+5) = DUM_AGRID(1)
         EDAT(EPNTK+6) = DUM_AGRID(2)
         EDAT(EPNTK+7) = DUM_AGRID(3)

      ELSE IF (TYPE == 'TETRA4  ') THEN

         DUM_AGRID(1)  = EDAT(EPNTK+2)
         DUM_AGRID(2)  = EDAT(EPNTK+3)
         DUM_AGRID(3)  = EDAT(EPNTK+4)
         DUM_AGRID(4)  = EDAT(EPNTK+5)
         EDAT(EPNTK+2) = DUM_AGRID(1)
         EDAT(EPNTK+3) = DUM_AGRID(3)
         EDAT(EPNTK+4) = DUM_AGRID(2)
         EDAT(EPNTK+5) = DUM_AGRID(4)

      ENDIF

! **********************************************************************************************************************************
 9902 FORMAT(' *INFORMATION: ',A,I8,' GRID ORDER HAS BEEN CHANGED TO: ',20I8)

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

      END SUBROUTINE EDAT_FIXUP

! ##################################################################################################################################

      SUBROUTINE EMG_QUIT

      IMPLICIT NONE

! **********************************************************************************************************************************
      WRITE(ERR,9999) NUM_EMG_FATAL_ERRS
      WRITE(F06,9999) NUM_EMG_FATAL_ERRS
      CALL OUTA_HERE ( 'Y' )

! **********************************************************************************************************************************
 9999 FORMAT(' Stopping processing due to ', I8, ' errors found in subr EMG')

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

      END SUBROUTINE EMG_QUIT

      END SUBROUTINE EMG


text:
      subroutine iso(iaer_prof,tamoy,trmoy,pizmoy,
     s            tamoyp,trmoyp,
     s            palt,nt,mu,rm,gb,
     s                     xf)

      include "paramdef.inc"
      double precision psl(-1:nqmax_p,-mu:mu)
c  dimension for gauss integration
      integer mu
      real rm(-mu:mu),gb(-mu:mu)
c  dimension for os computation
      real  xf(-1:1)
c array for sos computation
      double precision xpl(-mu:mu),bp(0:mu,-mu:mu),ch(0:nt)
      double precision xdel(0:nt),ydel(0:nt),h(0:nt),altc(0:nt)
      double precision i1(0:nt,-mu:mu),i2(0:nt,-mu:mu),i3(-mu:mu)
      double precision in(-mu:mu),inm1(-mu:mu),inm2(-mu:mu)
      double precision acu,acu2,ta,piz
      double precision tr,trp,tap,hr,ha,zx,yy,dd,ppp2,ppp1,ca
      double precision cr,ratio,taup,th,xt1,xt2,aaaa,ron,beta0,beta2
      double precision tavion0,tavion1,tavion2,tavion,zi1,xpk,ypk,x,y
      double precision z,xi1,xi2,bpjk,bpjmk,f,a,b,c,d,xx,a1,d1,g1
      double precision y1,xpj,xxx,ii1,ii2
      Real tamoy,trmoy,pizmoy
      Real tamoyp,trmoyp,palt
      Real delta,sigma
      integer snt,nt,iplane,ntp,j,it,itp,i,ig,k,index,iwr,m
      integer jj,l
      logical ier
      integer igmax,iaer_prof

      common/sixs_del/delta,sigma
      common/sixs_ier/iwr,ier
      common /multorder/ igmax

      snt=nt
      iplane=0
      acu=1.e-20
      acu2=1.e-3
      ta=tamoy
      piz=pizmoy
      tr=trmoy
      do 615 m=-1,1
 615  xf(m)=0.
c
c     molecular ratio within the layer
c     computations are performed assuming a scale of 8km for
c     molecules and 2km for aerosols
c
c the optical thickness above plane are recomputed to give o.t above pla
      trp=trmoy-trmoyp
      tap=tamoy-tamoyp
c     print *, 'tamoy,trmoy,pizmoy,tap,trp,palt,nt'
c     print *,tamoy,trmoy,pizmoy,tap,trp,palt,nt
      acu=1.e-20
c if plane observations recompute scale height for aerosol knowing:
c the aerosol optical depth as measure from the plane 	= tamoyp
c the rayleigh scale   height = 			= hr (8km)
c the rayleigh optical depth  at plane level 		= trmoyp
c the altitude of the plane 				= palt
c the rayleigh optical depth for total atmos		= trmoy
c the aerosol  optical depth for total atmos		= tamoy
c if not plane observations then ha is equal to 2.0km
c ntp local variable: if ntp=nt     no plane observation selected
c                        ntp=nt-1   plane observation selected

      hr=8.0
c     it's a mixing rayleigh+aerosol
      if(palt.le.900..and.palt.gt.0.0)then
      if (tap.gt.1.e-03) then
         ha=-palt/log(tap/ta)
         else
         ha=2.
         endif
      ntp=nt-1
      else
      ha=2.0
      ntp=nt
      endif
c
      ta=tamoy
      tr=trmoy
      piz=pizmoy
c
c compute mixing rayleigh, aerosol
c case 1: pure rayleigh
c case 2: pure aerosol
c case 3: mixing rayleigh-aerosol
c
      if((ta.le.acu2).and.(tr.gt.ta)) then
      do j=0,ntp
      h(j)=j*tr/ntp
      ydel(j)=1.0
      xdel(j)=0.0
      enddo
      endif
      if((tr.le.acu2).and.(ta.gt.tr)) then
      do j=0,ntp
      h(j)=j*ta/ntp
      ydel(j)=0.0
      xdel(j)=piz
      enddo
      endif
c
      if(tr.gt.acu2.and.ta.gt.acu2.and.iaer_prof.eq.0)then
      ydel(0)=1.0
      xdel(0)=0.0
      h(0)=0.
      altc(0)=300.
      zx=300.
      iplane=0
      do 14 it=0,ntp
      if (it.eq.0) then
         yy=0.
         dd=0.
         goto 111
      endif
      yy=h(it-1)
      dd=ydel(it-1)
 111  ppp2=300.0
      ppp1=0.0
      itp=it
      call discre(ta,ha,tr,hr,itp,ntp,yy,dd,ppp2,ppp1,
     s    zx)
      if(ier)return
      xxx=-zx/ha
      if (xxx.lt.-18) then
         ca=0.
         else
         ca=ta*dexp(xxx)
         endif
      xxx=-zx/hr
      cr=tr*dexp(xxx)
      h(it)=cr+ca
      altc(it)=zx
c     print *,it,cr,ca,h(it),zx
      cr=cr/hr
      ca=ca/ha
      ratio=cr/(cr+ca)
      xdel(it)=(1.e+00-ratio)*piz
      ydel(it)=ratio
  14  continue
      endif

      if(tr.gt.acu2.and.ta.gt.acu2.and.iaer_prof.eq.1)then
      call aero_prof(ta,piz,tr,hr,ntp,xmus,
     s   h,ch,ydel,xdel,altc)
      endif

c update plane layer if necessary
      if (ntp.eq.(nt-1)) then
c compute position of the plane layer
         taup=tap+trp
         iplane=-1
         do i=0,ntp
         if (taup.ge.h(i)) iplane=i
         enddo
c update the layer from the end to the position to update if necessary
         th=0.005
         xt1=abs(h(iplane)-taup)
         xt2=abs(h(iplane+1)-taup)
         if ((xt1.gt.th).and.(xt2.gt.th)) then
         do i=nt,iplane+1,-1
            xdel(i)=xdel(i-1)
            ydel(i)=ydel(i-1)
            h(i)=h(i-1)
            altc(i)=altc(i-1)
         enddo
         else
         nt=ntp
         if (xt2.lt.xt1) iplane=iplane+1
         endif
         h(iplane)=taup
         if ( tr.gt.acu2.and.ta.gt.acu2) then
         ca=ta*exp(-palt/ha)
         cr=tr*exp(-palt/hr)
         cr=cr/hr
         ca=ca/ha
         ratio=cr/(cr+ca)
         xdel(iplane)=(1.e+00-ratio)*piz
         ydel(iplane)=ratio
         altc(iplane)=palt
         endif
         if ( tr.gt.acu2.and.ta.le.acu2) then
         ydel(iplane)=1.
         xdel(iplane)=0.
         altc(iplane)=palt
         endif
         if ( tr.le.acu2.and.ta.gt.acu2) then
         ydel(iplane)=0.
         xdel(iplane)=1.*piz
         altc(iplane)=palt
         endif
      endif
c
c     print *,ha,hr,palt,tamoy,trmoy,tap,trp
c     do i=0,nt
c     print *,i,h(i),xdel(i),ydel(i),altc(i)
c     enddo
c
      aaaa=delta/(2-delta)
      ron=(1-aaaa)/(1+2*aaaa)
c
c     rayleigh phase function
c
      beta0=1.
      beta2=0.5*ron
c
c    primary scattering
c
      ig=1
      tavion0=0.
      tavion1=0.
      tavion2=0.
      tavion=0.
      do 16 j=-mu,mu
      i3(j)=0.
   16 continue
c
c     kernel computations
c
      call kernel(0,mu,rm,xpl,psl,bp)
      do 100 j=-mu,mu
      do 101 k=0,nt
      i2(k,j)=0.0000
  101 continue
  100 continue
c
c     vertical integration, primary upward radiation
c

      do 108 k=1,mu
      i1(nt,k)=1.0
      zi1=i1(nt,k)
      yy=rm(k)
      do 108 i=nt-1,0,-1
      i1(i,k)=exp(-(ta+tr-h(i))/yy)
  108 continue
c
c     vertical integration, primary downward radiation
c
      do 109 k=-mu,-1
      do 109 i=0,nt
      i1(i,k)=0.00
  109 continue
c
c     inm2 is inialized with scattering computed at n-2
c     i3 is inialized with primary scattering
c
      do 20 k=-mu,mu
      if(k) 21,20,23
   21 index=nt
      go to 25
   23 index=0
   25 continue
      inm1(k)=i1(index,k)
      inm2(k)=i1(index,k)
      i3(k)=i1(index,k)
   20 continue
      tavion=i1(iplane,mu)
      tavion2=i1(iplane,mu)
c
c     loop on successive order
c
  503 ig=ig+1
c     write(6,*) 'ig ',ig
c
c     successive orders
c
c     multiple scattering source function at every level within the laye
c
c
      do455 k=1,mu
      xpk=xpl(k)
      ypk=xpl(-k)
      do 455 i=0,nt
      ii1=0.
      ii2=0.
      x=xdel(i)
      y=ydel(i)
      do477 j=1,mu
      xpj=xpl(j)
      z=gb(j)
      xi1=i1(i,j)
      xi2=i1(i,-j)
      bpjk=bp(j,k)*x+y*(beta0+beta2*xpj*xpk)
      bpjmk=bp(j,-k)*x+y*(beta0+beta2*xpj*ypk)
      ii2=ii2+z*(xi1*bpjk+xi2*bpjmk)
      ii1=ii1+z*(xi1*bpjmk+xi2*bpjk)
 477  continue
      i2(i,k)=ii2
      i2(i,-k)=ii1
 455  continue
c
c     vertical integration, upward radiation
c
      do 48 k=1,mu
      i1(nt,k)=0.0
      zi1=i1(nt,k)
      yy=rm(k)
      do 48 i=nt-1,0,-1
      jj=i+1
      f=h(jj)-h(i)
      a=(i2(jj,k)-i2(i,k))/f
      b=i2(i,k)-a*h(i)
      c=exp(-f/yy)
      d=1.e+00-c
      xx=h(i)-h(jj)*c
      zi1=c*zi1+(d*(b+a*yy)+a*xx)*0.5e+00
      i1(i,k)=zi1
   48 continue
c
c     vertical integration, downward radiation
c
      do 50 k=-mu,-1
      i1(0,k)=0.
      zi1=i1(0,k)
      yy=rm(k)
      do 50 i=1,nt
      jj=i-1
      f=h(i)-h(jj)
      c=exp(f/yy)
      d=1.e+00-c
      a=(i2(i,k)-i2(jj,k))/f
      b=i2(i,k)-a*h(i)
      xx=h(i)-h(jj)*c
      zi1=c*zi1+(d*(b+a*yy)+a*xx)*0.5e+00
      i1(i,k)=zi1
   50 continue
c
c     in is the nieme scattering order
c
      do 30 k=-mu,mu
      if(k) 31,30,33
   31 index=nt
      go to 34
   33 index=0
   34 continue
      in(k)=i1(index,k)
   30 continue
      tavion0=i1(iplane,mu)
c
c   convergence test (geometrical serie)
c
      if(ig.gt.2) then
      z=0.
      a1=tavion2
      d1=tavion1
      g1=tavion0
      if (a1.ge.acu.and.d1.ge.acu.and.tavion.ge.acu) then
         y=((g1/d1-d1/a1)/((1.-g1/d1)**2)*(g1/tavion))
         y=abs(y)
         z=max(y,z)
      endif
      do 99 l=-mu,mu
      if (l.eq.0) goto 99
      a1=inm2(l)
      d1=inm1(l)
      g1=in(l)
      if(a1.eq.0.) go to 99
      if(d1.eq.0.) go to 99
      if(i3(l).eq.0.) go to 99
      y=((g1/d1-d1/a1)/((1-g1/d1)**2)*(g1/i3(l)))
      y=abs(y)
      z=max(y,z)
  99  continue
      if(z.lt.0.0001) then
c
c     successful test (geometrical serie)
c
      do 606 l=-mu,mu
      if (l.eq.0) goto 606
      y1=1.
      d1=inm1(l)
      g1=in(l)
      if(d1.eq.0.0) go to 606
      y1=1-g1/d1
      g1=g1/y1
      i3(l)=i3(l)+g1
  606 continue
      d1=tavion1
      g1=tavion0
      y1=1.
      if (d1.ge.acu) then
      if (abs(g1-d1).ge.acu) then
         y1=1.-g1/d1
         g1=g1/y1
      endif
      tavion=tavion+g1
      endif
      go to 505
      endif
c
c     inm2 is the (n-2)ieme scattering order
c
      do 26 k=-mu,mu
      inm2(k)=inm1(k)
   26 continue
      tavion2=tavion1
      endif
c
c     inm1 is the (n-1)ieme scattering order
c
      do 27 k=-mu,mu
      inm1(k)=in(k)
   27 continue
      tavion1=tavion0
c
c     sum of the n-1 orders
c
      do 610 l=-mu,mu
      i3(l)=i3(l)+in(l)
  610 continue
      tavion=tavion+tavion0
c
c     stop if order n is less than 1% of the sum
c
      z=0.
      do 611 l=-mu,mu
      if(i3(l).ne.0)then
      y=abs(in(l)/i3(l))
      z=max(z,y)
      endif
  611 continue
      if(z.lt.0.00001) go to 505
c
c      stop if order n is greater than 20 in any case
c
      if(ig-igmax) 503,503,505
  505 continue
c
c
      xf(1)=xf(1)+i3(mu)
      xf(-1)=tavion
      do k=1,mu
      xf(0)=xf(0)+rm(k)*gb(k)*i3(-k)
      enddo
      nt=snt
      return
      end