text:
!
! Copyright (C) 2001-2011 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!----------------------------------------------------------------------------
SUBROUTINE electrons()
  !----------------------------------------------------------------------------
  !
  ! ... This routine is a driver of the self-consistent cycle.
  ! ... It uses the routine c_bands for computing the bands at fixed
  ! ... Hamiltonian, the routine sum_band to compute the charge density,
  ! ... the routine v_of_rho to compute the new potential and the routine
  ! ... mix_rho to mix input and output charge densities.
  ! ... It prints on output the total energy and its decomposition in
  ! ... the separate contributions.
  !
  USE kinds,                ONLY : DP
  USE constants,            ONLY : eps8, pi
  USE io_global,            ONLY : stdout, ionode
  USE cell_base,            ONLY : at, bg, alat, omega, tpiba2
  USE ions_base,            ONLY : zv, nat, nsp, ityp, tau, compute_eextfor
  USE basis,                ONLY : starting_pot, starting_wfc
  USE fft_base,             ONLY : dfftp
  USE gvect,                ONLY : ngm, gstart, nl, nlm, g, gg, gcutm
  USE gvecs,                ONLY : doublegrid, ngms
  USE klist,                ONLY : xk, wk, nelec, ngk, nks, nkstot, lgauss
  USE lsda_mod,             ONLY : lsda, nspin, magtot, absmag, isk
  USE vlocal,               ONLY : strf
  USE wvfct,                ONLY : nbnd, et, npwx, ecutwfc
  USE ener,                 ONLY : etot, hwf_energy, eband, deband, ehart, &
                                   vtxc, etxc, etxcc, ewld, demet, epaw
  USE scf,                  ONLY : scf_type, scf_type_COPY, &
                                   create_scf_type, destroy_scf_type, &
                                   rho, rho_core, rhog_core, &
                                   v, vltot, vrs, kedtau, vnew
  USE control_flags,        ONLY : mixing_beta, tr2, ethr, niter, nmix, &
                                   iprint, istep, lscf, lmd, conv_elec, &
                                   restart, io_level, do_makov_payne,  &
                                   iverbosity, textfor
  USE io_files,             ONLY : iunwfc, iunocc, nwordwfc, output_drho, &
                                   iunefield, iunpaw
  USE buffers,              ONLY : save_buffer
  USE extfield,             ONLY : tefield, etotefield
  USE wavefunctions_module, ONLY : evc, psic
  USE spin_orb,             ONLY : domag
  USE control_flags,        ONLY : adapt_thr, tr2_init, tr2_multi
  USE mp_global,            ONLY : intra_bgrp_comm, nbgrp, mpime, &
                                   inter_bgrp_comm, my_bgrp_id
  USE mp,                   ONLY : mp_sum
  !
  !
  USE uspp_param,           ONLY : nh, nhm ! used for PAW
  USE dfunct,                 only : newd
  !
  !
  IMPLICIT NONE
  !
  ! ... a few local variables
  !
  real(dp) :: v_kin_r !subsitute dummy variable for v%kin_r
  real(dp) :: eth     !usbsitutue dummy variable for ldaU%eth
  REAL(DP) :: &
      dr2,          &! the norm of the diffence between potential
      charge,       &! the total charge
      deband_hwf,   &! deband for the Harris-Weinert-Foulkes functional
      mag           ! local magnetization
  INTEGER :: &
      i,            &! counter on polarization
      idum,         &! dummy counter on iterations
      iter,         &! counter on iterations
      ik_,          &! used to read ik from restart file
      kilobytes
  REAL(DP) :: &
      tr2_min,     &! estimated error on energy coming from diagonalization
      descf,       &! correction for variational energy
      en_el=0.0_DP,&! electric field contribution to the total energy
      eext=0.0_DP   ! external forces contribution to the total energy
  LOGICAL :: &
      first
  !
  ! ... auxiliary variables for calculating and storing temporary copies of
  ! ... the charge density and of the HXC-potential
  !
  type (scf_type), save :: rhoin ! used to store rho_in of current/next iteration
  !
  ! ... external functions
  !
  REAL(DP), EXTERNAL :: ewald, get_clock
  REAL(DP) :: tr2_final ! final threshold for exx minimization
                        ! when using adaptive thresholds.
  iter = 0
  ik_  = 0
  dr2  = 0.0_dp
  tr2_final = tr2
  !
  WRITE( stdout, 9000 ) get_clock( 'MiniDFT' )
  !
  CALL memstat( kilobytes )
  !
  IF ( kilobytes > 0 ) WRITE( stdout, 9001 ) kilobytes/1000.0
  !
  CALL flush_unit( stdout )
  !
  !
  CALL start_clock( 'electrons' )
  !
  !
  CALL flush_unit( stdout )
  !
  ! ... calculates the ewald contribution to total energy
  !
  ewld = ewald( alat, nat, nsp, ityp, zv, at, bg, tau, &
                omega, g, gg, ngm, gcutm, gstart, .false., strf )
  !
  call create_scf_type ( rhoin )
  !
10 CONTINUE
  !
  ! ... Convergence threshold for iterative diagonalization
  !
  ! ... for the first scf iteration of each ionic step (after the first),
  ! ... the threshold is fixed to a default value of 1.D-6
  !
  IF ( istep > 0 ) ethr = 1.D-6
  !
  WRITE( stdout, 9002 )
  !
  CALL flush_unit( stdout )

  !
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  !%%%%%%%%%%%%%%%%%%%%          iterate !          %%%%%%%%%%%%%%%%%%%%%
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  !
  DO idum = 1, niter
     !
     IF ( check_stop_now() ) RETURN
     !
     iter = iter + 1
     !
     WRITE( stdout, 9010 ) iter, ecutwfc, mixing_beta
     !
     CALL flush_unit( stdout )
     !
     ! ... Convergence threshold for iterative diagonalization is
     ! ... automatically updated during self consistency
     !
     IF ( iter > 1 .AND. ik_ == 0 ) THEN
        !
        IF ( iter == 2 ) ethr = 1.D-2
        ethr = MIN( ethr, 0.1D0*dr2 / MAX( 1.D0, nelec ) )
        ! ... do not allow convergence threshold to become too small:
        ! ... iterative diagonalization may become unstable
        ethr = MAX( ethr, 1.D-13 )
        !
     END IF
     !
     first = ( iter == 1 )
     !
     ! ... deband = - \sum_v <\psi_v | V_h + V_xc |\psi_v> is calculated a
     ! ... first time here using the input density and potential ( to be
     ! ... used to calculate the Harris-Weinert-Foulkes energy )
     !
     deband_hwf = delta_e()
     !
     ! save input current density in rhoin
     call scf_type_COPY( rho, rhoin )
     !
     scf_step: DO
        !
        ! ... tr2_min is set to an estimate of the error on the energy
        ! ... due to diagonalization - used only for the first scf iteration
        !
        tr2_min = 0.D0
        !
        IF ( first ) tr2_min = ethr*MAX( 1.D0, nelec )
        !
        ! ... diagonalization of the KS hamiltonian
        !
           !
           CALL c_bands( iter, ik_, dr2 )
           !
        !
        IF ( check_stop_now() ) RETURN
        !
        ! ... xk, wk, isk, et, wg are distributed across pools;
        ! ... the first node has a complete copy of xk, wk, isk,
        ! ... while eigenvalues et and weights wg must be
        ! ... explicitely collected to the first node
        ! ... this is done here for et, in sum_band for wg
        !
        CALL poolrecover( et, nbnd, nkstot, nks )
        !
        ! ... the new density is computed here
        CALL sum_band()
        ! PAW : sum_band computes new becsum (stored in uspp modules) and a
        ! subtly different copy in rho%bec (scf module)
        !
        ! ... the Harris-Weinert-Foulkes energy is computed here using only
        ! ... quantities obtained from the input density
        !
        hwf_energy = eband + deband_hwf + (etxc - etxcc) + ewld + ehart + demet
        !
        !
        ! ... calculate total and absolute magnetization
        !
        !
        ! ... eband  = \sum_v \epsilon_v    is calculated by sum_band
        ! ... deband = - \sum_v <\psi_v | V_h + V_xc |\psi_v>
        ! ... eband + deband = \sum_v <\psi_v | T + Vion |\psi_v>
        !
        deband = delta_e()
        !
        ! ... mix_rho mixes several quantities: rho in g-space, tauk (for meta-gga)
        ! ... ns (for lda+u) and becsum (for paw)
        !
        CALL mix_rho( rho, rhoin, mixing_beta, dr2, tr2_min, iter, nmix, conv_elec )
        !
        ! ... if convergence is achieved or if the self-consistency error
        ! ... (dr2) is smaller than the estimated error due to diagonalization
        ! ... (tr2_min), rhoin and rho are unchanged: rhoin contains the input
        ! ...  density and rho contains the output density
        ! ... In the other cases rhoin contains the mixed charge density
        ! ... (the new input density) while rho is unchanged
        !
        IF ( first .and. nat > 0) THEN
           !
           ! ... first scf iteration: check if the threshold on diagonalization
           ! ... (ethr) was small enough wrt the error in self-consistency (dr2)
           ! ... if not, perform a new diagonalization with reduced threshold
           !
           first = .FALSE.
           !
           IF ( dr2 < tr2_min ) THEN
              !
              WRITE( stdout, '(/,5X,"Threshold (ethr) on eigenvalues was ", &
                               &    "too large:",/,5X,                      &
                               & "Diagonalizing with lowered threshold",/)' )
              !
              ethr = 0.1D0*dr2 / MAX( 1.D0, nelec )
              !
              CYCLE scf_step
              !
           END IF
           !
        END IF
        !
        not_converged_electrons : &
        IF ( .NOT. conv_elec ) THEN
           ! ... no convergence yet: calculate new potential from mixed
           ! ... charge density (i.e. the new estimate)
           !
           CALL v_of_rho( rhoin, rho_core, rhog_core, &
                          ehart, etxc, vtxc, eth, etotefield, charge, v)
           !
           ! ... estimate correction needed to have variational energy:
           ! ... T + E_ion (eband + deband) are calculated in sum_band
           ! ... and delta_e using the output charge density rho;
           ! ... E_H (ehart) and E_xc (etxc) are calculated in v_of_rho
           ! ... above, using the mixed charge density rhoin%of_r.
           ! ... delta_escf corrects for this difference at first order
           !
           descf = delta_escf()
           !
           ! ... now copy the mixed charge density in R- and G-space in rho
           !
           CALL scf_type_COPY( rhoin, rho )
           !
           ! ... write the charge density to file
           ! ... also write ldaU ns coeffs and PAW becsum
           !call infomsg( "electrons.f90:275","skipping write_rho; not implemented")
           !CALL write_rho( rho, nspin )
           !
        ELSE not_converged_electrons
           !
           ! ... convergence reached:
           ! ... 1) the output HXC-potential is saved in vr
           ! ... 2) vnew contains V(out)-V(in) ( used to correct the forces ).
           !
           vnew%of_r(:,:) = v%of_r(:,:)
           !
           CALL v_of_rho( rho,rho_core,rhog_core, &
                          ehart, etxc, vtxc, eth, etotefield, charge, v)
           !

           !
           vnew%of_r(:,:) = v%of_r(:,:) - vnew%of_r(:,:)
           !
           ! ... note that rho is here the output, not mixed, charge density
           ! ... so correction for variational energy is no longer needed
           !
           descf = 0._dp
           !
        END IF not_converged_electrons
        !
        !
        ! ... if we didn't cycle before we can exit the do-loop
        !
        EXIT scf_step
        !
     END DO scf_step
     !
     !
     ! ... define the total local potential (external + scf)
     !
     CALL set_vrs( vrs, vltot, v%of_r, kedtau, v_kin_r, dfftp%nnr, nspin, doublegrid )
     !
     ! ... in the US case we have to recompute the self-consistent
     ! ... term in the nonlocal potential
     !
     ! ... PAW: newd contains PAW updates of NL coefficients
     !write(*,*) "electrons.f90:480 this call to newd() might be moot if .not.okpaw"
     CALL newd()
     !
     ! ... save converged wfc if they have not been written previously
     !
     !BMA: eliminiate filesystem activity for miniDFT
     !IF ( nks == 1 .AND. (io_level < 2) ) &
     !   CALL save_buffer ( evc, nwordwfc, iunwfc, nks )
     !
     ! ... calculate the polarization
     !
        en_el=0.d0
     !
     ! ... write recover file
     !
     CALL save_in_electrons( iter, dr2 )
     !
     !
     WRITE( stdout, 9000 ) get_clock( 'MiniDFT' )
     !
     IF ( conv_elec ) WRITE( stdout, 9101 )
     !
     IF ( conv_elec .OR. MOD( iter, iprint ) == 0 ) THEN
        !
        CALL print_ks_energies()
        !
     END IF
     !
     IF ( ABS( charge - nelec ) / charge > 1.D-7 ) THEN
        WRITE( stdout, 9050 ) charge, nelec
        IF ( ABS( charge - nelec ) / charge > 1.D-3 ) THEN
           IF (.not.lgauss) THEN
              CALL errore( 'electrons', 'charge is wrong: smearing is needed', 1 )
           ELSE
              CALL errore( 'electrons', 'charge is wrong', 1 )
           END IF
        END IF
     END IF
     !
     etot = eband + ( etxc - etxcc ) + ewld + ehart + deband + demet + descf +en_el
     IF( textfor ) THEN
        eext =  compute_eextfor()
        etot = etot + eext
     END IF
     !
     etot = etot
     hwf_energy = hwf_energy
     !
     IF ( tefield ) THEN
        etot = etot + etotefield
        hwf_energy = hwf_energy + etotefield
     END IF
     !
     !
     IF ( ( conv_elec .OR. MOD( iter, iprint ) == 0 ) .AND. .NOT. lmd ) THEN
        !
        IF ( dr2 > eps8 ) THEN
           WRITE( stdout, 9081 ) etot, hwf_energy, dr2
        ELSE
           WRITE( stdout, 9083 ) etot, hwf_energy, dr2
        END IF
        !
        WRITE( stdout, 9060 ) &
            ( eband + deband ), ehart, ( etxc - etxcc ), ewld
        !
        !
        IF ( textfor)             WRITE( stdout, &
            '(/5x,"Energy of the external Forces = ", F18.8)' ) eext
        IF ( tefield )            WRITE( stdout, 9061 ) etotefield
        IF ( ABS (descf) > eps8 ) WRITE( stdout, 9069 ) descf
        !
        ! ... With Fermi-Dirac population factor, etot is the electronic
        ! ... free energy F = E - TS , demet is the -TS contribution
        !
        IF ( lgauss ) WRITE( stdout, 9070 ) demet
        !
     ELSE IF ( conv_elec .AND. lmd ) THEN
        !
        IF ( dr2 > eps8 ) THEN
           WRITE( stdout, 9081 ) etot, hwf_energy, dr2
        ELSE
           WRITE( stdout, 9083 ) etot, hwf_energy, dr2
        END IF
        !
     ELSE
        !
        IF ( dr2 > eps8 ) THEN
           WRITE( stdout, 9080 ) etot, hwf_energy, dr2
        ELSE
           WRITE( stdout, 9082 ) etot, hwf_energy, dr2
        END IF
     END IF
     !
     !
     IF ( lsda ) WRITE( stdout, 9017 ) magtot, absmag
     !
     !
     !
     CALL flush_unit( stdout )
     !
     IF ( conv_elec ) THEN
        !
        !
        ! ... if system is charged add a Makov-Payne correction to the energy
        !
        !
        ! ... print out ESM potentials if desired
        !
        !
        WRITE( stdout, 9110 ) iter
        !
        ! ... jump to the end
        !
        IF ( output_drho /= ' ' ) CALL remove_atomic_rho()
        !
        CALL stop_clock( 'electrons' )
        !
        call destroy_scf_type ( rhoin )
        !
        RETURN
        !
     END IF
     !
     ! ... uncomment the following line if you wish to monitor the evolution
     ! ... of the force calculation during self-consistency
     !
     !CALL forces()
     !
  END DO
  !
  WRITE( stdout, 9101 )
  WRITE( stdout, 9120 ) iter
  !
  CALL flush_unit( stdout )
  !
  IF ( output_drho /= ' ' ) CALL remove_atomic_rho()
  !
  CALL stop_clock( 'electrons' )
  !
  RETURN
  !
  ! ... formats
  !
9000 FORMAT(/'     total cpu time spent up to now is ',F10.1,' secs' )
9001 FORMAT(/'     per-process dynamical memory: ',f7.1,' Mb' )
9002 FORMAT(/'     Self-consistent Calculation' )
9010 FORMAT(/'     iteration #',I3,'     ecut=', F9.2,' Ry',5X,'beta=',F4.2 )
9017 FORMAT(/'     total magnetization       =', F9.2,' Bohr mag/cell', &
            /'     absolute magnetization    =', F9.2,' Bohr mag/cell' )
9018 FORMAT(/'     total magnetization       =',3F9.2,' Bohr mag/cell' &
       &   ,/'     absolute magnetization    =', F9.2,' Bohr mag/cell' )
9050 FORMAT(/'     WARNING: integrated charge=',F15.8,', expected=',F15.8 )
9060 FORMAT(/'     The total energy is the sum of the following terms:',/,&
            /'     one-electron contribution =',F17.8,' Ry' &
            /'     hartree contribution      =',F17.8,' Ry' &
            /'     xc contribution           =',F17.8,' Ry' &
            /'     ewald contribution        =',F17.8,' Ry' )
9061 FORMAT( '     electric field correction =',F17.8,' Ry' )
9062 FORMAT( '     - averaged Fock potential =',F17.8,' Ry' )
9064 FORMAT( '     + Fock energy             =',F17.8,' Ry' )
9065 FORMAT( '     Hubbard energy            =',F17.8,' Ry' )
9066 FORMAT( '     est. exchange err (dexx)  =',F17.8,' Ry' )
9069 FORMAT( '     scf correction            =',F17.8,' Ry' )
9070 FORMAT( '     smearing contrib. (-TS)   =',F17.8,' Ry' )
9071 FORMAT( '     Magnetic field            =',3F12.7,' Ry' )
9072 FORMAT( '     Magnetic field            =',F12.7, ' Ry' )
9073 FORMAT( '     lambda                    =',F11.2,' Ry' )
9074 FORMAT( '     Dispersion Correction     =',F17.8,' Ry' )
9080 FORMAT(/'     total energy              =',0PF17.8,' Ry' &
            /'     Harris-Foulkes estimate   =',0PF17.8,' Ry' &
            /'     estimated scf accuracy    <',0PF17.8,' Ry' )
9081 FORMAT(/'!    total energy              =',0PF17.8,' Ry' &
            /'     Harris-Foulkes estimate   =',0PF17.8,' Ry' &
            /'     estimated scf accuracy    <',0PF17.8,' Ry' )
9082 FORMAT(/'     total energy              =',0PF17.8,' Ry' &
            /'     Harris-Foulkes estimate   =',0PF17.8,' Ry' &
            /'     estimated scf accuracy    <',1PE17.1,' Ry' )
9083 FORMAT(/'!    total energy              =',0PF17.8,' Ry' &
            /'     Harris-Foulkes estimate   =',0PF17.8,' Ry' &
            /'     estimated scf accuracy    <',1PE17.1,' Ry' )
9085 FORMAT(/'     total all-electron energy =',0PF17.6,' Ry' )
9101 FORMAT(/'     End of self-consistent calculation' )
9110 FORMAT(/'     convergence has been achieved in ',i3,' iterations' )
9120 FORMAT(/'     convergence NOT achieved after ',i3,' iterations: stopping' )
9121 FORMAT(/'     scf convergence threshold =',1PE17.1,' Ry' )
  !
  CONTAINS
     !
     !-----------------------------------------------------------------------
     !
     !-----------------------------------------------------------------------
     FUNCTION check_stop_now()
       !-----------------------------------------------------------------------
       !
       USE check_stop,    ONLY : global_check_stop_now => check_stop_now
       !
       IMPLICIT NONE
       !
       LOGICAL :: check_stop_now
       INTEGER :: unit
       !
       unit = stdout
       !
       check_stop_now = global_check_stop_now( unit )
       !
       IF ( check_stop_now ) conv_elec = .FALSE.
       !
       RETURN
       !
     END FUNCTION check_stop_now
     !
     !-----------------------------------------------------------------------
     FUNCTION delta_e()
       !-----------------------------------------------------------------------
       ! ... delta_e = - \int rho%of_r(r)  v%of_r(r)
       !               - \int rho%kin_r(r) v%kin_r(r) [for Meta-GGA]
       !               - \sum rho%ns       v%ns       [for LDA+U]
       !               - \sum becsum       D1_Hxc     [for PAW]
       IMPLICIT NONE
       REAL(DP) :: delta_e, delta_e_hub
       !
       delta_e = - SUM( rho%of_r(:,:)*v%of_r(:,:) )
       !
       !
       delta_e = omega * delta_e / ( dfftp%nr1*dfftp%nr2*dfftp%nr3 )
       !
       CALL mp_sum( delta_e, intra_bgrp_comm )
       !
       !
       !
       RETURN
       !
     END FUNCTION delta_e
     !
     !-----------------------------------------------------------------------
     FUNCTION delta_escf()
       !-----------------------------------------------------------------------
       !
       ! ... delta_escf = - \int \delta rho%of_r(r)  v%of_r(r)
       !                  - \int \delta rho%kin_r(r) v%kin_r(r) [for Meta-GGA]
       !                  - \sum \delta rho%ns       v%ns       [for LDA+U]
       !                  - \sum \delta becsum       D1         [for PAW]
       ! ... calculates the difference between the Hartree and XC energy
       ! ... at first order in the charge density difference \delta rho(r)
       IMPLICIT NONE
       !
       REAL(DP) :: delta_escf, delta_escf_hub
       !
       delta_escf = - SUM( ( rhoin%of_r(:,:)-rho%of_r(:,:) )*v%of_r(:,:) )
       !
       !
       delta_escf = omega * delta_escf / ( dfftp%nr1*dfftp%nr2*dfftp%nr3 )
       !
       CALL mp_sum( delta_escf, intra_bgrp_comm )
       !

       RETURN
       !
     END FUNCTION delta_escf
     !
     !-----------------------------------------------------------------------
     !
END SUBROUTINE electrons

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
c===========================================================================
c
c   This file is part of TISEAN
c
c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
c
c   TISEAN 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   TISEAN 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 TISEAN; if not, write to the Free Software
c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
c
c===========================================================================
c   c2g.f
c   Gaussian kernel correlation integral from c2
c   author T. Schreiber (1998)
c===========================================================================
      parameter(meps=1000)
      dimension e(meps), c(meps), lw(meps)
      character*72 file, fout, aline
      double precision g,gd,h,d,f,func,gg,ggd,err,dum1,dum2,a,b,dc,de
      external func, funcd
      common h,d,f
      data iverb/1/

      call whatido("Gaussian kernel correlation sum from c2",iverb)
      isout=igetout(fout,iverb)
      if(nstrings().eq.0) call usage()
      call nthstring(1,file)
      call infile(file,iunit,iverb)
      if(isout.eq.1) call addsuff(fout,file,"_g")
      call outfile(fout,iunit2,iverb)
 1    read(iunit,'(a)',end=999) aline
 4    if(aline(1:1).ne."#") goto 1
      if(aline(1:1).eq."#")
     .   read(aline(index(aline,"m=")+2:72),'(i20)',err=1) m
      me=0
 2    read(iunit,'(a)') aline
      if(aline(1:72).eq." ") goto 3
      me=me+1
      read(aline,*,err=999,end=999) ee, cc
      if(cc.le.0.) goto 3
      e(me)=log(ee)
      c(me)=log(cc)
      goto 2
 3    write(iunit2,'(4h#m= ,i5)') m
      call indexx(me,e,lw)
      call index2sort(me,e,lw)
      call index2sort(me,c,lw)
      do 10 j=1,me
         h=exp(e(j))
         g=0
         gd=0
         do 20 k=1,me-1
            f=exp((e(k+1)*c(k)-e(k)*c(k+1))/(e(k+1)-e(k)))
            d=(c(k+1)-c(k))/(e(k+1)-e(k))
            a=e(k)
            b=e(k+1)
            gg=0.
            ggd=0.
            if(b.ne.a) call dqk15(func,a,b,gg,err,dum1,dum2)
            if(b.ne.a) call dqk15(funcd,a,b,ggd,err,dum1,dum2)
            g=g+gg
 20         gd=gd+ggd
         dc=c(me)
         de=e(me)
         cgauss=g/(h**2)+exp(-exp(2*de)/(2*h**2))
         cgd=gd/(h**4)+(2+exp(2*de)/h**2)*exp(-exp(2*de)/(2*h**2))
 10      write(iunit2,*) h, cgauss, -2+cgd/cgauss
      write(iunit2,'()')
      write(iunit2,'()')
      goto 4
 999  stop
      end

      double precision function func(u)
      double precision h,d,f,u
      common h,d,f

      func=f*exp((2+d)*u-exp(2*u)/(2*h**2))
      end

      double precision function funcd(u)
      double precision h,d,f,u
      common h,d,f

      funcd=f*exp((4+d)*u-exp(2*u)/(2*h**2))
      end

      subroutine usage()
c usage message

      call whatineed(
     .   "[-o outfile -V# -h] file")
      call pout("file_g")
      call pall()
      stop
      end

ext: f

max_stars_repo_name: homiak/jabaws

max_issues_repo_name: homiak/jabaws

max_forks_repo_name: homiak/jabaws


text:
      SUBROUTINE INPTT1
C
C     READ DATA BLOCK(S) FROM A NASTRAN USER TAPE WHICH MUST BE SET UP.
C
C     CALL TO THIS MODULE IS
C
C     INPUTT1  /O1,O2,O3,O4,O5/V,N,P1/V,N,P2/V,N,P3/V,N,P4  $
C
C     PARAMETERS P1 AND P2 ARE INTEGER INPUT, P3 AND P4 ARE BCD
C
C               P1= 0, NO ACTION TAKEN BEFORE READ (DEFAULT)
C                 =+N, SKIP FORWARD N DATA BLOCKS BEFORE READ
C                 =-1, USER TAPE IS REWOUND BEFORE READ
C                 =-2, A NEW REEL IS MOUNTED BEFORE READ
C                 =-3, THE NAMES OF ALL DATA BLOCKS ON USER TAPE ARE
C                      PRINTED AND READ OCCURS AT BEGINNING OF TAPE
C                 =-4, AN OUTPUT TAPE IS TO BE DISMOUNTED
C                      AFTER AN END-OF-FILE MARK IS WRITTEN.
C                      A NEW INPUT REEL WILL THEN BE MOUNTED.
C                 =-5, SEARCH USER TAPE FOR FIRST VERSION OF DATA
C                      BLOCKS REQUESTED.
C                      IF ANY ARE NOT FOUND, A FATAL TERMINATION
C                      OCCURS.
C                 =-6, SEARCH USER TAPE FOR FINAL VERSION OF DATA
C                      BLOCKS REQUESTED.
C                      IF ANY ARE NOT FOUND, A FATAL TERMINATION
C                      OCCURS.
C                 =-7, SEARCH USER TAPE FOR FIRST VERSION OF DATA
C                      BLOCKS REQUESTED.
C                      IF ANY ARE NOT FOUND, A WARNING OCCURS.
C                 =-8, SEARCH USER TAPE FOR FINAL VERSION OF DATA
C                      BLOCKS REQUESTED.
C                      IF ANY ARE NOT FOUND, A WARNING OCCURS.
C                 =-9, REWIND AND UNLOAD USER TAPE
C
C               P2= 0, FILE NAME IS INPT
C                 = 1, FILE NAME IS INP1
C                 = 2, FILE NAME IS INP2
C                 = 3, FILE NAME IS INP3
C                 = 4, FILE NAME IS INP4
C                 = 5, FILE NAME IS INP5
C                 = 6, FILE NAME IS INP6
C                 = 7, FILE NAME IS INP7
C                 = 8, FILE NAME IS INP8
C                 = 9, FILE NAME IS INP9
C                 THE MPL DEFAULT VALUE FOR P2 IS 0
C
C               P3=    TAPE ID CODE FOR USER TAPE, AN ALPHANUMERIC
C                      VARIABLE WHOSE VALUE MUST MATCH A CORRESPONDING
C                      VALUE ON THE USER TAPE.
C                      THIS CHECK IS DEPENDENT ON THE VALUE OF
C                      P1 AS FOLLOWS..
C                       *P1*             *TAPE ID CHECKED*
C                        +N                     NO
C                         0                     NO
C                        -1                    YES
C                        -2                    YES (ON NEW REEL)
C                        -3                    YES (WARNING CHECK)
C                        -4                    YES (ON NEW REEL)
C                        -5                    YES
C                        -6                    YES
C                        -7                    YES
C                        -8                    YES
C                        -9                     NO
C                      THE MPL DEFAULT VALUE FOR P3 IS XXXXXXXX
C
C
      EXTERNAL        RSHIFT,ANDF
      LOGICAL         TAPEUP,TAPBIT
      INTEGER         OUBUF,OUTPUT,P1,P2,P3,P4,ZERO,RSHIFT,ANDF,NONE(2),
     1                TRL(7),NAME(2),SUBNAM(2),INN(10),OUT(5),NAMEX(2),
     2                IDHDR(7),IDHDRX(7),P3X(2),NT(5,3),DX(3),TAPCOD(2),
     3                BCDBIN(4)
      CHARACTER       UFM*23,UWM*25,UIM*29,SFM*25
      COMMON /XMSSG / UFM,UWM,UIM,SFM
      COMMON /MACHIN/ MACH
      COMMON /BLANK / P1,P2,P3(2),P4(2)
     1       /SYSTEM/ KSYSTM(65)
     2       /ZZZZZZ/ X(1)
      EQUIVALENCE     (KSYSTM(1),NB  ), (KSYSTM( 2),NOUT),
     1                (KSYSTM(9),NLPP), (KSYSTM(12),LINE)
      DATA    SUBNAM/ 4HINPT, 4HT1  / , MSC      / 4HMSC /
      DATA    OUT   / 201,202,203,204,205/, MASK / 65535 /
      DATA    ZERO  , MONE,MTWO,MTRE,MFOR/ 0,-1,-2,-3,-4 /,
     1        MFIV  , MSIX,METE,MNIN     /-5,-6,-8,-9    /
      DATA    INN   / 4HINPT,4HINP1,4HINP2,4HINP3,4HINP4 ,
     1                4HINP5,4HINP6,4HINP7,4HINP8,4HINP9 /
      DATA    IDHDR / 4HNAST,4HRAN ,4HUSER,4H TAP,4HE ID,4H COD,4HE - /
      DATA    BCDBIN/ 4HBCD ,4H    ,4HBINA,4HRY          /
      DATA    NONE  / 4H (NO,4HNE) /, IPT1,IPT4/ 1H1,1H4 /
C
C
      IPTX = IPT1
      IF (P4(1) .EQ. MSC) GO TO 20
      GO TO 100
C
C
      ENTRY INPUT1
C     ============
C
C     INPUT1 HANDELS MSC/OUTPUT1 DATA.
C     INPUT1 IS CALLED FROM INPTT1, WITH P4 = 'MSC', OR IT IS CALLED
C     FROM INPTT4
C
   20 IPTX = IPT4
      IF (P3(1).EQ.BCDBIN(1) .AND. P3(2).EQ.BCDBIN(2)) GO TO 9918
      IF (P3(1).EQ.BCDBIN(3) .AND. P3(2).EQ.BCDBIN(4)) GO TO 9918
      WRITE  (NOUT,30) UIM
   30 FORMAT (A29,'. INPUTT1 IS REQUESTED TO READ INPUT TAPE GENERATED',
     1       ' IN MSC/OUTPUT1 COMPATIBLE RECORDS')
C
  100 LCOR = KORSZ(X) - 2*NB
      IF (LCOR .LE. 0) CALL MESAGE (-8,LCOR,SUBNAM)
      INBUF = LCOR  + 1
      OUBUF = INBUF + NB
      TAPCOD(1) = P3(1)
      TAPCOD(2) = P3(2)
      IF (P2.LT.0 .OR. P2.GT.9) GO TO 9907
      IN = INN(P2+1)
      IF (IPTX .EQ. IPT4) WRITE (NOUT,110) UIM,NB,IN
  110 FORMAT (A29,', CURRENT NASTRAN BUFFER SIZE IS',I9,' WORDS', /5X,
     1       'SYNCHRONIZED BUFFSIZE IS REQUIRED IN CURRENT NASTRAN AND',
     2       ' THE VERSION THAT WROTE ',A4,' TAPE (OR FILE)', /5X,
     3       3(4H====),/)
      IFILE  = IN
      IF (MACH .GE. 5) GO TO 120
      TAPEUP = TAPBIT(IN)
      IF (.NOT.TAPEUP ) GO TO 9909
  120 IF (P1 .LT. MNIN) GO TO 9908
C
      IF (P1 .EQ. MNIN) GO TO 5000
      IF (P1 .LT. MFOR) GO TO 3000
      IF (P1 .EQ. MTRE) GO TO 2000
      IF (P1 .LE. ZERO) GO TO 150
C
      CALL OPEN (*9901,IN,X(INBUF),2)
      DO 130 I = 1,P1
      CALL READ (*9906,*9906,IN,NAMEX,2,0,NF)
  130 CALL SKPFIL (IN,1)
      GO TO 250
C
  150 IF (P1.NE.MTWO .AND. P1.NE.MFOR) GO TO 190
C
C     P1 = -2 OR P1 = -4 IS ACCEPTABLE ONLY ON IBM OR UNIVAC
C
      IF (MACH.NE.2 .AND. MACH.NE.3) GO TO 9908
C
      IOLD = -P1/2
      CALL OPEN (*9901,IN,X(INBUF),2)
      CALL TPSWIT (IN,IOLD,1,TAPCOD)
C
  190 IF (P1.NE.MONE .AND. P1.NE.MTWO .AND. P1.NE.MFOR) GO TO 230
C
C     OPEN USER TAPE TO READ WITH REWIND AHD TAPE ID CHECK
C
      IF (P1.NE.MONE .AND. P1.NE.MTWO .AND. P1.NE.MFOR .AND.
     1    IPTX.EQ.IPT4) GO TO 230
      CALL OPEN (*9901,IN,X(INBUF),0)
      CALL READ (*9911,*9912,IN,DX,3,0,NF)
      CALL READ (*9911,*9912,IN,IDHDRX,7,0,NF)
      DO 210 KF = 1,7
      IF (IDHDRX(KF) .NE. IDHDR(KF)) GO TO 9913
  210 CONTINUE
      CALL READ (*9911,*9912,IN,P3X,2,1,NF)
      IF (P3X(1).NE.P3(1) .OR. P3X(2).NE.P3(2)) GO TO 9910
      CALL SKPFIL (IN,1)
      GO TO 250
C
C     OPEN USER TAPE TO READ WITHOUT REWIND AND NO TAPE ID CHECK
C
  230 CALL OPEN (*9901,IN,X(INBUF),2)
      IF (IPTX .EQ. IPT4) CALL FWDREC (*9912,IX)
C
  250 DO 1000 I = 1,5
      OUTPUT = OUT(I)
      TRL(1) = OUTPUT
      CALL RDTRL (TRL)
      IF (TRL(1) .LE. 0) GO TO 1000
      CALL FNAME (OUTPUT,NAME)
      IF (NAME(1).EQ.NONE(1) .AND. NAME(2).EQ.NONE(2)) GO TO 1000
C
C     PASS FILE NAME HEADER RECORD
C
      CALL READ (*9904,*9905,IN,NAMEX,2,0,NF)
C
C     READ TRAILER RECORD, SIX WORDS (OR 3 WORDS, IPTX=4 ONLY)
C
      CALL READ (*9904,*300,IN,TRL(2),6,1,NF)
      GO TO 340
C
C     JUST A NOTE, FROM G.CHAN/UNISYS -
C     LEVEL 17.5 USED 2 RECORDS HERE FOR THE MATRIX NAME (2 BCD WORDS,
C     1ST RECORD) AND 7 TRAILER WORDS (2ND RECORD)
C
  300 IF (IPTX.NE.IPT4 .OR. NF.LT.3) GO TO 9905
      TRL(5) = TRL(2)
      TRL(6) = TRL(3)
      TRL(7) = TRL(4)
      DO 320 J = 2,7
      J1 = J/2 + 4
      J2 = MOD(J-1,2)*16
      TRL(J) = ANDF(RSHIFT(TRL(J1),J2),MASK)
  320 CONTINUE
C
C     OPEN OUTPUT DATA BLOCK TO WRITE WITH REWIND
C
  340 CALL OPEN (*9902,OUTPUT,X(OUBUF),1)
C
C     COPY CONTENTS OF USER TAPE ONTO OUTPUT DATA BLOCK, INCLUDING
C     FILE NAME IN RECORD 0
C
      CALL CPYFIL (IN,OUTPUT,X,LCOR,NF)
C
C     CLOSE OUTPUT DATA BLOCK WITH REWIND AND EOF
C
      CALL CLOSE (OUTPUT,1)
C
C     WRITE TRAILER
C
      TRL(1) = OUTPUT
      CALL WRTTRL (TRL)
      CALL PAGE2 (-3)
      WRITE  (NOUT,400) UIM,NAME,IN,NAMEX
  400 FORMAT (A29,' 4105,     DATA BLOCK ',2A4,' RETRIEVED FROM USER ',
     1       'TAPE',A4, /5X,'NAME OF DATA BLOCK WHEN PLACED ON USER ',
     2       'TAPE WAS ',2A4 )
C
 1000 CONTINUE
C
C     CLOSE NASTRAN USER TAPE WITHOUT REWIND
C
      CALL CLOSE (IN,2)
      RETURN
C
C     OBTAIN LIST OF DATA BLOCKS ON USER TAPE.
C
 2000 CALL OPEN (*9901,IN,X(INBUF),0)
      CALL READ (*9911,*9912,IN,DX,3,0,NF)
      CALL READ (*9911,*9912,IN,IDHDRX,7,0,NF)
      DO 2005 KF = 1,7
      IF (IDHDRX(KF) .NE. IDHDR(KF)) GO TO 9913
 2005 CONTINUE
      CALL READ (*9911,*9912,IN,P3X,2,1,NF)
      IF (P3X(1).NE.P3(1) .OR. P3X(2).NE.P3(2)) GO TO 9914
 2006 CALL SKPFIL (IN,1)
      KF = 0
 2007 CALL PAGE1
      LINE = LINE + 5
      WRITE  (NOUT,2010) IN
 2010 FORMAT (1H0,50X,A4,14H FILE CONTENTS ,/46X,4HFILE,18X,4HNAME/1H0)
 2020 CALL READ (*2050,*9915,IN,NAMEX,2,1,NF)
      CALL SKPFIL (IN,1)
      KF = KF + 1
      LINE = LINE + 1
      WRITE (NOUT,2030) KF,NAMEX
 2030 FORMAT (45X,I5,18X,2A4)
      IF (LINE - NLPP) 2020,2007,2007
 2050 CALL REWIND (IN)
      CALL SKPFIL (IN,1)
      GO TO 250
C
C
C     SEARCH MODE
C
 3000 CONTINUE
C
C     EXAMINE OUTPUT REQUESTS AND FILL NAME TABLE
C
      NNT = 0
      DO 3050 I = 1,5
      OUTPUT = OUT(I)
      TRL(1) = OUTPUT
      CALL RDTRL (TRL)
      IF (TRL(1) .LE. 0) GO TO 3020
      CALL FNAME (OUTPUT,NAME)
      IF (IPTX.EQ.IPT4 .AND. NAME(1).EQ.NONE(1) .AND. NAME(2).EQ.NONE(2)
     1   ) GO TO 3010
      NT(I,1) = 0
      NT(I,2) = NAME(1)
      NT(I,3) = NAME(2)
      NNT = NNT + 1
      GO TO 3050
 3010 NT(I,2) = NAME(1)
      NT(I,3) = NAME(2)
 3020 NT(I,1) = -1
 3050 CONTINUE
C
      IF (NNT .GT. 0) GO TO 3070
      CALL PAGE2 (-2)
      WRITE  (NOUT,3060) UWM,IPTX
 3060 FORMAT (A25,' 4137,  ALL OUTPUT DATA BLOCKS FOR INPUTT',A1,
     1       ' ARE PURGED.')
      RETURN
C
C     CHECK TAPE ID LABEL.
C
 3070 CALL OPEN (*9901,IN,X(INBUF),0)
      CALL READ (*9911,*9912,IN,DX,3,0,NF)
      CALL READ (*9911,*9912,IN,IDHDRX,7,0,NF)
      DO 3080 KF = 1,7
      IF (IDHDRX(KF) .NE. IDHDR(KF)) GO TO 9913
 3080 CONTINUE
      CALL READ (*9911,*9912,IN,P3X,2,1,NF)
      IF (P3X(1).NE.P3(1) .OR. P3X(2).NE.P3(2)) GO TO 9910
      CALL SKPFIL (IN,1)
C
C
C     BEGIN SEARCH OF TAPE.
C
      KF = 0
 3110 CALL READ (*3500,*9915,IN,NAMEX,2,0,NF)
      KF = KF + 1
C
      DO 3200 I = 1,5
      NAME(1) = NT(I,2)
      NAME(2) = NT(I,3)
      IF (NT(I,1) .LT. 0) GO TO 3200
      IF (NAME(1).NE.NAMEX(1) .OR. NAME(2).NE.NAMEX(2)) GO TO 3200
      NT(I,1) = NT(I,1) + 1
      IF (NT(I,1).EQ.1 .OR. P1.EQ.MSIX .OR. P1.EQ.METE) GO TO 3150
      CALL PAGE2 (-3)
      WRITE  (NOUT,3140) UWM,NAME,KF,IN
 3140 FORMAT (A25,' 4138,  DATA BLOCK ',2A4,' (DATA BLOCK COUNT =',I5,
     2        ') HAS PREVIOUSLY BEEN RETRIEVED FROM', /36X ,
     3        'USER TAPE ',A4,' AND WILL BE IGNORED.')
      GO TO 3205
 3150 CALL READ (*9904,*3160,IN,TRL(2),6,1,NF)
      GO TO 3180
 3160 IF (IPTX.NE.IPT4 .OR. NF.LT.3) GO TO 9905
      TRL(5) = TRL(2)
      TRL(6) = TRL(3)
      TRL(7) = TRL(4)
      DO 3170 J = 2,7
      J1 = J/2 + 4
      J2 = MOD(J-1,2)*16
      TRL(J) = ANDF(RSHIFT(TRL(J1),J2),MASK)
 3170 CONTINUE
 3180 OUTPUT = OUT(I)
      CALL OPEN (*9902,OUTPUT,X(OUBUF),1)
      CALL CPYFIL (IN,OUTPUT,X,LCOR,NF)
      CALL CLOSE (OUTPUT,1)
      TRL(1) = OUTPUT
      CALL WRTTRL (TRL)
      CALL PAGE2 (-2)
      WRITE  (NOUT,3185) UIM,NAME,IN,KF
 3185 FORMAT (A29,' 4139, DATA BLOCK ',2A4,' RETRIEVED FROM USER TAPE ',
     1       A4,' (DATA BLOCK COUNT =',I5,1H))
      IF (NT(I,1) .GT. 1) GO TO 3190
      NNT = NNT - 1
      GO TO 3210
 3190 WRITE  (NOUT,3195) UWM
 3195 FORMAT (A25,' 4140, SECONDARY VERSION OF DATA BLOCK HAS REPLACED',
     1        ' EARLIER ONE.')
      CALL PAGE2 (-2)
      GO TO 3210
 3200 CONTINUE
C
 3205 CALL SKPFIL (IN,1)
 3210 IF (NNT.GT.0 .OR. P1.EQ.MSIX .OR. P1.EQ.METE) GO TO 3110
      GO TO 3900
C
 3500 IF (NNT .LE. 0) GO TO 3900
      CALL PAGE2 (-7)
      IF (P1.EQ.MFIV .OR. P1.EQ.MSIX) GO TO 9916
      WRITE  (NOUT,3510) UWM
 3510 FORMAT (A25,' 4141, ONE OR MORE DATA BLOCKS NOT FOUND ON USER ',
     1        'TAPE.')
      DO 3530 I = 1,5
      IF (NT(I,1) .NE. 0) GO TO 3530
      WRITE  (NOUT,3520) NT(I,2),NT(I,3)
 3520 FORMAT (20X,21HNAME OF DATA BLOCK = ,2A4)
 3530 CONTINUE
      IF (P1.EQ.MFIV .OR. P1.EQ.MSIX) GO TO 9995
C
 3900 CONTINUE
      CALL SKPFIL (IN,-1)
      CALL CLOSE (IN,2)
      RETURN
C
 5000 CONTINUE
      CALL UNLOAD (IN)
      RETURN
C
C     ERRORS
C
 9901 WRITE  (NOUT,9951) SFM,IPTX,IN
 9951 FORMAT (A25,' 4107, MODULE INPTT',A1,' UNABLE TO OPEN NASTRAN ',
     1       'FILE ',A4,1H.)
      GO TO 9995
C
 9902 WRITE  (NOUT,9952) SFM,IPTX,OUTPUT
 9952 FORMAT (A25,' 4108, SUBROUTINE INPTT',A1,' UNABLE TO OPEN OUTPUT',
     1       ' DATA BLOCK',I5)
      GO TO 9995
C
 9904 CALL MESAGE (-2,IFILE,SUBNAM)
C
 9905 CALL MESAGE (-3,IFILE,SUBNAM)
C
 9906 WRITE  (NOUT,9956) UFM,IPTX,P1,IN,I
 9956 FORMAT (A22,' 4111, MODULE INPUTT',A1,' IS UNABLE TO SKIP FORWARD'
     1,       I10,' DATA BLOCKS ON PERMANENT NASTRAN FILE ',A4,1H., /5X,
     2        'NUMBER OF DATA BLOCKS SKIPPED =',I5)
      LINE = LINE + 1
      GO TO 9995
C
 9907 WRITE  (NOUT,9957) UFM,IPTX,P2
 9957 FORMAT (A23,' 4112, MODULE INPUTT',A1,' - ILLEGAL VALUE FOR ',
     1        'SECOND PARAMETER =',I20)
      GO TO 9995
C
 9908 WRITE  (NOUT,9958) UFM,IPTX,P1
 9958 FORMAT (A23,' 4113, MODULE INPUTT',A1,' - ILLEGAL VALUE FOR ',
     1       'FIRST PARAMETER =',I20)
      GO TO 9995
C
 9909 WRITE  (NOUT,9959) UFM,IN
 9959 FORMAT (A23,' 4127, USER TAPE ',A4,' NOT SET UP.')
      GO TO 9995
C
 9910 WRITE  (NOUT,9960) UFM,P3X,IPTX,P3
 9960 FORMAT (A23,' 4136, USER TAPE ID CODE -',2A4,'- DOES NOT MATCH ',
     1       'THIRD INPUTT',A1,' DMAP PARAMETER -',2A4,2H-.)
      GO TO 9995
C
 9911 WRITE  (NOUT,9961) UFM,IPTX,IN
 9961 FORMAT (A23,' 4132, MODULE INPUTT',A1,' - END-OF-FILE ENCOUNTERED'
     1,   ' WHILE ATTEMPTING TO READ TAPE ID CODE ON USER TAPE ',A4,1H.)
      GO TO 9995
C
 9912 WRITE  (NOUT,9962) UFM,IPTX,IN
 9962 FORMAT (A23,' 4133, MODULE INPUTT',A1,' - END-OF-RECORD ',
     1       'ENCOUNTERED WHILE ATTEMPTING TO READ TAPE ID CODE ON ',
     2       'USER TAPE ',A4,1H.)
      GO TO 9995
C
 9913 WRITE  (NOUT,9963) UFM,IPTX,IDHDRX
 9963 FORMAT (A23,' 4134, MODULE INPUTT',A1,
     1        ' - ILLEGAL TAPE CODE HEADER = ',7A4)
      GO TO 9995
C
 9914 WRITE  (NOUT,9964) UWM,P3X,P3
 9964 FORMAT (A25,' 4135, USER TAPE ID CODE -',2A4,'- DOES NOT MATCH ',
     1        'THIRD INPUTT1 DMAP PARAMETER -',2A4,2H-.)
      LINE = LINE + 2
      GO TO 2006
C
 9915 WRITE  (NOUT,9965) SFM,IPTX
 9965 FORMAT (A25,' 4106, MODULE INPUTT',A1,' - SHORT RECORD.')
      GO TO 9995
C
 9916 WRITE  (NOUT,9966) UFM
 9966 FORMAT (A23,' 4142, ONE OR MORE DATA BLOCKS NOT FOUND ON USER ',
     1       'TAPE',/)
      DO 9917 I = 1,5
      IF (NT(I,1) .NE. 0) GO TO 9917
      WRITE (NOUT,9967) NT(I,2),NT(I,3)
      LINE = LINE + 1
 9917 CONTINUE
 9967 FORMAT (20X,'NAME OF DATA BLOCK = ',2A4)
      GO TO 9995
C
 9918 WRITE  (NOUT,9968) UFM,P3
 9968 FORMAT (A23,', ILLEGAL TAPE LABEL NAME -',2A4,'-  POSSIBLY ',
     1       'THE 4TH PARAMETER OF INPTT4 IS IN ERROR')
C
C
 9995 LINE = LINE + 2
      CALL MESAGE (-61,0,0)
      RETURN
C
      END


text:
      function sigz2(r)

      include 'commonblocks'

      psizh=pot(r,zdisk)
      psi00=pot(r,0.0)
cc
      truesigz2=(psizh-psi00)/log(0.419974)
cc
      call splintd(rr(0),fszrat(0),szrat2(0),nrspl+1,r,fcor)
      sigz2=truesigz2*fcor
      if(sigz2.lt.0.) sigz2 = 1.e-10

      return
      end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!
! !MODULE: ENS_Cpl_ESMFMod --- Sub module of the ESMF coupler gridded component
!                              of the EARTH ensemble system.
!
! !DESCRIPTION: ENS_Cpl_ESMFMod --- Routines to run the superstructure of the
!                                   EARTH Coupler Grid Component.
! !REVISION HISTORY:
!
!  April 2006     Weiyu Yang Initial code.
!
! !INTERFACE:
!

 MODULE ENS_Cpl_ESMFMod

!!USES:
!------

 USE ENS_Cpl_Run_ESMFMod
 USE ENS_CplState_ESMFMod

 IMPLICIT none

 END MODULE ENS_Cpl_ESMFMod


text:
c
c-----------------------------------------------------------------------
c subroutine: fftsym
c compute dft for real, symmetric, n-point sequence x(m) using
c n/2-point fft
c symmetric sequence means x(m)=x(n-m), m=1,...,n/2-1
c note: index m is sequence index--not fortran index
c-----------------------------------------------------------------------
c
      subroutine fftsym(x, n, y)
      dimension x(1), y(1)
c
c x = real array which on input contains the n/2+1 points of the
c     input sequence (symmetrical)
c     on output x contains the n/2+1 real points of the transform of
c     the input--i.e. the zero valued imaginary parts are not returned
c n = true size of input
c y = scratch array of size n/2+2
c
c
c for n = 2, compute dft directly
c
      if (n.gt.2) go to 10
      t = x(1) + x(2)
      x(2) = x(1) - x(2)
      x(1) = t
      return
  10  twopi = 8.*atan(1.0)
c
c first compute b0 term, where b0=sum of odd values of x(m)
c
      no2 = n/2
      no4 = n/4
      nind = no2 + 1
      b0 = 0.
      do 20 i=2,nind,2
        b0 = b0 + x(i)
  20  continue
      b0 = b0*2.
c
c for n = 4 skip recursion loop
c
      if (n.eq.4) go to 40
c
c form new sequence, y(m)=x(2*m)+(x(2*m+1)-x(2*m-1))
c
      do 30 i=2,no4
        ind = 2*i
        t1 = x(ind) - x(ind-2)
        y(i) = x(ind-1) + t1
        ind1 = no2 + 2 - i
        y(ind1) = x(ind-1) - t1
  30  continue
  40  y(1) = x(1)
      y(no4+1) = x(no2+1)
c
c take n/2 point (real) fft of y
c
      call fast(y, no2)
c
c form original dft by unscrambling y(k)
c use recursion to give sin(tpn*i) multiplier
c
      tpn = twopi/float(n)
      cosi = 2.*cos(tpn)
      sini = 2.*sin(tpn)
      cosd = cosi/2.
      sind = sini/2.
      nind = no4 + 1
      do 50 i=2,nind
        ind = 2*i
        bk = y(ind)/sini
        ak = y(ind-1)
        x(i) = ak + bk
        nind1 = n/2 + 2 - i
        x(nind1) = ak - bk
        temp = cosi*cosd - sini*sind
        sini = cosi*sind + sini*cosd
        cosi = temp
  50  continue
      x(1) = b0 + y(1)
      x(no2+1) = y(1) - b0
      return
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: ysBach/irafdocgen

max_forks_repo_name: ysBach/irafdocgen


text:
module blas
  use blas_s1
  use blas_s2
  use blas_s3
  use blas_d1
  use blas_d2
  use blas_d3
  use blas_c1
  use blas_c2
  use blas_c3
  use blas_z1
  use blas_z2
  use blas_z3
  implicit none
  public
end module blas

max_stars_repo_name: awvwgk/blas-interface

max_issues_repo_name: awvwgk/blas-interface

max_forks_repo_name: awvwgk/blas-interface


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 APATIT (PMSA  , FL    , IPOINT, INCREM, NOSEG ,
     +                   NOFLUX, IEXPNT, IKNMRK, NOQ1  , NOQ2  ,
     +                   NOQ3  , NOQ4  )

!
!     Description of the module :
!     precipitation and dissolution of an apatite-like P-mineral
!
!
! Name    T   L I/O   Description                                   Units
! ----    --- -  -    -------------------                            ----
! CPHD    R*4 1 I     concentration dissolved phosphate            [gP/m3]
! CPHDE   R*4 1 I     saturation concentration dissolved phosphate [gP/m3]
! CPHPR   R*4 1 I     concentration apatite                        [gP/m3]
! DELT    R*4 1 I     timestep                                         [d]
! FPRC    R*4 1 O     precipitation flux                         [gP/m3/d]
! FRR     R*4 1 I     ratio of apatite and vivianite preciptation rates[-]
! FSOL    R*4 1 O     dissolution flux                           [gP/m3/d]
! KPRC    R*4 1 I     precipitation rate                             [1/d]
! KSOL    R*4 1 I     dissolution rate                           [m3/gP/d]
! POROS   R*4 1 I     porosity                                         [-]
! TCPRC   R*4 1 I     temperature coefficient of precipitation         [-]
! TCSOL   R*4 1 I     temperature coefficient of dissolution           [-]
! TEMP    R*4 1 I     temperature                                     [oC]
! TMPPRC  R*4 1 -     temperature function for precipitation           [-]
! TMPSOL  R*4 1 -     temperature function for dissolution             [-]
!
!     Logical Units : -

!     Modules called : -

!     Name     Type   Library
!     ------   -----  ------------
!
      IMPLICIT REAL (A-H,J-Z)

      INTEGER  NOSEG , NOFLUX, NOQ1  , NOQ2  , NOQ3  ,  NOQ4
      INTEGER  IPOINT(*)       , INCREM(*),
     +         IEXPNT(4,*)     , IKNMRK(*)
      REAL     PMSA(*)         , FL(*)

      REAL     KSOL , KPRC   , FSOL   , FPRC  , FRR   ,
     +         TEMP , TMPSOL , TMPPRC , TCSOL , TCPRC ,
     +         CPHD , CPHPR  , CPHDE  , POROS ,
     +         DELT

      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)
!
      IFLUX = 0

      DO 9000 ISEG = 1 , NOSEG

      CALL DHKMRK(1,IKNMRK(ISEG),IKMRK1)

      IF (IKMRK1.EQ.1.OR.IKMRK1.EQ.3) THEN

         CPHD    = MAX(PMSA(IP1),0.0)
         CPHPR   = MAX(PMSA(IP2),0.0)
         CPHDE   = PMSA(IP3)
         KPRC    = PMSA(IP4)
         TCPRC   = PMSA(IP5)
         KSOL    = PMSA(IP6)
         TCSOL   = PMSA(IP7)
         FRR     = PMSA(IP8)
         TEMP    = PMSA(IP9)
         POROS   = PMSA(IP10)
         DELT    = PMSA(IP11)
!
!     Calculation of the precipitation or dissolution flux
!
         FPRC = 0.0
         FSOL = 0.0
         TMPPRC = TCPRC**(TEMP - 20.0)
         TMPSOL = TCSOL**(TEMP - 20.0)
!
         FPRC = FRR * KPRC * TMPPRC * ( CPHD / POROS - CPHDE ) * POROS
         FSOL = KSOL * TMPSOL * CPHPR * ( CPHDE - CPHD / POROS )
!
         IF ( FPRC .LT. 0.0) FPRC = 0.0
         IF ( FSOL .LT. 0.0) FSOL = 0.0
         IF ( FSOL*DELT .GE. CPHPR) FSOL = 0.5 * CPHPR / DELT
!
!     Output of module
!
         FL(1+IFLUX) = FPRC
         FL(2+IFLUX) = FSOL
!
!     End active cells block
!
      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 )
!
 9000 CONTINUE
!
      RETURN
!
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
Founded in 2003, the Law Office of MaryAlice Coleman is a law firm that handles employment matters on behalf of private and public sector employees seeking redress for unlawful and improper conduct in their workplaces. The firm has successfully litigated cases against the stats of California, the University of California, School Districts, corporations and other employers. The firm’s attorneys also have a broad range of experience representing clients in general litigation matters including personal injury, business disputes, and defending professional licenses in disciplinary matters.

  wiki:calbar:98365 CalBar Profile  Admitted to the bar on 198163.

Up the creek without a paddle?  Check out the Legal Services page.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
SUBROUTINE invar(stress,sigm,dsbar,theta)
!
! This subroutine forms the stress invariants in 2- or 3-d.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::stress(:)
 REAL(iwp),INTENT(OUT),OPTIONAL::sigm,dsbar,theta
 REAL(iwp)::sx,sy,sz,txy,dx,dy,dz,xj3,sine,s1,s2,s3,s4,s5,s6,ds1,ds2,ds3, &
   d2,d3,sq3,zero=0.0_iwp,small=1.e-10_iwp,one=1.0_iwp,two=2.0_iwp,       &
   three=3.0_iwp,six=6.0_iwp,thpt5=13.5_iwp
 INTEGER::nst
 nst=UBOUND(stress,1)
 SELECT CASE(nst)
 CASE(4)
   sx=stress(1)
   sy=stress(2)
   txy=stress(3)
   sz=stress(4)
   sigm=(sx+sy+sz)/three
   dsbar=SQRT((sx-sy)**2+(sy-sz)**2+(sz-sx)**2+six*txy**2)/SQRT(two)
   IF(dsbar<small)THEN
     theta=zero
   ELSE
     dx=(two*sx-sy-sz)/three
     dy=(two*sy-sz-sx)/three
     dz=(two*sz-sx-sy)/three
     xj3=dx*dy*dz-dz*txy**2
     sine=-thpt5*xj3/dsbar**3
     IF(sine>=one)sine=one
     IF(sine<-one)sine=-one
     theta=ASIN(sine)/three
   END IF
 CASE(6)
   sq3=SQRT(three)
   s1=stress(1)
   s2=stress(2)
   s3=stress(3)
   s4=stress(4)
   s5=stress(5)
   s6=stress(6)
   sigm=(s1+s2+s3)/three
   d2=((s1-s2)**2+(s2-s3)**2+(s3-s1)**2)/six+s4*s4+s5*s5+s6*s6
   ds1=s1-sigm
   ds2=s2-sigm
   ds3=s3-sigm
   d3=ds1*ds2*ds3-ds1*s5*s5-ds2*s6*s6-ds3*s4*s4+two*s4*s5*s6
   dsbar=sq3*SQRT(d2)
   IF(dsbar<small)THEN
     theta=zero
   ELSE
     sine=-three*sq3*d3/(two*SQRT(d2)**3)
     IF(sine>=one)sine=one
     IF(sine<-one)sine=-one
     theta=ASIN(sine)/three
   END IF
 CASE DEFAULT
   WRITE(*,*)"wrong size for nst in invar"
 END SELECT
RETURN
END SUBROUTINE invar

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


  !! ----------------------------------------------------------------------- !!
  !>
  !! Bathymetry depth hh(:,:) ( in m-unit ) must be filled here.
  !<
  block

    integer :: i, j
    !! --

    open(21,file="example/bath.dat", status='unknown')

    do i=1,nx
      read(21,*)(hh(i,j), j=1,ny)
    enddo
    close(21)
    !-------------------

  end block


text:
      program t
      implicit none
C     open stmt null (default) in the blank specifier
      integer::output
      open (87, status='old', file='infile', blank='null')
      read (87, '(i10)') output
      print *,output

      read (87, '(i10)') output
      print *,output
      endprogram t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
      REAL*4 function scamaxf(n,zx,incx)
c
c     finds the index of element having max. absolute value.
c     jack dongarra, 1/15/85.
c     modified 3/93 to return if incx .le. 0.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      COMPLEX*8 zx(*)
      integer i,incx,ix,n
      REAL*4 scabs1
