text:
        PROGRAM EEZMPA
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C       PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C       PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C       PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C       PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
        PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
C Define the array that holds the area map.  It is put in a labeled
C common block only because, on some machines, having a local array
C that large causes problems.
C
        COMMON /CLRCMN/ IAM(400000)
C
C Dimension the arrays needed by ARSCAM for edges.
C
        DIMENSION XCS(10000),YCS(10000)
C
C Dimension the arrays needed by ARSCAM and ARDRLN for area and group
C ids.
C
        DIMENSION IAI(10),IAG(10)
C
C Define an array for RGB triples.
C
        DIMENSION RGB(3,14)
C
C Define an array for indices determining the ordering of the colors.
C
        DIMENSION IOC(14)
C
C Define an array for aspect source flags.
C
        DIMENSION IF(13)
C
C Declare the routine which will color the areas and the one which
C will draw the lines of latitude and longitude over water.
C
        EXTERNAL COLRAM,COLRLN
C
C Define the required RGB triples and indices.
C
        DATA RGB / 0.70,0.70,0.70 , 0.75,0.50,1.00 , 0.50,0.00,1.00 ,
     +             0.00,0.00,1.00 , 0.00,0.50,1.00 , 0.00,1.00,1.00 ,
     +             0.00,1.00,0.60 , 0.00,1.00,0.00 , 0.70,1.00,0.00 ,
     +             1.00,1.00,0.00 , 1.00,0.75,0.00 , 1.00,0.38,0.38 ,
     +             1.00,0.00,0.38 , 1.00,0.00,0.00 /
C
        DATA IOC / 6,2,5,12,10,11,1,3,4,8,9,7,13,14 /
C
C Open GKS.
C
        CALL GOPKS (IERRF, ISZDM)
        CALL GOPWK (IWKID, LUNIT, IWTYPE)
        CALL GACWK (IWKID)
C
C Re-set certain aspect source flags to "individual".
C
        CALL GQASF (IE,IF)
        IF(11)=1
        IF(12)=1
        CALL GSASF (IF)
C
C Force solid fill.
C
        CALL GSFAIS (1)
C
C Define 15 different color indices.  The first 14 are spaced through
C the color spectrum and the final one is black.
C
        DO 101 J=1,14
          I=IOC(J)
          CALL GSCR(IWKID,J,RGB(1,I),RGB(2,I),RGB(3,I))
  101   CONTINUE
C
        CALL GSCR(IWKID,15,0.,0.,0.)
C
C Set up EZMAP, but don't draw anything.
C
        CALL MAPSTC ('OU','PO')
        CALL MAPROJ ('ME',0.,0.,0.)
        CALL MAPSET ('MA',0.,0.,0.,0.)
C
C Set the number of vertical strips and the group identifiers to
C be used by MAPBLA.
C
        CALL MAPSTI ('VS',150)
        CALL MAPSTI ('G1',1)
        CALL MAPSTI ('G2',2)
C
C Initialize EZMAP.
C
        CALL MAPINT
C
C Initialize the area map.
C
        CALL ARINAM (IAM,400000)
C
C Add edges to the area map.
C
        CALL MAPBLA (IAM)
C
C Pre-process the area map.
C
        CALL ARPRAM (IAM,0,0,0)
C
C Compute and print the amount of space used in the area map.
C
        ISU=400000-(IAM(6)-IAM(5)-1)
        PRINT * , 'SPACE USED IN AREA MAP IS ',ISU
C
C Set the background color.
C
        CALL GSCR (IWKID,0,1.,1.,1.)
C
C Color the map.
C
        CALL ARSCAM (IAM,XCS,YCS,10000,IAI,IAG,10,COLRAM)
C
C In black, draw a perimeter and outline all the countries.  We turn
C off the labels (since they seem to detract from the appearance of
C the plot) and we reduce the minimum vector length so as to include
C all of the points in the boundaries.
C
C Flush PLOTIT's buffers and set polyline color index to black.
C
        CALL PLOTIT(0,0,0)
        CALL GSPLCI(15)
C
        CALL MAPSTI ('LA',0)
        CALL MAPSTI ('MV',1)
        CALL MAPLBL
        CALL MAPLOT
C
C Draw lines of latitude and longitude over water.  They will be in
C black because of the GSPLCI call above.
C
        CALL MAPGRM (IAM,XCS,YCS,10000,IAI,IAG,10,COLRLN)
C
C Advance the frame.
C
        CALL FRAME
C
C Close GKS.
C
        CALL GDAWK (IWKID)
        CALL GCLWK (IWKID)
        CALL GCLKS
C
C Done.
C
        STOP
C
      END
C
      SUBROUTINE COLRAM (XCS,YCS,NCS,IAI,IAG,NAI)
      DIMENSION XCS(*),YCS(*),IAI(*),IAG(*)
C
C For each area, one gets a set of points (using normalized device
C coordinates), two group identifiers and their associated area
C identifiers.  If both of the area identifiers are zero or negative,
C the area need not be color-filled; otherwise, it is filled with
C a color obtained from MAPACI.  If the area is defined by more than
C 150 points, we'd like to know about it.  (I'm assuming that the
C device being used won't handle polygons defined by more points than
C that.)
C
        IF (IAI(1).GE.0.AND.IAI(2).GE.0) THEN
          ITM=MAX(IAI(1),IAI(2))
          IF (ITM.GT.0) THEN
            IF (NCS.GT.150) PRINT * , 'COLRAM - NCS TOO BIG - ',NCS
C
C Set area fill color index.
C
            CALL GSFACI(MAPACI(ITM))
C
            CALL GFA (NCS-1,XCS,YCS)
          END IF
        END IF
C
C Done.
C
        RETURN
C
      END
C
      SUBROUTINE COLRLN (XCS,YCS,NCS,IAI,IAG,NAI)
      DIMENSION XCS(*),YCS(*),IAI(*),IAG(*)
C
C For each line segment, one gets a set of points (using normalized
C device coordinates), two group identifiers and their associated
C area identifiers.  If both of the area identifiers are zero or
C negative, the segment is not drawn; otherwise, we use MAPACI to
C see if the segment is over water and, if so, we draw the segment.
C If the segment is defined by more than 150 points, we'd like to
C know about it.
C
        IF (IAI(1).GE.0.AND.IAI(2).GE.0) THEN
          ITM=MAX(IAI(1),IAI(2))
          IF (MAPACI(ITM).EQ.1) THEN
            IF (NCS.GT.150) PRINT * , 'COLRLN - NCS TOO BIG - ',NCS
            CALL GPL (NCS,XCS,YCS)
          END IF
        END IF
C
C Done.
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
c
c
c =========================================================
      subroutine out1(maxmx,meqn,mbc,mx,xlower,dx,q,t,iframe,aux,maux)
c =========================================================
c
c
c     # Output the results for a general system of conservation laws
c     # in 1 dimension to a hdf file as a scientific data set.
c     # See http://hdf.ncsa.uiuc.edu/ for more info on HDF.
c     # The results in this output file can be plotted in MATLAB
c     # using the "plotclaw1" script.
c
c     # Revised 2003 by Peter Blossey
c     # Adapted code written by Sorin Mitran to standard F77 CLAWPACK style
c
      implicit double precision (a-h,o-z)
c
      parameter   (nDim = 1)
      dimension    q(1-mbc:maxmx+mbc, meqn)
      dimension aux(1-mbc:maxmx+mbc, maux)
      character*14 fname
      character*13 qname2
      character*9  qname
c
c     # HDF: Declare variables that describe datasets and HDF files.
c
      integer    sd_id, sds_id, sds_rank
      integer    sds_dims, sds_start, sds_edges, sds_stride
      dimension  sds_dims(nDim), sds_start(nDim)
      dimension  sds_edges(nDim), sds_stride(nDim)
      dimension  qbuf(21), qout(mx)
c
c     # HDF: Declare external HDF functions
c
      integer  sfstart, sfcreate, sfwdata, sfscompress, sfendacc, sfend
      external sfstart, sfcreate, sfwdata, sfscompress, sfendacc, sfend
c
c     # HDF: Set up HDF constants
c
      integer 	DFACC_READ, DFACC_WRITE, DFACC_CREATE
      parameter(DFACC_READ = 1, DFACC_WRITE = 2, DFACC_CREATE = 4)

      integer   DFNT_FLOAT64, DFNT_INT32
      parameter(DFNT_FLOAT64 = 6, DFNT_INT32 = 24)

      integer   SUCCEED, FAIL
      parameter(SUCCEED = 0, FAIL = -1)
c
c     # HDF: Set up compression constants for HDF file.
c
      integer   COMP_CODE_DEFLATE, DEFLATE_LEVEL
      parameter (COMP_CODE_DEFLATE = 4, DEFLATE_LEVEL = 6)
c
c     # Write the results to the file fort.q<iframe>
c     # Use format required by matlab script  plotclaw2.m
c     # The same format is used by the amrclaw package.
c     # Here it's adapted to output just the single grid.
c
c     # first create the file name and open file
c
      fname = 'fort.q'
     &     // char(ichar('0') + mod(iframe/1000,10))
     &     // char(ichar('0') + mod(iframe/100,10))
     &     // char(ichar('0') + mod(iframe/10,10))
     &     // char(ichar('0') + mod(iframe,10))
     &     // '.hdf'
c
c     # Specify grid number and create a string which will describe this
c     # grid in HDF file.  This could be an input for simulations with
c     # multiple grids, as in AMRCLAW.
c
      ngrids_out=1
      qname = 'grid_'
     &     // char(ichar('0') + mod(ngrids_out/1000,10))
     &     // char(ichar('0') + mod(ngrids_out/100,10))
     &     // char(ichar('0') + mod(ngrids_out/10,10))
     &     // char(ichar('0') + mod(ngrids_out,10))
c
c     # HDF: create hdf file.
c
      sd_id = sfstart(fname, DFACC_CREATE)
      if (sd_id.eq.FAIL) THEN
         WRITE(*,*) 'Failed to create HDF file (call to sfstart)'
         STOP
      end if
c
c     # HDF: create a data set for parameters describing q in HDF file.
c
      sds_rank = 1
      sds_dims(1) = 21

      sds_id = sfcreate(sd_id,qname,DFNT_FLOAT64,sds_rank,sds_dims)
      if (sds_id.eq.FAIL) THEN
         WRITE(*,*) 'Failed to create scientific data set in HDF file'
         STOP
      end if
c
c     # HDF: set up parameters describing data set.
c
      sds_start(1)  = 0
      sds_edges(1)  = sds_dims(1)
      sds_stride(1) = 1
      qbuf(1) = ngrids_out
      qbuf(2) = nDim
      qbuf(3) = t
      qbuf(4) = meqn
      qbuf(5) = 1.
      qbuf(6) = mx
      qbuf(7) = 0.
      qbuf(8) = 0.
      qbuf(9) = 0.
      qbuf(10) = xlower
      qbuf(11) = 0.
      qbuf(12) = 0.
      qbuf(13) = 0.
      qbuf(14) = xlower+mx*dx
      qbuf(15) = 0.
      qbuf(16) = 0.
      qbuf(17) = 0.
      qbuf(18) = dx
      qbuf(19) = 0.
      qbuf(20) = 0.
      qbuf(21) = 0.
      istat = sfwdata(sds_id,sds_start,sds_stride,sds_edges,qbuf)
      istat = sfendacc(sds_id)
c
c     # Loop over fields in q
c
      do m = 1,meqn
c
c     # HDF: create a data set for parameters describing q in HDF file.
c
         qname2 = qname // '_'
     &        // char(ichar('0') + mod(m/100,10))
     &        // char(ichar('0') + mod(m/10,10))
     &        // char(ichar('0') + mod(m,10))
c
         sds_rank = nDim
         sds_dims(1) = mx
c
         sds_id = sfcreate(sd_id,qname,DFNT_FLOAT64,sds_rank,sds_dims)
         if (sds_id.eq.FAIL) THEN
            WRITE(*,*) 'Failed to create data set in HDF file'
            STOP
         end if
c
c     # HDF: set up parameters describing data set.
c
         sds_start(1)  = 0
         sds_edges(1)  = sds_dims(1)
         sds_stride(1) = 1

c     # Copy current field of q into qout.
c
         do i = 1,mx
            qout(i) = q(i,m)
         end do
c
c     # HDF: set compression mode and write data to hdf file.
c
         istat=sfscompress(sds_id,COMP_CODE_DEFLATE,DEFLATE_LEVEL)
         istat = sfwdata(sds_id,sds_start,sds_stride,sds_edges,qout)
         if (istat.eq.FAIL) then
            WRITE(*,*) 'Failed to write SDS (call to sfwdata)'
            STOP
         end if
c
c     # HDF: Close the data set
c
         istat = sfendacc(sds_id)
         if (istat.eq.FAIL) then
            WRITE(*,*) 'Failed to close SDS (call to sfendacc)'
            STOP
         end if
      end do
c
c     # HDF: Close HDF file.
c
      istat = sfend(sd_id)
      if (istat.eq.FAIL) then
         WRITE(*,*) 'Failed to close ', fname, ' (call to sfend)'
         STOP
      end if

      return
      end

