text:
C***********************************************************************
C    Module:  vortex.f
C
C    Copyright (C) 2011 Mark Drela
C
C    This program is free software; you can redistribute it and/or modify
C    it under the terms of the GNU General Public License as published by
C    the Free Software Foundation; either version 2 of the License, or
C    (at your option) any later version.
C
C    This program is distributed in the hope that it will be useful,
C    but WITHOUT ANY WARRANTY; without even the implied warranty of
C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C    GNU General Public License for more details.
C
C    You should have received a copy of the GNU General Public License
C    along with this program; if not, write to the Free Software
C    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
C***********************************************************************

      SUBROUTINE VRTXCO(IMAX,II, NBLDS, LDUCT, RAKE,
     &                  XI,XV,GAM, ADW,VIND_GAM,VIND_ADW)
C
      PARAMETER (NTDIM=5000)
      DIMENSION XI(IMAX), XV(IMAX), GAM(IMAX)
      DIMENSION VIND_ADW(3,IMAX), VIND_GAM(3,IMAX,IMAX)
C
      DIMENSION A(3), B(3), UVW(3), UVW_A(3,3), UVW_B(3,3)
      DIMENSION VSUM(3), VADW(3)
      DIMENSION THETSPC(5000)
C
      LOGICAL LDUCT
C-----------------------------------------
C     Calculates "Vortex Momentum"
C     Gamma-swirl influence coefficients
C
C     Input:
C       IMAX         array dimension
C       II           number of radial points on blade (circulation stations)
C       NN           number of Fourier harmonics
C       NBLDS        number of blades
C       LDUCT        T for duct outer BC
C       RAKE         blade rake angle from Y axis in XY plane
C       XI(i)        r/R control point coordinate array
C       XV(i)        r/R vortex  leg   coordinate array
C       GAM(i)       circulation array
C       ADW          wake advance ratio  V/wR
C
C     Output:
C       VT_GAM(i,j)  sensitivity of swirl velocity to circulation
C       VT_ADW(i)    sensitivity of swirl velocity to wake advance ratio
C-----------------------------------------
      BLDS = FLOAT(NBLDS)
C
      PI = 4.0*ATAN(1.0)
C
      XI0   = XV(1)
      XITIP = XV(II+1)
      TANRAK = TAN(RAKE)

      DO I=1, II
        DO J=1, II
          VIND_GAM(1,I,J) = 0.
          VIND_GAM(2,I,J) = 0.
          VIND_GAM(3,I,J) = 0.
        ENDDO
        VIND_ADW(1,I) = 0.
        VIND_ADW(2,I) = 0.
        VIND_ADW(3,I) = 0.
      ENDDO
C
C--- Set up variable theta spacing for near, intermediate and far field
      DTH1  = PI / 60.
      RAD1  = 2.0
      THET1 = RAD1/ADW
      DTH2  = PI / 20.
      RAD2  = 4.0
      THET2 = RAD2/ADW
      DTH3  = PI / 8.
      RAD3  = 50.0
      THET3 = RAD3/ADW
C
      N1 = IFIX(2.0*(THET1      )/(DTH2+DTH1))
      N2 = IFIX(2.0*(THET2-THET1)/(DTH3+DTH2))
      DDTH1 = (DTH2-DTH1)/FLOAT(N1-1)
      DDTH2 = (DTH3-DTH2)/FLOAT(N2-1)
C
      THET = 0.0
      DTH = DTH1
      DO I = 1, NTDIM
        IF(THET.LT.THET1) THEN
          THETSPC(I) = THET
          THET = THET + DTH
          DTH = DTH + DDTH1
         ELSEIF(THET.LT.THET2) THEN
          THETSPC(I) = THET
          THET = THET + DTH
          DTH = DTH + DDTH2
         ELSEIF(THET.LT.THET3) THEN
          THETSPC(I) = THET
          DTH = DTH3
          THET = THET + DTH
         ELSE
          NTHET = I-1
          GO TO 100
        ENDIF
      END DO
      WRITE(*,*) 'Too many vortex segments for spacing array'
      NTHET = NTDIM
 100  CONTINUE
C
      IF(LDUCT) THEN
C----- use simple mean swirl to get swirl at blade
       DO I=1, II
         DO J=1, II
           VIND_GAM(1,I,J) = 0.
           VIND_GAM(2,I,J) = 0.
           VIND_GAM(3,I,J) = 0.
         ENDDO
         VIND_GAM(3,I,I) =  BLDS/(4.0*PI*XI(I))
         VIND_ADW(3,I)   =  0.0
         VIND_ADW(2,I)   =  0.0
         VIND_GAM(1,I,I) =  VIND_GAM(3,I,I)*XI(I) /ADW
         VIND_ADW(1,I)   = -VIND_GAM(1,I,I)*GAM(I)/ADW
       ENDDO
C
      ELSE
C
C----- Do a discrete vortex integration of slipstream vortices
       DTBLD = 2.0*PI/FLOAT(NBLDS)
       write(*,20) nblds*nthet
 20    format(/'Vortex points/radial station = ',I6)
C
C--- velocity influences for point - R0
       DO I = 1, II
         R0X = XI(I)*TANRAK
         R0Y = XI(I)
         R0Z = 0.0
c         write(*,*) 'CP ',r0x,r0y,r0z
C
         VIND_ADW(1,I) = 0.0
         VIND_ADW(2,I) = 0.0
         VIND_ADW(3,I) = 0.0
C
C----- For each vortex trailing leg (II+1 legs starting at XV(J))
         DO J = 1, II+1
C
           VSUM(1) = 0.0
           VSUM(2) = 0.0
           VSUM(3) = 0.0
           VADW(1) = 0.0
           VADW(2) = 0.0
           VADW(3) = 0.0
C
           RV   = XV(J)
           XXV  = RV*TANRAK
C
C----- For each blade
           THETOFF  = 0.0
           DO N = 1, NBLDS
C
C----- For angles around helix to the far-field
             THET1 = THETSPC(1)
             R1X = XXV + THET1*XITIP*ADW
             R1Y = RV*COS(THET1+THETOFF)
             R1Z = RV*SIN(THET1+THETOFF)
             R1_ADW = THET1*XITIP
C
             DO L = 1, NTHET-1
C
               THET2 = THETSPC(L+1)
               R2X = XXV + THET2*XITIP*ADW
               R2Y = RV*COS(THET2+THETOFF)
               R2Z = RV*SIN(THET2+THETOFF)
               R2_ADW = THET2*XITIP
C
               A(1) = R1X - R0X
               A(2) = R1Y - R0Y
               A(3) = R1Z - R0Z
               B(1) = R2X - R0X
               B(2) = R2Y - R0Y
               B(3) = R2Z - R0Z
               CALL VORSEGVEL(A, B, UVW, UVW_A, UVW_B )
C
               VSUM(1) = VSUM(1) + UVW(1)
               VSUM(2) = VSUM(2) + UVW(2)
               VSUM(3) = VSUM(3) + UVW(3)
C
               VADW(1) = VADW(1) + UVW_A(1,1)*R1_ADW
     &                           + UVW_B(1,1)*R2_ADW
               VADW(2) = VADW(2) + UVW_A(2,1)*R1_ADW
     &                           + UVW_B(2,1)*R2_ADW
               VADW(3) = VADW(3) + UVW_A(3,1)*R1_ADW
     &                           + UVW_B(3,1)*R2_ADW
C
               THET1 = THET2
               R1X = R2X
               R1Y = R2Y
               R1Z = R2Z
               R1_ADW = R2_ADW
C
c               if(i.eq.1) then
c                write(88,10) r1x,r1y,r1z
c                write(*,10) 'r1 ',i,j,n,l,r2x,r2y,r2z
c               endif
 10            format(3f10.3)
 12            format(A,2I6,2f10.3)

             END DO ! L loop
C
             THETOFF = THETOFF + DTBLD
C
           END DO ! N loop
C
           VSUM(3) = -VSUM(3)
           VADW(3) = -VADW(3)
C---- alternate + and - influence for each vortex line into velocity
C     influence matrix
C
C---Open wake, interdigitate all vortex lines
           IF(J.LE.II) THEN
             VIND_GAM(1,I,J) = -VSUM(1)
             VIND_GAM(2,I,J) = -VSUM(2)
             VIND_GAM(3,I,J) = -VSUM(3)
             VIND_ADW(1,I) = VIND_ADW(1,I) - VADW(1)*GAM(J)
             VIND_ADW(2,I) = VIND_ADW(2,I) - VADW(2)*GAM(J)
             VIND_ADW(3,I) = VIND_ADW(3,I) - VADW(3)*GAM(J)
           ENDIF
           IF(J.GT.1) THEN
             VIND_GAM(1,I,J-1) = VIND_GAM(1,I,J-1) + VSUM(1)
             VIND_GAM(2,I,J-1) = VIND_GAM(2,I,J-1) + VSUM(2)
             VIND_GAM(3,I,J-1) = VIND_GAM(3,I,J-1) + VSUM(3)
             VIND_ADW(1,I) = VIND_ADW(1,I) + VADW(1)*GAM(J-1)
             VIND_ADW(2,I) = VIND_ADW(2,I) + VADW(2)*GAM(J-1)
             VIND_ADW(3,I) = VIND_ADW(3,I) + VADW(3)*GAM(J-1)
           ENDIF
C
         END DO ! J loop
C
       END DO ! I loop
C
      ENDIF
C
      RETURN
      END ! VRTXCO

      SUBROUTINE VORSEGVEL(A, B, UVW, UVW_A, UVW_B )
C-------------------------------------------------------------------
C     Calculates the velocity induced by a vortex segment
C     of unit strength, with no core radius.
C
C     The point where the velocity is calculated is at 0,0,0.
C
C     Positive circulation is by righthand rule from A to B.
C
C  Input:
C     A(3)    coordinates of vertex #1 of the vortex
C     B(3)    coordinates of vertex #2 of the vortex
C
C  Output:
C     UVW(3)      induced velocity
C     UVW_A(3,3)  dUVW/dA  sensitivity
C     UVW_B(3,3)  dUVW/dB  sensitivity
C
C-------------------------------------------------------------------
      DIMENSION A(3), B(3), UVW(3), UVW_A(3,3), UVW_B(3,3)
C
C
      DIMENSION AXB(3), AXB_A(3,3), AXB_B(3,3)
C
      DATA PI4  / 12.56637062 /
C
      ASQ = A(1)**2 + A(2)**2 + A(3)**2
      BSQ = B(1)**2 + B(2)**2 + B(3)**2
C
      AMAG = SQRT(ASQ)
      BMAG = SQRT(BSQ)
C
      DO K = 1, 3
        UVW(K) = 0.
        DO L = 1, 3
          UVW_A(K,L) = 0.
          UVW_B(K,L) = 0.
        ENDDO
      ENDDO
C
C---- contribution from the vortex leg
      IF (AMAG*BMAG .NE. 0.0) THEN
        AXB(1) = A(2)*B(3) - A(3)*B(2)
        AXB(2) = A(3)*B(1) - A(1)*B(3)
        AXB(3) = A(1)*B(2) - A(2)*B(1)
C
        AXB_A(1,1) =  0.0
        AXB_A(1,2) =  B(3)
        AXB_A(1,3) = -B(2)
        AXB_A(2,1) = -B(3)
        AXB_A(2,2) =  0.0
        AXB_A(2,3) =  B(1)
        AXB_A(3,1) =  B(2)
        AXB_A(3,2) = -B(1)
        AXB_A(3,3) =  0.0
C
        AXB_B(1,1) =  0.0
        AXB_B(1,2) = -A(3)
        AXB_B(1,3) =  A(2)
        AXB_B(2,1) =  A(3)
        AXB_B(2,2) =  0.0
        AXB_B(2,3) = -A(1)
        AXB_B(3,1) = -A(2)
        AXB_B(3,2) =  A(1)
        AXB_B(3,3) =  0.0
C
        ADB = A(1)*B(1) + A(2)*B(2) + A(3)*B(3)
C
        DEN     =     AMAG*BMAG + ADB
        DEN_ASQ = 0.5*BMAG/AMAG
        DEN_BSQ = 0.5*AMAG/BMAG
C
        T = (1.0/AMAG + 1.0/BMAG) / DEN
C
        T_ADB = -T/DEN
        T_ASQ = -T/DEN*DEN_ASQ - 0.5/(DEN*AMAG*ASQ)
        T_BSQ = -T/DEN*DEN_BSQ - 0.5/(DEN*BMAG*BSQ)
C
        DO K = 1, 3
          UVW(K) = AXB(K)*T
C
          DO L = 1, 3
            UVW_A(K,L) = AXB(K)*T_ASQ  *  A(L)*2.0
     &                 + AXB(K)*T_ADB  *  B(L)
     &                 + AXB_A(K,L)*T
            UVW_B(K,L) = AXB(K)*T_BSQ  *  B(L)*2.0
     &                 + AXB(K)*T_ADB  *  A(L)
     &                 + AXB_B(K,L)*T
          ENDDO
        ENDDO
      ENDIF
C
      DO K = 1, 3
        UVW(K) = UVW(K)/PI4
        DO L = 1, 3
          UVW_A(K,L) = UVW_A(K,L)/PI4
          UVW_B(K,L) = UVW_B(K,L)/PI4
        ENDDO
      ENDDO
C
      RETURN
      END ! VORVEL


text:
      SUBROUTINE pyrmd4(dvol,x,xcen)
!---------------------------------------------------
!***  calculates volume of a pyramid with a quad base (1GP scheme)
!---------------------------------------------------

      IMPLICIT NONE

      REAL(kind=8):: dvol,x(3,4),xcen(3)

      REAL(kind=8) :: x1(3),x2(3),x3(3),xp(3),xcq(3)

      x1(:) = -x(:,1)+x(:,2)+x(:,3)-x(:,4)
      x2(:) = -x(:,1)-x(:,2)+x(:,3)+x(:,4)
      xcq(:) = .25d0 * (x(:,1)+x(:,2)+x(:,3)+x(:,4))
      x3(:) = xcq(:)-xcen(:)
      CALL vecpro(x1,x2,xp)
      dvol = DOT_PRODUCT(x3,xp)
      dvol = dvol/12.d0

      END SUBROUTINE pyrmd4

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


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

C=======================================================================
      SUBROUTINE TRNXYZ (XN, YN, XN3, YN3, ZN3, IXNP, NRNP, ZCORD)
C=======================================================================

C   --*** NEWXYZ *** (GEN3D) Calculate 3D coordinates
C   --   Written by Amy Gilkey - revised 05/09/88
C   --   Modified by Greg Sjaardema - 02/06/89
C   --       Added Warp Function
C   --       Added Gradient to Rotations (not for center blocks)
C   --
C   --NEWXYZ calculates the coordinate array for the 3D database.
C   --
C   --Parameters:
C   --   XN, YN - IN - the 2D coordinates, destroyed
C   --   XN3, YN3, ZN3 - OUT - the 3D coordinates
C   --   IXNP - IN - the new index for each node
C   --   NRNP - IN - the number of new nodes generated for each node
C   --   NPCEN - IN - the node numbers of the center nodes by column and row
C   --   ZCORD - SCRATCH - size = NNREPL, holds z coordinate for transformations
C   --
C   --Common Variables:
C   --   Uses NDIM, NUMNP of /DBNUMS/
C   --   Uses NDIM3, NUMNP3 of /DBNUM3/
C   --   Uses DOTRAN, NNREPL, DIM3, NRTRAN, D3TRAN, ZGRAD,
C   --      CENTER, NUMCOL, NUMROW of /PARAMS/
C   --   Uses XOFFS, YOFFS, ZOFFS of /XYZOFF/
C   --   Uses ROT3D, ROTMAT of /XYZROT/

      INCLUDE 'g3_dbnums.blk'
      INCLUDE 'g3_dbnum3.blk'
      INCLUDE 'g3_params.blk'

      REAL XN(NUMNP), YN(NUMNP),
     &   XN3(NUMNP3), YN3(NUMNP3), ZN3(NUMNP3)
      INTEGER IXNP(*), NRNP(*)
      REAL ZCORD(NNREPL)

C   --Copy Y coordinate from original

         DO 520 INP = 1, NUMNP
            JNP0 = IXNP(INP) - 1
            DO 510 NR = 1, NRNP(INP)
               YN3(JNP0+NR) = YN(INP)
  510       CONTINUE
  520    CONTINUE

C      --For translations, repeat the X coordinate and add a Z coordinate

      IBLK = 0
      NXTNR = 1

      ZEND = 0.0
    1    CONTINUE
         IBLK = IBLK + 1
         IF (NRTRAN(IBLK) .GT. 0) THEN
            ZBEG = ZEND
            ZEND = ZBEG + D3TRAN(IBLK)
            CALL INIGRD (-ZBEG, -ZEND, ZGRAD(IBLK),
     *         NRTRAN(IBLK), NRTRAN(IBLK)+1, ZCORD(NXTNR) )
            NXTNR = NXTNR + NRTRAN(IBLK)
            IF (IBLK .LT. MAXINT) GO TO 1
         END IF

C      --Repeat the X coordinate and add the calculated Z coordinate

      DO 30 INP = 1, NUMNP
         JNP = IXNP(INP) - 1
         DO 20 NR = 1, NNREPL
            XN3(JNP+NR) = XN(INP)
            ZN3(JNP+NR) = ZCORD(NR)
   20    CONTINUE
   30 CONTINUE
      RETURN
      END

ext: f

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
       subroutine iff_save(string)
c
c//////////////////////////////////////////////////////////////////////
c Copyright (c) 1997--2000 Matthew Newville, The University of Chicago
c Copyright (c) 1992--1996 Matthew Newville, University of Washington
c
c Permission to use and redistribute the source code or binary forms of
c this software and its documentation, with or without modification is
c hereby granted provided that the above notice of copyright, these
c terms of use, and the disclaimer of warranty below appear in the
c source code and documentation, and that none of the names of The
c University of Chicago, The University of Washington, or the authors
c appear in advertising or endorsement of works derived from this
c software without specific prior written permission from all parties.
c
c THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
c EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
c MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
c IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
c CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
c TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
c SOFTWARE OR THE USE OR OTHER DEALINGS IN THIS SOFTWARE.
c//////////////////////////////////////////////////////////////////////
c
c  save session to a  PAD save file
       implicit none
       include 'consts.h'
       include 'keywrd.h'
       include 'arrays.h'
       include 'encod.h'
       save
       integer       max_groups, igr
       parameter (max_groups=256)
       integer        i, k, istrln,  ndfkey, ilen, inn, get_array
       integer        npack, iex, ier, lun, isn, isf, iff_eval_in
       integer        j,jl
       character*512  string*(*), file, s2, build
       character*256  groups(max_groups)
       character*10   stat, spad*20, vers, defkey(2)*64
       logical s_str, s_sca, s_arr, s_sys, s_psc, dsave
       double precision getsca
       parameter (stat = 'unknown' , vers = '1.02')
       external istrln, iff_eval_in, get_array, getsca

       if (int(getsca('&sync_level',0)).ge.1)  call iff_sync
       file  = 'ifeffit.sav'
       call  gettxt('&build', build)
       npack = 8
       s_str = .true.
       s_sca = .true.
       s_arr = .true.
       s_psc = .false.
       s_sys = .false.
       igr   = 0
       do 10 i = 1, max_groups
          groups(i) = undef
 10    continue
c  interpret any and all keyword/value pairs for setting options
       call bkeys(string, mkeys, keys, values, nkeys)
       ndfkey    = 1
       defkey(1) = 'file'
       do 100 i = 1, nkeys
          k = istrln( keys(i))
          if ((values(i).eq.undef).and.(i.le.ndfkey)) then
             values(i) = keys(i)
             keys(i)   = defkey(i)
          end if
          if (keys(i).eq.'file')  then
             file = values(i)
          elseif (keys(i).eq.'group')  then
             igr = igr + 1
             if (igr.le.max_groups) groups(igr) = values(i)
          elseif ((keys(i).eq.'npad').or.(keys(i).eq.'npack')) then
             ier = iff_eval_in(values(i),npack)
          elseif (keys(i).eq.'no_strings')  then
             s_str = .false.
          elseif (keys(i).eq.'with_strings')  then
             s_str = .true.
          elseif (keys(i).eq.'no_scalars')  then
             s_sca = .false.
          elseif (keys(i).eq.'with_scalars')  then
             s_sca = .true.
          elseif (keys(i).eq.'pad_scalars')  then
             s_psc = .true.
          elseif (keys(i).eq.'no_pad_scalars')  then
             s_psc = .false.
          elseif (keys(i).eq.'no_arrays')  then
             s_arr = .false.
          elseif (keys(i).eq.'with_arrays')  then
             s_arr = .true.
          elseif (keys(i).eq.'no_sys')  then
             s_sys = .false.
          elseif (keys(i).eq.'with_sys')  then
             s_sys = .true.
          elseif (keys(i).eq.'arrays_only')  then
             s_sca = .false.
             s_str = .false.
             s_arr = .true.
          elseif (keys(i).eq.'scalars_only')  then
             s_sca = .true.
             s_str = .false.
             s_arr = .false.
          elseif (keys(i).eq.'strings_only')  then
             s_sca = .false.
             s_str = .true.
             s_arr = .false.
          else
             messg = keys(i)(1:k)//' " will be ignored'
             call warn(1,' *** save: unknown keyword " '//messg)
          end if
 100   continue
       call undels(file)
       iex  = 0
       ier  = 0
       lun  = 3
       call openfl(lun, file, stat, iex, ier)
       if ((ier.lt.0).or.(iex.lt.0)) then
          call warn(2,' *** save: error opening file ')
          return
       end if
       write(lun,'(3a,i3)') '#:IFEFFIT SAVE File: v',vers,
     $      '  npad =',npack
       ilen = istrln(build)
       write(lun,'(2a)') '#:build ', build(1:ilen)
       write(lun,'(a)')  '#:strings'
       if (s_str) then
          do 210 i = 1, maxtxt
             if ((txtnam(i).ne.' ').and. (text(i).ne.' ') .and.
     $            (txtnam(i)(1:1).ne.'&')) then
                isn = max(1,istrln(txtnam(i)))
                isf = max(1,istrln(text(i)))
                write(lun,'(4a)') '#= ',txtnam(i)(1:isn),
     $               '= ',text(i)(1:isf)
             endif
 210      continue
       end if
c
       write(lun,'(a)') '#:system'
       if (s_sys) then
          do 220 i = 1, maxsca
             if ((scanam(i).ne.' ').and.(scanam(i)(1:1).eq.'&')) then
                call pad(scalar(i),npack,spad)
                isn = max(1,istrln(scanam(i)))
                s2  = scafrm(i)
                isf = istrln(s2)
                if (icdsca(1,i).eq.-1) then
                   write(lun,'(4a)') '#? ', scanam(i)(1:isn), ' ',
     $                  spad(1:npack)
                elseif ((s2.eq.' ').or.(isf.eq.0)) then
                   write(lun,'(4a)') '#% ', scanam(i)(1:isn), ' ',
     $                  spad(1:npack)
                else
                   write(lun,'(4a)') '#= ', scanam(i)(1:isn), ' = ',
     $                  s2(1:isf)
                end if
             end if
 220      continue
       end if

       if (s_psc) then
          write(lun,'(a)') '#:scalars (pad)'
       else
          write(lun,'(a)') '#:scalars (nopad)'
       endif
       if (s_sca) then
          do 230 i = 1, maxsca
             if ((scanam(i).ne.' ').and.(scanam(i)(1:1).ne.'&')) then
                call pad(scalar(i),npack,spad)
                spad(npack+1:)= ' '
                if (.not.s_psc) write(spad, '(g19.13)') scalar(i)
                isn = max(1,istrln(scanam(i)))
                s2  = scafrm(i)
                isf = istrln(s2)
                if (icdsca(1,i).eq.-1) then
                   write(lun,'(4a)') '#? ', scanam(i)(1:isn),' ',spad
                elseif ((s2.eq.' ').or.(isf.eq.0)) then
                   write(lun,'(4a)') '#% ', scanam(i)(1:isn),' ',spad
                else
                   write(lun,'(4a)') '#= ', scanam(i)(1:isn), ' = ',
     $                  s2(1:isf)
                end if
             end if
 230      continue
       end if
c arrays
       write(lun,'(a)') '#:arrays'
       if (s_arr) then
          do 340 i = 1, maxarr-1
             if ((arrnam(i).ne.blank) .and. (arrnam(i).ne.undef)) then
                dsave = .true.
                if (igr.ge.1) then
                   dsave = .false.
                   j     = index(arrnam(i),'.')
                   do 305 k = 1, igr
                      jl = istrln(groups(k))
                      if (arrnam(i)(1:j-1).eq.groups(k)(1:jl)) then
                         dsave = .true.
                         goto 307
                      endif
 305               continue
 307               continue
                endif
                if (dsave) then
                   isn = max(1,istrln(arrnam(i)))
                   s2  = arrfrm(i)
                   isf = istrln(arrfrm(i))
                   if ((s2.eq.' ').or.(isf.eq.0)) then
                      inn  = get_array(arrnam(i),'',0,tmparr)
                      write(lun,'(a,i6,2a)') '#% [', inn, '] ',
     $                     arrnam(i)(1:isn)
                      call wrpadd(lun,npack,tmparr,inn)
                   else
                      write(lun,'(a,i6,4a)') '#= [', narray(i), '] ',
     $                     arrnam(i)(1:isn), ' = ', s2(1:isf)
                   end if
                end if
             endif
 340      continue
       end if
       if (lun.gt.0) close(lun)
       return
       end
c
       subroutine iff_restore(string)
c
c  restore session from a  PAD save file
       implicit none
       include 'consts.h'
       include 'keywrd.h'
       include 'arrays.h'
       include 'encod.h'
       save
       integer   i, j, k, istrln,  ndfkey, ierr, in1, in2
       integer   npack, iex, ier, lun, nwords, ntmp, iread_ky
       integer   iofstr, iofsca, ilen, i_sync
       character*512  str, file, s1, s2, string*(*)
       character*256  defkey(2)*64
       character*10   stat, pre*2, vers, words(2)*(20)
       double precision unpad, xvers, getsca
       logical      sc_pad
       parameter (stat = 'old')
       external unpad, istrln, iofstr, iofsca, iread_ky, getsca
       file = 'ifeffit.sav'
c  interpret any and all keyword/value pairs for setting options
       call bkeys(string, mkeys, keys, values, nkeys)
cc       print*, ' iff_restore : ', string(1:70)
cc       print*, ' iff_restore : ', file(1:40)
       ndfkey    = 1
       defkey(1) = 'file'
       do 100 i = 1, nkeys
          k = istrln( keys(i))
          if ((values(i).eq.undef).and.(i.le.ndfkey)) then
             values(i) = keys(i)
             keys(i)   = defkey(i)
          end if
          if (keys(i).eq.'file')  then
             file = values(i)
          else
             messg = keys(i)(1:k)//' " will be ignored'
             call warn(1,' *** restore: unknown keyword " '//messg)
          end if
 100   continue
c
c note that the sync_level is stored here (otherwise we'd "restore" it!)
       i_sync = int(getsca('&sync_level',0))
       if (i_sync.ge.3)  call iff_sync

       call undels(file)
       iex  = 0
       ier  = 0
       lun  = -1
       call openfl(lun, file, stat, iex, ier)
       if ((ier.lt.0).or.(iex.lt.0).or.(lun.le.0)) then
          ilen  = istrln(file)
          messg = ' *** restore: cannot open file '//file(1:ilen)
          call warn(2,messg)
          if (lun.gt.0) close(lun)
          return
       end if
c top line
       ilen = iread_ky(lun, pre, str)
       vers  = 'null'
       npack = 0
       if ((pre.eq.'#:').and.(str(1:20).eq.'IFEFFIT SAVE File: v')) then
          vers = str(21:25)
          call triml(vers)
          call lower(vers)
          call str2dp(vers,xvers,ier)
          s2   = str(26:)
          call lower(s2)
          nwords = 2
          call bwords(s2,nwords,words)
          if (words(1).eq.'npad') call str2in(words(2),npack,ierr)
       endif
       if ((npack.le.0).or.(npack.ge.16).or.(vers.eq.'null')) goto 1010
c next line must be 'build'
       if (xvers.ge.1.01) then
          ilen = iread_ky(lun, pre, str)
          if ((pre.ne.'#:').or.(str(1:5).ne.'build')) goto 1010
       endif
c next line must be 'strings'
       ilen = iread_ky(lun, pre, str)
       if ((pre.ne.'#:').or.(str(1:3).ne.'str')) goto 1010
c read strings
cc       print*,' strings:'
 150   continue
         ilen = iread_ky(lun, pre, str)
         if ((pre.eq.'#:').and.(str(1:3).eq.'sys')) goto 200
         if (pre.ne.'#=') goto 1010
         nwords   = 2
         words(1) = ' '
         words(2) = ' '
         call bwords(str,nwords,words)
         i         = iofstr(words(1),1)
         txtnam(i) = words(1)
         call strclp(str, words(1), words(2), text(i))
       go to 150
c scalars
 200   continue

 250   continue
       nwords = 2
       words(1) = ' '
       words(2) = ' '
       ilen = iread_ky(lun, pre, str)
cc       print*, ' pre, str = ', pre, ',  ', str(1:30)
       if ((pre.eq.'#:').and.(str(1:3).eq.'sca')) goto 300
       call bwords(str, nwords,words)
       i         = iofsca( words(1), 1)
cc       print*, ' iofsca = ', i
       scanam(i) = words(1)
       scafrm(i) = ' '
       if (pre.eq.'#?') then
          scalar(i)   = unpad(words(2), npack)
          icdsca(1,i) = -1
       elseif (pre.eq.'#%') then
          scalar(i)   = unpad(words(2), npack)
          icdsca(1,i) = i + jscale
       elseif (pre.eq.'#=') then
          call iff_set('def',str,.false.)
       else
          goto 1010
       endif
       goto 250
 300   continue
       sc_pad = .true.
       if ( (xvers.ge.1.015)  .and.
     $      (index(str(4:),'(nopad)').ne.0) ) then
          sc_pad = .false.
       endif
 350   continue
       nwords = 2
       words(1) = ' '
       words(2) = ' '
       ilen = iread_ky(lun, pre, str)
       if ((pre.eq.'#:').and.(str(1:3).eq.'arr')) goto 400
       call bwords(str, nwords,words)
       i         = iofsca( words(1), 1)
cc       print*, ' iofsca = ', i
       scanam(i) = words(1)
       scafrm(i) = ' '
       if (pre.eq.'#?') then
          if (sc_pad) then
             scalar(i)   = unpad(words(2), npack)
          else
             call str2dp(words(2), scalar(i), ier)
          endif
          icdsca(1,i) = -1
       elseif (pre.eq.'#%') then
          if (sc_pad) then
             scalar(i)   = unpad(words(2), npack)
          else
             call str2dp(words(2), scalar(i), ier)
          endif
          icdsca(1,i) = i + jscale
       elseif (pre.eq.'#=') then
          call iff_set('def',str,.false.)
       else
          goto 1010
       endif
       goto 350
 400   continue
c arrays
cc       print*,' arrays:'
 650   continue
       ilen = iread_ky(lun, pre,str)
cc       print*,' pre, ilen ', pre, ' ', ilen
       if (ilen.lt.0) goto 800
       in1 = index(str,'[')
       in2 = index(str,']')
       call str2in(str(in1+1:in2-1),ntmp, ierr)
       str = str(in2+1:)
       call triml(str)
       if (pre.eq.'#=') then
c define array by formula
          call iff_set('def',str,.false.)
       elseif (pre.eq.'#%') then
c
c read array data
          j  = index(str,'.')
          s1 = str(1:j-1)
          s2 = str(j+1:)
          call triml(s1)
          call triml(s2)
          call rdpadd(lun,npack,tmparr, ntmp)
          call set_array(s2,s1,tmparr,ntmp,1)
       else
          goto 1010
       endif
       go to 650

 800   continue
cc       print*, ' line 800'
c
c arrays
       if (lun.gt.0) close(lun)
       if (i_sync.ge.1)  call iff_sync
       return
 1010  continue
       call warn(2,' *** restore: not a valid save file')
       return
       end

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
!##############################################################################
!# ****************************************************************************
!# <name> flagship_signals </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains routines for signal handling
!#
!# The following routines are available:
!#
!# 1.) flagship_SIGINT
!#     -> Performs signal handling for SIGINT
!#
!# 2.) flagship_SIGQUT
!#     -> Performs signal handling for SIGQUIT
!# </purpose>
!##############################################################################

module flagship_signals

#include "flagship.h"

!$ use omp_lib
  use fsystem
  use genoutput
  use signals

  implicit none

  private
  public :: flagship_SIGINT
  public :: flagship_SIGQUIT

contains

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

!<function>

  function flagship_SIGINT(signum) result(sigcount)

!<description>
    ! This subroutine performs signal handling for SIGINT. In essence,
    ! it counts the number if SIGINT`s received and terminates if user
    ! sent SIGINT more than three times.
!</description>

!<input>
    integer, intent(in) :: signum
!</input>

!<result>
    ! signal
    integer :: sigcount
!</result>
!</function>

    ! local variables
    integer, save :: icount = 0

    sigcount = icount

    if (signum .eq. -1) then

      ! Reset counter
      icount = 0

    elseif (signum .eq. SIGINT) then

      ! Increase counter
      icount = icount+1
      if (icount .ge. 3) then
        call output_line('Simulation terminated due to user interruption (SIGINT)')
        call sys_halt()
      end if

    end if
  end function flagship_SIGINT

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

!<function>

  function flagship_SIGQUIT(signum) result(sigcount)

!<description>
    ! This subroutine performs signal handling for SIGQUIT.
!</description>

!<input>
    integer, intent(in) :: signum
!</input>

!<result>
    ! signal
    integer :: sigcount
!</result>
!</function>

    sigcount = 0
    call output_line('Simulation terminated due to user interruption (SIGQUIT)')
    call sys_halt()
  end function flagship_SIGQUIT

end module flagship_signals

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                     C
C                Copyright (C)  1995                                  C
C        University Corporation for Atmospheric Research              C
C                All Rights Reserved                                  C
C                                                                     C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C      File:           tx07f.f
C
C      Author:         Fred Clare (converted to Fortran by David
C                      Younghans)
C                      National Center for Atmospheric Research
C                      PO 3000, Boulder, Colorado
C                      (translated to Fortran by David Younghans)
C
C      Date:           Tue Nov  27 16:00:58 MST 1995
C
C      Description:    Demonstrates TextItem text spacings and
C                      aspect ratios.
C
      program tx07f
      implicit none
      external NhlFAppClass
      external NhlFNcgmWorkstationClass
      external NhlFCairoWindowWorkstationClass
      external NhlFPSWorkstationClass
      external NhlFPDFWorkstationClass
      external NhlFCairoPSPDFWorkstationClass
      external NhlFCairoImageWorkstationClass
      external NhlFTextItemClass

      integer   i,ierr,appid,wid,exid,labid,rlist
      real ypos,aspect
      real bkg_color(3),spacings(3)
      character*25  label
      character*7  wks_type

      data bkg_color/1.0,1.0,1.0/
      data spacings/0.0,1.5,0.6/
C
C Set the display.  Default is to display output to an X workstation.
C
      wks_type = "x11"
C
C  Initialize.
C
      call NhlFInitialize
      call NhlFRLCreate(rlist,'SETRL')

      call NhlFRLClear(rlist)
      call NhlFRLSetString(rlist,'appUsrDir','./',ierr)
      call NhlFRLSetString(rlist,'appDefaultParent','True',ierr)
      call NhlFCreate(appid,'tx07',NhlFAppClass,0,rlist,ierr)

      if (wks_type.eq."ncgm".or.wks_type.eq."NCGM") then
C
C Create a meta file workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkMetaName',
     &        './tx07f.ncgm',ierr)
         call NhlFRLSetFloatArray(rlist,'wkBackgroundColor',
     &        bkg_color,3,ierr)
         call NhlFCreate(wid,'tx07Work',NhlFncgmWorkstationClass,
     &        0,rlist,ierr)

      else if (wks_type.eq."x11".or.wks_type.eq."X11") then
C
C Create an X11 workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetInteger(rlist,'wkPause','True',ierr)
         call NhlFRLSetFloatArray(rlist,'wkBackgroundColor',
     &        bkg_color,3,ierr)
         call NhlFCreate(wid,'tx07Work',
     &        NhlFCairoWindowWorkstationClass,0, rlist,ierr)

      else if (wks_type.eq."oldps".or.wks_type.eq."OLDPS") then
C
C Create an older-style PS workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkPSFileName',
     &        './tx07f.ps',ierr)
         call NhlFRLSetFloatArray(rlist,'wkBackgroundColor',
     &        bkg_color,3,ierr)
         call NhlFCreate(wid,'tx07Work',NhlFpsWorkstationClass,
     &        0,rlist,ierr)
      else if (wks_type.eq."oldpdf".or.wks_type.eq."OLDPDF") then

C
C Create an older-style PDF workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkPDFFileName',
     &        './tx07f.pdf',ierr)
         call NhlFRLSetFloatArray(rlist,'wkBackgroundColor',
     &        bkg_color,3,ierr)
         call NhlFCreate(wid,'tx07Work',NhlFpdfWorkstationClass,
     &        0,rlist,ierr)
      else if (wks_type.eq."pdf".or.wks_type.eq."PDF".or.
     +         wks_type.eq."ps".or.wks_type.eq."PS") then

C
C Create a cairo PS/PDF workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkFileName',
     &        './tx07f',ierr)
         call NhlFRLSetString(rlist,'wkFormat',wks_type,ierr)
         call NhlFRLSetFloatArray(rlist,'wkBackgroundColor',
     &        bkg_color,3,ierr)
         call NhlFCreate(wid,'tx07Work',
     &        NhlFCairoPSpdfWorkstationClass,0,rlist,ierr)
      else if (wks_type.eq."png".or.wks_type.eq."PNG") then

C
C Create a cairo PNG workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkFileName',
     &        './tx07f',ierr)
         call NhlFRLSetString(rlist,'wkFormat',wks_type,ierr)
         call NhlFRLSetFloatArray(rlist,'wkBackgroundColor',
     &        bkg_color,3,ierr)
         call NhlFCreate(wid,'tx07Work',
     &        NhlFCairoImageWorkstationClass,0,rlist,ierr)
      endif
C
C Create two TextItem objects.
C
      call NhlFSetColor(wid,1,0.0, 0.0, 1.0,ierr)
      call NhlFSetColor(wid,2,0.4, 0.0, 0.4,ierr)
      call NhlFSetColor(wid,3,1.0, 0.0, 0.0,ierr)

      call NhlFCreate(exid,'Example String',
     &     NhlFtextItemClass,wid,0,ierr)
      call NhlFCreate(labid,'Label String',NhlFtextItemClass,wid,0,ierr)
C
C  Set up example string.
C
      call NhlFRLClear(rlist)
      call NhlFRLSetString(rlist,'txString', 'NCAR Graphics',ierr)
      call NhlFRLSetInteger(rlist,'txFontColor', 1,ierr)
      call NhlFRLSetInteger(rlist,'txFont', 25,ierr)
      call NhlFRLSetFloat(rlist,'txPosXF',  0.5,ierr)
      call NhlFRLSetFloat(rlist,'txFontHeightF', .05,ierr)
      call NhlFSetValues(exid,rlist,ierr)
C
C  Set up label string.
C
      call NhlFRLClear(rlist)
      call NhlFRLSetInteger(rlist,'txFontColor', 2,ierr)
      call NhlFRLSetInteger(rlist,'txFont', 21,ierr)
      call NhlFRLSetFloat(rlist,'txPosXF',  0.5,ierr)
      call NhlFRLSetFloat(rlist,'txFontHeightF', .04,ierr)
      call NhlFSetValues(labid,rlist,ierr)
C
C  Spacings
C
      do 100 i=0,2
         ypos = 0.83 - 0.17 * real(i)
         call NhlFRLClear(rlist)

         write (label,900) spacings(i+1)
         call NhlFRLSetString(rlist,'txString',label,ierr)
         call NhlFRLSetFloat(rlist,'txPosYF',ypos,ierr)
         call NhlFSetValues(labid,rlist,ierr)
         call NhlFDraw(labid,ierr)

         call NhlFRLClear(rlist)
         call NhlFRLSetFloat(rlist,'txConstantSpacingF',
     &        spacings(i+1),ierr)
         call NhlFRLSetFloat(rlist,'txPosYF',ypos - 0.07,ierr)
         call NhlFSetValues(exid,rlist,ierr)
         call NhlFDraw(exid,ierr)
 100  continue
 900  format ('txConstantSpacingF = ',f3.1)
C
C  Aspect ratios.
C
      do 200 i=0,1
         ypos = 0.32 - 0.17 * real(i)
         aspect = 1.3 * real(i) + 0.7
         call NhlFRLClear(rlist)

         write (label,910) aspect
         call NhlFRLSetString(rlist,'txString',label,ierr)
         call NhlFRLSetFloat(rlist,'txPosYF',ypos,ierr)
         call NhlFSetValues(labid,rlist,ierr)
         call NhlFDraw(labid,ierr)

         call NhlFRLClear(rlist)
         call NhlFRLSetFloat(rlist,'txConstantSpacingF',0.0,ierr)
         call NhlFRLSetFloat(rlist,'txPosYF',ypos - 0.07,ierr)
         call NhlFRLSetFloat(rlist,'txFontAspectF',aspect,ierr)
         call NhlFSetValues(exid,rlist,ierr)
         call NhlFDraw(exid,ierr)
 200  continue
 910  format ('txFontAspectF = ',f3.1)
C
C  Plot title.
C
      call NhlFRLClear(rlist)
      call NhlFRLSetInteger(rlist,'txFontColor', 3,ierr)
      call NhlFRLSetInteger(rlist,'txFont', 25,ierr)
      call NhlFRLSetFloat(rlist,'txPosYF',  0.93,ierr)
      call NhlFRLSetFloat(rlist,'txFontHeightF', .045,ierr)
      call NhlFRLSetString(rlist,'txString',
     &     'Text Spacings and Aspect Ratios',ierr)
      call NhlFSetValues(labid,rlist,ierr)
      call NhlFDraw(labid,ierr)

      call NhlFFrame(wid,ierr)
C
C  Close things down.
C
      call NhlFDestroy(exid,ierr)
      call NhlFDestroy(labid,ierr)
      call NhlFDestroy(wid,ierr)
      call NhlFClose

      stop
      end

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


!<ERROR> *fman* standard demo code format not found for ibset


text:
      SUBROUTINE DDRVRF1( NOUT, NN, NVAL, THRESH, A, LDA, ARF, WORK )
*
*  -- LAPACK test routine (version 3.2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2008
*
*     .. Scalar Arguments ..
      INTEGER            LDA, NN, NOUT
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      INTEGER            NVAL( NN )
      DOUBLE PRECISION   A( LDA, * ), ARF( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  DDRVRF1 tests the LAPACK RFP routines:
*      DLANSF
*
*  Arguments
*  =========
*
*  NOUT          (input) INTEGER
*                The unit number for output.
*
*  NN            (input) INTEGER
*                The number of values of N contained in the vector NVAL.
*
*  NVAL          (input) INTEGER array, dimension (NN)
*                The values of the matrix dimension N.
*
*  THRESH        (input) DOUBLE PRECISION
*                The threshold value for the test ratios.  A result is
*                included in the output file if RESULT >= THRESH.  To have
*                every test ratio printed, use THRESH = 0.
*
*  A             (workspace) DOUBLE PRECISION array, dimension (LDA,NMAX)
*
*  LDA           (input) INTEGER
*                The leading dimension of the array A.  LDA >= max(1,NMAX).
*
*  ARF           (workspace) DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2).
*
*  WORK          (workspace) DOUBLE PRECISION array, dimension ( NMAX )
*
*  =====================================================================
*     ..
*     .. Parameters ..
      DOUBLE PRECISION   ONE
      PARAMETER          ( ONE = 1.0D+0 )
      INTEGER            NTESTS
      PARAMETER          ( NTESTS = 1 )
*     ..
*     .. Local Scalars ..
      CHARACTER          UPLO, CFORM, NORM
      INTEGER            I, IFORM, IIN, IIT, INFO, INORM, IUPLO, J, N,
     +                   NERRS, NFAIL, NRUN
      DOUBLE PRECISION   EPS, LARGE, NORMA, NORMARF, SMALL
*     ..
*     .. Local Arrays ..
      CHARACTER          UPLOS( 2 ), FORMS( 2 ), NORMS( 4 )
      INTEGER            ISEED( 4 ), ISEEDY( 4 )
      DOUBLE PRECISION   RESULT( NTESTS )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLANSY, DLANSF, DLARND
      EXTERNAL           DLAMCH, DLANSY, DLANSF, DLARND
*     ..
*     .. External Subroutines ..
      EXTERNAL           DTRTTF
*     ..
*     .. Scalars in Common ..
      CHARACTER*32       SRNAMT
*     ..
*     .. Common blocks ..
      COMMON             / SRNAMC / SRNAMT
*     ..
*     .. Data statements ..
      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
      DATA               UPLOS / 'U', 'L' /
      DATA               FORMS / 'N', 'T' /
      DATA               NORMS / 'M', '1', 'I', 'F' /
*     ..
*     .. Executable Statements ..
*
*     Initialize constants and the random number seed.
*
      NRUN = 0
      NFAIL = 0
      NERRS = 0
      INFO = 0
      DO 10 I = 1, 4
         ISEED( I ) = ISEEDY( I )
   10 CONTINUE
*
      EPS = DLAMCH( 'Precision' )
      SMALL = DLAMCH( 'Safe minimum' )
      LARGE = ONE / SMALL
      SMALL = SMALL * LDA * LDA
      LARGE = LARGE / LDA / LDA
*
      DO 130 IIN = 1, NN
*
         N = NVAL( IIN )
*
         DO 120 IIT = 1, 3
*
*           IIT = 1 : random matrix
*           IIT = 2 : random matrix scaled near underflow
*           IIT = 3 : random matrix scaled near overflow
*
            DO J = 1, N
               DO I = 1, N
                  A( I, J) = DLARND( 2, ISEED )
               END DO
            END DO
*
            IF ( IIT.EQ.2 ) THEN
               DO J = 1, N
                  DO I = 1, N
                     A( I, J) = A( I, J ) * LARGE
                  END DO
               END DO
            END IF
*
            IF ( IIT.EQ.3 ) THEN
               DO J = 1, N
                  DO I = 1, N
                     A( I, J) = A( I, J) * SMALL
                  END DO
               END DO
            END IF
*
*           Do first for UPLO = 'U', then for UPLO = 'L'
*
            DO 110 IUPLO = 1, 2
*
               UPLO = UPLOS( IUPLO )
*
*              Do first for CFORM = 'N', then for CFORM = 'C'
*
               DO 100 IFORM = 1, 2
*
                  CFORM = FORMS( IFORM )
*
                  SRNAMT = 'DTRTTF'
                  CALL DTRTTF( CFORM, UPLO, N, A, LDA, ARF, INFO )
*
*                 Check error code from DTRTTF
*
                  IF( INFO.NE.0 ) THEN
                     IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) THEN
                        WRITE( NOUT, * )
                        WRITE( NOUT, FMT = 9999 )
                     END IF
                     WRITE( NOUT, FMT = 9998 ) SRNAMT, UPLO, CFORM, N
                     NERRS = NERRS + 1
                     GO TO 100
                  END IF
*
                  DO 90 INORM = 1, 4
*
*                    Check all four norms: 'M', '1', 'I', 'F'
*
                     NORM = NORMS( INORM )
                     NORMARF = DLANSF( NORM, CFORM, UPLO, N, ARF, WORK )
                     NORMA = DLANSY( NORM, UPLO, N, A, LDA, WORK )
*
                     RESULT(1) = ( NORMA - NORMARF ) / NORMA / EPS
                     NRUN = NRUN + 1
*
                     IF( RESULT(1).GE.THRESH ) THEN
                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) THEN
                           WRITE( NOUT, * )
                           WRITE( NOUT, FMT = 9999 )
                        END IF
                        WRITE( NOUT, FMT = 9997 ) 'DLANSF',
     +                      N, IIT, UPLO, CFORM, NORM, RESULT(1)
                        NFAIL = NFAIL + 1
                     END IF
   90             CONTINUE
  100          CONTINUE
  110       CONTINUE
  120    CONTINUE
  130 CONTINUE
*
*     Print a summary of the results.
*
      IF ( NFAIL.EQ.0 ) THEN
         WRITE( NOUT, FMT = 9996 ) 'DLANSF', NRUN
      ELSE
         WRITE( NOUT, FMT = 9995 ) 'DLANSF', NFAIL, NRUN
      END IF
      IF ( NERRS.NE.0 ) THEN
         WRITE( NOUT, FMT = 9994 ) NERRS, 'DLANSF'
      END IF
*
 9999 FORMAT( 1X, ' *** Error(s) or Failure(s) while testing DLANSF
     +         ***')
 9998 FORMAT( 1X, '     Error in ',A6,' with UPLO=''',A1,''', FORM=''',
     +        A1,''', N=',I5)
 9997 FORMAT( 1X, '     Failure in ',A6,' N=',I5,' TYPE=',I5,' UPLO=''',
     +        A1, ''', FORM =''',A1,''', NORM=''',A1,''', test=',G12.5)
 9996 FORMAT( 1X, 'All tests for ',A6,' auxiliary routine passed the ',
     +        'threshold (',I5,' tests run)')
 9995 FORMAT( 1X, A6, ' auxiliary routine:',I5,' out of ',I5,
     +        ' tests failed to pass the threshold')
 9994 FORMAT( 26X, I5,' error message recorded (',A6,')')
*
      RETURN
*
*     End of DDRVRF1
*
      END


text:
*------------------   Begin FS2Write  ----------------------------------
      subroutine FS2Write(f_fs2, nc_f_fs2, fas,
     :                    char_head, int_head, real_head, comments,
     :                    f_smc, nc_f_smc,
     :                    nacc_start, nacc_stop, npwr2,
     :                    nfreq, flow, fhigh, df, df_smooth, ftaper,
     :                    string,nc_string,datetime_l,
     :                    tskip, tlength_in, df_intrp, f_intrp_low,
     :                    f_intrp_high, signnpw2, dc_remove, nzpad)

*  Writes Fourier amplitude spectra into SMC format

* If datetime_l = .true., write date/time into the comments
* If nc_string .ne. 0, write "string" into the comments

* The arrays char_head, int_head, real_head must be dimensioned
* by the calling program as follows:
*
*      character char_head(11)*80, comments(N)*80 ! see "N" below
*      integer int_head(48)
*      real real_head(50)

*
* default values for int_head = -32768
*                    real_head = 1.7e+38

* Dates: 06/15/01 - written by D. Boore, patterned after smcwrite
*        07/10/01 - Added more processing parameters to comments in output file
*        11/11/01 - Write flow, df, fhigh, nfreq, npw2, dt to comments
*        12/18/01 - Store 0.0 if value is less than what can be written in the
*                   e format.
*        06/25/09 - Use get_date, get_time rather than datetime
!        12/29/10 - Include more precision in the output freq and fas.
!        08/01/12 - I reset int_head(47) to 8 (high precision) AFTER
!                   writing out int_head!  I corrected this.

      real real_head(*)
      integer int_head(*)
      character*80 char_head(*), comments(*)
      character f_fs2*(*), f_smc*(*), string*(*)
      real fas(*)
      logical datetime_l, dc_remove
      character date_c*10, time_c*11

* Open output file:
      call get_lun(nu_out)
      open(unit=nu_out, file=f_fs2(1:nc_f_fs2), status='unknown')

* Write headers:

* First write the 11 lines of text headers:

      char_head(1) = ' '
      char_head(1) = '6 FOURIER AMPLITUDE SPECTRUM'

      do i = 1, 11
        write(nu_out,'(a)') char_head(i)
      end do

* Fill in the int and real headers:

      int_head(42) = nacc_start
      int_head(43) = nacc_stop
      int_head(44) = npwr2
      int_head(45) = nfreq

      real_head(42) = flow
      real_head(43) = fhigh
      real_head(44) = df
      real_head(45) = df_smooth
      real_head(46) = ftaper

* Write to comments:

      ncomments = int_head(16)

      ncomments = ncomments + 1
      comments(ncomments) = '|'
      ncomments = ncomments + 1
      comments(ncomments) =
     :  '| Fourier amplitude spectrum computed for file:'
      ncomments = ncomments + 1
      comments(ncomments) = '|   '//f_smc(1:nc_f_smc)

      if(nc_string .ne. 0) then
        ncomments = ncomments + 1
        comments(ncomments) = string(1:nc_string)
      end if

      if(datetime_l) then
        iflag = 0
        call get_date(date_c)
        call get_time(time_c)
        ncomments = ncomments + 1
        comments(ncomments) =
     :  '|   file written on '//date_c//' at time '//time_c
      end if

      ncomments = ncomments + 1
      comments(ncomments) =
     :  '|   tskip, tlength_in, ftaper = '
      write(comments(ncomments)(42:80),
     :  '(1x,f7.2,1x,f8.2,1x,f4.1)')
     :    tskip, tlength_in, ftaper
      ncomments = ncomments + 1
      comments(ncomments) =
     :  '|   df_intrp(0=FFT freqs), f_intrp_low, f_intrp_high = '
      write(comments(ncomments)(60:80),
     :  '(1x,f3.1,1x,f7.3,1x,f7.3)')
     :       df_intrp, f_intrp_low, f_intrp_high

      ncomments = ncomments + 1
      comments(ncomments) =
!        1234567890123456789012345678901234567890
     :  '|   df_smooth(Hz;0=none), '//
     :  'signnpw2(<0= backup for npw2, no zpad) = '
      write(comments(ncomments)(66:80),
     :  '(1x,f5.2,1x,f4.1)')
     :    df_smooth, signnpw2

      ncomments = ncomments + 1
      comments(ncomments) =
!        1234567890123456789012345678901234567890
     :  '|   npw2, nzpad, dt = '
      write(comments(ncomments)(23:80),
     :  '(1x,i6,1x,i5, 1x,f8.5)')
     :       int_head(44), nzpad, 1.0/real_head(2)

      ncomments = ncomments + 1
      comments(ncomments) =
     :  '|   dc_remove? = '
      write(comments(ncomments)(20:80),
     :  '(1x,l1)')
     :    dc_remove

      ncomments = ncomments + 1
      comments(ncomments) =
     :  '|   flow (rhead(42)), df(rhead(44)), fhigh(rhead(43)), '//
     :  'nfreq(ihead(45)) = '
      ncomments = ncomments + 1
      comments(ncomments)(1:1) = '|'
      write(comments(ncomments)(2:80),
     :  '(6x,f8.4,9x,f8.6,7x,f8.4,13x,i5)')
     :    real_head(42), real_head(44), real_head(43), int_head(45)

      ncomments = ncomments + 1
      comments(ncomments) = '|'

      int_head(16) = ncomments

      int_head(47) = 8

      write(nu_out, '(8I10)') (int_head(i), i=1, 48)

      write(nu_out, '(5e15.7)') (real_head(i), i = 1, 50)

      if (int_head(16) .gt. 0) then
        do i = 1, int_head(16)
          write(nu_out,'(a)') comments(i)
        end do
      end if

* Check for a small number.  There is probably a better way of doing this
* involving writing one line at a time rather than replacing a value of FAS
      do i = 1, int_head(45)
        if (abs(fas(i)) .le. 1.0e-09) then
          write(*,'(a)') '  replacing fas with 0.0; i, fas = '
          write(*,*) i, fas(i)
          fas(i) = 0.0
        end if
      end do

      write(nu_out, '(5(1pe14.7e1))') (fas(i), i = 1, int_head(45))

* That should be it; close the output file.

      close(unit=nu_out)

      return
      end
*------------------   End FS2Write   ----------------------------------

ext: for

max_stars_repo_name: ZhangHCFJEA/bbp

max_issues_repo_name: ZhangHCFJEA/bbp

max_forks_repo_name: ZhangHCFJEA/bbp