c
      scamaxf = 0.
      if( n.lt.1 .or. incx.le.0 )return
      scamaxf = scabs1(zx(1))
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      scamaxf = scabs1(zx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(scabs1(zx(ix)).le.scamaxf) go to 5
         scamaxf = i
         scamaxf = scabs1(zx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
   20 scamaxf = scabs1(zx(1))
      do 30 i = 2,n
         if(scabs1(zx(i)).le.scamaxf) go to 30
         scamaxf = i
         scamaxf = scabs1(zx(i))
   30 continue
      return
      end

max_stars_repo_name: drhpc/OpenBLAS

max_issues_repo_name: drhpc/OpenBLAS

max_forks_repo_name: drhpc/OpenBLAS


text:
!-----subroutine---------------------------------------------g.guidati--
!
     subroutine DETFIELD
!
!-----------------------------------------------------------------------
!
!     franco guidati    IAG
!                       Uni Stuttgart
!
!
!     scope             determine sigma distribution
!
!.......................................................................
!               declarations
!.......................................................................
      USE TINoiseGeneric
      USE TINoiseInput
      USE TINoiseDDD
      USE TINoiseRHSin
      USE TIPrecision

IMPLICIT                        NONE

   ! Local variables.

INTEGER(4)                   :: klo,khi,i,n_in,nvar,nj,j1,j2,ivar,nok,nbad
REAL(DbKi)                   :: kk2, Isumwell, Isum, Itot,eps,h1,hmin
REAL(DbKi)                   :: ys1,ys2,poverall,phipot
REAL(DbKi)                   :: bigben(2),bigben2(2),dtime,etime
REAL(DbKi)                   :: phif(10)

COMPLEX(DbKi)                :: value,dval1,dval2,dval11,dval12,dval22
COMPLEX(DbKi)                :: btrans,btrans1,btrans2,btrans11,btrans12,btrans22
COMPLEX(DbKi)                :: bbb,bbb1,bbb2

EXTERNAL         :: SIGSUM,rkqs , RHSI

      nvar    =  10
      eps     =  1.0D-04
      h1      =  0.1d0
      hmin    =  0.0d0
      lderiv = .false.

!.......................................................................
!
!.......................................................................

!	write(79,*) 'ZONE I=
 !     write(*,*) etime(bigben2),dtime(bigben)
      ys1 = -200.0
      ys2 = 100.0

      poverall = 0.0
      nj = 101
!      write(722,*) 'variables = x1,x2,b,b1,b2,b11,b12,b22'
!      write(722,*) 'ZONE I=',nj,', J=',nj

      do j1=1,nj
	write(*,*) j1
        do j2=1,nj

	  x1 = (  -1.0 + 3.0 * dble(j1) / dble(nj)  )
	  x2 = -1.500324 + 3.0 * dble(j2) / dble(nj)

	  do ivar=1,nvar
	    phif(ivar) = 0.0
	  end do

 	  do ipath=1,npath
!            if(ipath.eq.npath/2+1) then
            call RHSINT(value,dval1,dval2,dval11,dval12,dval22)

            phif( 1) = phif( 1) + dble(value)
            phif( 2) = phif( 2) + dble(-imag*value)
            phif( 3) = phif( 3) + dble(dval1)
            phif( 4) = phif( 4) + dble(-imag*dval1)
            phif( 5) = phif( 5) + dble(dval2)
            phif( 6) = phif( 6) + dble(-imag*dval2)
            phif( 7) = phif( 7) + dble(dval11)
            phif( 8) = phif( 8) + dble(-imag*dval11)
            phif( 9) = phif( 9) + dble(dval12)
            phif(10) = phif(10) + dble(-imag*dval12)
!            endif
	  end do
          s1 = 0.0d0
          s2 = sworkat(na+1)
          call odeint(phif,nvar,s1,s2,eps,h1,hmin,nok,nbad,SIGSUM,rkqs)

          btrans   = dcmplx(phif(1),phif(2))
          btrans1  = dcmplx(phif(3),phif(4))
          btrans2  = dcmplx(phif(5),phif(6))
          btrans11 = dcmplx(phif(7),phif(8))
          btrans12 = dcmplx(phif(9),phif(10))
	  btrans22 = -btrans11 - kwave2 * btrans
	  phipot = x1
	  bbb = btrans * exp(-imag*kwave*mach_ti*phipot)
          bbb1 = (btrans1-btrans*imag*kwave*mach_ti)* exp(-imag*kwave*mach_ti*x1)
          bbb2 = btrans2 * exp(-imag*kwave*mach_ti*x1)
!	  write(722,*) real(x1),real(x2),real(bbb),real(-imag*bbb)
!     >                 ,real(btrans1),real(btrans2),
!     >                 real(btrans11),real(btrans12),real(btrans22)
	end do
      end do
!      call flush(722)

    !  write(*,*) etime(bigben2),dtime(bigben)

!.......................................................................
!               end of subroutine
!.......................................................................
      return
!***********************************************************************
      end

max_stars_repo_name: CarloSucameli/Openfast-noise

max_issues_repo_name: CarloSucameli/Openfast-noise

max_forks_repo_name: CarloSucameli/Openfast-noise


text:
      SUBROUTINE DT_RNDMIN(Uin,Cin,Cdin,Cmin,Iin,Jin)

      IMPLICIT NONE
      DOUBLE PRECISION Cdin , Cin , Cmin , Uin
      INTEGER Iin , Jin , kkk
      SAVE

C random number generator
      INCLUDE 'inc/dtrand'

      DIMENSION Uin(97)

      DO kkk = 1 , 97
         U(kkk) = Uin(kkk)
      END DO
      C = Cin
      CD = Cdin
      CM = Cmin
      I = Iin
      J = Jin

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
	real r (180), t (180)
c
	do 10 i = 1, 180
		t (i) = 2.*i - 2.
		r (i) = 1. + sin (t (i) * .01745329)
10              continue
	call setpen (2)
c       call setmtf ('card.mtf', 8)
	call page (30., 22., 'CARDIOIDE', 9, 1)
	call polreg (15., 11., 0., 10., 0., 360., 0, 0, 0., 0)
	call setpen (1)
	call limits (0., 2.2, 0., 360.)
	call poline (r, t, 180, -1, 1, 0, 1)
	call limits (0., 2.2, 360., 0.)
	call poline (r, t, 180, -1, 1, 0, 1)
	call setpen (6)
	call raxes ('RADIUS', 6, 0.5, 0.5, 1, 1, 0.)
	call setpen (3)
	call thaxes ('ANGLE', -5, 0.5, 30., 1, 1, 2.1)
	call setpen (5)
	call endpg ('kard')
	end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


! { dg-do compile }
! { dg-options "-O2 -fdump-tree-profile_estimate" }

subroutine test(block, array)
integer :: i, block(9), array(2)

do i = array(1), array(2), 2
    block(i) = i
end do

do i = array(1), array(2), -2
    block(i) = block(i) + i
end do

end subroutine test

! { dg-final { scan-tree-dump-times "Fortran loop preheader heuristics of edge\[^:\]*: 1.0%" 2 "profile_estimate" } }


text:
file tb.v
file ../common/ahbl_slave_assumptions.v
list $HDL/hazard3.f
include .


text:
      subroutine itransform(p,tp,x,ip,jp,kp)
************************************************************************
*     Author: R.K. Ellis                                               *
*     July, 1998.                                                      *
*     Given p ((n+1)-phase space) produce tp (n-phase space)           *
*     by creating jp                                                   *
************************************************************************
      implicit none
      include 'constants.f'
      include 'npart.f'
      double precision p(mxpart,4),tp(mxpart,4),k(4),kt(4),ks(4),
     . kDp(3:mxpart),ksDp(3:mxpart),kDk,ksDks,x
      integer ip,kp,j,nu,jp

      do nu=1,4
      tp(ip,nu)=x*p(ip,nu)
      tp(kp,nu)=p(kp,nu)
c---just so it is non-zero
      tp(jp,nu)=p(jp,nu)

      k(nu) =-p(ip,nu)-p(kp,nu)-p(jp,nu)
      kt(nu) =-tp(ip,nu)-tp(kp,nu)
      ks(nu)=k(nu)+kt(nu)
      enddo

      kDk=k(4)**2-k(1)**2-k(2)**2-k(3)**2
      ksDks=ks(4)**2-ks(1)**2-ks(2)**2-ks(3)**2

      do j=3,npart+2
      if (j .eq. jp) goto 20
      kDp(j)=k(4)*p(j,4)-k(1)*p(j,1)-k(2)*p(j,2)-k(3)*p(j,3)
      ksDp(j)=ks(4)*p(j,4)-ks(1)*p(j,1)-ks(2)*p(j,2)-ks(3)*p(j,3)

      do nu=1,4
      tp(j,nu)=p(j,nu)-two*ksDp(j)*ks(nu)/ksDks+two*kDp(j)*kt(nu)/kDk
      enddo
 20   continue
      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
program example4
  ! array index dependencies. Only iteration i of the first loop must
  ! complete before iteration i of the second loop.
  implicit none
  integer :: i
  real :: a(4)
  a = 0.0
  do i=1,4
     !$omp task out(a(i))
     call my_work1(i,a)
     !$omp end task
  end do

  do i=1,4
     !$omp task in(a(i))
     call my_work2(i,a)
     !$omp end task
  end do
end program example4

subroutine my_work1(i,a)
  implicit none
  real :: a(4)
  integer :: i
  print *, "hello", i
  a(i)=i
end subroutine my_work1

subroutine my_work2(i,a)
  implicit none
  real :: a(i)
  integer :: i
  print *, "world", i,a(i)
end subroutine my_work2

max_stars_repo_name: stfc/PSycloneBench

max_issues_repo_name: stfc/PSycloneBench

max_forks_repo_name: stfc/PSycloneBench