ext: f

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
! { dg-do compile }
!
! Check whether BESXY functions take scalars and
! arrays as arguments (PR31760).
!
PROGRAM test_erf
  REAL :: r = 0.0, ra(2) = (/ 0.0, 1.0 /)

  r  = BESJ0(r)
  r  = BESJ1(r)
  r  = BESJN(0, r)

  r  = BESY0(r)
  r  = BESY1(r)
  r  = BESYN(0, r)

  ra = BESJ0(ra)
  ra = BESJ1(ra)
  ra = BESJN(0, ra)

  ra = BESY0(ra)
  ra = BESY1(ra)
  ra = BESYN(0, ra)

  r  = BESSEL_J0(r)
  r  = BESSEL_J1(r)
  r  = BESSEL_JN(0, r)

  r  = BESSEL_Y0(r)
  r  = BESSEL_Y1(r)
  r  = BESSEL_YN(0, r)

  ra = BESSEL_J0(ra)
  ra = BESSEL_J1(ra)
  ra = BESSEL_JN(0, ra)

  ra = BESSEL_Y0(ra)
  ra = BESSEL_Y1(ra)
  ra = BESSEL_YN(0, ra)

END PROGRAM

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Shell Phone 1 on the North Davis Phones node

Phone Number: (530)unknown

Location: on the side of the Davis Shell gas station, near the entrance to the car wash.

Description: It comes in a pair. It is the one on the right. With a few more stickers than the other.

Charge: 50 cents for 15 minutes of local calls. 25 cents for long distance calls in the US. One dollar for a three minute call to Mexico.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C***********************************************************************
C                                                                      *
C            NACHOS II - A Finite Element Computer Program             *
C                        for Incompressible Flow Problems              *
C                                                                      *
C     Copyright (c) 1986,2019   National Technology & Engineering      *
C                               Solutions of Sandia, LLC (NTESS)       *
C                                                                      *
C                            All rights reserved.                      *
C                                                                      *
C     This software is distributed under the BSD 3-Clause License.     *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE PTFUNC (KIND,IVAR,S,T,VAR,V)
C
C     ******************************************************************
C
C     SUBROUTINE TO EVALUATE A DEPENDENT VARIABLE AT A GIVEN POINT (S,T)
C     WITHIN AN ELEMENT
C
C     ******************************************************************
C
      DIMENSION VAR(9)
C
C     ******************************************************************
C
      GO TO (10, 10, 20, 20, 30, 30), KIND
C
C     TRIANGULAR ELEMENTS (TRI6/3,TRI6/6)
C
   10 CONTINUE
      IF (IVAR.EQ.3) THEN
      V=S*VAR(1)+T*VAR(2)+(1.-S-T)*VAR(3)
       ELSE
      V=(2.*S*S-S)*VAR(1)+(2.*T*T-T)*VAR(2)+(2.*(S+T)**2-3.*(S+T)+1.)*
     1VAR(3)+(4.*S*T)*VAR(4)+(-4.*T*T-4.*S*T+4.*T)*VAR(5)+(-4.*S*S-4.*S*
     2T+4.*S)*VAR(6)
      END IF
      RETURN
C
C     QUADRILATERAL ELEMENTS (QUAD8/4,QUAD8/8)
C
   20 CONTINUE
      IF (IVAR.EQ.3) THEN
      V=.25*(1.-S)*(1.-T)*VAR(1)+.25*(1.+S)*(1.-T)*VAR(2)+.25*(1.+S)*
     1(1.+T)*VAR(3)+.25*(1.-S)*(1.+T)*VAR(4)
       ELSE
      V=.25*(1.-S)*(1.-T)*(-S-T-1.)*VAR(1)+.25*(1.+S)*(1.-T)*(S-T-1.)*
     1VAR(2)+.25*(1.+S)*(1.+T)*(S+T-1.)*VAR(3)+.25*(1.-S)*(1.+T)*
     2(-S+T-1.)*VAR(4)+.50*(1.-S*S)*(1.-T)*VAR(5)+.50*(1.+S)*(1.-T*T)*
     3VAR(6)+.50*(1.-S*S)*(1.+T)*VAR(7)+.50*(1.-S)*(1.-T*T)*VAR(8)
      END IF
      RETURN
C
C     QUADRILATERAL ELEMENTS (QUAD9/4,QUAD9/9)
C
   30 CONTINUE
      IF (IVAR.EQ.3) THEN
      V=.25*(1.-S)*(1.-T)*VAR(1)+.25*(1.+S)*(1.-T)*VAR(2)+.25*(1.+S)*
     1(1.+T)*VAR(3)+.25*(1.-S)*(1.+T)*VAR(4)
       ELSE
      V=.25*(1.-S)*(1.-T)*(S*T)*VAR(1)+.25*(1.+S)*(1.-T)*(-S*T)*VAR(2)+
     1.25*(1.+S)*(1.+T)*(S*T)*VAR(3)+.25*(1.-S)*(1.+T)*(-S*T)*VAR(4)+
     2.50*(1.-S*S)*(1.-T)*(-T)*VAR(5)+.50*(1.+S)*(1.-T*T)*(S)*VAR(6)+
     3.50*(1.-S*S)*(1.+T)*(T)*VAR(7)+.50*(1.-S)*(1.-T*T)*(-S)*VAR(8)+
     4(1.-S*S)*(1.-T*T)*VAR(9)
      END IF
C
      RETURN
      END

max_stars_repo_name: nachos-ii/nachos-ii

max_issues_repo_name: nachos-ii/nachos-ii

max_forks_repo_name: nachos-ii/nachos-ii


text:
!! R627 alloc-opt
!    or SOURCE = source-expr
!   - F2008 addition
!   - SOURCE= in ALLOCATE statement can give array bounds as well as the value

program alloc_opt
  real :: old(2,4)
  real, allocatable :: new(:,:)

  old(1,1:4) = [11,12,13,14]
  old(2,1:4) = [21,22,23,24]

  ALLOCATE (new, SOURCE=old(1:2,1:3))

end program

max_stars_repo_name: ouankou/rose

max_issues_repo_name: WildeGeist/rose

max_forks_repo_name: WildeGeist/rose


text:
      SUBROUTINE swchem(nw,wl,nz,tlev,airden,
     $     j,sq,jlabel,tpflag)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  Load various "weighting functions" (products of cross section and        =*
*=  quantum yield at each altitude and each wavelength).  The altitude       =*
*=  dependence is necessary to ensure the consideration of pressure and      =*
*=  temperature dependence of the cross sections or quantum yields.          =*
*=  The actual reading, evaluation and interpolation is done in separate     =*
*=  subroutines for ease of management and manipulation.  Please refer to    =*
*=  the inline documentation of the specific subroutines for detail          =*
*=  information.                                                             =*
*-----------------------------------------------------------------------------*
*=  PARAMETERS:                                                              =*
*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
*=           wavelength grid                                                 =*
*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
*=           working wavelength grid                                         =*
*=  NZ     - INTEGER, number of altitude levels in working altitude grid  (I)=*
*=  TLEV   - REAL, temperature (K) at each specified altitude level       (I)=*
*=  AIRDEN - REAL, air density (molec/cc) at each altitude level          (I)=*
*=  J      - INTEGER, counter for number of weighting functions defined  (IO)=*
*=  SQ     - REAL, cross section * quantum yield (cm^2) for each          (O)=*
*=           photolysis reaction defined, at each defined wavelength and     =*
*=           at each defined altitude level                                  =*
*=  JLABEL - CHARACTER*50, string identifier for each photolysis reaction (O)=*
*=           defined                                                         =*
*-----------------------------------------------------------------------------*

      IMPLICIT NONE
      INCLUDE 'params'

* input

      INTEGER nw
      REAL wl(kw)

      INTEGER nz
      REAL tlev(kz)
      REAL airden(kz)

* weighting functions

      CHARACTER*50 jlabel(kj)
      REAL sq(kj,kz,kw)
      INTEGER tpflag(kj)

* input/output:
      INTEGER j

* local:
      REAL wc(kw)
      INTEGER iw
*_______________________________________________________________________

* complete wavelength grid

      DO 5, iw = 1, nw - 1
         wc(iw) = (wl(iw) + wl(iw+1))/2.
 5    CONTINUE

*____________________________________________________________________________

******** Ox Photochemistry

* A1.  O2 + hv -> O + O
* reserve first position.  Cross section parameterization in Schumman-Runge and
* Lyman-alpha regions are zenith-angle dependent, will be written in
* subroutine seto2.f.
* declare temperature dependence, tpflag = 1

      j = 1
      jlabel(j) = 'O2 -> O + O'
      tpflag(j) = 1

*A2.  O3 + hv ->  (both channels)
      CALL r01(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

******** HOx Photochemistry

*B1. HO2 + hv -> OH + O
      CALL r39(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*B3. H2O2 + hv -> 2 OH
      CALL r08(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

******** NOx Photochemistry

*C1.  NO2 + hv -> NO + O(3P)
      CALL r02(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*C2.  NO3 + hv ->  (both channels)
      CALL r03(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*C3.  N2O + hv -> N2 + O(1D)
      CALL r44(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*C5.  N2O5 + hv -> (both channels)
      CALL r04(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*C6.  HNO2 + hv -> OH + NO
      CALL r05(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*C7.  HNO3 + hv -> OH + NO2
      CALL r06(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*C8.  HNO4 + hv -> HO2 + NO2
      CALL r07(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

* NO3-(aq) + hv -> NO2 + O-     (for snow)
* NO3-(aq) + hv -> NO2- + O(3P) (for snow)
      CALL r118(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

******** Organic Photochemistry

*D1.  CH2O + hv -> (both channels)
c      CALL r10(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)
      CALL pxCH2O(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D2.  CH3CHO + hv -> (all three channels)
      CALL r11(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D3.  C2H5CHO + hv -> C2H5 + HCO
      CALL r12(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D9.  CH3(OOH) + hv -> CH3O + OH
      CALL r16(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D10. HOCH2OOH -> HOCH2O. + OH
      CALL r121(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D12. CH3(ONO2) + hv -> CH3O + NO2
      CALL r17(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D13. CH3(OONO2) -> CH3(OO) + NO2
      call r134(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     CH3CH2(ONO2) -> CH3CH2O + NO2
      CALL r106(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     C2H5(ONO2) -> C2H5O + NO2
      call r141(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     n-C3H7ONO2 -> n-C3H7O + NO2
      call r142(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     1-C4H9ONO2 -> 1-C4H9O + NO2
      call r143(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     2-C4H9ONO2 -> 2-C4H9O + NO2
      call r144(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     CH3CH(ONO2)CH3 -> CH3CHOCH3 + NO2
      CALL r107(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     CH2(OH)CH2(ONO2) -> CH2(OH)CH2(O.) + NO2
      CALL r108(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     CH3COCH2(ONO2) -> CH3COCH2(O.) + NO2
      CALL r109(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     C(CH3)3(ONO2) -> C(CH3)3(O.) + NO2
      CALL r110(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     C(CH3)3(ONO) -> C(CH3)3(O) + NO
      call r135(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D14. PAN + hv -> CH3CO(OO) + NO2
*     PAN + hv -> CH3CO(O) + NO3
      CALL r18(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D15. CH3CH2COO2NO2 -> CH3CH2CO(OO) + NO2
*     CH3CH2COO2NO2 -> CH2CH2CO(O) + NO3
      CALL r120(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D16. CH2=CHCHO + hv -> Products
      CALL r122(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D17. CH2=C(CH3)CHO + hv -> Products
      CALL r104(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D18. CH3COCH=CH2 + hv -> Products
      CALL r103(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D19. CH2(OH)CHO + hv -> Products
      CALL r101(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D20. CH3COCH3 + hv -> Products
      CALL r15(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     CH3COCH2CH3 -> CH3CO + CH2CH3
      CALL r119(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D21. CH2(OH)COCH3 -> CH3CO + CH2(OH)
*     CH2(OH)COCH3 -> CH2(OH)CO + CH3
      CALL r112(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D22. CHOCHO + hv -> Products
      CALL r13(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D23. CH3COCHO + hv -> Products
      CALL r14(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*     CH3COCOCH3 + hv -> Products
      CALL r102(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D25. CH3CO(OH) + hv -> Products
      CALL r138(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D26. CH3CO(OOH) + hv -> Products
      CALL r123(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*D28. CH3COCO(OH) + hv -> Products
      CALL r105(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

* (CH3)2NNO -> products
      call r124(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

******** FOx Photochemistry

*E12. CF2O + hv -> Products
      CALL r22(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

******** ClOx Photochemistry

*F1.  Cl2 + hv -> Cl + Cl

      CALL r47(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F2.  ClO -> Cl + O
      call r125(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F3.  ClOO + hv -> Products
      CALL r31(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F4.  OCLO -> Products
      call r132(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F7.  ClOOCl -> Cl + ClOO
      CALL r111(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F13. HCl -> H + Cl
      CALL r137(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F14. HOCl -> HO + Cl
      call r130(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F15. NOCl -> NO + Cl
      call r131(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F16. ClNO2 -> Cl + NO2
      call r126(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F17. ClONO -> Cl + NO2
      call r136(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F18. ClONO2 + hv -> Products
      CALL r45(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F19. CCl4 + hv -> Products
      CALL r20(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F20. CH3OCl + hv -> Cl + CH3O
      CALL r139(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F21. CHCl3 -> Products
      CALL r140(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F23. CH3Cl + hv -> Products
      CALL r30(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F24. CH3CCl3 + hv -> Products
      CALL r29(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F30. CCl2O + hv -> Products
      CALL r19(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F32. CClFO + hv -> Products
      CALL r21(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F33. CCl3F (CFC-11) + hv -> Products
      CALL r26(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F34. CCl2F2 (CFC-12) + hv -> Products
      CALL r27(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F36. CF2ClCFCl2 (CFC-113) + hv -> Products
      CALL r23(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F37. CF2ClCF2Cl (CFC-114) + hv -> Products
      CALL r24(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F38. CF3CF2Cl (CFC-115) + hv -> Products
      CALL r25(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F40. CHClF2 (HCFC-22) + hv -> Products
      CALL r38(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F42. CF3CHCl2 (HCFC-123) + hv -> Products
      CALL r32(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F43. CF3CHFCl (HCFC-124) + hv -> Products
      CALL r33(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F45. CH3CFCl2 (HCFC-141b) + hv -> Products
      CALL r34(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F46. CH3CF2Cl (HCFC-142b) + hv -> Products
      CALL r35(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F56. CF3CF2CHCl2 (HCFC-225ca) + hv -> Products
      CALL r36(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*F57. CF2ClCF2CHFCl (HCFC-225cb) + hv -> Products
      CALL r37(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

**** BrOx Photochemistry

*G1.  Br2 -> Br + Br
      CALL r115(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G3.  BrO -> Br + O
      CALL r114(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G6.  HOBr -> OH + Br
      CALL r113(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G7.  BrNO -> Br + NO
      call r127(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G8.  BrONO -> Br + NO2
*     BrONO -> BrO + NO
      call r129(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G9.  BrNO2 -> Br + NO2
      call r128(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G10. BrONO2 + hv -> Products
      CALL r46(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G11. BrCl -> Br + Cl
      call r133(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G13. CH3Br + hv -> Products
      CALL r28(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G15. CHBr3 + hv -> Products
      CALL r09(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G25. CF2Br2 (Halon-1202) + hv -> Products
      CALL r40(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G26. CF2BrCl (Halon-1211) + hv -> Products
      CALL r41(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G27. CF3Br (Halon-1301) + hv -> Products
      CALL r42(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*G35. CF2BrCF2Br (Halon-2402) + hv -> Products
      CALL r43(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

**** IOx Photochemistry

*H01. I2 -> I + I
      CALL r146(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*H02. IO -> I + O
      CALL r147(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*H05. IOH -> I + OH
      CALL r148(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

*H24. perfluoro n-iodo propane -> products
      CALL r145(nw,wl,wc,nz,tlev,airden,j,sq,jlabel,tpflag)

****************************************************************

      IF (j .GT. kj) STOP '1002'
      RETURN
      END

max_stars_repo_name: SeregaOsipov/TUV

max_issues_repo_name: SeregaOsipov/TUV

max_forks_repo_name: SeregaOsipov/TUV


text:
      SUBROUTINE MIEV0( XX, CREFIN, PERFCT, MIMCUT, ANYANG, NUMANG, XMU,
     &                  NMOM, IPOLZN, MOMDIM, PRNT, QEXT, QSCA, GQSC,
     &                  PMOM, SFORW, SBACK, S1, S2, TFORW, TBACK,
     &                  SPIKE )

c       NoPMOM version:  saves storage by deleting everything
c       having to do with PMOM calculation:

c       * deletes modules LPCoef, LPCO1T, LPCO2T

c       * NMOM,IPOLZN,PMOM,MOMDIM,CalcMo,NPQuan deleted from argument
c         lists and everywhere else they occur internally, although the
c         first four are retained in the MIEV0 argument list for
c         compatibility with complete version of MIEV0

c       * MIEV0:   shrinks size of LitA,LitB arrays to 2;  deletes call
c                  to LPCoef;  saving of An,Bn in LitA,LitB arrays
c                  deleted

c       * MiPrnt:  PMOM print deleted

c       * CkInMi:  NMOM test changed to allow only zero value (nonzero
c                  value would have no impact but indicates user
c                  conceptual error or negligence)

c       * TestMi:  saving and restoring of NMOM,IPOLZN deleted;
c                  30-loop deleted

c    Computes Mie scattering and extinction efficiencies; asymmetry
c    factor;  forward- and backscatter amplitude;  scattering
c    amplitudes vs. scattering angle for incident polarization parallel
c    and perpendicular to the plane of scattering;
c    some quantities needed in polarized radiative transfer;  and
c    information about whether or not a resonance has been hit.

c    Input and output arguments are described in file MIEV.doc;
c    internal variables are described in MIEV0.f file.

c    CALLING TREE:

c        MIEV0
c            TESTMI
c                TSTBAD
c                MIPRNT
c                ERRMSG
c            CKINMI
c                WRTBAD
c                WRTDIM
c                ERRMSG
c            SMALL1
c            SMALL2
c            ERRMSG
c            BIGA
c                CONFRA
c                    ERRMSG
c            MIPRNT

      IMPLICIT  NONE

c ----------------------------------------------------------------------
c --------  I / O SPECIFICATIONS FOR SUBROUTINE MIEV0  -----------------
c ----------------------------------------------------------------------
      LOGICAL  ANYANG, PERFCT, PRNT(*)
      INTEGER  IPOLZN, MOMDIM, NUMANG, NMOM
      REAL     GQSC, MIMCUT, PMOM( 0:MOMDIM, * ), QEXT, QSCA, SPIKE,
     $         XMU(*), XX
      COMPLEX  CREFIN, SFORW, SBACK, S1(*), S2(*), TFORW(*), TBACK(*)
c ----------------------------------------------------------------------

c                                  ** NOTE --  MaxTrm = 10100  is neces-
c                                  ** sary to do some of the test probs,
c                                  ** but 1100 is sufficient for most
c                                  ** conceivable applications
c     .. Parameters ..

      INTEGER   MAXANG, MXANG2
      PARAMETER ( MAXANG = 501, MXANG2 = MAXANG / 2 + 1 )
      INTEGER   MAXTRM
      PARAMETER ( MAXTRM = 10100 )
      REAL      ONETHR
      PARAMETER ( ONETHR = 1. / 3. )
c     ..
c     .. Local Scalars ..

      LOGICAL   NOABS, PASS1, YESANG
      INTEGER   I, J, N, NANGD2, NTRM
      REAL      CHIN, CHINM1, COEFF, DENAN, DENBN, FN, MIM, MM, MRE,
     &          NP1DN, PSIN, PSINM1, RATIO, RIORIV, RN, RTMP, TAUN,
     &          TCOEF, TWONP1, XINV
      COMPLEX   AN, ANM1, ANP, ANPM, BN, BNM1, BNP, BNPM, CDENAN,
     &          CDENBN, CIOR, CIORIV, CTMP, ZET, ZETN, ZETNM1
c     ..
c     .. Local Arrays ..

      REAL      PIN( MAXANG ), PINM1( MAXANG ), RBIGA( MAXTRM )
      COMPLEX   CBIGA( MAXTRM ), LITA( 2 ), LITB( 2 ), SM( MAXANG ),
     &          SMS( MXANG2 ), SP( MAXANG ), SPS( MXANG2 )
c     ..
c     .. External Subroutines ..

      EXTERNAL  BIGA, CKINMI, ERRMSG, MIPRNT, SMALL1, SMALL2, TESTMI
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, AIMAG, CMPLX, CONJG, COS, MAX, MIN, REAL, SIN
c     ..
      SAVE      PASS1

c     .. Statement Functions ..

      REAL      SQ
c     ..
c     .. Statement Function definitions ..

      SQ( CTMP ) = REAL( CTMP )**2 + AIMAG( CTMP )**2
c     ..
      DATA      PASS1 /.TRUE./

c                    ** Save some input variables and replace them
c                    ** with values needed to do the self-test

      IF( PASS1 ) CALL TESTMI( .FALSE., XX, CREFIN, MIMCUT, PERFCT,
     &                         ANYANG, NUMANG, XMU, QEXT, QSCA, GQSC,
     &                         SFORW, SBACK, S1, S2, TFORW, TBACK )

c                                        ** Check input and calculate
c                                        ** certain variables from input

   10 CONTINUE
      CALL CKINMI( NUMANG, MAXANG, XX, PERFCT, CREFIN, NMOM, ANYANG,
     &             XMU )

      IF( PERFCT .AND. XX.LE.0.1 ) THEN
c                                            ** Use totally-reflecting
c                                            ** small-particle limit

         CALL SMALL1( XX, NUMANG, XMU, QEXT, QSCA, GQSC, SFORW, SBACK,
     &                S1, S2, TFORW, TBACK, LITA, LITB )

         NTRM = 2
         GO TO  100

      END IF

      NOABS = .TRUE.

      IF( .NOT.PERFCT ) THEN

         CIOR = CREFIN
         IF( AIMAG( CIOR ).GT.0.0 ) CIOR   = CONJG( CIOR )

         MRE    =   REAL( CIOR )
         MIM    = - AIMAG( CIOR )
         NOABS  = MIM.LE.MIMCUT
         CIORIV = 1.0 / CIOR
         RIORIV = 1.0 / MRE

         IF( XX*MAX( 1.0, ABS(CIOR) ) .LE. 0.1 ) THEN

c                                    ** Use general-refractive-index
c                                    ** small-particle limit
c                                    ** ( Ref. 2, p. 1508 )

            CALL SMALL2( XX, CIOR, .NOT.NOABS, NUMANG, XMU, QEXT, QSCA,
     &                   GQSC, SFORW, SBACK, S1, S2, TFORW, TBACK, LITA,
     &                   LITB )

            NTRM = 2
            GO TO  100

         END IF

      END IF

      NANGD2 = ( NUMANG + 1 ) / 2
      YESANG = NUMANG.GT.0
c                              ** Estimate number of terms in Mie series
c                              ** ( Ref. 2, p. 1508 )
      IF( XX.LE.8.0 ) THEN

         NTRM   = XX + 4.*XX**ONETHR + 1.

      ELSE IF( XX.LT.4200. ) THEN

         NTRM   = XX + 4.05*XX**ONETHR + 2.

      ELSE

         NTRM   = XX + 4.*XX**ONETHR + 2.

      END IF

      IF( NTRM + 1.GT.MAXTRM )
     &    CALL ERRMSG('MIEV0--PARAMETER MaxTrm TOO SMALL',.TRUE.)

c                            ** Calculate logarithmic derivatives of
c                            ** J-Bessel-fcn., A-sub-(1 to NTrm)

      IF( .NOT.PERFCT ) CALL BIGA( CIOR, XX, NTRM, NOABS, YESANG, RBIGA,
     &                             CBIGA )

c                            ** Initialize Ricatti-Bessel functions
c                            ** (psi,chi,zeta)-sub-(0,1) for upward
c                            ** recurrence ( Ref. 1, Eq. 19 )
      XINV   = 1.0 / XX
      PSINM1 = SIN( XX )
      CHINM1 = COS( XX )
      PSIN   = PSINM1*XINV - CHINM1
      CHIN   = CHINM1*XINV + PSINM1
      ZETNM1 = CMPLX( PSINM1, CHINM1 )
      ZETN   = CMPLX( PSIN, CHIN )
c                                     ** Initialize previous coeffi-
c                                     ** cients for -GQSC- series
      ANM1 = ( 0.0, 0.0 )
      BNM1 = ( 0.0, 0.0 )
c                             ** Initialize angular function  pi
c                             ** and sums for S+, S- ( Ref. 2, p. 1507 )
      IF( ANYANG ) THEN

         DO 20 J = 1, NUMANG
            PINM1( J ) = 0.0
            PIN( J )   = 1.0
            SP( J )    = ( 0.0, 0.0 )
            SM( J )    = ( 0.0, 0.0 )
   20    CONTINUE

      ELSE

         DO 30 J = 1, NANGD2
            PINM1( J ) = 0.0
            PIN( J )   = 1.0
            SP( J )    = ( 0.0, 0.0 )
            SM( J )    = ( 0.0, 0.0 )
            SPS( J )   = ( 0.0, 0.0 )
            SMS( J )   = ( 0.0, 0.0 )
   30    CONTINUE

      END IF

c                         ** Initialize Mie sums for efficiencies, etc.
      QSCA   = 0.0
      GQSC   = 0.0
      SFORW      = ( 0., 0.)
      SBACK       = ( 0., 0.)
      TFORW( 1 ) = ( 0., 0.)
      TBACK( 1 ) = ( 0., 0.)

c ---------  LOOP TO SUM MIE SERIES  -----------------------------------

      MM     = +1.0
      SPIKE  = 1.0

      DO 60  N = 1, NTRM
c                           ** Compute various numerical coefficients
         FN     = N
         RN     = 1.0 / FN
         NP1DN  = 1.0 + RN
         TWONP1 = 2*N + 1
         COEFF  = TWONP1 / ( FN*( N + 1 ) )
         TCOEF  = TWONP1*( FN*( N + 1 ) )

c                              ** Calculate Mie series coefficients
         IF( PERFCT ) THEN
c                                   ** Totally-reflecting case

            AN     = ( ( FN*XINV )*PSIN - PSINM1 ) /
     &               ( ( FN*XINV )*ZETN - ZETNM1 )
            BN     = PSIN / ZETN

         ELSE IF( NOABS ) THEN
c                                      ** No-absorption case

            CDENAN = ( RIORIV*RBIGA(N) + (FN*XINV) ) * ZETN - ZETNM1
            AN =   ( ( RIORIV*RBIGA(N) + (FN*XINV) ) * PSIN - PSINM1 )
     $             / CDENAN
            CDENBN = (  MRE * RBIGA(N) + (FN*XINV) ) * ZETN - ZETNM1
            BN =   ( (  MRE * RBIGA(N) + (FN*XINV) ) * PSIN - PSINM1 )
     $             / CDENBN

         ELSE
c                                       ** Absorptive case

            CDENAN = ( CIORIV * CBIGA(N) + (FN*XINV) ) * ZETN - ZETNM1
            CDENBN = (   CIOR * CBIGA(N) + (FN*XINV) ) * ZETN - ZETNM1
            AN =   ( ( CIORIV * CBIGA(N) + (FN*XINV) ) * PSIN - PSINM1 )
     $             / CDENAN
            BN =   ( (   CIOR * CBIGA(N) + (FN*XINV) ) * PSIN - PSINM1 )
     $             / CDENBN

            QSCA = QSCA + TWONP1 * ( SQ( AN ) + SQ( BN ) )

         END IF

         IF( .NOT.PERFCT .AND. N.GT.XX ) THEN
c                                               ** Flag resonance spikes
            DENAN  = ABS( CDENAN )
            DENBN  = ABS( CDENBN )
            RATIO  = DENAN / DENBN

            IF( RATIO.LE.0.2 .OR. RATIO.GE.5.0 )
     &          SPIKE  = MIN( SPIKE, DENAN, DENBN )

         END IF
c                                  ** Increment Mie sums for non-angle-
c                                  ** dependent quantities

         SFORW      = SFORW      + TWONP1 * ( AN + BN )
         TFORW( 1 ) = TFORW( 1 ) + TCOEF  * ( AN - BN )
         SBACK      = SBACK      + ( MM * TWONP1 ) * ( AN - BN )
         TBACK( 1 ) = TBACK( 1 ) + ( MM * TCOEF )  * ( AN + BN )
         GQSC = GQSC + ( FN - RN ) * REAL( ANM1 * CONJG( AN )
     $                                   + BNM1 * CONJG( BN ) )
     $          + COEFF * REAL( AN * CONJG( BN ) )

         IF( YESANG ) THEN
c                                      ** Put Mie coefficients in form
c                                      ** needed for computing S+, S-
c                                      ** ( Ref. 2, p. 1507 )
            ANP = COEFF * ( AN + BN )
            BNP = COEFF * ( AN - BN )
c                                      ** Increment Mie sums for S+, S-
c                                      ** while upward recursing
c                                      ** angular functions pi and tau
            IF( ANYANG ) THEN
c                                         ** Arbitrary angles

c                                              ** vectorizable loop
               DO 40 J = 1, NUMANG
                  RTMP = ( XMU( J ) * PIN( J ) ) - PINM1( J )
                  TAUN =  FN * RTMP - PINM1( J )
                  SP( J )  = SP( J ) + ANP * ( PIN( J ) + TAUN )
                  SM( J )  = SM( J ) + BNP * ( PIN( J ) - TAUN )
                  PINM1( J ) = PIN( J )
                  PIN( J ) = ( XMU( J ) * PIN( J ) ) + NP1DN * RTMP
   40          CONTINUE

            ELSE
c                                  ** Angles symmetric about 90 degrees
               ANPM   = MM * ANP
               BNPM   = MM * BNP
c                                          ** vectorizable loop
               DO 50 J = 1, NANGD2
                  RTMP = ( XMU( J ) * PIN( J ) ) - PINM1( J )
                  TAUN =  FN * RTMP - PINM1( J )
                  SP ( J ) = SP ( J ) +  ANP * ( PIN( J ) + TAUN )
                  SMS( J ) = SMS( J ) + BNPM * ( PIN( J ) + TAUN )
                  SM ( J ) = SM ( J ) +  BNP * ( PIN( J ) - TAUN )
                  SPS( J ) = SPS( J ) + ANPM * ( PIN( J ) - TAUN )
                  PINM1( J ) = PIN( J )
                  PIN( J ) = ( XMU( J ) * PIN( J ) ) + NP1DN * RTMP
   50          CONTINUE

            END IF

         END IF
c                          ** Update relevant quantities for next
c                          ** pass through loop
         MM     = -MM
         ANM1   = AN
         BNM1   = BN
c                           ** Upward recurrence for Ricatti-Bessel
c                           ** functions ( Ref. 1, Eq. 17 )

         ZET    = ( TWONP1*XINV )*ZETN - ZETNM1
         ZETNM1 = ZETN
         ZETN   = ZET
         PSINM1 = PSIN
         PSIN   = REAL( ZETN )

   60 CONTINUE

c ---------- END LOOP TO SUM MIE SERIES --------------------------------

      QEXT   = 2./ XX**2 * REAL( SFORW )

      IF( PERFCT .OR. NOABS ) THEN

         QSCA = QEXT

      ELSE

         QSCA = 2./ XX**2 * QSCA

      END IF

      GQSC   = 4./ XX**2 * GQSC
      SFORW = 0.5 * SFORW
      SBACK = 0.5 * SBACK
      TFORW( 2 ) = 0.5 * (   SFORW + 0.25 * TFORW( 1 ) )
      TFORW( 1 ) = 0.5 * (   SFORW - 0.25 * TFORW( 1 ) )
      TBACK( 2 ) = 0.5 * (   SBACK + 0.25 * TBACK( 1 ) )
      TBACK( 1 ) = 0.5 * ( - SBACK + 0.25 * TBACK( 1 ) )

      IF( YESANG ) THEN
c                                ** Recover scattering amplitudes
c                                ** from S+, S- ( Ref. 1, Eq. 11 )
         IF( ANYANG ) THEN
c                                         ** vectorizable loop
            DO 70 J = 1, NUMANG
               S1( J ) = 0.5*( SP( J ) + SM( J ) )
               S2( J ) = 0.5*( SP( J ) - SM( J ) )
   70       CONTINUE

         ELSE
c                                         ** vectorizable loop
            DO 80 J = 1, NANGD2
               S1( J ) = 0.5*( SP( J ) + SM( J ) )
               S2( J ) = 0.5*( SP( J ) - SM( J ) )
   80       CONTINUE
c                                         ** vectorizable loop
            DO 90 J = 1, NANGD2
               S1( NUMANG + 1 - J ) = 0.5*( SPS( J ) + SMS( J ) )
               S2( NUMANG + 1 - J ) = 0.5*( SPS( J ) - SMS( J ) )
   90       CONTINUE

         END IF

      END IF

  100 CONTINUE
      IF( AIMAG( CREFIN ).GT.0.0 ) THEN
c                                         ** Take complex conjugates
c                                         ** of scattering amplitudes
         SFORW  = CONJG( SFORW )
         SBACK  = CONJG( SBACK )

         DO 110 I = 1, 2
            TFORW( I ) = CONJG( TFORW( I ) )
            TBACK( I ) = CONJG( TBACK( I ) )
  110    CONTINUE

         DO 120 J = 1, NUMANG
            S1( J ) = CONJG( S1( J ) )
            S2( J ) = CONJG( S2( J ) )
  120    CONTINUE

      END IF

      IF( PASS1 ) THEN
c                           ** Compare test case results with
c                           ** correct answers and abort if bad;
c                           ** otherwise restore user input and proceed

         CALL TESTMI( .TRUE., XX, CREFIN, MIMCUT, PERFCT, ANYANG,
     &                NUMANG, XMU, QEXT, QSCA, GQSC, SFORW, SBACK, S1,
     &                S2, TFORW, TBACK )

         PASS1  = .FALSE.
         GO TO  10

      END IF

      IF ( PRNT(1) .OR. PRNT(2) )
     $   CALL  MIPRNT( PRNT, XX, PERFCT, CREFIN, NUMANG, XMU, QEXT,
     $                 QSCA, GQSC, SFORW, SBACK, TFORW, TBACK, S1, S2 )

      RETURN

      END

      SUBROUTINE CKINMI( NUMANG, MAXANG, XX, PERFCT, CREFIN, NMOM,
     &                   ANYANG, XMU )

c        Check for bad input to MIEV0
c        (NoPMOM version)

c     Routines called :  ERRMSG, WRTBAD, WRTDIM

      IMPLICIT  NONE

c     .. Scalar Arguments ..

      LOGICAL   ANYANG, PERFCT
      INTEGER   MAXANG, NMOM, NUMANG
      REAL      XX
      COMPLEX   CREFIN
c     ..
c     .. Array Arguments ..

      REAL      XMU( * )
c     ..
c     .. Local Scalars ..

      LOGICAL   INPERR
      INTEGER   I
c     ..
c     .. External Functions ..

      LOGICAL   WRTBAD, WRTDIM
      EXTERNAL  WRTBAD, WRTDIM
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, AIMAG, REAL
c     ..

      INPERR = .FALSE.

      IF( NUMANG.GT.MAXANG ) INPERR = WRTDIM( 'MaxAng', NUMANG )
      IF( NUMANG.LT.0 ) INPERR = WRTBAD( 'NUMANG' )

      IF( XX.LT.0.) INPERR = WRTBAD( 'XX' )

      IF( .NOT.PERFCT .AND. REAL( CREFIN ).LE.0. )
     &    INPERR = WRTBAD( 'CREFIN' )

      IF( NMOM.NE.0 ) THEN

         INPERR = WRTBAD( 'NMOM' )
         WRITE( *, '(A)' )
     &      ' For nonzero NMOM, use complete version of MIEV0'
      END IF

      IF( ANYANG ) THEN
c                                ** Allow for slight imperfections in
c                                ** computation of cosine
         DO 10 I = 1, NUMANG

             IF ( XMU(I).LT.-1.00001 .OR. XMU(I).GT.1.00001 )
     $            INPERR = WRTBAD( 'XMU' )

   10    CONTINUE

      ELSE

         DO 20 I = 1, ( NUMANG + 1 ) / 2

             IF ( XMU(I).LT.-0.00001 .OR. XMU(I).GT.1.00001 )
     $            INPERR = WRTBAD( 'XMU' )

   20    CONTINUE

      END IF

      IF( INPERR ) CALL ERRMSG( 'MIEV0--INPUT ERROR(S).  Aborting...',
     &                          .TRUE.)

      IF( XX.GT.20000.0 .OR. REAL( CREFIN ).GT.10.0 .OR.
     &    ABS( AIMAG( CREFIN ) ).GT.10.0 )
     &    CALL ERRMSG( 'MIEV0--XX or CREFIN outside tested range',
     &    .FALSE. )

      RETURN
      END

      SUBROUTINE BIGA( CIOR, XX, NTRM, NOABS, YESANG, RBIGA, CBIGA )

c        Calculate logarithmic derivatives of J-Bessel-function

c     Input :  CIOR, XX, NTRM, NOABS, YESANG  (defined in MIEV0)

c    Output :  RBIGA or CBIGA  (defined in MIEV0)

c    Routines called :  CONFRA

c    INTERNAL VARIABLES :

c       CONFRA     Value of Lentz continued fraction for cBigA(NTrm),
c                     used to initialize downward recurrence.
c       DOWN       = True, use down-recurrence.  False, do not.
c       F1,F2,F3   Arithmetic statement functions used in determining
c                     whether to use up-  or down-recurrence
c                     ( Ref. 2, Eqs. 6-8 )
c       MRE        Real refractive index
c       MIM        Imaginary refractive index
c       REZINV     1 / ( MRE * XX ); temporary variable for recurrence
c       ZINV       1 / ( CIOR * XX ); temporary variable for recurrence

      IMPLICIT  NONE

c     .. Scalar Arguments ..

      LOGICAL   NOABS, YESANG
      INTEGER   NTRM
      REAL      XX
      COMPLEX   CIOR
c     ..
c     .. Array Arguments ..

      REAL      RBIGA( * )
      COMPLEX   CBIGA( * )
c     ..
c     .. Local Scalars ..

      LOGICAL   DOWN
      INTEGER   N
      REAL      MIM, MRE, REZINV, RTMP
      COMPLEX   CTMP, ZINV
c     ..
c     .. External Functions ..

      COMPLEX   CONFRA
      EXTERNAL  CONFRA
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, AIMAG, COS, EXP, REAL, SIN
c     ..
c     .. Statement Functions ..

      REAL      F1, F2, F3
c     ..
c     .. Statement Function definitions ..

      F1( MRE ) = -8.0 + MRE**2*( 26.22 +
     &            MRE*( -0.4474 + MRE**3*( 0.00204 - 0.000175*MRE ) ) )

      F2( MRE ) = 3.9 + MRE*( -10.8 + 13.78*MRE )

      F3( MRE ) = -15.04 + MRE*( 8.42 + 16.35*MRE )
c     ..

c                                  ** Decide whether BigA can be
c                                  ** calculated by up-recurrence
      MRE = REAL( CIOR )
      MIM = ABS( AIMAG( CIOR ) )

      IF( MRE.LT.1.0 .OR. MRE.GT.10.0 .OR. MIM.GT.10.0 ) THEN

         DOWN = .TRUE.

      ELSE IF( YESANG ) THEN

         DOWN = .TRUE.
         IF( MIM*XX .LT. F2( MRE ) ) DOWN   = .FALSE.

      ELSE

         DOWN = .TRUE.
         IF( MIM*XX .LT. F1( MRE ) ) DOWN   = .FALSE.

      END IF

      ZINV   = 1.0 / ( CIOR*XX )
      REZINV = 1.0 / ( MRE*XX )

      IF( DOWN ) THEN
c                          ** Compute initial high-order BigA using
c                          ** Lentz method ( Ref. 1, pp. 17-20 )

         CTMP = CONFRA( NTRM, ZINV )

c                                   *** Downward recurrence for BigA
c                                   *** ( Ref. 1, Eq. 22 )
         IF( NOABS ) THEN
c                                            ** No-absorption case
            RBIGA( NTRM ) = REAL( CTMP )

            DO 10 N = NTRM, 2, -1
               RBIGA( N - 1 ) = ( N*REZINV ) -
     &                          1.0 / ( ( N*REZINV ) + RBIGA( N ) )
   10       CONTINUE

         ELSE
c                                            ** Absorptive case
            CBIGA( NTRM ) = CTMP

            DO 20 N = NTRM, 2, -1
               CBIGA( N-1 ) = (N*ZINV) - 1.0 / ( (N*ZINV) + CBIGA( N ) )
   20       CONTINUE

         END IF

      ELSE
c                              *** Upward recurrence for BigA
c                              *** ( Ref. 1, Eqs. 20-21 )
         IF( NOABS ) THEN
c                                            ** No-absorption case
            RTMP   = SIN( MRE*XX )
            RBIGA( 1 ) = -REZINV + RTMP /
     &                   ( RTMP*REZINV - COS( MRE*XX ) )

            DO 30 N = 2, NTRM
               RBIGA( N ) = -( N*REZINV ) +
     &                      1.0 / ( ( N*REZINV ) - RBIGA( N - 1 ) )
   30       CONTINUE

         ELSE
c                                                ** Absorptive case
            CTMP = EXP( - (0.,2.) * CIOR * XX )
            CBIGA( 1 ) = - ZINV + (1.-CTMP) /
     $                          ( ZINV * (1.-CTMP) - (0.,1.)*(1.+CTMP) )

            DO 40 N = 2, NTRM
               CBIGA( N ) = - (N*ZINV) + 1.0 / ((N*ZINV) - CBIGA( N-1 ))
   40       CONTINUE

         END IF

      END IF

      RETURN
      END

      COMPLEX FUNCTION CONFRA( N, ZINV )

c         Compute Bessel function ratio A-sub-N from its
c         continued fraction using Lentz method ( Ref. 1, pp. 17-20 )

c         ZINV = Reciprocal of argument of A

c    I N T E R N A L    V A R I A B L E S
c    ------------------------------------

c    CAK      Term in continued fraction expansion of A (Ref. 1, Eq. 25)
c    CAPT     Factor used in Lentz iteration for A (Ref. 1, Eq. 27)
c    CDENOM   Denominator in capT  ( Ref. 1, Eq. 28B )
c    CNUMER   Numerator   in capT  ( Ref. 1, Eq. 28A )
c    CDTD     Product of two successive denominators of capT factors
c                 ( Ref. 1, Eq. 34C )
c    CNTN     Product of two successive numerators of capT factors
c                 ( Ref. 1, Eq. 34B )
c    EPS1     Ill-conditioning criterion
c    EPS2     Convergence criterion
c    KK       Subscript k of cAk  ( Ref. 1, Eq. 25B )
c    KOUNT    Iteration counter ( used only to prevent runaway )
c    MAXIT    Max. allowed no. of iterations
c    MM       + 1  and - 1, alternately

      IMPLICIT  NONE

c     .. Scalar Arguments ..

      INTEGER   N
      COMPLEX   ZINV
c     ..
c     .. Local Scalars ..

      INTEGER   KK, KOUNT, MAXIT, MM
      REAL      EPS1, EPS2
      COMPLEX   CAK, CAPT, CDENOM, CDTD, CNTN, CNUMER
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, AIMAG, REAL
c     ..
      DATA      EPS1 / 1.E-2 / , EPS2 / 1.E-8 /
      DATA      MAXIT / 10000 /

c                                      *** Ref. 1, Eqs. 25a, 27
      CONFRA = ( N + 1 ) * ZINV
      MM     = - 1
      KK     = 2*N + 3
      CAK    = ( MM*KK ) * ZINV
      CDENOM = CAK
      CNUMER = CDENOM + 1.0 / CONFRA
      KOUNT  = 1

   10 CONTINUE
      KOUNT = KOUNT + 1

      IF( KOUNT.GT.MAXIT )
     &    CALL ERRMSG('ConFra--Iteration failed to converge',.TRUE.)

c                                         *** Ref. 2, Eq. 25b
      MM  = - MM
      KK  = KK + 2
      CAK = ( MM*KK ) * ZINV
c                                         *** Ref. 2, Eq. 32
      IF( ABS( CNUMER / CAK ).LE.EPS1 .OR.
     &    ABS( CDENOM / CAK ).LE.EPS1 ) THEN

c                                  ** Ill-conditioned case -- stride
c                                  ** two terms instead of one

c                                         *** Ref. 2, Eqs. 34
         CNTN   = CAK * CNUMER + 1.0
         CDTD   = CAK * CDENOM + 1.0
         CONFRA = ( CNTN / CDTD ) * CONFRA
c                                             *** Ref. 2, Eq. 25b
         MM  = - MM
         KK  = KK + 2
         CAK = ( MM*KK ) * ZINV
c                                         *** Ref. 2, Eqs. 35
         CNUMER = CAK + CNUMER / CNTN
         CDENOM = CAK + CDENOM / CDTD
         KOUNT  = KOUNT + 1
         GO TO  10

      ELSE
c                                ** Well-conditioned case

c                                        *** Ref. 2, Eqs. 26, 27
         CAPT   = CNUMER / CDENOM
         CONFRA = CAPT * CONFRA
c                                    ** Check for convergence
c                                    ** ( Ref. 2, Eq. 31 )

         IF (      ABS( REAL (CAPT) - 1.0 ).GE.EPS2
     $        .OR. ABS( AIMAG(CAPT) )      .GE.EPS2 )  THEN

c                                        *** Ref. 2, Eqs. 30A-B
            CNUMER = CAK + 1.0 / CNUMER
            CDENOM = CAK + 1.0 / CDENOM
            GO TO  10

         END IF

      END IF

      RETURN

      END

      SUBROUTINE MIPRNT( PRNT, XX, PERFCT, CREFIN, NUMANG, XMU, QEXT,
     &                   QSCA, GQSC, SFORW, SBACK, TFORW, TBACK, S1,
     &                   S2 )

c         Print scattering quantities of a single particle
c         (NoPMOM version)

      IMPLICIT  NONE

c     .. Scalar Arguments ..

      LOGICAL   PERFCT
      INTEGER   NUMANG
      REAL      GQSC, QEXT, QSCA, XX
      COMPLEX   CREFIN, SBACK, SFORW
c     ..
c     .. Array Arguments ..

      LOGICAL   PRNT( * )
      REAL      XMU( * )
      COMPLEX   S1( * ), S2( * ), TBACK( * ), TFORW( * )
c     ..
c     .. Local Scalars ..

      INTEGER   I
      REAL      I1, I2
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC AIMAG, CONJG, REAL
c     ..

      IF( PERFCT ) WRITE( *, '(''1'',10X,A,1P,E11.4)' )
     &    'Perfectly Conducting Case, size parameter =', XX

      IF( .NOT.PERFCT ) WRITE( *, '(''1'',10X,3(A,1P,E11.4))' )
     &    'Refractive Index:  Real ', REAL( CREFIN ), '  Imag ',
     &    AIMAG( CREFIN ), ',   Size Parameter =', XX

      IF( PRNT( 1 ) .AND. NUMANG.GT.0 ) THEN

         WRITE( *, '(/,A)' )
     &      '    cos(angle)  ------- S1 ---------  ------- S2 ---------'
     &      // '  --- S1*conjg(S2) ---   i1=S1**2   i2=S2**2  (i1+i2)/2'
     &      // '  DEG POLZN'

         DO 10 I = 1, NUMANG
            I1     = REAL( S1( I ) )**2 + AIMAG( S1( I ) )**2
            I2     = REAL( S2( I ) )**2 + AIMAG( S2( I ) )**2
            WRITE( *, '( I4, F10.6, 1P,10E11.3 )'   )
     $              I, XMU(I), S1(I), S2(I), S1(I)*CONJG(S2(I)),
     $              I1, I2, 0.5*(I1+I2), (I2-I1)/(I2+I1)
   10    CONTINUE

      END IF

      IF( PRNT( 2 ) ) THEN

         WRITE ( *, '(/,A,9X,A,17X,A,17X,A,/,(0P,F7.2, 1P,6E12.3) )' )
     $           '  Angle', 'S-sub-1', 'T-sub-1', 'T-sub-2',
     $               0.0,     SFORW,    TFORW(1),  TFORW(2),
     $              180.,     SBACK,    TBACK(1),  TBACK(2)
         WRITE ( *, '(/,4(A,1P,E11.4))' )
     $           ' Efficiency Factors,  extinction:', QEXT,
     $                              '   scattering:', QSCA,
     $                              '   absorption:', QEXT-QSCA,
     $                           '   rad. pressure:', QEXT-GQSC

      END IF

      RETURN

      END

      SUBROUTINE SMALL1( XX, NUMANG, XMU, QEXT, QSCA, GQSC, SFORW,
     &                   SBACK, S1, S2, TFORW, TBACK, A, B )

c       Small-particle limit of Mie quantities in totally reflecting
c       limit ( Mie series truncated after 2 terms )

c        A,B       First two Mie coefficients, with numerator and
c                  denominator expanded in powers of XX ( a factor
c                  of XX**3 is missing but is restored before return
c                  to calling program )  ( Ref. 2, p. 1508 )

      IMPLICIT  NONE

c     .. Parameters ..

      REAL      TWOTHR, FIVTHR, FIVNIN
      PARAMETER ( TWOTHR = 2. / 3., FIVTHR = 5. / 3., FIVNIN = 5. / 9. )
c     ..
c     .. Scalar Arguments ..

      INTEGER   NUMANG
      REAL      GQSC, QEXT, QSCA, XX
      COMPLEX   SBACK, SFORW
c     ..
c     .. Array Arguments ..

      REAL      XMU( * )
      COMPLEX   A( * ), B( * ), S1( * ), S2( * ), TBACK( * ), TFORW( * )
c     ..
c     .. Local Scalars ..

      INTEGER   J
      REAL      RTMP
      COMPLEX   CTMP
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC AIMAG, CMPLX, CONJG, REAL
c     ..
c     .. Statement Functions ..

      REAL      SQ
c     ..
c     .. Statement Function definitions ..

      SQ( CTMP ) = REAL( CTMP )**2 + AIMAG( CTMP )**2
c     ..

      A( 1 ) = CMPLX( 0., TWOTHR*( 1. - 0.2*XX**2 ) ) /
     &         CMPLX( 1. - 0.5*XX**2, TWOTHR*XX**3 )

      B( 1 ) = CMPLX( 0., -( 1. - 0.1*XX**2 ) / 3. ) /
     &         CMPLX( 1. + 0.5*XX**2, -XX**3 / 3. )

      A( 2 ) = CMPLX( 0.,   XX**2 / 30. )
      B( 2 ) = CMPLX( 0., - XX**2 / 45. )

      QSCA   = 6.* XX**4 *( SQ( A(1) ) + SQ( B(1) ) +
     &                      FIVTHR*( SQ( A(2) ) + SQ( B(2) ) ) )
      QEXT   = QSCA
      GQSC   = 6.* XX**4 * REAL( A(1)*CONJG( A(2) + B(1) ) +
     &                          ( B(1) + FIVNIN*A(2) )*CONJG( B(2) ) )

      RTMP   = 1.5 * XX**3
      SFORW  = RTMP * ( A(1) + B(1) + FIVTHR * ( A(2) + B(2) ) )
      SBACK  = RTMP * ( A(1) - B(1) - FIVTHR * ( A(2) - B(2) ) )
      TFORW( 1 ) = RTMP*( B(1) + FIVTHR * ( 2.*B(2) - A(2) ) )
      TFORW( 2 ) = RTMP*( A(1) + FIVTHR * ( 2.*A(2) - B(2) ) )
      TBACK( 1 ) = RTMP*( B(1) - FIVTHR * ( 2.*B(2) + A(2) ) )
      TBACK( 2 ) = RTMP*( A(1) - FIVTHR * ( 2.*A(2) + B(2) ) )

      DO 10 J = 1, NUMANG
         S1( J ) = RTMP*( A(1) + B(1)*XMU( J ) +
     &                    FIVTHR*( A(2)*XMU( J ) +
     &                             B(2)*( 2.*XMU( J )**2 - 1.) ) )
         S2( J ) = RTMP*( B(1) + A(1)*XMU( J ) +
     &                    FIVTHR*( B(2)*XMU( J ) +
     &                             A(2)*( 2.*XMU( J )**2 - 1.) ) )
   10 CONTINUE

c                                     ** Recover actual Mie coefficients
      A( 1 ) = XX**3 * A(1)
      A( 2 ) = XX**3 * A(2)
      B( 1 ) = XX**3 * B(1)
      B( 2 ) = XX**3 * B(2)

      RETURN
      END

      SUBROUTINE SMALL2( XX, CIOR, CALCQE, NUMANG, XMU, QEXT, QSCA,
     &                   GQSC, SFORW, SBACK, S1, S2, TFORW, TBACK, A,
     &                   B )

c       Small-particle limit of Mie quantities for general refractive
c       index ( Mie series truncated after 2 terms )

c        A,B       First two Mie coefficients, with numerator and
c                  denominator expanded in powers of XX ( a factor
c                  of XX**3 is missing but is restored before return
c                  to calling program )  ( Ref. 2, p. 1508 )

c        CIORSQ    Square of refractive index

      IMPLICIT  NONE

c     .. Parameters ..

      REAL      TWOTHR, FIVTHR
      PARAMETER  ( TWOTHR = 2./3., FIVTHR = 5./3. )
c     ..
c     .. Scalar Arguments ..

      LOGICAL   CALCQE
      INTEGER   NUMANG
      REAL      GQSC, QEXT, QSCA, XX
      COMPLEX   CIOR, SBACK, SFORW
c     ..
c     .. Array Arguments ..

      REAL      XMU( * )
      COMPLEX   A( * ), B( * ), S1( * ), S2( * ), TBACK( * ), TFORW( * )
c     ..
c     .. Local Scalars ..

      INTEGER   J
      REAL      RTMP
      COMPLEX   CIORSQ, CTMP
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC AIMAG, CMPLX, CONJG, REAL
c     ..
c     .. Statement Functions ..

      REAL      SQ
c     ..
c     .. Statement Function definitions ..

      SQ( CTMP ) = REAL( CTMP )**2 + AIMAG( CTMP )**2
c     ..

      CIORSQ = CIOR**2
      CTMP   = CMPLX( 0., TWOTHR )*( CIORSQ - 1.0 )
      A( 1 ) = CTMP*( 1. - 0.1*XX**2 +
     &         ( CIORSQ / 350.+ 1./ 280.)*XX**4 ) /
     &         ( CIORSQ + 2.+ ( 1.- 0.7*CIORSQ )*XX**2 -
     &         ( CIORSQ**2 / 175.- 0.275*CIORSQ + 0.25 )*XX**4 +
     &         XX**3 * CTMP*( 1.- 0.1*XX**2 ) )

      B( 1 ) = ( XX**2 / 30.) * CTMP * ( 1.+
     &         ( CIORSQ / 35.- 1./ 14.)*XX**2 ) /
     &         ( 1.- ( CIORSQ / 15.- 1./ 6.)*XX**2 )

      A( 2 ) = ( 0.1*XX**2 )*CTMP*( 1.- XX**2 / 14.) /
     &         ( 2.*CIORSQ + 3.- ( CIORSQ / 7.- 0.5 )*XX**2 )

      QSCA   = 6.* XX**4 * ( SQ( A(1) ) + SQ( B(1) ) +
     &                       FIVTHR*SQ( A(2) ) )
      GQSC   = 6.* XX**4 * REAL( A(1)*CONJG( A(2) + B(1) ) )
      QEXT   = QSCA
      IF( CALCQE ) QEXT   = 6.*XX*REAL( A(1) + B(1) +
     &                      FIVTHR*A(2) )

      RTMP   = 1.5 * XX**3
      SFORW  = RTMP*( A(1) + B(1) + FIVTHR*A(2) )
      SBACK  = RTMP*( A(1) - B(1) - FIVTHR*A(2) )
      TFORW( 1 ) = RTMP*( B(1) - FIVTHR*A(2) )
      TFORW( 2 ) = RTMP*( A(1) + 2.*FIVTHR*A(2) )
      TBACK( 1 ) = TFORW(1)
      TBACK( 2 ) = RTMP*( A(1) - 2.*FIVTHR*A(2) )

      DO 10 J = 1, NUMANG
         S1( J ) = RTMP*( A(1) + ( B(1) + FIVTHR*A(2) )*XMU( J ) )
         S2( J ) = RTMP*( B(1) + A(1)*XMU( J ) +
     &                    FIVTHR * A(2)*( 2.*XMU( J )**2 - 1.) )
   10 CONTINUE

c                                     ** Recover actual Mie coefficients
      A( 1 ) = XX**3 * A(1)
      A( 2 ) = XX**3 * A(2)
      B( 1 ) = XX**3 * B(1)
      B( 2 ) = ( 0., 0.)

      RETURN
      END

      SUBROUTINE TESTMI( COMPAR, XX, CREFIN, MIMCUT, PERFCT, ANYANG,
     &                   NUMANG, XMU, QEXT, QSCA, GQSC, SFORW, SBACK,
     &                   S1, S2, TFORW, TBACK )

c         Set up to run test case when  COMPAR = False;  when  = True,
c         compare Mie code test case results with correct answers
c         and abort if even one result is inaccurate.
c         (NoPMOM version)

c         The test case is :  Mie size parameter = 10
c                             refractive index   = 1.5 - 0.1 i
c                             scattering angle = 140 degrees
c                             1 Sekera moment

c         Results for this case may be found among the test cases
c         at the end of reference (1).

c         *** NOTE *** When running on some computers, esp. in single
c         precision, the Accur criterion below may have to be relaxed.
c         However, if Accur must be set larger than 10**-3 for some
c         size parameters, your computer is probably not accurate
c         enough to do Mie computations for those size parameters.

c     Routines called :  ERRMSG, MIPRNT, TSTBAD

      IMPLICIT  NONE

c     .. Scalar Arguments ..

      LOGICAL   ANYANG, COMPAR, PERFCT
      INTEGER   NUMANG
      REAL      GQSC, MIMCUT, QEXT, QSCA, XX
      COMPLEX   CREFIN, SBACK, SFORW
c     ..
c     .. Array Arguments ..

      REAL      XMU( * )
      COMPLEX   S1( * ), S2( * ), TBACK( * ), TFORW( * )
c     ..
c     .. Local Scalars ..

      LOGICAL   ANYSAV, OK, PERSAV
      INTEGER   N, NUMSAV
      REAL      ACCUR, CALC, EXACT, MIMSAV, TESTGQ, TESTQE, TESTQS,
     &          XMUSAV, XXSAV
      COMPLEX   CRESAV, TESTS1, TESTS2, TESTSB, TESTSF
c     ..
c     .. Local Arrays ..

      LOGICAL   PRNT( 2 )
      COMPLEX   TESTTB( 2 ), TESTTF( 2 )
c     ..
c     .. External Functions ..

      LOGICAL   TSTBAD
      EXTERNAL  TSTBAD
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG, MIPRNT
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, AIMAG, REAL
c     ..
c     .. Statement Functions ..

      LOGICAL   WRONG
c     ..
      SAVE      XXSAV, CRESAV, MIMSAV, PERSAV, ANYSAV, NUMSAV, XMUSAV

      DATA   TESTQE / 2.459791 /,
     $       TESTQS / 1.235144 /,
     $       TESTGQ / 1.139235 /,
     $       TESTSF / ( 61.49476, -3.177994 ) /,
     $       TESTSB / ( 1.493434, 0.2963657 ) /,
     $       TESTS1 / ( -0.1548380, -1.128972) /,
     $       TESTS2 / ( 0.05669755, 0.5425681) /,
     $       TESTTF / ( 12.95238, -136.6436 ), ( 48.54238, 133.4656 ) /,
     $       TESTTB / ( 41.88414, -15.57833 ), ( 43.37758, -15.28196 )/

      DATA      ACCUR / 1.E-4 /
c     ..
c     .. Statement Function definitions ..

      WRONG( CALC, EXACT ) = ABS( ( CALC - EXACT ) / EXACT ).GT.ACCUR
c     ..

      IF( .NOT.COMPAR ) THEN
c                                   ** Save certain user input values
         XXSAV  = XX
         CRESAV = CREFIN
         MIMSAV = MIMCUT
         PERSAV = PERFCT
         ANYSAV = ANYANG
         NUMSAV = NUMANG
         XMUSAV = XMU( 1 )
c                                   ** Reset input values for test case
         XX     = 10.0
         CREFIN = ( 1.5, -0.1 )
         MIMCUT = 0.0
         PERFCT = .FALSE.
         ANYANG = .TRUE.
         NUMANG = 1
         XMU( 1 ) = -0.7660444

      ELSE
c                                    ** Compare test case results with
c                                    ** correct answers and abort if bad
         OK     = .TRUE.

        IF ( WRONG( QEXT,TESTQE ) )
     $       OK =  TSTBAD( 'QEXT', ABS((QEXT - TESTQE) / TESTQE) )

        IF ( WRONG( QSCA,TESTQS ) )
     $       OK =  TSTBAD( 'QSCA', ABS((QSCA - TESTQS) / TESTQS) )

        IF ( WRONG( GQSC,TESTGQ ) )
     $       OK =  TSTBAD( 'GQSC', ABS((GQSC - TESTGQ) / TESTGQ) )

         IF( WRONG( REAL( SFORW ),REAL( TESTSF ) ) .OR.
     &       WRONG( AIMAG( SFORW ),AIMAG( TESTSF ) ) )
     &       OK = TSTBAD( 'SFORW', ABS( ( SFORW - TESTSF ) / TESTSF ) )

         IF( WRONG( REAL( SBACK ),REAL( TESTSB ) ) .OR.
     &       WRONG( AIMAG( SBACK ),AIMAG( TESTSB ) ) )
     &       OK = TSTBAD( 'SBACK', ABS( ( SBACK - TESTSB ) / TESTSB ) )

         IF( WRONG( REAL( S1( 1 ) ),REAL( TESTS1 ) ) .OR.
     &       WRONG( AIMAG( S1( 1 ) ),AIMAG( TESTS1 ) ) )
     &       OK = TSTBAD( 'S1', ABS( ( S1( 1 ) - TESTS1 ) / TESTS1 ) )

         IF( WRONG( REAL( S2( 1 ) ),REAL( TESTS2 ) ) .OR.
     &       WRONG( AIMAG( S2( 1 ) ),AIMAG( TESTS2 ) ) )
     &       OK = TSTBAD( 'S2', ABS( ( S2( 1 ) - TESTS2 ) / TESTS2 ) )

         DO 10 N = 1, 2

           IF ( WRONG(  REAL(TFORW(N)),  REAL(TESTTF(N)) ) .OR.
     $          WRONG( AIMAG(TFORW(N)), AIMAG(TESTTF(N)) ) )
     $          OK =  TSTBAD( 'TFORW', ABS( (TFORW(N) - TESTTF(N)) /
     $                                       TESTTF(N) ) )
           IF ( WRONG(  REAL(TBACK(N)),  REAL(TESTTB(N)) ) .OR.
     $          WRONG( AIMAG(TBACK(N)), AIMAG(TESTTB(N)) ) )
     $          OK =  TSTBAD( 'TBACK', ABS( (TBACK(N) - TESTTB(N)) /
     $                                        TESTTB(N) ) )

   10    CONTINUE

         IF( .NOT.OK ) THEN

            PRNT( 1 ) = .TRUE.
            PRNT( 2 ) = .TRUE.

            CALL  MIPRNT( PRNT, XX, PERFCT, CREFIN, NUMANG, XMU, QEXT,
     $                    QSCA, GQSC, SFORW, SBACK, TFORW, TBACK, S1,S2)

            CALL ERRMSG( 'MIEV0 -- Self-test failed', .TRUE.)

         END IF
c                                       ** Restore user input values
         XX     = XXSAV
         CREFIN = CRESAV
         MIMCUT = MIMSAV
         PERFCT = PERSAV
         ANYANG = ANYSAV
         NUMANG = NUMSAV
         XMU( 1 ) = XMUSAV

      END IF

      RETURN
      END

max_stars_repo_name: peterkuma/alcf

max_issues_repo_name: peterkuma/alcf

max_forks_repo_name: peterkuma/alcf


text:
      SUBROUTINE PSGEEQU( M, N, A, IA, JA, DESCA, R, C, ROWCND, COLCND,
     $                    AMAX, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      INTEGER            IA, INFO, JA, M, N
      REAL               AMAX, COLCND, ROWCND
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * )
      REAL               A( * ), C( * ), R( * )
*     ..
*
*  Purpose
*  =======
*
*  PSGEEQU computes row and column scalings intended to equilibrate an
*  M-by-N distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA:JA+N-1) and
*  reduce its condition number.  R returns the row scale factors and C
*  the column scale factors, chosen to try to make the largest entry in
*  each row and column of the distributed matrix B with elements
*  B(i,j) = R(i) * A(i,j) * C(j) have absolute value 1.
*
*  R(i) and C(j) are restricted to be between SMLNUM = smallest safe
*  number and BIGNUM = largest safe number.  Use of these scaling
*  factors is not guaranteed to reduce the condition number of
*  sub( A ) but works well in practice.
*
*  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
*  =========
*
*  M       (global input) INTEGER
*          The number of rows to be operated on i.e the number of rows
*          of the distributed submatrix sub( A ). M >= 0.
*
*  N       (global input) INTEGER
*          The number of columns to be operated on i.e the number of
*          columns of the distributed submatrix sub( A ). N >= 0.
*
*  A       (local input) REAL pointer into the local memory
*          to an array of dimension ( LLD_A, LOCc(JA+N-1) ), the
*          local pieces of the M-by-N distributed matrix whose
*          equilibration factors are to be computed.
*
*  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.
*
*  R       (local output) REAL array, dimension LOCr(M_A)
*          If INFO = 0 or INFO > IA+M-1, R(IA:IA+M-1) contains the row
*          scale factors for sub( A ). R is aligned with the distributed
*          matrix A, and replicated across every process column. R is
*          tied to the distributed matrix A.
*
*  C       (local output) REAL array, dimension LOCc(N_A)
*          If INFO = 0,  C(JA:JA+N-1) contains the column scale factors
*          for sub( A ). C is aligned with the distributed matrix A, and
*          replicated down every process row. C is tied to the distri-
*          buted matrix A.
*
*  ROWCND  (global output) REAL
*          If INFO = 0 or INFO > IA+M-1, ROWCND contains the ratio of
*          the smallest R(i) to the largest R(i) (IA <= i <= IA+M-1).
*          If ROWCND >= 0.1 and AMAX is neither too large nor too small,
*          it is not worth scaling by R(IA:IA+M-1).
*
*  COLCND  (global output) REAL
*          If INFO = 0, COLCND contains the ratio of the smallest C(j)
*          to the largest C(j) (JA <= j <= JA+N-1). If COLCND >= 0.1, it
*          is not worth scaling by C(JA:JA+N-1).
*
*  AMAX    (global output) REAL
*          Absolute value of largest distributed matrix element.  If
*          AMAX is very close to overflow or very close to underflow,
*          the matrix should be scaled.
*
*  INFO    (global output) INTEGER
*          = 0:  successful exit
*          < 0:  If the i-th argument is an array and the j-entry had
*                an illegal value, then INFO = -(i*100+j), if the i-th
*                argument is a scalar and had an illegal value, then
*                INFO = -i.
*          > 0:  If INFO = i,  and i is
*                <= M:  the i-th row of the distributed matrix sub( A )
*                       is exactly zero,
*                >  M:  the (i-M)-th column of the distributed
*                       matrix sub( A ) is exactly zero.
*
*  =====================================================================
*
*     .. 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 )
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      CHARACTER          COLCTOP, ROWCTOP
      INTEGER            I, IACOL, IAROW, ICOFF, ICTXT, IDUMM, IIA,
     $                   IOFFA, IROFF, J, JJA, LDA, MP, MYCOL, MYROW,
     $                   NPCOL, NPROW, NQ
      REAL               BIGNUM, RCMAX, RCMIN, SMLNUM
*     ..
*     .. Local Arrays ..
      INTEGER            DESCC( DLEN_ ), DESCR( DLEN_ )
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, CHK1MAT, DESCSET, IGAMX2D,
     $                   INFOG2L, PCHK1MAT, PB_TOPGET, PXERBLA, SGAMN2D,
     $                   SGAMX2D
*     ..
*     .. External Functions ..
      INTEGER            INDXL2G, NUMROC
      REAL               PSLAMCH
      EXTERNAL           INDXL2G, NUMROC, PSLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN, MOD
*     ..
*     .. Executable Statements ..
*
*     Get grid parameters
*
      ICTXT = DESCA( CTXT_ )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*     Test the input parameters.
*
      INFO = 0
      IF( NPROW.EQ.-1 ) THEN
         INFO = -(600+CTXT_)
      ELSE
         CALL CHK1MAT( M, 1, N, 2, IA, JA, DESCA, 6, INFO )
         CALL PCHK1MAT( M, 1, N, 2, IA, JA, DESCA, 6, 0, IDUMM, IDUMM,
     $                  INFO )
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( ICTXT, 'PSGEEQU', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
         ROWCND = ONE
         COLCND = ONE
         AMAX = ZERO
         RETURN
      END IF
*
      CALL PB_TOPGET( ICTXT, 'Combine', 'Rowwise', ROWCTOP )
      CALL PB_TOPGET( ICTXT, 'Combine', 'Columnwise', COLCTOP )
*
*     Get machine constants and local indexes.
*
      SMLNUM = PSLAMCH( ICTXT, 'S' )
      BIGNUM = ONE / SMLNUM
      CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, IIA, JJA,
     $              IAROW, IACOL )
      IROFF = MOD( IA-1, DESCA( MB_ ) )
      ICOFF = MOD( JA-1, DESCA( NB_ ) )
      MP = NUMROC( M+IROFF, DESCA( MB_ ), MYROW, IAROW, NPROW )
      NQ = NUMROC( N+ICOFF, DESCA( NB_ ), MYCOL, IACOL, NPCOL )
      IF( MYROW.EQ.IAROW )
     $   MP = MP - IROFF
      IF( MYCOL.EQ.IACOL )
     $   NQ = NQ - ICOFF
      LDA = DESCA( LLD_ )
*
*     Assign descriptors for R and C arrays
*
      CALL DESCSET( DESCR, M, 1, DESCA( MB_ ), 1, 0, 0, ICTXT,
     $               MAX( 1, MP ) )
      CALL DESCSET( DESCC, 1, N, 1, DESCA( NB_ ), 0, 0, ICTXT, 1 )
*
*     Compute row scale factors.
*
      DO 10 I = IIA, IIA+MP-1
         R( I ) = ZERO
   10 CONTINUE
*
*     Find the maximum element in each row.
*
      IOFFA = (JJA-1)*LDA
      DO 30 J = JJA, JJA+NQ-1
         DO 20 I = IIA, IIA+MP-1
            R( I ) = MAX( R( I ), ABS( A( IOFFA + I ) ) )
   20    CONTINUE
         IOFFA = IOFFA + LDA
   30 CONTINUE
      CALL SGAMX2D( ICTXT, 'Rowwise', ROWCTOP, MP, 1, R( IIA ),
     $              MAX( 1, MP ), IDUMM, IDUMM, -1, -1, MYCOL )
*
*     Find the maximum and minimum scale factors.
*
      RCMIN = BIGNUM
      RCMAX = ZERO
      DO 40 I = IIA, IIA+MP-1
         RCMAX = MAX( RCMAX, R( I ) )
         RCMIN = MIN( RCMIN, R( I ) )
   40 CONTINUE
      CALL SGAMX2D( ICTXT, 'Columnwise', COLCTOP, 1, 1, RCMAX, 1, IDUMM,
     $              IDUMM, -1, -1, MYCOL )
      CALL SGAMN2D( ICTXT, 'Columnwise', COLCTOP, 1, 1, RCMIN, 1, IDUMM,
     $              IDUMM, -1, -1, MYCOL )
      AMAX = RCMAX
*
      IF( RCMIN.EQ.ZERO ) THEN
*
*        Find the first zero scale factor and return an error code.
*
         DO 50 I = IIA, IIA+MP-1
            IF( R( I ).EQ.ZERO .AND. INFO.EQ.0 )
     $         INFO = INDXL2G( I, DESCA( MB_ ), MYROW, DESCA( RSRC_ ),
     $                NPROW ) - IA + 1
   50    CONTINUE
         CALL IGAMX2D( ICTXT, 'Columnwise', COLCTOP, 1, 1, INFO, 1,
     $                 IDUMM, IDUMM, -1, -1, MYCOL )
         IF( INFO.NE.0 )
     $      RETURN
      ELSE
*
*        Invert the scale factors.
*
         DO 60 I = IIA, IIA+MP-1
            R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
   60    CONTINUE
*
*        Compute ROWCND = min(R(I)) / max(R(I))
*
         ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
*
      END IF
*
*     Compute column scale factors
*
      DO 70 J = JJA, JJA+NQ-1
         C( J ) = ZERO
   70 CONTINUE
*
*     Find the maximum element in each column,
*     assuming the row scaling computed above.
*
      IOFFA = (JJA-1)*LDA
      DO 90 J = JJA, JJA+NQ-1
         DO 80 I = IIA, IIA+MP-1
            C( J ) = MAX( C( J ), ABS( A( IOFFA + I ) )*R( I ) )
   80    CONTINUE
         IOFFA = IOFFA + LDA
   90 CONTINUE
      CALL SGAMX2D( ICTXT, 'Columnwise', COLCTOP, 1, NQ, C( JJA ),
     $              1, IDUMM, IDUMM, -1, -1, MYCOL )
*
*     Find the maximum and minimum scale factors.
*
      RCMIN = BIGNUM
      RCMAX = ZERO
      DO 100 J = JJA, JJA+NQ-1
         RCMIN = MIN( RCMIN, C( J ) )
         RCMAX = MAX( RCMAX, C( J ) )
  100 CONTINUE
      CALL SGAMX2D( ICTXT, 'Columnwise', COLCTOP, 1, 1, RCMAX, 1, IDUMM,
     $              IDUMM, -1, -1, MYCOL )
      CALL SGAMN2D( ICTXT, 'Columnwise', COLCTOP, 1, 1, RCMIN, 1, IDUMM,
     $              IDUMM, -1, -1, MYCOL )
*
      IF( RCMIN.EQ.ZERO ) THEN
*
*        Find the first zero scale factor and return an error code.
*
         DO 110 J = JJA, JJA+NQ-1
            IF( C( J ).EQ.ZERO .AND. INFO.EQ.0 )
     $         INFO = M + INDXL2G( J, DESCA( NB_ ), MYCOL,
     $                DESCA( CSRC_ ), NPCOL ) - JA + 1
  110    CONTINUE
         CALL IGAMX2D( ICTXT, 'Columnwise', COLCTOP, 1, 1, INFO, 1,
     $                 IDUMM, IDUMM, -1, -1, MYCOL )
         IF( INFO.NE.0 )
     $      RETURN
      ELSE
*
*        Invert the scale factors.
*
         DO 120 J = JJA, JJA+NQ-1
            C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
  120    CONTINUE
*
*        Compute COLCND = min(C(J)) / max(C(J))
*
         COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
*
      END IF
*
      RETURN
*
*     End of PSGEEQU
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details
C=======================================================================
      SUBROUTINE RDEB1 (NDB, IDELB, NUMELB, NUMLNK, NUMATR,
     &   LINK, ATRIB, ATRNM, NAMLEN)
C=======================================================================

C   --*** RDEB1 *** (EXPLORE) Read database element block misc.
C   --
C   --RDEB1 reads the element block connectivity and attribute information
C   --from the database.  An error message is displayed if the end of file
C   --is read.
C   --
C   --Parameters:
C   --   NDB - IN - the database file
C   --   IDELB - IN - the element block number (for errors)
C   --   NUMELB - IN - the number of elements in the block
C   --   NUMLNK - IN - the number of nodes per element
C   --   NUMATR - IN - the number of attributes
C   --   LINK - OUT - the element connectivity for this block
C   --   ATRIB - OUT - the attributes for this block
C   --
      include 'exodusII.inc'

      INTEGER LINK(*)
      REAL ATRIB(*)

      CHARACTER*80 ERRMSG
      CHARACTER*(NAMLEN) ATRNM(*)

      IF (NUMELB .GT. 0 .AND. NUMLNK .GT. 0) THEN
        CALL EXGELC (NDB, IDELB, LINK, IERR)
        IF (IERR .NE. 0) GO TO 100
      END IF

      IF (NUMELB .GT. 0 .AND. NUMATR .GT. 0) THEN
        CALL EXGEAT (NDB, IDELB, ATRIB, IERR)
        IF (IERR .NE. 0) GO TO 110

        CALL EXGEAN (NDB, IDELB, NUMATR, ATRNM, IERR)
        IF (IERR .NE. 0) GO TO 115
      END IF

      RETURN

  100 CONTINUE
      WRITE (ERRMSG, 10000, IOSTAT=IDUM)
     *  'CONNECTIVITY for block', IDELB
      GOTO 120
  110 CONTINUE
      WRITE (ERRMSG, 10000, IOSTAT=IDUM)
     *  'ATTRIBUTES for block', IDELB
      GOTO 120
 115  CONTINUE
      WRITE (ERRMSG, 10000, IOSTAT=IDUM)
     *  'ATTRIBUTE NAMES for block', IDELB
      GOTO 120
  120 CONTINUE
      CALL WDBERR (IERR, ERRMSG)

      RETURN

10000  FORMAT (5 (A, I12))
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
      SUBROUTINE DGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
     $                   WORK, INFO )
*
*  -- LAPACK driver routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDB, M, N, NRHS, RANK
      DOUBLE PRECISION   RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            JPVT( * )
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  This routine is deprecated and has been replaced by routine DGELSY.
*
*  DGELSX computes the minimum-norm solution to a real linear least
*  squares problem:
*      minimize || A * X - B ||
*  using a complete orthogonal factorization of A.  A is an M-by-N
*  matrix which may be rank-deficient.
*
*  Several right hand side vectors b and solution vectors x can be
*  handled in a single call; they are stored as the columns of the
*  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
*  matrix X.
*
*  The routine first computes a QR factorization with column pivoting:
*      A * P = Q * [ R11 R12 ]
*                  [  0  R22 ]
*  with R11 defined as the largest leading submatrix whose estimated
*  condition number is less than 1/RCOND.  The order of R11, RANK,
*  is the effective rank of A.
*
*  Then, R22 is considered to be negligible, and R12 is annihilated
*  by orthogonal transformations from the right, arriving at the
*  complete orthogonal factorization:
*     A * P = Q * [ T11 0 ] * Z
*                 [  0  0 ]
*  The minimum-norm solution is then
*     X = P * Z' [ inv(T11)*Q1'*B ]
*                [        0       ]
*  where Q1 consists of the first RANK columns of Q.
*
*  Arguments
*  =========
*
*  M       (input) INTEGER
*          The number of rows of the matrix A.  M >= 0.
*
*  N       (input) INTEGER
*          The number of columns of the matrix A.  N >= 0.
*
*  NRHS    (input) INTEGER
*          The number of right hand sides, i.e., the number of
*          columns of matrices B and X. NRHS >= 0.
*
*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
*          On entry, the M-by-N matrix A.
*          On exit, A has been overwritten by details of its
*          complete orthogonal factorization.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,M).
*
*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
*          On entry, the M-by-NRHS right hand side matrix B.
*          On exit, the N-by-NRHS solution matrix X.
*          If m >= n and RANK = n, the residual sum-of-squares for
*          the solution in the i-th column is given by the sum of
*          squares of elements N+1:M in that column.
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B. LDB >= max(1,M,N).
*
*  JPVT    (input/output) INTEGER array, dimension (N)
*          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
*          initial column, otherwise it is a free column.  Before
*          the QR factorization of A, all initial columns are
*          permuted to the leading positions; only the remaining
*          free columns are moved as a result of column pivoting
*          during the factorization.
*          On exit, if JPVT(i) = k, then the i-th column of A*P
*          was the k-th column of A.
*
*  RCOND   (input) DOUBLE PRECISION
*          RCOND is used to determine the effective rank of A, which
*          is defined as the order of the largest leading triangular
*          submatrix R11 in the QR factorization with pivoting of A,
*          whose estimated condition number < 1/RCOND.
*
*  RANK    (output) INTEGER
*          The effective rank of A, i.e., the order of the submatrix
*          R11.  This is the same as the order of the submatrix T11
*          in the complete orthogonal factorization of A.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension
*                      (max( min(M,N)+3*N, 2*min(M,N)+NRHS )),
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            IMAX, IMIN
      PARAMETER          ( IMAX = 1, IMIN = 2 )
      DOUBLE PRECISION   ZERO, ONE, DONE, NTDONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, DONE = ZERO,
     $                   NTDONE = ONE )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
     $                   SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLANGE
      EXTERNAL           DLAMCH, DLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           DGEQPF, DLAIC1, DLASCL, DLASET, DLATZM, DORM2R,
     $                   DTRSM, DTZRQF, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN
*     ..
*     .. Executable Statements ..
*
      MN = MIN( M, N )
      ISMIN = MN + 1
      ISMAX = 2*MN + 1
*
*     Test the input arguments.
*
      INFO = 0
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
         INFO = -7
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGELSX', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( MIN( M, N, NRHS ).EQ.0 ) THEN
         RANK = 0
         RETURN
      END IF
*
*     Get machine parameters
*
      SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
      BIGNUM = ONE / SMLNUM
      CALL DLABAD( SMLNUM, BIGNUM )
*
*     Scale A, B if max elements outside range [SMLNUM,BIGNUM]
*
      ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
      IASCL = 0
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
*
*        Scale matrix norm up to SMLNUM
*
         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
         IASCL = 1
      ELSE IF( ANRM.GT.BIGNUM ) THEN
*
*        Scale matrix norm down to BIGNUM
*
         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
         IASCL = 2
      ELSE IF( ANRM.EQ.ZERO ) THEN
*
*        Matrix all zero. Return zero solution.
*
         CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
         RANK = 0
         GO TO 100
      END IF
*
      BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
      IBSCL = 0
      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
*
*        Scale matrix norm up to SMLNUM
*
         CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
         IBSCL = 1
      ELSE IF( BNRM.GT.BIGNUM ) THEN
*
*        Scale matrix norm down to BIGNUM
*
         CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
         IBSCL = 2
      END IF
*
*     Compute QR factorization with column pivoting of A:
*        A * P = Q * R
*
      CALL DGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), INFO )
*
*     workspace 3*N. Details of Householder rotations stored
*     in WORK(1:MN).
*
*     Determine RANK using incremental condition estimation
*
      WORK( ISMIN ) = ONE
      WORK( ISMAX ) = ONE
      SMAX = ABS( A( 1, 1 ) )
      SMIN = SMAX
      IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
         RANK = 0
         CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
         GO TO 100
      ELSE
         RANK = 1
      END IF
*
   10 CONTINUE
      IF( RANK.LT.MN ) THEN
         I = RANK + 1
         CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
     $                A( I, I ), SMINPR, S1, C1 )
         CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
     $                A( I, I ), SMAXPR, S2, C2 )
*
         IF( SMAXPR*RCOND.LE.SMINPR ) THEN
            DO 20 I = 1, RANK
               WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
               WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
   20       CONTINUE
            WORK( ISMIN+RANK ) = C1
            WORK( ISMAX+RANK ) = C2
            SMIN = SMINPR
            SMAX = SMAXPR
            RANK = RANK + 1
            GO TO 10
         END IF
      END IF
*
*     Logically partition R = [ R11 R12 ]
*                             [  0  R22 ]
*     where R11 = R(1:RANK,1:RANK)
*
*     [R11,R12] = [ T11, 0 ] * Y
*
      IF( RANK.LT.N )
     $   CALL DTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO )
*
*     Details of Householder rotations stored in WORK(MN+1:2*MN)
*
*     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
*
      CALL DORM2R( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ),
     $             B, LDB, WORK( 2*MN+1 ), INFO )
*
*     workspace NRHS
*
*     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
*
      CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
     $            NRHS, ONE, A, LDA, B, LDB )
*
      DO 40 I = RANK + 1, N
         DO 30 J = 1, NRHS
            B( I, J ) = ZERO
   30    CONTINUE
   40 CONTINUE
*
*     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
*
      IF( RANK.LT.N ) THEN
         DO 50 I = 1, RANK
            CALL DLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA,
     $                   WORK( MN+I ), B( I, 1 ), B( RANK+1, 1 ), LDB,
     $                   WORK( 2*MN+1 ) )
   50    CONTINUE
      END IF
*
*     workspace NRHS
*
*     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
*
      DO 90 J = 1, NRHS
         DO 60 I = 1, N
            WORK( 2*MN+I ) = NTDONE
   60    CONTINUE
         DO 80 I = 1, N
            IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN
               IF( JPVT( I ).NE.I ) THEN
                  K = I
                  T1 = B( K, J )
                  T2 = B( JPVT( K ), J )
   70             CONTINUE
                  B( JPVT( K ), J ) = T1
                  WORK( 2*MN+K ) = DONE
                  T1 = T2
                  K = JPVT( K )
                  T2 = B( JPVT( K ), J )
                  IF( JPVT( K ).NE.I )
     $               GO TO 70
                  B( I, J ) = T1
                  WORK( 2*MN+K ) = DONE
               END IF
            END IF
   80    CONTINUE
   90 CONTINUE
*
*     Undo scaling
*
      IF( IASCL.EQ.1 ) THEN
         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
         CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
     $                INFO )
      ELSE IF( IASCL.EQ.2 ) THEN
         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
         CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
     $                INFO )
      END IF
      IF( IBSCL.EQ.1 ) THEN
         CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
      ELSE IF( IBSCL.EQ.2 ) THEN
         CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
      END IF
*
  100 CONTINUE
*
      RETURN
*
*     End of DGELSX
*
      END

max_stars_repo_name: pkeuchel/oomph-lib

max_issues_repo_name: pkeuchel/oomph-lib

max_forks_repo_name: pkeuchel/oomph-lib


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998 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.
!
      subroutine checkjac(cotet,node,pnew,kontet,c1,jflag,
     &     iedtet,iedgmid,ipoeled,ieled,iedge,ibadnodes,nbadnodes,
     &     iwrite)
!
!     check the Jacobian in all integration points of all elements
!     belonging to midnode "node" if the coordinates of "node" are
!     changed from cotet(1..3,node) to pnew(1..3); if the Jacobian
!     is strictly positive, the coordinates are changed, if not, they
!     are left unchanged
!     "node" belongs to edge "iedge"
!
      implicit none
!
      integer node,kontet(4,*),i,ielem,indexe,iflag,jflag,iwrite,
     &     j,n,iedtet(6,*),iedgmid(*),ipoeled(*),ieled(2,*),
     &     iedge,k,m,ibadnodes(*),nbadnodes,id,listed,neighedge
!
      real*8 cotet(3,*),pnew(3),pold(3),xsj,c1,c2,xi,et,ze,xl(3,10),
     &     shp(4,10)
!
      include "gauss.f"
!
      iflag=2
!
!     backup of the old coordinates
!
      do i=1,3
        pold(i)=cotet(i,node)
      enddo
!
!     storing the new coordinates in cotet
!
      c2=c1
      do i=1,3
        cotet(i,node)=c2*pnew(i)+(1.d0-c2)*pold(i)
      enddo
!
!     trial loops; if bad elements the projection is reduced and the
!     elements are rechecked
!
      n=3
      loop1: do j=1,n
      indexe=ipoeled(iedge)
!
!       loop over all elements belonging to edge "iedge"
!
        loop2: do
          if(indexe.eq.0) exit loop1
          ielem=ieled(1,indexe)
!
!         vertex nodes
!
          do k=1,4
            do m=1,3
              xl(m,k)=cotet(m,kontet(k,ielem))
            enddo
          enddo
!
!         middle nodes
!
          do k=1,6
            do m=1,3
              xl(m,k+4)=cotet(m,iedgmid(iedtet(k,ielem)))
            enddo
          enddo
!
          do k=1,4
            xi=gauss3d5(1,k)
            et=gauss3d5(2,k)
            ze=gauss3d5(3,k)
            call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
            if(xsj.le.0.d0) then
c              write(*,*) 'checkjac bad element ',ielem,' node ',node
              exit loop2
            endif
          enddo
          indexe=ieled(2,indexe)
        enddo loop2
!
!       last loop: revert to non-projected coordinates
!
        if(j.eq.n) then
          do i=1,3
            cotet(i,node)=pold(i)
          enddo
          exit
        endif
!
!       not last loop: reduce the projection and retry
!
        c2=c2/2.d0
        do i=1,3
          cotet(i,node)=c2*pnew(i)+(1.d0-c2)*pold(i)
        enddo
      enddo loop1
!
!     if node was not completely projected: list all edges of
!     neighboring elements (corresponds 1-to-1 to midnodes)
!     -> application of fminsirefine for all neighboring edges/midnodes
!        which are not on the free surface
!
      if(j.ne.1) then
        indexe=ipoeled(iedge)
        do
          if(indexe.eq.0) exit
!
!         adjacent element
!
          ielem=ieled(1,indexe)
          do k=1,6
            neighedge=iedtet(k,ielem)
!
!           edge of adjacent element
!
            listed=0
            call nident(ibadnodes,neighedge,nbadnodes,id)
            if(id.gt.0) then
              if(ibadnodes(id).eq.neighedge) then
                listed=1
              endif
            endif
            if(listed.eq.0) then
              nbadnodes=nbadnodes+1
              do m=nbadnodes,id+2,-1
                ibadnodes(m)=ibadnodes(m-1)
              enddo
              ibadnodes(id+1)=neighedge
            endif
          enddo
          indexe=ieled(2,indexe)
        enddo
      endif
!
!     if jflag=1 and j>1: last chance for a complete projection
!     did not work
!
      if((jflag.eq.1).and.(j.ne.1)) then
        write(*,*) '*WARNING in checkjac: projection of midnode ',node
        write(*,*) '         had to be reduced to keep the adjacent'
        write(*,*) '         elements regular'
        write(*,*)
        write(40,*) node
        iwrite=1
      endif
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop