   module fire_types

   implicit none

   integer, parameter :: namsize = 256
   integer, parameter :: linsize = 512

   TYPE mdl_poly_type
     integer :: nVtx
     real    :: cntr_lon
     real    :: cntr_lat
     real    :: area
     real, allocatable :: vtx_lon(:)
     real, allocatable :: vtx_lat(:)
     logical :: active
     logical :: x180
   END TYPE mdl_poly_type

   TYPE model_grid_type
     integer :: nlons
     integer :: nlats
     integer :: nlevs
     integer :: maxPolyVtx
     integer :: nPolygons
     integer :: lon_shift = 0
     integer :: PolarCapNdx(2) = 0
     real    :: dx
     real, allocatable :: mdl_lons(:,:)
     real, allocatable :: mdl_lats(:,:)
     real, allocatable :: mdl_cell_area(:,:)
     real, allocatable :: z_at_w(:,:,:)                    ! height above ground at interfaces (meters)
     real, allocatable :: landmask(:,:)                    ! land/water mask (1/0)
     character(len=namsize) :: filename                    ! name of input file for CAM
     character(len=5 ) :: mdlType                          ! 'WRF' || 'CAMFV' || 'CAMSE' || 'MPAS'
     logical :: mdl_is_WRF
     logical :: mdl_is_CAM, mdl_is_CAMFV, mdl_is_CAMSE
     logical :: mdl_is_MPAS = .false.
     TYPE(mdl_poly_type), allocatable :: mdl_poly(:)
   END TYPE model_grid_type

   end module fire_types


text:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Extract camera serial number from Voyager flight label.

      subroutine vgrlab_cam(lab,cam)
      implicit none
      character*7200 lab	!input Voyager flight label
      integer*4 cam		!output camera serial number (4-7)

      integer*4 i
      integer*4 sc_id,camera_id
      integer*4 camera_sn(2,2)/7,6,	!VGR1 NA,WA
     &			       5,4/	!VGR2 NA,WA

      call vgrlab_sc_id(lab,sc_id)
      call vgrlab_camera_id(lab,camera_id)

      if (sc_id.eq.-999 .or. camera_id.eq.-999) then
         cam = -999
      else
         cam = camera_sn(camera_id,sc_id)
      endif
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Extract spacecraft ID from Voyager flight label.

      subroutine vgrlab_sc_id(lab,sc_id)
      implicit none
      character*7200 lab	!input Voyager flight label
      integer*4 sc_id		! 1=VGR-1, 2=VGR-2

      integer*4 i

      sc_id = -999

      i = index(lab(1:),'VGR-1')
      if (i.ne.0) then
         sc_id = 1		!Voyager 1
      else
         i = index(lab(1:),'VGR-2')
         if (i.ne.0) sc_id=2
      endif
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Extract camera id from Voyager flight label.

      subroutine vgrlab_camera_id(lab,camera_id)
      implicit none
      character*7200 lab	!input Voyager flight label
      integer*4 camera_id	! 1=NA, 2=WA

      integer*4 i

      camera_id = -999

      i = index(lab(1:),'NA CAMERA')
      if (i.ne.0) then
         camera_id = 1		!narrow angle
      else
         i = index(lab(1:),'WA CAMERA')
         if (i.ne.0) camera_id=2
      endif
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Extract FDS count from Voyager flight label.
c
      subroutine vgrlab_fds(lab,fds)
      implicit none
      character*7200 lab	!Voyager flight label (input)
      integer*4 fds		!FDS count (output)

      integer*4 i,mod16,mod60

c note: the FDS is returned as an integer as 100*mod16 + mod60
c Example: FDS 16368.32 is returened as 1636832

       i = index(lab(1:),' FDS ')
       if (i.ne.0 .and. lab(i+5:i+5).ne.'*') then
          read(lab(i+5:),'(bn,i5)') mod16
          read(lab(i+11:),'(bn,i2)') mod60
          fds = 100*mod16 + mod60
       else
          fds = -999
       endif
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Extract filter, gain, scan rate, operating mode, and exposure time
c from Voyager flight label.

      subroutine vgrlab_camparams(lab,filter,gain,scan,mode,t)
      implicit none
c     ...input argument
      character*7200 lab	!Voyager flight label
c     ...output arguments
      integer*4 filter		!filter position (0-7)
      integer*4 gain		!0=low gain, 1=high gain
      integer*4 scan		!scan rate: 1, 2, 3, 5, 10
      integer*4 mode		!operating mode
      real*4 t			!exposure time in milliseconds

      integer*4 i

      filter = -999
      gain = -999
      scan = -999
      mode = -999
      t = -999.

      i = index(lab(1:),' FILT ')	!example: FILT 0(CLEAR)
      if (i.ne.0 .and. lab(i+6:i+6).ne.'*')
     &	    read(lab(i+6:),'(bn,i1)') filter

      if (index(lab(1:),' LO GAIN').gt.0) then
         gain = 0
      elseif (index(lab(1:),' HI GAIN ').gt.0) then
         gain = 1
      endif

      i = index(lab(1:),'  SCAN RATE')	!example: SCAN RATE 10:1
      if (i.ne.0 .and. lab(i+12:i+12) .ne. '*')
     & 	   read(lab(i+12:),'(bn,i2)') scan

      i = index(lab(1:),' MODE')	!example: MODE 4(BOTALT)
      if (i.ne.0 .and. lab(i+6:i+6).ne.'*')
     &	   read (lab(i+6:),'(bn,i1)') mode

      i = index(lab(1:),' EXP ')	!example: EXP 00360.0 MSEC
      if (i.ne.0 .and. lab(i+5:i+5).ne.'*')
     &	   read(lab(i+5:),'(bn,f7.0)') t
      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Extract PICNO from Voyager flight label.
c
      subroutine vgrlab_picno(lab,picno,ind1)
      implicit none
      character*7200 lab	!Voyager flight label (input)
      character*10 picno	!picno='xxxxPS+DDD' (output)
      integer*4 ind1

      integer*4 i

      i = index(lab(1:),' PICNO ')	!Example: PICNO 0548J1-001
      if (i.gt.0 .and. lab(i+7:i+7).ne.'X'
     +          .and. lab(i+7:i+7).ne.'*') then
         picno = lab(i+7:i+16)
      else
         ind1 = -1
      endif
      return
      end

ext: f

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
! RUN: bbc -emit-fir -use-alloc-runtime %s -o - | FileCheck %s

! Test lowering of allocatables using runtime for allocate/deallocate statements.
! CHECK-LABEL: _QPpointer_runtime(
subroutine pointer_runtime(n)
    integer :: n
    character(:), pointer :: scalar, array(:)
    ! CHECK-DAG: %[[sBoxAddr:.*]] = fir.alloca !fir.box<!fir.ptr<!fir.char<1,?>>> {{{.*}}uniq_name = "_QFpointer_runtimeEscalar"}
    ! CHECK-DAG: %[[sNullAddr:.*]] = fir.zero_bits !fir.ptr<!fir.char<1,?>>
    ! CHECK-DAG: %[[sInitBox:.*]] = fir.embox %[[sNullAddr]] typeparams %c0{{.*}} : (!fir.ptr<!fir.char<1,?>>, index) -> !fir.box<!fir.ptr<!fir.char<1,?>>>
    ! CHECK-DAG: fir.store %[[sInitBox]] to %[[sBoxAddr]] : !fir.ref<!fir.box<!fir.ptr<!fir.char<1,?>>>>

    ! CHECK-DAG: %[[aBoxAddr:.*]] = fir.alloca !fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>> {{{.*}}uniq_name = "_QFpointer_runtimeEarray"}
    ! CHECK-DAG: %[[aNullAddr:.*]] = fir.zero_bits !fir.ptr<!fir.array<?x!fir.char<1,?>>>
    ! CHECK-DAG: %[[aNullShape:.*]] = fir.shape %c0{{.*}} : (index) -> !fir.shape<1>
    ! CHECK-DAG: %[[aInitBox:.*]] = fir.embox %[[aNullAddr]](%[[aNullShape]]) typeparams %c0{{.*}} : (!fir.ptr<!fir.array<?x!fir.char<1,?>>>, !fir.shape<1>, index) -> !fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>
    ! CHECK-DAG: fir.store %[[aInitBox]] to %[[aBoxAddr]] : !fir.ref<!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>>

    allocate(character(10):: scalar, array(30))
    ! CHECK-DAG: %[[sBoxCast1:.*]] = fir.convert %[[sBoxAddr]] : (!fir.ref<!fir.box<!fir.ptr<!fir.char<1,?>>>>) -> !fir.ref<!fir.box<none>>
    ! CHECK-DAG: %[[ten1:.*]] = fir.convert %c10{{.*}} : (i32) -> i64
    ! CHECK: fir.call @{{.*}}PointerNullifyCharacter(%[[sBoxCast1]], %[[ten1]], %c1{{.*}}, %c0{{.*}}, %c0{{.*}})
    ! CHECK-NOT: PointerSetBounds
    ! CHECK: %[[sBoxCast2:.*]] = fir.convert %[[sBoxAddr]] : (!fir.ref<!fir.box<!fir.ptr<!fir.char<1,?>>>>) -> !fir.ref<!fir.box<none>>
    ! CHECK: fir.call @{{.*}}PointerAllocate(%[[sBoxCast2]]

    ! CHECK-DAG: %[[aBoxCast1:.*]] = fir.convert %[[aBoxAddr]] : (!fir.ref<!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>>) -> !fir.ref<!fir.box<none>>
    ! CHECK-DAG: %[[ten2:.*]] = fir.convert %c10{{.*}} : (i32) -> i64
    ! CHECK: fir.call @{{.*}}PointerNullifyCharacter(%[[aBoxCast1]], %[[ten2]], %c1{{.*}}, %c1{{.*}}, %c0{{.*}})
    ! CHECK: %[[aBoxCast2:.*]] = fir.convert %[[aBoxAddr]] : (!fir.ref<!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>>) -> !fir.ref<!fir.box<none>>
    ! CHECK: fir.call @{{.*}}PointerSetBounds(%[[aBoxCast2]]
    ! CHECK: %[[aBoxCast3:.*]] = fir.convert %[[aBoxAddr]] : (!fir.ref<!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>>) -> !fir.ref<!fir.box<none>>
    ! CHECK: fir.call @{{.*}}PointerAllocate(%[[aBoxCast3]]

    deallocate(scalar, array)
    ! CHECK: %[[sBoxCast3:.*]] = fir.convert %[[sBoxAddr]] : (!fir.ref<!fir.box<!fir.ptr<!fir.char<1,?>>>>) -> !fir.ref<!fir.box<none>>
    ! CHECK: fir.call @{{.*}}PointerDeallocate(%[[sBoxCast3]]
    ! CHECK: %[[aBoxCast4:.*]] = fir.convert %[[aBoxAddr]] : (!fir.ref<!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>>) -> !fir.ref<!fir.box<none>>
    ! CHECK: fir.call @{{.*}}PointerDeallocate(%[[aBoxCast4]]

    ! only testing that the correct length is set in the descriptor.
    allocate(character(n):: scalar, array(40))
    ! CHECK: %[[n:.*]] = fir.load %arg0 : !fir.ref<i32>
    ! CHECK-DAG: %[[ncast1:.*]] = fir.convert %[[n]] : (i32) -> i64
    ! CHECK-DAG: %[[sBoxCast4:.*]] = fir.convert %[[sBoxAddr]] : (!fir.ref<!fir.box<!fir.ptr<!fir.char<1,?>>>>) -> !fir.ref<!fir.box<none>>
    ! CHECK: fir.call @{{.*}}PointerNullifyCharacter(%[[sBoxCast4]], %[[ncast1]], %c1{{.*}}, %c0{{.*}}, %c0{{.*}})
    ! CHECK-DAG: %[[ncast2:.*]] = fir.convert %[[n]] : (i32) -> i64
    ! CHECK-DAG: %[[aBoxCast5:.*]] = fir.convert %[[aBoxAddr]] : (!fir.ref<!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>>) -> !fir.ref<!fir.box<none>>
    ! CHECK: fir.call @{{.*}}PointerNullifyCharacter(%[[aBoxCast5]], %[[ncast2]], %c1{{.*}}, %c1{{.*}}, %c0{{.*}})
  end subroutine

max_stars_repo_name: colemancda/llvm-project

max_issues_repo_name: colemancda/llvm-project

max_forks_repo_name: colemancda/llvm-project


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine adspo4 ( pmsa   , fl     , ipoint , increm , noseg  ,
     &                    noflux , iexpnt , iknmrk , noq1   , noq2   ,
     &                    noq3   , noq4   )
!>\file
!>       P-ad/desorption to particulate inorganic matter. 3 options for sorption formulation.

!
!     Description of the module :
!     P-adsorption onto particulate inorganic matter. 3 options for
!     sorption formulation.
!
!        ----- old version -----
! Name    T   L I/O   Description                                   Units
! ----    --- -  -    -------------------                            ----
! AAP     R*4 1 I AAP concentration                                [gP/m3]
! SWAdsP  R*4 1 I swithc for adsorption option to use                  [-]
! DELT    R*4 1 I DELWAQ timestep                                    [scu]
! EQAAP   R*4 1 L calculated equlibrium AAP concentration          [gP/m3]
! FL (1)  R*4 1 O adsorption/desorption flux                     [gP/m3/d]
! TIM     R*4 1 I TIM   concentration                             [gDM/m3]
! KD      R*4 1 I partition coefficent PO4-AAP             [-] or [m3/gDM]
! MAXADS  R*4 1 I maximum adsorption (capacity)                   [gP/gDM]
! PMSA(9) R*4 1 O calculated equlibrium AAP concentration          [gP/m3]
! PO4     R*4 1 I PO4 concentration                                [gP/m3]
! RCADS   R*4 1 I first order rate constatnt adsorption              [1/d]
!
!        ----- new version -----
! Name    T   L I/O   Description                                   Units
! ----    --- -  -    -------------------                            ----
! AAP     R*4 1 I AAP concentration                                [gP/m3]
! PO4     R*4 1 I PO4 concentration                                [gP/m3]
! EQAAP   R*4 1 - calculated equlibrium AAP concentration          [gP/m3]
! KDADS   R*4 1 I distribution coeff. [-], or adsorpt. constant    [m3/gP]
! KADS20  R*4 1 I molar adsorption constant at 20 oC    [mole(a-1).l(a-1)]
! KADS    R*4 1 - molar adsorption constant             [mole(a-1).l(a-1)]
! KSORP   R*4 1 I first order rate constant for adsorption           [1/d]
! FCAP    R*4 1 I adsorption capacity for Fe                      [gP/gFe]
! FRA     R*4 1 - correction factor for oxidised iron                  [-]
! FADS    R*4 1 O adsorption/desorption flux                     [gP/m3/d]
! CADST   R*4 1 - molar total concentration adsorption sites    [moleFe/l]
! CADS    R*4 1 - molar concentration free adsorption sites     [moleFe/l]
! IM1     R*4 1 I IM1 concentration                               [gDM/m3]
! IM2     R*4 1 I IM2 concentration                               [gDM/m3]
! IM3     R*4 1 I IM3 concentration                               [gDM/m3]
! TFE     R*4 1 - total Fe concentration                          [gFe/m3]
! FRFE1   R*4 1 I fraction of Fe in IM1                          [gFe/gDW]
! FRFE2   R*4 1 I fraction of Fe in IM2                          [gFe/gDW]
! FRFE3   R*4 1 I fraction of Fe in IM3                          [gFe/gDW]
! FRFEOX  R*4 1 I fraction of oxidised Fe                              [-]
! FIM1    R*4 1 O fraction of adsorbed P in IM1                        [-]
! FIM2    R*4 1 O fraction of adsorbed P in IM2                        [-]
! FIM3    R*4 1 O fraction of adsorbed P in IM3                        [-]
! PH      R*4 1 I pH                                                   [-]
! OH      R*4 1 - hydroxyl concentration                          [mole/l]
! AOH     R*4 1 I reaction constant for hydroxyl                       [-]
! TEMP    R*4 1 I temperature                                         [oC]
! TC      R*4 1 I temperature coeffcient adsorption constant           [-]
! OXY     R*4 1 I dissolved oxygen concentration                  [gO2/m3]
! CROXY   R*4 1 I critical dissolved oxygen concentration         [gO2/m3]
! POROS   R*4 1 I prorosity                                            [-]
! DELT    R*4 1 I DELWAQ timestep                                    [scu]
! SWADSP  R*4 1 I switch for selection of the adsorption option        [-]
!
!     Logical Units : -

!     Modules called : -

!     Name     Type   Library
!     ------   -----  ------------
!
      IMPLICIT REAL (A-H,J-Z)
!
      REAL     PMSA  ( * ) , FL    (*)
      INTEGER  IPOINT( * ) , INCREM(*) , NOSEG , NOFLUX,
     +         IEXPNT(4,*) , IKNMRK(*) , NOQ1, NOQ2, NOQ3, NOQ4
!
      INTEGER  IVERSN
      REAL     AAP   , PO4   , EQAAP , KDADS , KADS20 , KADS   , KSORP ,
     +         FCAP  , FRA   , FADS  , CADST , CADS   , IM1    , IM2   ,
     +         IM3   , TFE   , FRFE1 , FRFE2 , FRFE3  , FRFEOX , FIM1  ,
     +         FIM2  , FIM3  , PH    , OH    , AOH    , TEMP   , TC    ,
     +         OXY   , CROXY , POROS , DELT  , SWADSP , QIM1   , QIM2  ,
     +         QIM3  , EQAAPM
      INTEGER  NR_MES, ILUMON
      SAVE     NR_MES
      DATA     NR_MES / 0 /
!

      IP1  = IPOINT( 1)
      IP2  = IPOINT( 2)
      IP3  = IPOINT( 3)
      IP4  = IPOINT( 4)
      IP5  = IPOINT( 5)
      IP6  = IPOINT( 6)
      IP7  = IPOINT( 7)
      IP8  = IPOINT( 8)
      IP9  = IPOINT( 9)
      IP10 = IPOINT(10)
      IP11 = IPOINT(11)
      IP12 = IPOINT(12)
      IP13 = IPOINT(13)
      IP14 = IPOINT(14)
      IP15 = IPOINT(15)
      IP16 = IPOINT(16)
      IP17 = IPOINT(17)
      IP18 = IPOINT(18)
      IP19 = IPOINT(19)
      IP20 = IPOINT(20)
      IP21 = IPOINT(21)
      IP22 = IPOINT(22)
      IP23 = IPOINT(23)
      IP24 = IPOINT(24)
      IP25 = IPOINT(25)
      IP26 = IPOINT(26)
      IP27 = IPOINT(27)
      IP28 = IPOINT(28)
      IP29 = IPOINT(29)
      IP30 = IPOINT(30)
      IP31 = IPOINT(31)
      IP32 = IPOINT(32)
      IP33 = IPOINT(33)
      IP34 = IPOINT(34)
!
      IFLUX = 0
!
      IVERSN = NINT( PMSA( IP23) )
!
!     Use the old version when IVERSN=0
!
      DO 9000 ISEG = 1 , NOSEG

!!    CALL DHKMRK(1,IKNMRK(ISEG),IKMRK1)
!!    IF (IKMRK1.EQ.1.OR.IKMRK1.EQ.3) THEN
      IF (BTEST(IKNMRK(ISEG),0)) THEN
!
      IF (IVERSN .EQ. 0) THEN
!
        dPAds = 0.0
        EQAAP = 0.0
        AtotP = 0.0
        Kads  = 0.0
!
        SWAdsP = PMSA(IP1 )
        PO4    = MAX( PMSA(IP2 ), 0.0)
        AAP    = MAX( PMSA(IP3 ), 0.0)
        IM1    = MAX( PMSA(IP4 ) , 0.0)
        IM2    = MAX( PMSA(IP5 ) , 0.0)
        IM3    = MAX( PMSA(IP6 ) , 0.0)
        RCADS  = PMSA(IP12)
        KD     = PMSA(IP7 )
        MAXADS = PMSA(IP8 )
        DELT   = PMSA(IP9 )
!
        TIM    = IM1 + IM2 + IM3
        QIM1 = 0.0
        QIM2 = 0.0
        QIM3 = 0.0
        IF ( TIM .LT. 1E-10) THEN
           FIM1 = 0.0
           FIM2 = 0.0
           FIM3 = 0.0
        ELSE
           FIM1 = IM1 / TIM
           FIM2 = IM2 / TIM
           FIM3 = IM3 / TIM
           IF ( IM1 .GT. 1E-10 ) QIM1 = AAP*FIM1/IM1
           IF ( IM2 .GT. 1E-10 ) QIM2 = AAP*FIM2/IM2
           IF ( IM3 .GT. 1E-10 ) QIM3 = AAP*FIM3/IM3
        ENDIF
!
!------ Error messages
!
        IF (KD.LT.0.0) CALL ERRSYS ('KD in ADSPO4 lower then zero', 1 )
        IF (((ABS(KD) .LT. 1E-20) .OR. (ABS(PO4).LT.1E-20))
     &      .AND. (ABS(MAXADS) .LT. 1E-20) )
     &     CALL ERRSYS
     &              ('(KD or PO4) and MAXADS equal zero in ADSPO4', 1 )
!
!------ (1) Instantanaeous equilibrium partitioning
!           SWAdsP = 0
!
        IF (NINT(SWAdsP) .EQ. 0) THEN
            dPAds =(((AAP + PO4) / (1.+ KD)) - AAP) /  DELT
            EQAAP = -1.0
        ENDIF
!
!------ (2) Kinetic Langmuir sorption
!           SWAdsP = 1
!
        IF (NINT(SWAdsP) .EQ. 1) THEN
          IF ( (MAXADS .LT. 1.E-10) .OR. (KD .LT. 1.E-10)) THEN
             EQAAP = 0.0
             dPAds = 0.0
          ELSE
             EQAAP = TIM * (KD * PO4 * MAXADS) / (KD * PO4 + MAXADS)
             dPAds = RCADS * (EQAAP - AAP)
          ENDIF
        ENDIF
!
!------ (3) GEM formulation for sorption
!           SWAdsP = 2
!
        IF (NINT(SWAdsP) .EQ. 2) THEN
           KAds20   = PMSA(IP10)
           TCKAds   = PMSA(IP11)
           RCadsP   = PMSA(IP25)
           aOHPO4   = PMSA(IP13)
           fr_Fe    = PMSA(IP24)
           frFeox   = PMSA(IP17)
           OXY      = PMSA(IP18)
           CrOXY    = PMSA(IP19)
           pH       = PMSA(IP20)
           Temp     = PMSA(IP21)
           poros    = PMSA(IP22)
!
!--------- (3a) Eqs 6.38, 6.37 of GEM report
           OH = 10.0**(pH-14.0)
           IF (OXY .GE. CrOXY*POROS) THEN
              fOxSor = 1.0
           ELSE
              fOxSor = frFeox
           ENDIF
!
!--------- (3b) Eqn 6.36 of GEM report : total sorption capacity (mol/l)
           AtotP = fOxSor * fr_Fe * TIM/(56000.0 * poros)
!
!--------- (3c) Eqn 6.35 of GEM report : free sorption capacity (mol/l)
           Afree = AtotP - (AAP/(31000.0 *poros))
!
!--------- (3d) Eqn 6.34 of GEM report : temp.correction of Kads
           Kads = Kads20 * TCKads**(Temp-20.0)
!
!--------- (3e) Eqn 6.33 of GEM report : equilibrium conc of Pads (gP/m3)
           EQAAP = (AAP + PO4) *
     &             (1.0 - 1.0/(Kads*Afree*(OH**aOHPO4) + 1.0))
!
!--------- (3f) Eqn 6.31 of GEM report : the adsorption flux (gP/m3/d)
           dPads = RCAdsP * (EQAAP - AAP)
!
        ENDIF
!
!---- Output of module
!
      FL(1+IFLUX) = dPAds
      PMSA(IP26)  = EQAAP
      PMSA(IP27)  = AtotP
      PMSA(IP28)  = Kads
      PMSA(IP29)  = FIM1
      PMSA(IP30)  = FIM2
      PMSA(IP31)  = FIM3
      PMSA(IP32)  = QIM1
      PMSA(IP33)  = QIM2
      PMSA(IP34)  = QIM3
!
!---- End active cells block
!
      ELSE
!
!     Use the new version when IVERSN=1
!
        SWADSP = PMSA(IP1 )
        PO4    = MAX( PMSA(IP2 ) , 0.0)
        AAP    = MAX( PMSA(IP3 ) , 0.0)
        IM1    = MAX( PMSA(IP4 ) , 0.0)
        IM2    = MAX( PMSA(IP5 ) , 0.0)
        IM3    = MAX( PMSA(IP6 ) , 0.0)
        KDADS  = PMSA(IP7 )
        FCAP   = PMSA(IP8 )
        DELT   = PMSA(IP9 )
        KSORP  = PMSA(IP12)
        FRFE1  = PMSA(IP14)
        FRFE2  = PMSA(IP15)
        FRFE3  = PMSA(IP16)
!
!     Calculation of the total concentration of iron and
!     the fractions of adsorbed phosphate in the inorganic matter
!     fractions IM1-3
!
        TFE  = FRFE1 * IM1 + FRFE2 * IM2 + FRFE3 * IM3
        QIM1 = 0.0
        QIM2 = 0.0
        QIM3 = 0.0
        IF ( TFE .LT. 1E-10) THEN
           FIM1 = 0.0
           FIM2 = 0.0
           FIM3 = 0.0
        ELSE
           FIM1 = FRFE1 * IM1 / TFE
           FIM2 = FRFE2 * IM2 / TFE
           FIM3 = FRFE3 * IM3 / TFE
           IF ( IM1 .GT. 1E-10 ) QIM1 = AAP*FIM1/IM1
           IF ( IM2 .GT. 1E-10 ) QIM2 = AAP*FIM2/IM2
           IF ( IM3 .GT. 1E-10 ) QIM3 = AAP*FIM3/IM3
        ENDIF
!
!     Error messages
!
        IF (KDADS .LT. 0.0) CALL ERRSYS ('KDADS in ADSPO4 negative', 1 )
        IF (((ABS(KDADS) .LT. 1E-20) .OR. (ABS(PO4) .LT. 1E-20))
     +     .AND. (ABS(FCAP) .LT. 1E-20) )
     +     CALL ERRSYS
     +     ('(KDADS or PO4) and FCAP equal zero in ADSPO4', 1 )
!
!     Start the calculation of the sorption flux
!     Use one of three options
!
        CADS  = 0.0
        CADST = 0.0
        EQAAP = 0.0
        FADS  = 0.0
        KADS  = 0.0
!
!     SWADSP = 0 : Instantaneous equilibrium partitioning
!
        IF (NINT(SWADSP) .EQ. 0) THEN
            FADS  =(((AAP + PO4) / (1.0 + KDADS)) - AAP) / DELT
            EQAAP = -1.0
        ENDIF
!
!     SWADSP = 1 : Langmuir sorption
!
        IF (NINT(SWADSP) .EQ. 1) THEN
          IF ( (FCAP .LT. 1E-10) .OR. (KDADS .LT. 1E-10)) THEN
             EQAAP = 0.0
             FADS  = 0.0
          ELSE
             CADST = FCAP * TFE
             EQAAP = (KDADS * CADST * PO4) / (KDADS * PO4 + 1.0)
             FADS  = KSORP * (EQAAP - AAP)
          ENDIF
        ENDIF
!
!     SWADSP = 2 : pH dependent Langmuir sorption
!
        IF (NINT(SWADSP) .EQ. 2) THEN
           KADS20   = PMSA(IP10)
           TC       = PMSA(IP11)
           AOH      = PMSA(IP13)
           FRFEOX   = PMSA(IP17)
           OXY      = PMSA(IP18)
           CROXY    = PMSA(IP19)
           PH       = PMSA(IP20)
           TEMP     = PMSA(IP21)
           POROS    = PMSA(IP22)

           IF (POROS .LT. 1E-10) CALL ERRSYS
     +        ('POROS in ADSPO4 equals zero', 1 )
           IF ( TFE .LT. 1.E-10 ) THEN
              IF ( NR_MES .LT. 25 ) THEN
                 CALL GETMLU(ILUMON)
                 NR_MES = NR_MES + 1
                 WRITE ( ILUMON , * ) 'WARNING :zero TFE in ADSPO4',
     +                                ' segment=',ISEG,' TFE=',TFE
              ENDIF
              IF ( NR_MES .EQ. 25 ) THEN
                 CALL GETMLU(ILUMON)
                 NR_MES = NR_MES + 1
                 WRITE(ILUMON,*) ' 25 WARNINGS on zero TFE'
                 WRITE(ILUMON,*) ' Further messages on algae surpressed'
              ENDIF
           ENDIF
!
!     Calculate pH dependency (hydroxyl) and factor for redox potential
!
           OH = 10.0**(PH-14.0)
           IF (OXY .GE. (CROXY*POROS) ) THEN
              FRA = 1.0
           ELSE
              FRA = FRFEOX
           ENDIF
!
!     Calculate total sorption capacity (mol/l)
!
           CADST = FRA * TFE / (56000.0 * POROS)
!
!     Calculate free sorption capacity (mol/l)
!
           CADS  = CADST - (AAP / (31000.0 * POROS) )
!
!     Calculate temperature corrected KADS
!
           KADS  = KADS20 * TC**(TEMP-20.0)
!
!     Calculate equilibrium concentration of adsorbed P (gP/m3)
!
           IF ( ABS(CADS) .LT. 1.E-20 ) THEN
              EQAAP = 0.0
           ELSE
              EQAAP = (AAP + PO4)/(1.0 + OH**AOH/(KADS*CADS))
           ENDIF
!
!     Maximize EQAAP on equivalent CADST
!
           IF ( CADS .LT. 0.0 ) THEN
              EQAAPM = 0.9 * CADST * (31000.0 * POROS)
              EQAAP = EQAAPM
           ENDIF
!
!     Calculate the adsorption flux (gP/m3/d)
!
           FADS  = KSORP * (EQAAP - AAP)

        ENDIF
!
!     Output of module
!
      FL(1+IFLUX) = FADS
      PMSA(IP26)  = EQAAP
      PMSA(IP27)  = CADST
      PMSA(IP28)  = KADS
      PMSA(IP29)  = FIM1
      PMSA(IP30)  = FIM2
      PMSA(IP31)  = FIM3
      PMSA(IP32)  = QIM1
      PMSA(IP33)  = QIM2
      PMSA(IP34)  = QIM3
!
!     End active cells block
!
      ENDIF
!
      ENDIF
!
      IFLUX = IFLUX + NOFLUX
      IP1   = IP1   + INCREM (  1 )
      IP2   = IP2   + INCREM (  2 )
      IP3   = IP3   + INCREM (  3 )
      IP4   = IP4   + INCREM (  4 )
      IP5   = IP5   + INCREM (  5 )
      IP6   = IP6   + INCREM (  6 )
      IP7   = IP7   + INCREM (  7 )
      IP8   = IP8   + INCREM (  8 )
      IP9   = IP9   + INCREM (  9 )
      IP10  = IP10  + INCREM ( 10 )
      IP11  = IP11  + INCREM ( 11 )
      IP12  = IP12  + INCREM ( 12 )
      IP13  = IP13  + INCREM ( 13 )
      IP14  = IP14  + INCREM ( 14 )
      IP15  = IP15  + INCREM ( 15 )
      IP16  = IP16  + INCREM ( 16 )
      IP17  = IP17  + INCREM ( 17 )
      IP18  = IP18  + INCREM ( 18 )
      IP19  = IP19  + INCREM ( 19 )
      IP20  = IP20  + INCREM ( 20 )
      IP21  = IP21  + INCREM ( 21 )
      IP22  = IP22  + INCREM ( 22 )
      IP23  = IP23  + INCREM ( 23 )
      IP24  = IP24  + INCREM ( 24 )
      IP25  = IP25  + INCREM ( 25 )
      IP26  = IP26  + INCREM ( 26 )
      IP27  = IP27  + INCREM ( 27 )
      IP28  = IP28  + INCREM ( 28 )
      IP29  = IP29  + INCREM ( 29 )
      IP30  = IP30  + INCREM ( 30 )
      IP31  = IP31  + INCREM ( 31 )
      IP32  = IP32  + INCREM ( 32 )
      IP33  = IP33  + INCREM ( 33 )
      IP34  = IP34  + INCREM ( 34 )
!
 9000 CONTINUE
!
      RETURN
!
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
This user is BannedGroup banned from Davis Wiki for making threats against other individuals and creating sockpuppet accounts.

If you are the individual associated with this account, this is notice that you are no longer permitted to edit the Davis Wiki until this ban is lifted.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
MODULE meshmo_db
  USE param_db,ONLY: msts, mnam
  IMPLICIT NONE
  SAVE

  ! GENERAL VARIABLES
  INTEGER(kind=4):: r_nrfset !>0 if remeshing is asked for a number of sets

  INTEGER(kind=4)::nrmstra=0 ! Strategy remesh counter

  INTEGER(kind=4)::ntrmstra=0 ! Total remesh counter

 ! REAL (kind=8) :: p_elm_size, & !present element size
 !                  ini_dtime     !initial time increment to compare

! VARIABLES ASSOCIATED TO REMESHING

  CHARACTER(len=mnam),POINTER:: r_refset(:)   !label of the element set to modify
  LOGICAL,POINTER:: r_meshmo(:)               ! .TRUE. if remesh is done over a set

  LOGICAL :: r_last = .FALSE. !to know if in previous step it was performed
  !  LOGICAL,POINTER:: r_last(:)  !to know if in previous step it was performed
  LOGICAL :: r_elm_zone = .FALSE., & !true for zone remeshing estrategy
             r_z_actv                !true for distorted element

  REAL (kind=8) :: r_elm_size,  & !uniform element size for the mesher
                   r_elm_ratio, & !ratio to modify elm_size between remeshings
                   r_min_size,  & !minimum element size to trigger remeshing
  !                 r_min_dt_r,  & !minimum time incr. ratio to trigger remeshing
                   r_t_freq,    & !time frequency to check remeshing criterion
                   r_t_star,    & !time start to check remeshing time freq.
                   r_l_dist,    & !limit element distortion
                   r_b_dist,    & !band width of element distortion
                   scalef = 1000d0  !to avoid very small elems in GID

  INTEGER (kind=4) :: r_s_freq, & !step frequency to check remeshing criterion
                      r_s_star, & !step start to check remeshing time freq.
                      r_size,   & !size of r_times array
                      n_first,  & !number of first point in the new mesh
                                  !mesh (to avoid conflict with added nodes)
                      itransf=1,& !transfer of GP variables,
                                  ! = 1, Superconvergent Patch Recovery
                                  ! = 2, based on averaged nodal vars (old)
                      ltype  =1,& !Line Type for boundary interpolation (in 2D)
                                  ! = 1, PolyLine
                                  ! = 2, NurbsLine
                      r_crit      !element distortion measure criteria

  REAL (kind=8), POINTER :: r_times(:)

  ! scratch variable for remeshing
  INTEGER(kind=4):: nelbnd             !Number of element of boundary
  INTEGER(kind=4),POINTER:: lnbnd(:)   !Boundary conectivities
  INTEGER(kind=4),POINTER:: nlbnd(:)   !Number of continuos lines in boundary
  INTEGER(kind=4),POINTER:: nodlb(:,:) !Node labels in boundary line
  INTEGER(kind=4),POINTER:: fside(:)   !Free side array
  REAL(kind=8),POINTER:: cdbnd(:,:)    !Boundary coordinates

  CHARACTER(len=mnam),POINTER:: strnam(:)  !surface names list

  ! ==========================================================================
  ! work (scratch database) - no need to dump for restart

  INTEGER (kind=4), POINTER :: nodset(:),  &   ! (numnp), nodes of the eset
                               l_old(:,:), &   ! (nnode,nelem) old elements
                               elemp(:),   &   ! (numnp)element where projects
                               l_new(:,:)      ! (nnode,ne_new) new elements

  REAL (kind=8), POINTER ::  &
                 strnd(:,:), &  !(nstre,npoin) nodal stresses (old mesh)
                 sn_new(:,:),&  !(nstre,numnp) nodal stresses (new mesh)
                 sg_old(:,:,:),&  !(nvarg,ngaus2,nelem), old gaussian variables
                 sg_new(:,:) ,&  !(nvarg,ngaus2*ne_new) gaus.vars (new mesh)
                 v_new(:,:), &  !(ndofn,numnp) nodal velocities (new mesh)
                 p_new(:)       !(numnp) nodal pressure (new mesh), if split

  ! Derived type for POLYLINE data - used for contour definition
  TYPE pline
    INTEGER (kind=4) :: nfirst, &  ! first node (local numbering on the contour)
                        nlast,  &  ! last node
                        ifix(7)    ! 7=MAX(nd1) fixity code for the polyline
    LOGICAL          :: fside      ! .TRUE. if a free-side
    TYPE (pline), POINTER :: next
  END TYPE pline
  TYPE (pline), POINTER :: headp=>NULL(), tailp=>NULL()  !first and last lines
  INTEGER (kind=4) :: nplins  ! number of lines defining a contour
  INTEGER (kind=4) :: nstart  ! starting node number (nfirst in the first line)
  INTEGER(kind=4),POINTER:: lstrt(:)   !starting node number in each contour line
  INTEGER (kind=4) :: criter  ! ??

  ! new mesh
  INTEGER (kind=4) :: numpo,      &  ! number of points in old mesh set
                      numpn,      &  ! number of points in new mesh set
                      maxnn,      &  ! maximum node number after deleting nodes
                      ne_new         ! number of elements in the new mesh
  REAL   (kind=8), ALLOCATABLE :: &
                      coorn(:,:), &  ! (ndime,numnp) coordinates of a new mesh
                      coran(:,:)     ! (ndime,numnp) coordinates of a new mesh

CONTAINS

  SUBROUTINE dump_meshmo()
    USE ctrl_db,ONLY: npoin
    ! dumps mesh modifications data
    IMPLICIT NONE

    !Local variables

    WRITE(50,ERR=9999) r_nrfset
    WRITE(50,ERR=9999) nrmstra
    WRITE(50,ERR=9999) ntrmstra

    ! variables associated to remeshing
    IF (r_nrfset > 0) THEN
      WRITE(50,ERR=9999) r_refset(1:r_nrfset)
      WRITE(50,ERR=9999) r_meshmo(1:r_nrfset)   !Remesh status
      WRITE(50,ERR=9999) r_last,r_elm_zone, r_z_actv
      WRITE(50,ERR=9999) r_elm_size, r_elm_ratio, r_min_size, r_t_freq, r_t_star,r_l_dist,r_b_dist, &
                scalef, r_s_freq, r_s_star, r_size, n_first, itransf, ltype, r_crit
!      WRITE(50,ERR=9999) r_elm_size, r_elm_ratio, r_min_size, r_min_dt_r, r_t_freq, r_t_star,     &
!                scalef, r_s_freq, r_s_star, r_size, n_first, itransf, ltype
      IF( r_size > 0 )WRITE(50,ERR=9999) r_times(1:r_size)
    END IF

  RETURN
  9999 CALL runen2('')
  END SUBROUTINE dump_meshmo

  SUBROUTINE rest_meshmo()
    USE ctrl_db,ONLY: npoin
    ! reads mesh modifications data for a restart
    IMPLICIT NONE

    !Local variables

    READ(51) r_nrfset
    READ(51) nrmstra
    READ(51) ntrmstra

    ! variables associated to remeshing
    IF (r_nrfset > 0) THEN
      ALLOCATE(r_meshmo(r_nrfset),r_refset(r_nrfset))
      READ(51) r_refset(1:r_nrfset)
      READ(51) r_meshmo(1:r_nrfset)   !Remesh status
      READ(51) r_last,r_elm_zone, r_z_actv
      READ(51) r_elm_size, r_elm_ratio, r_min_size, r_t_freq, r_t_star,r_l_dist,r_b_dist, &
               scalef, r_s_freq, r_s_star, r_size, n_first, itransf, ltype, r_crit
!      READ(51) r_elm_size, r_elm_ratio, r_min_size, r_min_dt_r, r_t_freq, r_t_star,      &
!               scalef, r_s_freq, r_s_star, r_size, n_first, itransf, ltype
      IF( r_size > 0 )THEN
        ALLOCATE( r_times(r_size) )
        READ(51)  r_times(1:r_size)
      END IF
    END IF

    RETURN
  END SUBROUTINE rest_meshmo

  !===================== routines to handle POLYLINES ======================
  SUBROUTINE new_pline(elm)
  !initialize a list of polylines
  IMPLICIT NONE

    !--- Dummy arguments
    TYPE(pline),POINTER:: elm

    ALLOCATE(elm)
    elm%nfirst = 0
    elm%nlast = 0
    elm%fside = .FALSE.
    elm%ifix(1:7) = 0
    NULLIFY(elm%next)

  RETURN
  END SUBROUTINE new_pline

  SUBROUTINE add_pline (new, head, tail)
    !This subroutine adds data to the end of the list
    !Dummy arguments
    TYPE (pline), POINTER :: new, head, tail

    !Check if a list is empty
    IF (.NOT. ASSOCIATED (head)) THEN
      !list is empty, start it
      head => new
      tail => new
      NULLIFY (tail%next)

    ELSE
      !add data to the list
      tail%next => new
      NULLIFY (new%next)
      tail => new

    ENDIF
  END SUBROUTINE add_pline

  SUBROUTINE delete_pline(head,tail)
  !This subroutine adds data to the end of the list
  IMPLICIT NONE
    !--- Dummy arguments
    TYPE(pline),POINTER:: head, tail
    !--- Local variables
    TYPE(pline),POINTER:: elm

    NULLIFY(tail)
    DO
      IF (.NOT.ASSOCIATED(head)) EXIT
      elm => head%next
      DEALLOCATE(head)
      head => elm
    END DO

  RETURN
  END SUBROUTINE delete_pline

  !===================== routines to handle REFSET names====================
  SUBROUTINE search_refset(name,list,found)
  !search for name in a list of REFSETs previously stored
  IMPLICIT NONE
    !--- Dummy arguments
    CHARACTER(len=mnam), INTENT(IN) :: name    ! refset NAME
    CHARACTER(len=mnam), POINTER    :: list(:) ! LIST of previous stored REFSETs
    LOGICAL, INTENT(OUT) :: found  ! .TRUE. when NAME is stored in LIST
    !--- local variables
    INTEGER(kind=4) :: i,ndim

    found =.FALSE. !initialize
    IF (ASSOCIATED(list)) THEN

      ndim = SIZE(list) !Determine LIST vector dimension
      DO i = 1,ndim  ! loop in list of REFSETs stored
        IF ( TRIM(list(i)) == TRIM(name) ) THEN
          found = .TRUE. ! if refset is previuosly stored
          EXIT
        END IF
      END DO

    END IF

  RETURN
  END SUBROUTINE search_refset

END MODULE meshmo_db

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:
      subroutine tpraa(alamaa,atheta,iaapr,ipbtmx,jpfcmx,naapr,
     $ natmax,ncvaa,nerr,noutpt,npx3mx,nthdt,nttyo,nwarn,pcvaa,
     $ qpdaa,uaqsp,uthdtr)
c
c     Test and process the Pitzer data for aa' (anion, different
c     anion) pairs read from the DATA0 file. Find and flag errors,
c     such as duplication of data (e.g., two data blocks for the same
c     aa' pair). Calculate the conventional primitive Pitzer parameters
c     from the observable compound parameters read from the data file:
c
c       theta(aa') -> lambda(aa')
c
c     Check the coverage of entered Pitzer data against all possible
c     anion-distinct anion pairs that can be composed of the aqueous
c     anions present on the data file.
c
c     This subroutine is complementary to tprcc.f, which does the
c     same thing for cc' (cation, different cation) pairs.
c
c     This subroutine is called by:
c
c       EQPT/eqpt.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       nat    = the number of aqueous species
c       uaqsp  = array of names of aqueous species
c
c     Principal output:
c
c       nerr   = cumulative error counter
c       nwarn  = cumulative warning counter
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ipbtmx,jpfcmx,naapr,natmax,npx3mx
c
      integer noutpt,nttyo
c
      integer iaapr(2,naapr)
c
      integer ncvaa,nerr,nthdt,nwarn
c
      logical qpdaa(naapr)
c
      character(len=24) uaqsp(natmax),uthdtr(3,npx3mx)
c
      real*8 alamaa(jpfcmx,0:ipbtmx,naapr),atheta(jpfcmx,npx3mx)
c
      real*8 pcvaa
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,j,jth,jthpr,j2,j3,j4,j5,n,ncount,ndupl,nlistl,
     $ nn,nodatc
c
      integer ilnobl
c
      logical qmatch,qzerov
c
      character*56 ustr56
      character*24 unam1,unam2
      character*8 ux8
c
c-----------------------------------------------------------------------
c
c     Limit on the list of pairs for which no data were found.
c
      data nlistl / 20 /
c
c-----------------------------------------------------------------------
c
c     Initialize the data arrays.
c
      do n = 1,naapr
        qpdaa(n) = .false.
        do j = 1,jpfcmx
          do i = 0,ipbtmx
            alamaa(j,i,n) = 0.
          enddo
        enddo
      enddo
c
c     Check the entered data for aa' pairs.
c
      nodatc = 0
c
      do n = 1,naapr
        i = iaapr(1,n)
        j = iaapr(2,n)
        unam1 = uaqsp(i)
        unam2 = uaqsp(j)
c
        jthpr = 0
        ndupl = 0
c
c       Search for unam1, unam2 in the species triplet blocks.
c       This is actually done by searching the uthdtr array,
c       not the utripl array.
c
        do jth = 1,nthdt
          if (unam1(1:24) .eq. uthdtr(1,jth)(1:24)) then
            if (unam2(1:24) .eq. uthdtr(2,jth)(1:24)) then
              jthpr = jth
              go to 110
            endif
          endif
        enddo
  110   continue
c
        if (jthpr .gt. 0) then
c
c         Have found an entry in the triplet blocks.
c
          qpdaa(n) = .true.
c
c         Store the data.
c
          do j = 1,jpfcmx
            alamaa(j,0,n) = atheta(j,jthpr)
            do i = 1,ipbtmx
              alamaa(j,i,n) = 0.
            enddo
          enddo
c
c         Search for duplicates in the species triplets blocks.
c         Ignore duplications in the triplet blocks if the
c         values are all zeros (take this to mean "no data input")
c         or if the values all match the first-encountered
c         values (exact duplication).
c
          do jth = jthpr + 1,nthdt
            if (unam1(1:24) .eq. uthdtr(1,jth)(1:24)) then
              if (unam2(1:24) .eq. uthdtr(2,jth)(1:24)) then
                qzerov = .true.
                do j = 1,jpfcmx
                  if (atheta(j,jth) .ne. 0.) qzerov = .false.
                enddo
                qmatch = .true.
                do j = 1,jpfcmx
                  if ((atheta(j,jth) - atheta(j,jthpr)) .gt. 1.e-6)
     $            qmatch = .false.
                enddo
                if (.not.qzerov .and. .not.qmatch)
     $          ndupl = ndupl + 1
              endif
            endif
          enddo
c
          if (ndupl .gt. 0) then
            j2 = ilnobl(unam1)
            j3 = ilnobl(unam2)
            if (ndupl .eq. 1) then
              write (noutpt,1010) unam1(1:j2),unam2(1:j3)
              write (nttyo,1010) unam1(1:j2),unam2(1:j3)
 1010         format(/' * Error - (EQPT/tpraa) Have found a',
     $        ' duplicate data block on the DATA0 file',
     $        /7x,"for the aa' pair ",a,', ',a,'.')
            else
              ux8 = ' '
              write (ux8,'(i5)') ndupl
              call lejust(ux8)
              j5 = ilnobl(ux8)
              write (noutpt,1020) ux8(1:j5),unam1(1:j2),unam2(1:j3)
              write (nttyo,1020) ux8(1:j5),unam1(1:j2),unam2(1:j3)
 1020         format(/' * Error - (EQPT/tpraa) Have found ',a,
     $        ' duplicate data blocks on the DATA0 file',
     $        /7x,"for the aa' pair ",a,', ',a,'.')
            endif
            nerr = nerr + ndupl
          endif
        endif
c
        if (jthpr.le.0) then
c
c         No data block was found on the DATA0 file.
c         Note that qpdaa(n) is left with a value of .false.
c
          nodatc = nodatc + 1
        endif
c
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (nodatc .gt. 0) then
        write (noutpt,1040)
        write (nttyo,1040)
 1040   format(/' * Warning - (EQPT/tpraa) Did not find a data',
     $  ' block on the DATA0 file',/7x,'for any of the following',
     $  " aa' pairs:",/)
c
        ncount = 0
        do n = 1,naapr
          if (.not.qpdaa(n)) then
            ncount = ncount + 1
            i = iaapr(1,n)
            j = iaapr(2,n)
            unam1 = uaqsp(i)
            unam2 = uaqsp(j)
            j2 = ilnobl(unam1)
            j3 = ilnobl(unam2)
            ustr56 = unam1(1:j2) // ', ' // unam2(1:j3)
            j4 = ilnobl(ustr56)
            write (noutpt,1050) ustr56(1:j4)
            write (nttyo,1050) ustr56(1:j4)
 1050       format(9x,a)
            if (ncount .eq. nlistl) go to 200
          endif
        enddo
  200   continue
c
        nn = nodatc - ncount
        if (nn .gt. 0) then
          write (ux8,'(i5)') nn
          call lejust(ux8)
          j3 = ilnobl(ux8)
          write (noutpt,1060) ux8(1:j3)
          write (nttyo,1060) ux8(1:j3)
 1060     format(/9x,'plus ',a,' others')
        endif
        write (noutpt,1070)
        write (nttyo,1070)
 1070   format(1x)
        nwarn = nwarn + 1
      endif
c
      ncvaa = naapr - nodatc
      pcvaa = (100.*ncvaa)/float(naapr)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
  999 continue
c
      end


text:
A marketing campaign put on by the Campus Crusade for Christ (CCC) and other Christian groups in early 2002. The simple slogan I Found It was printed on flyers, silkscreened on bright orange shirts, and scribbled on chalkboards. Students wondered what the I Found It slogan meant and what it was trying to sell. After weeks of masking the campus with their slogan and keeping their interest a secret, the Christian organization revealed that what they found was Jesus Christ.

Whether or not the campaign managed to convert students is hard to say, but it kept students talking and interested in the I Found It phenomenon. The campaign was spoofed by the humorous AS Papers publication The Sneeze. There was also a shirt by an anonymous designer that simply said I Smoked It.

A couple of the parodies that were spotted include:
     I lost it. (in the insane sense)
     I lost it at UC Davis. (in the virginity sense)
     I smoked it. (In the smoking marijuana sense)

Those who observed the campaign may have also noticed that when the 2002 campaign was over, the campus was littered with bright orange scraps of paper and broken balloons. SPAC received several complaints about this, and reportedly held a meeting with organizers about it.

Two years before the I Found It campaign, the CCC had a nearly identical campaign slogan, I agree with Jessica: Black text on lime green. This campaign culminated with a similar event, where the speaker was a woman, well, named Jessica. It was not as ubiquitous as the I Found It campaign.

According to the CCCs http://www.demossnewspond.com/ccci/presskit/execprofiles.htm website newsroom, I Found It is a registered trademark. According to the http://www.uspto.gov US Patent and Trademark Office, this http://tess2.uspto.gov/bin/showfield?fdoc&statef9tmim.2.4 registration is http://tess2.uspto.gov/bin/gate.exe?fhelp&statef9tmim.2.4#FAQ_DEAD Dead.
Sadly due to policy you will need to run your own search as they expire.~Dave

 2007 I found it too campaign

http://ifoundittoo.com/ Website

Media Coverage: The California Aggie Aggie http://www.californiaaggie.com/news/2007/05/07/CampusNews/i.Found.It.Too.Campaign.Provokes.Controversy.Debate2896327.shtml article

 2007 Parodies
I Found it too... NOT! Sponsored by AGASA.
Some of the discussion regarding AGASAs parody can be seen on the http://community.livejournal.com/ucdavis/3763244.html davis communitys livejournal

ALSO, some facebook.com groups like http://ucdavis.facebook.com/group.php?gid2301308591 Pipe Down, You Didnt Find Shit have developed to counter this campaign.

 Sponsoring organizations

  Campus Crusade for Christ CCC http://www.kenandamy.com/prayerletters/March%2042007.pdf Newsletter for March 31, 2007.
  Davis Christian Fellowship

I saw the bumper sticker (black on yellow) from the original campaign many times around the Portland, OR area.  Someone came up with a great variation (rejoinder) on the theme:  I Found It, and underneath (in small letters) and now my finger stinks

The I Found It campaign is not a UC Davis invention, nor is it new...  A Christian group called Campus Crusade for Christ Campus Crusade ran a huge media campaign in 1976 and published the slogan on bulletin boards, bumper stickers, and pins...  I remember seeing the bumper sticker I Lost It when I was a kid, and asking my mom what it meant.  I always thought it meant losing the Christian religion, not losing your mind.   Users/SummerSong

My high school did a similar ploy when I was there; the shirts said I agree with Jason.  Jason was a local youth leader who came to speak on the last day of the campaign.  Problem was, the whole school thought the campaign had something to do with our rivalry game at the end of the week, and that whoever Jason was, he had something against the game. Users/AllisonEriksen

This campaign deserves a high shadyness rating. It is deceitful. Obscuring the truth is deceit and so too is using the pronoun it to intentially hide the meaning jesus. This tactic is OK for commercial advertising campaigns because they dont sell a way of life. However, selling a good lifestyle shouldnt involve obfuscation and propaganda. In my opinion, hiding the truth shouldnt mesh with christianity, yet with this campaign it does. Why not be honest and say what is meant, I found Jesus.  Shadyness.

As an aside note, after the ebay it campaign doesnt this liken jesus to a cheap plastic figurine purchased on ebay? no?
 A little, but honestly I rather have it than a religious faith. Personally I dont see why people were not upfront about this, anyone want to flashmob the next event they have, just to complain about not getting enough of it, or wanting more of it, or how it tastes so good, feels so fine, etc. ~Dave

I wikied it. ~Dave

20070430 17:23:00 nbsp I Photographed It? Users/JasonAller

 ~This doesnt show on the US Patent and Trademark Office searches. ~Dave
  ~ David Grundler doesnt show up on US Patent and Trademark Office http://tess2.uspto.gov/bin/showfield?ftoc&statef9tmim.1.1&p_searchsearchss&p_L50&BackReference&p_pluralyes&p_s_PARA1&p_tagrepl%7E%3APARA1%24LD&exprPARA1+AND+PARA2&p_s_PARA2David+Grundler&p_tagrepl%7E%3APARA2%24ALL&p_op_ALLADJ&a_defaultsearch&a_searchSubmit+Query&a_searchSubmit+Query searches.
   ~Thank goodness!  How much would that suck to have your name trademarked by someone else?
   ~Dang, I need to trademark the Tilde, but I am sure Coke^tm^ has it already. ~Dave

20070430 18:06:41 nbsp I popped it. Users/GregWebb

20070501 22:21:30 nbsp As to the question about this campaigns success, I can tell you that it is very successful. Our point is not to convert people, we are trying to introduce people to something that has changed our lives. We are not forcing it on anyone, we are inviting people to ask questions. As to the shadyness claim... no, we are not trying to deceive anyone. The whole point is for people to ask us about it. If we wore shirts that said I found Jesus, would any of you listen to anything we had to say? So if youre open to serious conversation about who Jesus is (or isnt) and what all this means, stop by the table on the quad.  Users/lauren

20070502 01:35:29 nbsp Someone in a lime green shirt asked me today, on the MU patio, what was going on while the UCD May 1, 2007: Day of Action May Day rally was occurring. I told him the issues that were being protested and he nodded along and seemed rather apathetic or nonsupportive. I mentioned that, since green and white were the suggested colors to were in the May Day event that the lime green shirt folks should join us, especially on issues of illegal immigration, as Leviticus 3334 states, When an alien lives with you in your land, do not mistreat him. The alien living with you must be treated as one of your nativeborn. Love him as yourself. The guy I talked to was not so enthused to take a biblical stance that countered conservative values and walked off.  Users/RobRoy

20070502 01:38:42 nbsp I understand the desire for church. It is a melding of community and fear of death. But this Ifoundittoo mumbo gumbo is typical Christian trickery of hiding the Jesus to lure in those that shy away from religion. In my childhood to teenage years I went to church over a thousand times, I know how they operate. I support your right to have the crutch of religion. But being forthright about what the real message is would be more respectful to the campus community. This is just like the Porn Nation event from early 2006. It was a Christian event but it was not marketed as such, they wanted you to show up, lured in by sex, and then be forced to listen to someones sermon of finding Jesus. No thanks; crusades are for insecure, accepting the world as the place that God made and converting those to be like them as a way of selfjustification. Id much rather march in a revolution, because do you really want to worship a God that created a world with George Bush as its leader? Users/RobRoy
  I kind of sympathize with this kind of campaign.  Its hard to start a conversation about faith, and they mainly use this as a way to start a dialog with people.  Its easier to start talking about my Discordian faith than my Christian faith, because people have less preconceptions.  Among those preconceptions that God appoints the leaders of this world.  Come on Rob, the founders of this country got over the whole divine right of kings thing, and decided that humans appoint their own leaders.  Im sure you can come to that conclusion as well.  They should of course support our rally too, as they are already correctly dressed.   Users/BrentLaabs

20070502 03:19:26 nbsp Well, I know that God does not choose our leaders  but that is because I do not believe in God. But, the Bible can contradict this  as it is ALL Gods plan. Did God give us free will? The first few chapters of Genesis would lead one to believe so but the Bible is full of rhetoric stating that god plans everything. I dont want to get into a theological debate. All religions are the same to me  suspensions of logical disbelief to get us thru our fear of death. But if that is what someone needs to get thru the day, whatevs. It is the bright green papers, the ganglike congregation of Tshirts and the coy advertising that annoys me. I know that I have annoyed many folks, knocking on their doors asking for a vote  but when I do, I say my cause forthright. I am not ambiguous, but with this campaign it makes me wonder, is Christ something to be ashamed of, because if so, the rooster crows three times for the lime shirted Peters. Users/RobRoy

20070502 03:39:30 nbsp I have this burning desire to film semimalicious acts on some of their folk, A la Messin with Sasquatch... but I always stop short... Users/CarlosOverstreet
   While there are some Pentacostals in CCC, they dont even approach a majority.  Users/WilliamLewis wl
   Yeah, I just thought it sounded better than Evangelical in that phrase. Users/CarlosOverstreet

20070507 21:04:53 nbsp Theres nothing to find. Users/DanXie

20070508 18:25:23 nbsp I dislike the it conotation, I could have sworn that jesus was male... Users/StevenDaubert

20070509 19:04:47 nbsp Im pretty sure it refers to a personal relationship with God, thus it isnt used as a pronoun referring to Jesus, but the noun relationship. Users/EvLiu

20090825 10:59:22 nbsp I was in middle school in Denver in 1976 when this campaign first launched.  Had great fun calling the 800 number and signing up all my Jewish frieds to receive the campaign kit, phone calls and other junk. Hilarious! Users/Judas

20100914 23:14:06 nbsp Actually, I Found It started even before 76 if memory serves because my aunt and I were recruited and it was before I graduated HS so it had to be before that. I also remember they were into record burning the whole you can hear the devils message when you play it enter rock groups name here backwards, subliminal messages in just about EVERYthing... blah, blah, blah ad nauseum... Users/aleabethlanaykra


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:35:01 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE RAININI__genmod
          INTERFACE
            SUBROUTINE RAININI
            END SUBROUTINE RAININI
          END INTERFACE
        END MODULE RAININI__genmod


text:
! Implementation code for vector_analysis module
submodule (vector_analysis) vector_analysis_sm

contains

    ! Calculate the gradient
    module function grad(p, phi) result(r)
        real(8) :: r(3)
        real(8), intent(in) :: p(3)
        procedure(f)        :: phi
        type(auto_var) :: q(3), s
        integer        :: i

        do i=1,3
            call make_var(p, q, i)
            s = phi(q)
            r(i) = s%get_derivative()
        end do

    end function grad

    ! Calculate the divergence
    module function div(p, phi)
        real(8) :: div
        real(8), intent(in) :: p(3)
        procedure(f)        :: phi
        real(8) :: q(3)

        q = grad(p, phi)

        div = sum(q)
    end function div

    !! Calculate the del_squared (== div(grad))
    !module function del_squared(x, y, phi)
    !    real(8) :: del_squared
    !    real(8), intent(in) :: x, y
    !    procedure(f)        :: phi
    !    real(8) :: grad(2)
    !    type(auto_var) :: xa, ya, ra
    !
    !    ! x-component varies, y-component constant
    !    call xa%set(x)
    !    call ya%set_constant(y)
    !    ra = phi(xa, ya)
    !    grad(1) = ra%get_derivative()
    !
    !    ! y-component varies, x-component constant
    !    call xa%set_constant(x)
    !    call ya%set(y)
    !    ra = phi(xa, ya)
    !    grad(2) = ra%get_derivative()
    !
    !    del_squared = sum(grad**2)
    !end function del_squared

    ! This is actually the z-component of curl
    !module function curl(x, y, phi)
    !    real(8) :: curl
    !    real(8), intent(in) :: x, y
    !    procedure(f)        :: phi
    !    real(8) :: grad(2)
    !    type(auto_var) :: xa, ya, ra
    !
    !    ! x-component varies, y-component constant
    !    call xa%set(x)
    !    call ya%set_constant(y)
    !    ra = phi(xa, ya)
    !    grad(1) = ra%get_derivative()
    !
    !    ! y-component varies, x-component constant
    !    call xa%set_constant(x)
    !    call ya%set(y)
    !    ra = phi(xa, ya)
    !    grad(2) = ra%get_derivative()
    !
    !    curl = grad(1) - grad(2)
    !end function curl

    ! Evaluate the given function
    module function evaluate(p, phi)
        real(8) :: evaluate
        real(8), intent(in) :: p(3)
        procedure(f)        :: phi
        type(auto_var) :: q(3), ra
        integer :: i

        ! Make parameters constant - no derivatives
        do i=1,3
            call q(i)%set_constant(p(i))
        end do
        ra = phi(q)
        evaluate = ra%get_value()
    end function evaluate

    ! Cross product of two vectors
    module pure function cross(a, b) result(r)
        real(8), intent(in) :: a(3), b(3)
        real(8) :: r(3)
        integer :: i, j, k

        do i=1,3
            j = merge(1,i+1,i==3)
            k = merge(1,j+1,j==3)
            r(i) = a(j)*b(k)-a(k)*b(j)
        end do
    end function cross

end submodule vector_analysis_sm

max_stars_repo_name: sgeard/code

max_issues_repo_name: sgeard/code

max_forks_repo_name: sgeard/code


text:
C$Procedure      ZZHLP014 ( private help text )

      SUBROUTINE ZZHLP014 ( BEGIN, FINISH, TEXT )

C$ Abstract
C
C     Fill out a portion of the help text needed by percy.
C
C     Private routine intended solely for the support of Inspekt
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     PRIVATE
C
C$ Declarations

      IMPLICIT NONE
      INTEGER               BEGIN ( * )
      INTEGER               FINISH( * )
      CHARACTER*(*)         TEXT  ( * )

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     BEGIN      O   Indexes of begins of text help
C     FINISH     O   Indexes of ends of text help
C     TEXT       O   A block of text help.
C
C$ Exceptions
C
C     Error free.
C
C$ Particulars
C
C     This routine simply fills begin and end markers as well
C     as actual text for a block of help text for percy.
C
C$ Examples
C
C     None.
C
C$ Restrictions
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber      (JPL)
C
C$ Literature_References
C
C     None.
C
C$ Version
C
C-    Inspekt Version 1.0.0, 1-AUG-1997 (WLT)
C
C
C-&

      INTEGER               I
      INTEGER               J

      J = FINISH ( 1 )
      I = BEGIN  ( 1 )

      FINISH(1) = J
      BEGIN (1) = I

      CALL REPMC ( TEXT(1), '*', '*', TEXT(1) )
      TEXT( 1308 ) = '[NOT] column <relation> value'
      TEXT( 1309 ) = '|endliteral'
      TEXT( 1310 ) = 'If'
      TEXT( 1311 ) = 'the column is a character or time colum'
     .//             'n, the value must be enclosed in either'
     .//             ' s'
      TEXT( 1312 ) = 'Allowed relations are EQ NE LT LE GT GE'
     .//             ' and LIKE (used for pattern'
      TEXT( 1313 ) = 'matching).'
      TEXT( 1314 ) = ' '
      TEXT( 1315 ) = '@@Looking at Data    --- SELECT'
      TEXT( 1316 ) = 'Quit Help'
      TEXT( 1317 ) = 'Help'
      TEXT( 1318 ) = 'Column and Table Abbreviations'
      TEXT( 1319 ) = 'Select Clause'
      TEXT( 1320 ) = 'From Clause'
      TEXT( 1321 ) = 'Where Clause'
      TEXT( 1322 ) = 'Order By'
      TEXT( 1323 ) = 'Combining Tables'
      TEXT( 1324 ) = 'Reports'
      TEXT( 1325 ) = 'Getting Too Much Data'
      TEXT( 1326 ) = ' '
      FINISH( 31 ) = 1326

      BEGIN ( 32 ) = 1327
      TEXT( 1327 ) = 'When you select a help topic for which '
     .//             'there is some text that'
      TEXT( 1328 ) = 'should be displayed, the help system be'
     .//             'gins sending this text'
      TEXT( 1329 ) = 'to your display.  If there is a lot of '
     .//             'text, some of it may'
      TEXT( 1330 ) = 'scroll by before you have a chance to r'
     .//             'ead it.  There'
      TEXT( 1331 ) = 'are two ways to deal with this.'
      TEXT( 1332 ) = '@newlist'
      TEXT( 1333 ) = '@numitem You can hit CTRL-S on your key'
     .//             'pad to cause output'
      TEXT( 1334 ) = 'to your display to be temporarily disab'
     .//             'led. Hit CTRL-Q to restart'
      TEXT( 1335 ) = 'the output.  This works on most termina'
     .//             'ls and terminal emulators.'
      TEXT( 1336 ) = ' '
      TEXT( 1337 ) = '@numitem You can tell Inspekt to wait o'
     .//             'nce it finishes displaying'
      TEXT( 1338 ) = 'a page full of text.'
      TEXT( 1339 ) = ' '
      TEXT( 1340 ) = 'To do this at the Inspekt prompt type:'
      TEXT( 1341 ) = '@literal'
      TEXT( 1342 ) = 'Inspekt> SET HELP WAIT;'
      TEXT( 1343 ) = '|endliteral'
      TEXT( 1344 ) = 'The Inspekt help system will now pause '
     .//             'after each page of'
      TEXT( 1345 ) = 'text it displays and wait for you to hi'
     .//             't a carriage return before'
      TEXT( 1346 ) = 'it displays the next page or related to'
     .//             'pics menu.'
      TEXT( 1347 ) = ' '
      TEXT( 1348 ) = 'To return to the original help system b'
     .//             'ehaviour, type the command'
      TEXT( 1349 ) = '@literal'
      TEXT( 1350 ) = 'Inspekt> SET HELP NO WAIT;'
      TEXT( 1351 ) = '|endliteral'
      TEXT( 1352 ) = ' '
      TEXT( 1353 ) = 'Note that the size of a page is also un'
     .//             'der your control.  You may'
      TEXT( 1354 ) = 'set the number of lines that will fit o'
     .//             'n a page by using the'
      TEXT( 1355 ) = 'command "SET PAGE HEIGHT".  This comman'
     .//             'd is described in the'
      TEXT( 1356 ) = '"SET PAGE ..." help topic.'
      TEXT( 1357 ) = ' '
      TEXT( 1358 ) = '@@Making Help Wait'
      TEXT( 1359 ) = 'Quit Help'
      TEXT( 1360 ) = 'Help'
      TEXT( 1361 ) = 'SET PAGE   ...'
      FINISH( 32 ) = 1361

      BEGIN ( 33 ) = 1362
      TEXT( 1362 ) = 'A numeric format is specified by creati'
     .//             'ng a picture of the'
      TEXT( 1363 ) = 'format.  For example to specify that a '
     .//             'number should start'
      TEXT( 1364 ) = 'with 3 digits and be displayed to 3 dec'
     .//             'imal places use'
      TEXT( 1365 ) = 'a picture such as this:'
      TEXT( 1366 ) = '@literal'
      TEXT( 1367 ) = '###.###'
      TEXT( 1368 ) = '|endliteral'
      TEXT( 1369 ) = ' '
      TEXT( 1370 ) = 'If the first character of the picture i'
     .//             's a minus sign,'
      TEXT( 1371 ) = 'the first character in the output strin'
     .//             'g will be'
      TEXT( 1372 ) = 'a blank if the number is non-negative, '
     .//             'a minus sign'
      TEXT( 1373 ) = 'if the number is negative.'
      TEXT( 1374 ) = ' '
      TEXT( 1375 ) = 'If the first character of the picture i'
     .//             's a plus sign,'
      TEXT( 1376 ) = 'the first character of the output strin'
     .//             'g will be a'
      TEXT( 1377 ) = 'plus if the number is positive, a blank'
     .//             ' if the number'
      TEXT( 1378 ) = 'is zero, and a minus sign if the number'
     .//             ' is negative.'
      TEXT( 1379 ) = ' '
      TEXT( 1380 ) = 'If the first character of the string is'
     .//             ' NOT a sign'
      TEXT( 1381 ) = '(plus or minus) the first character of '
     .//             'the output'
      TEXT( 1382 ) = 'string will be a minus sign if the numb'
     .//             'er is negative'
      TEXT( 1383 ) = 'and will be the first character of the '
     .//             'integer part'
      TEXT( 1384 ) = 'of the number otherwise.'
      TEXT( 1385 ) = ' '
      TEXT( 1386 ) = 'The integer portion of STRING will cont'
     .//             'ain at least'
      TEXT( 1387 ) = 'as many characters as appear before the'
     .//             ' decimal point'
      TEXT( 1388 ) = '(or last character if there is no decim'
     .//             'al point) but'
      TEXT( 1389 ) = 'after a leading + or -. There will ALWA'
     .//             'YS be at least'
      TEXT( 1390 ) = 'one digit output in integer portion of '
     .//             'STRING.'
      TEXT( 1391 ) = ' '
      TEXT( 1392 ) = 'If the picture begins with a any of the'
     .//             ' following'
      TEXT( 1393 ) = '@literal'
      TEXT( 1394 ) = '   ''+0'', ''-0'', or ''0'''
      TEXT( 1395 ) = '|endliteral'
      TEXT( 1396 ) = 'it is said to have a leading zero.  If '
     .//             'a picture has'
      TEXT( 1397 ) = 'a leading zero and the integer portion '
     .//             'is not large'
      TEXT( 1398 ) = 'enough to fill up the integer space spe'
     .//             'cified by'
      TEXT( 1399 ) = 'the picture, the output will be zero pa'
     .//             'dded from the sign (if'
      TEXT( 1400 ) = 'one is required) up to the first charac'
     .//             'ter of the'
      TEXT( 1401 ) = 'integer part of the number.'

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
      SUBROUTINE TG01DD( COMPZ, L, N, P, A, LDA, E, LDE, C, LDC, Z, LDZ,
     $                   DWORK, LDWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To reduce the descriptor system pair (C,A-lambda E) to the
C     RQ-coordinate form by computing an orthogonal transformation
C     matrix Z such that the transformed descriptor system pair
C     (C*Z,A*Z-lambda E*Z) has the descriptor matrix E*Z in an upper
C     trapezoidal form.
C     The right orthogonal transformations performed to reduce E can
C     be optionally accumulated.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     COMPZ   CHARACTER*1
C             = 'N':  do not compute Z;
C             = 'I':  Z is initialized to the unit matrix, and the
C                     orthogonal matrix Z is returned;
C             = 'U':  Z must contain an orthogonal matrix Z1 on entry,
C                     and the product Z1*Z is returned.
C
C     Input/Output Parameters
C
C     L       (input) INTEGER
C             The number of rows of matrices A and E.  L >= 0.
C
C     N       (input) INTEGER
C             The number of columns of matrices A, E, and C.  N >= 0.
C
C     P       (input) INTEGER
C             The number of rows of matrix C.  P >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading L-by-N part of this array must
C             contain the state dynamics matrix A.
C             On exit, the leading L-by-N part of this array contains
C             the transformed matrix A*Z.
C
C     LDA     INTEGER
C             The leading dimension of array A.  LDA >= MAX(1,L).
C
C     E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
C             On entry, the leading L-by-N part of this array must
C             contain the descriptor matrix E.
C             On exit, the leading L-by-N part of this array contains
C             the transformed matrix E*Z in upper trapezoidal form,
C             i.e.
C
C                      ( E11 )
C                E*Z = (     ) ,  if L >= N ,
C                      (  R  )
C             or
C
C                E*Z = ( 0  R ),  if L < N ,
C
C             where R is an MIN(L,N)-by-MIN(L,N) upper triangular
C             matrix.
C
C     LDE     INTEGER
C             The leading dimension of array E.  LDE >= MAX(1,L).
C
C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
C             On entry, the leading P-by-N part of this array must
C             contain the state/output matrix C.
C             On exit, the leading P-by-N part of this array contains
C             the transformed matrix C*Z.
C
C     LDC     INTEGER
C             The leading dimension of array C.  LDC >= MAX(1,P).
C
C     Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
C             If COMPZ = 'N':  Z is not referenced.
C             If COMPZ = 'I':  on entry, Z need not be set;
C                              on exit, the leading N-by-N part of this
C                              array contains the orthogonal matrix Z,
C                              which is the product of Householder
C                              transformations applied to A, E, and C
C                              on the right.
C             If COMPZ = 'U':  on entry, the leading N-by-N part of this
C                              array must contain an orthogonal matrix
C                              Z1;
C                              on exit, the leading N-by-N part of this
C                              array contains the orthogonal matrix
C                              Z1*Z.
C
C     LDZ     INTEGER
C             The leading dimension of array Z.
C             LDZ >= 1,        if COMPZ = 'N';
C             LDZ >= MAX(1,N), if COMPZ = 'U' or 'I'.
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= MAX(1, MIN(L,N) + MAX(L,N,P)).
C             For optimum performance
C             LWORK >= MAX(1, MIN(L,N) + MAX(L,N,P)*NB),
C             where NB is the optimal blocksize.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The routine computes the RQ factorization of E to reduce it
C     the upper trapezoidal form.
C
C     The transformations are also applied to the rest of system
C     matrices
C
C         A <- A * Z,  C <- C * Z.
C
C     NUMERICAL ASPECTS
C
C     The algorithm is numerically backward stable and requires
C     0( L*N*N )  floating point operations.
C
C     CONTRIBUTOR
C
C     A. Varga, German Aerospace Center, DLR Oberpfaffenhofen.
C     March 1999. Based on the RASP routine RPDSRQ.
C
C     REVISIONS
C
C     July 1999, V. Sima, Research Institute for Informatics, Bucharest.
C
C     KEYWORDS
C
C     Descriptor system, matrix algebra, matrix operations,
C     orthogonal transformation.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
C     .. Scalar Arguments ..
      CHARACTER          COMPZ
      INTEGER            INFO, L, LDA, LDC, LDE, LDWORK, LDZ, N, P
C     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), DWORK( * ),
     $                   E( LDE, * ), Z( LDZ, * )
C     .. Local Scalars ..
      LOGICAL            ILZ
      INTEGER            ICOMPZ, LN, WRKOPT
C     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
C     .. External Subroutines ..
      EXTERNAL           DGERQF, DLASET, DORMRQ, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          INT, MAX, MIN
C
C     .. Executable Statements ..
C
C     Decode COMPZ.
C
      IF( LSAME( COMPZ, 'N' ) ) THEN
         ILZ = .FALSE.
         ICOMPZ = 1
      ELSE IF( LSAME( COMPZ, 'U' ) ) THEN
         ILZ = .TRUE.
         ICOMPZ = 2
      ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
         ILZ = .TRUE.
         ICOMPZ = 3
      ELSE
         ICOMPZ = 0
      END IF
C
C     Test the input parameters.
C
      INFO = 0
      WRKOPT = MAX( 1, MIN( L, N ) + MAX( L, N, P ) )
      IF( ICOMPZ.EQ.0 ) THEN
         INFO = -1
      ELSE IF( L.LT.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( P.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, L ) ) THEN
         INFO = -6
      ELSE IF( LDE.LT.MAX( 1, L ) ) THEN
         INFO = -8
      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
         INFO = -10
      ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
         INFO = -12
      ELSE IF( LDWORK.LT.WRKOPT ) THEN
         INFO = -14
      END IF
      IF( INFO .NE. 0 ) THEN
         CALL XERBLA( 'TG01DD', -INFO )
         RETURN
      END IF
C
C     Initialize Q if necessary.
C
      IF( ICOMPZ.EQ.3 )
     $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
C
C     Quick return if possible.
C
      IF( L.EQ.0 .OR. N.EQ.0 ) THEN
         DWORK( 1 ) = ONE
         RETURN
      END IF
C
      LN = MIN( L, N )
C
C     Compute the RQ decomposition of E, E = R*Z.
C
C     Workspace: need   MIN(L,N) + L;
C                prefer MIN(L,N) + L*NB.
C
      CALL DGERQF( L, N, E, LDE, DWORK, DWORK( LN+1 ), LDWORK-LN, INFO )
      WRKOPT = MAX( WRKOPT, INT( DWORK( LN+1 ) ) + LN )
C
C     Apply transformation on the rest of matrices.
C
C     A <--  A * Z'.
C     Workspace: need   MIN(L,N) + L;
C                prefer MIN(L,N) + L*NB.
C
      CALL DORMRQ( 'Right', 'Transpose', L, N, LN, E( L-LN+1,1 ), LDE,
     $              DWORK, A, LDA, DWORK( LN+1 ), LDWORK-LN, INFO )
      WRKOPT = MAX( WRKOPT, INT( DWORK( LN+1 ) ) + LN )
C
C     C <-- C * Z'.
C     Workspace: need   MIN(L,N) + P;
C                prefer MIN(L,N) + P*NB.
C
      CALL DORMRQ( 'Right', 'Transpose', P, N, LN, E( L-LN+1,1 ), LDE,
     $              DWORK, C, LDC, DWORK( LN+1 ), LDWORK-LN, INFO )
      WRKOPT = MAX( WRKOPT, INT( DWORK( LN+1 ) ) + LN )
C
C     Z <-- Z1 * Z'.
C     Workspace: need   MIN(L,N) + N;
C                prefer MIN(L,N) + N*NB.
C
      IF( ILZ ) THEN
         CALL DORMRQ( 'Right', 'Transpose', N, N, LN, E( L-LN+1,1 ),
     $                LDE, DWORK, Z, LDZ, DWORK( LN+1 ), LDWORK-LN,
     $                INFO )
         WRKOPT = MAX( WRKOPT, INT( DWORK( LN+1 ) ) + LN )
      END IF
C
C     Set lower triangle of E to zero.
C
      IF( L.LT.N ) THEN
         CALL DLASET( 'Full', L, N-L, ZERO, ZERO, E, LDE )
         IF( L.GE.2 ) CALL DLASET( 'Lower', L-1, L, ZERO, ZERO,
     $                             E( 2, N-L+1 ), LDE )
      ELSE
         IF( N.GE.2 ) CALL DLASET( 'Lower', N-1, N, ZERO, ZERO,
     $                             E( L-N+2, 1 ), LDE )
      END IF
C
      DWORK(1) = WRKOPT
C
      RETURN
C *** Last line of TG01DD ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
      SUBROUTINE POWSQ(XPARAM, NVAR, FUNCT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION XPARAM(*)
      COMMON /MESAGE/ IFLEPO,ISCF
**********************************************************************
*
*   POWSQ OPTIMIZES THE GEOMETRY BY MINIMISING THE GRADIENT NORM.
*         THUS BOTH GROUND AND TRANSITION STATE GEOMETRIES CAN BE
*         CALCULATED. IT IS ROUGHLY EQUIVALENT TO FLEPO, FLEPO MINIMIZES
*         THE ENERGY, POWSQ MINIMIZES THE GRADIENT NORM.
*
*  ON ENTRY XPARAM = VALUES OF PARAMETERS TO BE OPTIMIZED.
*           NVAR   = NUMBER OF PARAMETERS TO BE OPTIMIZED.
*
*  ON EXIT  XPARAM = OPTIMIZED PARAMETERS.
*           FUNCT  = HEAT OF FORMATION IN KCALS.
*
**********************************************************************
C        *****  ROUTINE PERFORMS  A LEAST SQUARES MINIMIZATION  *****
C        *****  OF A FUNCTION WHICH IS A SUM OF SQUARES.        *****
C        *****  INITIALLY WRITTEN BY J.W. MCIVER JR. AT SUNY/   *****
C        *****  BUFFALO, SUMMER 1971.  REWRITTEN AND MODIFIED   *****
C        *****  BY A.K. AT SUNY BUFFALO AND THE UNIVERSITY OF   *****
C        *****  TEXAS.  DECEMBER 1973                           *****
C
      COMMON /GEOVAR/ NDUM,LOC(2,MAXPAR), IDUMY, XARAM(MAXPAR)
      COMMON /GEOM  / GEO(3,NUMATM)
      COMMON /LAST  / LAST
      COMMON /KEYWRD/ KEYWRD
      COMMON /TIME  / TIME0
      COMMON /NUMSCF/ NSCF
      COMMON /GEOSYM/ NDEP, LOCPAR(MAXPAR), IDEPFN(MAXPAR),
     1                 LOCDEP(MAXPAR)
      COMMON /GRADNT/ GRAD(MAXPAR),GNFINA
      COMMON /TIMDMP/ TLEFT, TDUMP
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
     1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /NUMCAL/ NUMCAL
      COMMON /SIGMA1/ GNEXT, AMIN, ANEXT
      COMMON /SIGMA2/ GNEXT1(MAXPAR), GMIN1(MAXPAR)
      COMMON /NLLCOM/ HESS(MAXPAR,MAXPAR),BMAT(MAXPAR,MAXPAR),
     1PMAT(MAXPAR*MAXPAR)
      COMMON /SCRACH/ PVEC
      DIMENSION IPOW(9), SIG(MAXPAR),
     1          E1(MAXPAR), E2(MAXPAR),
     2          P(MAXPAR), WORK(MAXPAR),
     3          PVEC(MAXPAR*MAXPAR), EIG(MAXPAR), Q(MAXPAR)
      LOGICAL DEBUG, RESTRT, TIMES, OKF, SCF1, RESFIL, LOG
      CHARACTER*241 KEYWRD
      SAVE ICALCN
      DATA  ICALCN /0/
      IF(ICALCN.NE.NUMCAL) THEN
         ICALCN=NUMCAL
         RESTRT=(INDEX(KEYWRD,'RESTART') .NE. 0)
         LOG=(INDEX(KEYWRD,'NOLOG') .EQ. 0)
         SCF1=(INDEX(KEYWRD,'1SCF') .NE. 0)
         TIME1=SECOND()
         TIME2=TIME1
         ICYC=0
         TIMES=(INDEX(KEYWRD,'TIME') .NE. 0)
         TLAST=TLEFT
         RESFIL=.FALSE.
         LAST=0
         ILOOP=1
         XINC=0.00529167D0
         RHO2=1.D-4
         TOL2=4.D-1
         IF(INDEX(KEYWRD,'PREC') .NE. 0) TOL2=1.D-2
         IF(INDEX(KEYWRD,'GNORM') .NE. 0) THEN
            TOL2=READA(KEYWRD,INDEX(KEYWRD,'GNORM'))
            IF(TOL2.LT.0.01.AND.INDEX(KEYWRD,' LET').EQ.0)THEN
               WRITE(6,'(/,A)')'  GNORM HAS BEEN SET TOO LOW, RESET TO 0
     1.01'
               TOL2=0.01D0
            ENDIF
         ENDIF
         DEBUG = (INDEX(KEYWRD,'POWSQ') .NE. 0)
         IF(RESTRT) THEN
C
C   RESTORE STORED DATA
C
            IPOW(9)=0
            CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,ILOOP,BMAT,IPOW)
            IF(SCF1) GOTO 390
            NSCF=IPOW(8)
            DO 10 I=1,NVAR
               GRAD(I)=GMIN1(I)
   10       GNEXT1(I)=GMIN1(I)
            WRITE(6,'('' XPARAM'',6F10.6)')(XPARAM(I),I=1,NVAR)
            IF(ILOOP .GT. 0) THEN
C#               ILOOP=ILOOP+1
               WRITE(6,'(//10X,'' RESTARTING AT POINT'',I3)')ILOOP
            ELSE
               WRITE(6,'(//10X,''RESTARTING IN OPTIMISATION'',
     1         '' ROUTINES'')')
            ENDIF
         ENDIF
*
*   DEFINITIONS:   NVAR   = NUMBER OF GEOMETRIC VARIABLES = 3*NUMAT-6
*
      ENDIF
      NVAR=ABS(NVAR)
      IF(DEBUG) THEN
         WRITE(6,'('' XPARAM'')')
         WRITE(6,'(5(2I3,F10.4))')(LOC(1,I),LOC(2,I),XPARAM(I),I=1,NVAR)
      ENDIF
      IF( .NOT. RESTRT) THEN
         DO 20 I=1,NVAR
   20    GRAD(I)=0.D0
         CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
      ENDIF
      IF(DEBUG) THEN
         WRITE(6,'('' STARTING GRADIENTS'')')
         WRITE(6,'(3X,8F9.4)')(GRAD(I),I=1,NVAR)
      ENDIF
      GMIN=SQRT(DOT(GRAD,GRAD,NVAR))
      DO 30 I=1,NVAR
         GNEXT1(I)=GRAD(I)
         GMIN1(I)=GNEXT1(I)
   30 CONTINUE
C
C    NOW TO CALCULATE THE HESSIAN MATRIX.
C
      IF(ILOOP.LT.0) GOTO 140
C
C   CHECK THAT HESSIAN HAS NOT ALREADY BEEN CALCULATED.
C
      ILPR=ILOOP
      DO 50 ILOOP=ILPR,NVAR
         TIME1=SECOND()
         XPARAM(ILOOP)=XPARAM(ILOOP) + XINC
         CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
         IF(SCF1) GOTO 390
         IF(DEBUG)WRITE(6,'(I3,12(8F9.4,/3X))')
     1    ILOOP,(GRAD(IF),IF=1,NVAR)
         GRAD(ILOOP)=GRAD(ILOOP)+1.D-5
         XPARAM(ILOOP)=XPARAM(ILOOP) - XINC
         DO 40 J=1,NVAR
   40    HESS(ILOOP,J)=-(GRAD(J)-GNEXT1(J))/XINC
         TIME2=SECOND()
         TSTEP=TIME2-TIME1
         IF(TIMES)WRITE(6,'('' TIME FOR STEP:'',F8.2,'' LEFT'',F8.2)')
     1    TSTEP, TLEFT
         IF(TLAST-TLEFT.GT.TDUMP)THEN
            TLAST=TLEFT
            RESFIL=.TRUE.
            IPOW(9)=2
            I=ILOOP
            IPOW(8)=NSCF
            CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW)
         ENDIF
         IF( TLEFT .LT. TSTEP*2.D0) THEN
C
C  STORE RESULTS TO DATE.
C
            IPOW(9)=1
            I=ILOOP
            IPOW(8)=NSCF
            CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW)
            STOP
         ENDIF
   50 CONTINUE
C        *****  SCALE -HESSIAN- MATRIX                           *****
      IF( DEBUG) THEN
         WRITE(6,'(//10X,''UN-NORMALIZED HESSIAN MATRIX'')')
         DO 60 I=1,NVAR
   60    WRITE(6,'(8F10.4)')(HESS(J,I),J=1,NVAR)
      ENDIF
      DO 80 I=1,NVAR
         SUM = 0.0D0
         DO 70 J=1,NVAR
   70    SUM = SUM+HESS(I,J)**2
   80 WORK(I) = 1.0D0/SQRT(SUM)
      DO 90 I=1,NVAR
         DO 90 J=1,NVAR
   90 HESS(I,J) = HESS(I,J)*WORK(I)
      IF( DEBUG) THEN
         WRITE(6,'(//10X,''HESSIAN MATRIX'')')
         DO 100 I=1,NVAR
  100    WRITE(6,'(8F10.4)')(HESS(J,I),J=1,NVAR)
      ENDIF
C        *****  INITIALIZE B MATIRX                        *****
      DO 120 I=1,NVAR
         DO 110 J=1,NVAR
  110    BMAT(I,J) = 0.0D0
  120 BMAT(I,I) = WORK(I)*2.D0
************************************************************************
*
*  THIS IS THE START OF THE BIG LOOP TO OPTIMIZE THE GEOMETRY
*
************************************************************************
      ILOOP=-99
      TSTEP=TSTEP*4
  130 CONTINUE
      IF(TLAST-TLEFT.GT.TDUMP)THEN
         TLAST=TLEFT
         RESFIL=.TRUE.
         IPOW(9)=2
         I=ILOOP
         IPOW(8)=NSCF
         CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW)
      ENDIF
      IF( TLEFT .LT. TSTEP*2.D0) THEN
C
C  STORE RESULTS TO DATE.
C
         IPOW(9)=1
         I=ILOOP
         IPOW(8)=NSCF
         CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW)
         IFLEPO=-1
         RETURN
      ENDIF
  140 CONTINUE
C        *****  FORM-A- DAGGER-A- IN PA SLONG WITH -P-     *****
      IJ=0
      DO 160 J=1,NVAR
         DO 160 I=1,J
            IJ=IJ+1
            SUM = 0.0D0
            DO 150 K=1,NVAR
  150       SUM = SUM + HESS(I,K)*HESS(J,K)
  160 PMAT(IJ) = SUM
      DO 180 I=1,NVAR
         SUM = 0.0D0
         DO 170 K=1,NVAR
  170    SUM = SUM-HESS(I,K)*GMIN1(K)
  180 P(I) = -SUM
      L=0
      IF(DEBUG) THEN
         WRITE(6,'(/10X,''P MATRIX IN POWSQ'')')
         CALL VECPRT(PMAT,NVAR)
      ENDIF
      CALL RSP(PMAT,NVAR,NVAR,EIG,PVEC)
C        *****  CHECK FOR ZERO EIGENVALUE                  *****
C#      WRITE(6,'(''  EIGS IN POWSQ:'')')
C#      WRITE(6,'(6F13.8)')(EIG(I),I=1,NVAR)
      IF(EIG(1).LT.RHO2) GO TO 240
C        *****  IF MATRIX IS NOT SINGULAR FORM INVERSE     *****
C        *****  BY BACK TRANSFORMING THE EIGENVECTORS      *****
      IJ=0
      DO 200 I=1,NVAR
         DO 200 J=1,I
            IJ=IJ+1
            SUM = 0.0D0
            DO 190 K=1,NVAR
  190       SUM = SUM+PVEC((K-1)*NVAR+J)*PVEC((K-1)*NVAR+I)/EIG(K)
  200 PMAT(IJ) = SUM
C        *****  FIND -Q- VECTOR                            *****
      L=0
      IL=L+1
      L=IL+I-1
      DO 230 I=1,NVAR
         SUM = 0.0D0
         DO 210 K=1,I
            IK=(I*(I-1))/2+K
  210    SUM = SUM+PMAT(IK)*P(K)
         IP1=I+1
         DO 220 K=IP1,NVAR
            IK=(K*(K-1))/2+I
  220    SUM=SUM+PMAT(IK)*P(K)
  230 Q(I) = SUM
      GO TO 260
  240 CONTINUE
C        *****  TAKE  -Q- VECTOR AS EIGENVECTOR OF ZERO     *****
C        *****  EIGENVALUE                                 *****
      DO 250 I=1,NVAR
  250 Q(I) = PVEC(I)
  260 CONTINUE
C        *****  FIND SEARCH DIRECTION                      *****
      DO 270 I=1,NVAR
         SIG(I) = 0.0D0
         DO 270 J=1,NVAR
  270 SIG(I) = SIG(I) + Q(J)*BMAT(I,J)
C        *****  DO A ONE DIMENSIONAL SEARCH                *****
      IF (DEBUG) THEN
         WRITE(6,'('' SEARCH VECTOR'')')
         WRITE(6,'(8F10.5)')(SIG(I),I=1,NVAR)
      ENDIF
      CALL SEARCH(XPARAM, ALPHA, SIG, NVAR, GMIN, OKF, FUNCT)
      IF( NVAR .EQ. 1) GOTO 390
C
C  FIRST WE ATTEMPT TO OPTIMIZE GEOMETRY USING SEARCH.
C  IF THIS DOES NOT WORK, THEN SWITCH TO LINMIN, WHICH ALWAYS WORKS,
C  BUT IS TWICE AS SLOW AS SEARCH.
C
      RMX = 0.0D0
      DO 280 K=1,NVAR
         RT = ABS(GMIN1(K))
         IF(RT.GT.RMX)RMX = RT
  280 CONTINUE
      IF(RMX.LT.TOL2) GO TO 390
C        *****  TWO STEP ESTIMATION OF DERIVATIVES         *****
      DO 290 K=1,NVAR
  290 E1(K) = (GMIN1(K)-GNEXT1(K))/(AMIN-ANEXT)
      RMU = DOT(E1,GMIN1,NVAR)/DOT(GMIN1,GMIN1,NVAR)
      DO 300 K=1,NVAR
  300 E2(K) = E1(K) - RMU*GMIN1(K)
C        *****  SCALE -E2- AND -SIG-                       *****
      SK = 1.0D0/SQRT(DOT(E2,E2,NVAR))
      DO 310 K=1,NVAR
  310 SIG(K) = SK*SIG(K)
      DO 320 K=1,NVAR
  320 E2(K) = SK*E2(K)
C        *****  FIND INDEX OF REPLACEMENT DIRECTION        *****
      PMAX = -1.0D+20
      DO 330 I=1,NVAR
         IF(ABS(P(I)*Q(I)).LE.PMAX) GO TO 330
         PMAX = ABS(P(I)*Q(I))
         ID = I
  330 CONTINUE
C        *****  REPLACE APPROPRIATE DIRECTION AND DERIVATIVE ***
      DO 340 K=1,NVAR
  340 HESS(ID,K) = -E2(K)
C        *****  REPLACE STARTING POINT                     *****
      DO 350 K=1,NVAR
  350 BMAT(K,ID) = SIG(K)/0.529167D0
      DO 360 K=1,NVAR
  360 GNEXT1(K) = GMIN1(K)
      TIME1=TIME2
      TIME2=SECOND()
      TLEFT=TLEFT-TIME2+TIME0
      TSTEP=TIME2-TIME1
      ICYC=ICYC+1
      IF(RESFIL)THEN
         WRITE(6,370)MIN(TLEFT,9999999.9D0),
     1MIN(GMIN,999999.999D0),FUNCT
         IF(LOG)WRITE(11,370)MIN(TLEFT,9999999.9D0),
     1MIN(GMIN,999999.999D0),FUNCT
  370    FORMAT('  RESTART FILE WRITTEN,  TIME LEFT:',F9.1,
     1' GRAD.:',F10.3,' HEAT:',G14.7)
         RESFIL=.FALSE.
      ELSE
         WRITE(6,380)ICYC,MIN(TSTEP,9999.99D0),
     1MIN(TLEFT,9999999.9D0),MIN(GMIN,999999.999D0),FUNCT
         IF(LOG)WRITE(11,380)ICYC,MIN(TSTEP,9999.99D0),
     1MIN(TLEFT,9999999.9D0),MIN(GMIN,999999.999D0),FUNCT
  380    FORMAT(' CYCLE:',I5,' TIME:',F6.1,' TIME LEFT:',F9.1,
     1' GRAD.:',F10.3,' HEAT:',G14.7)
      ENDIF
      IF(TIMES)WRITE(6,'('' TIME FOR STEP:'',F8.2,'' LEFT'',F8.2)')
     1TSTEP, TLEFT
      GO TO 130
  390 CONTINUE
      DO 400 I=1,NVAR
  400 GRAD(I)=0.D0
      LAST=1
      CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
      DO 410 I=1,NVAR
  410 GRAD(I)=GMIN1(I)
      GNFINA=SQRT(DOT(GRAD,GRAD,NVAR))
      IFLEPO=11
      IF(SCF1)IFLEPO=13
      RETURN
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive