text:
module m_Thomas
  implicit none
  contains

    subroutine InnerProductForm(a)
      implicit none
      real(8),intent(inout) :: a(:,:)
      real(8) :: dajk
      integer :: i,j,k,n

      n=size(a(1,:))
      do k=1,n
        do j=1,k-1
          dajk = a(j,k)
          do i=j+1,n
            a(i,k) = a(i,k)-a(i,j)*dajk
          enddo
        enddo
        a(k,k) = 1.0d0/a(k,k)
        do i=k+1,n
          a(i,k) = a(i,k)*a(k,k)
        enddo
      enddo
      return
    end subroutine InnerProductForm

    subroutine Substitution(a,b)
      implicit none
      real(8),intent(inout) :: a(:,:),b(:,:)
      real(8),allocatable :: l(:,:),u(:,:),y(:,:)
      integer ::  i,j,k,n,m
      n = size(a(1,:))
      m = size(b(1,:))
      allocate(l(n,n),u(n,n),y(n,n))
      !L,Uを生成
      do j = 1,n
        do i = 1,j-1
          u(i,j) = a(i,j)
        enddo
        u(j,j) = 1.0d0/a(j,j)
        l(j,j) = 1
        do i = j+1,n
          l(i,j) = a(i,j)
        enddo
      enddo
      !ynを導出
      do k = 1,m
        y(1,k) = b(1,k)
        do i = 2,n
            y(i,k) = b(i,k)
            do j=1,i-1
              y(i,k) = y(i,k)-l(i,j)*y(j,k)
            enddo
        enddo
      enddo
      !xnを導出 (=b')
      do k = 1,m
        b(n,k) = y(n,k)/u(n,n)
        do i = n-1,1,-1
          b(i,k) = y(i,k)
          do j = i+1,n
            b(i,k) = b(i,k)-u(i,j)*b(j,k)
          enddo
          b(i,k) = b(i,k)/u(i,i)
        enddo
      enddo
      !write(*,*) b,"b"
      return
    end subroutine Substitution

    subroutine InverseMatrix(a,b)
      !B'=A^-1*Bを返す
      implicit none
      real(8),intent(in) :: a(:,:)
      real(8),intent(inout) :: b(:,:)
      real(8),allocatable :: IPF(:,:)
      integer :: n
      n = size(a(1,:))
      IPF = a
      call InnerProductForm(IPF)
      call Substitution(IPF,b)
      return
    end subroutine InverseMatrix

    subroutine n_Thomas(c,a,b,r)
      implicit none
      real(8),intent(inout) :: c(:,:,:),a(:,:,:),b(:,:,:),r(:,:,:)
      !row,column,index
      integer :: i,n

      n = size(r(1,1,:))
      call InverseMatrix(a(:,:,1),b(:,:,1))
      do i = 2,n
        a(:,:,i) = a(:,:,i)-MATMUL(c(:,:,i),b(:,:,i-1))
        call InverseMatrix(a(:,:,i),b(:,:,i))
      enddo
      call InverseMatrix(a(:,:,1),r(:,:,1))
      do i = 2,n
        r(:,:,i) = r(:,:,i) - MATMUL(c(:,:,i),r(:,:,i-1))
        call InverseMatrix(a(:,:,i),r(:,:,i))
      enddo
      do i = n-1,1,-1
        r(:,:,i) = r(:,:,i)- MATMUL(b(:,:,i),r(:,:,i+1))
      enddo
      return
    end subroutine n_Thomas

    subroutine Thomas(c,a,b,r)
      implicit none
      real(8),intent(inout) :: c(:),a(:),b(:),r(:)
      integer :: i,n

      n = size(r)
      b(1) = b(1)/a(1)
      do i = 2,n
        a(i) = a(i)-c(i)*b(i-1)
        b(i) = b(i)/a(i)
      enddo
      r(1) = r(1)/a(1)
      do i = 2,n
        r(i) = (r(i)-c(i)*r(i-1))/a(i)
      enddo
      do i = n-1,1,-1
        r(i) = r(i)-b(i)*r(i+1)
      enddo
      return
    end subroutine Thomas

end module m_Thomas

max_stars_repo_name: ShunsukeKobayashi/fuzzy-broccoli

max_issues_repo_name: ShunsukeKobayashi/fuzzy-broccoli

max_forks_repo_name: ShunsukeKobayashi/fuzzy-broccoli


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

module cov_cutoff_mod

use     types_mod, only : r8
use utilities_mod, only : register_module, error_handler, E_ERR, E_MSG, &
                          do_output, do_nml_file, do_nml_term, nmlfileunit, &
                          find_namelist_in_file, check_namelist_read
use location_mod,  only : location_type

implicit none
private

public :: comp_cov_factor

! version controlled file description for error handling, do not edit
character(len=256), parameter :: source   = &
   "$URL$"
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate  = "$Date$"

!============================================================================

!---- namelist with default values
logical :: namelist_initialized = .false.

integer :: select_localization = 1
! Value 1 selects default Gaspari-Cohn cutoff
! Value 2 selects boxcar
! Value 3 selects ramped boxcar

namelist / cov_cutoff_nml / select_localization

!============================================================================

contains

!======================================================================

function comp_cov_factor(z_in, c, obs_loc, obs_type, target_loc, target_kind, &
   localization_override)
!----------------------------------------------------------------------
! function comp_cov_factor(z_in, c)
!
! Computes a covariance cutoff function from Gaspari and Cohn
! QJRMS, 125, 723-757.  (their eqn. 4.10)
!
! z_in is the distance while c is the cutoff distance.
! For distances greater than 2c, the cov_factor returned goes to 0.

! Other ramping shapes are also available and can be selected by a namelist
! parameter. At present, these include a boxcar with the given halfwidth
! and a ramp in which the weight is set at 1.0 up to the half-width
! distance and then decreases linearly to 0 at twice the half-width
! distance.

! Additional information is passed in about the location and specific type of the
! observation and the location and generic kind of the variable being targeted for
! increments. These can be used for more refined algorithms that want to
! make the cutoff a function of these additional arguments.

implicit none

real(r8),                      intent(in) :: z_in, c
type(location_type), optional, intent(in) :: obs_loc, target_loc
integer,             optional, intent(in) :: obs_type, target_kind
integer,             optional, intent(in) :: localization_override
real(r8)                                  :: comp_cov_factor

real(r8) :: z, r
integer  :: iunit, io
integer  :: localization_type

!--------------------------------------------------------
! Initialize namelist if not already done
if(.not. namelist_initialized) then

   call register_module(source, revision, revdate)

   namelist_initialized = .true.

   ! Read the namelist entry
   call find_namelist_in_file("input.nml", "cov_cutoff_nml", iunit)
   read(iunit, nml = cov_cutoff_nml, iostat = io)
   call check_namelist_read(iunit, io, "cov_cutoff_nml")

   if (do_nml_file()) write(nmlfileunit,nml=cov_cutoff_nml)
   if (do_nml_term()) write(     *     ,nml=cov_cutoff_nml)

   if (do_output()) then
      select case (select_localization)
         case (1)
            call error_handler(E_MSG,'comp_cov_factor:', &
               'Standard Gaspari Cohn localization selected')
         case (2)
            call error_handler(E_MSG,'comp_cov_factor:', &
               'Boxcar localization selected')
         case (3)
            call error_handler(E_MSG,'comp_cov_factor:', &
               'Ramped localization selected')
         case default
            call error_handler(E_ERR,'comp_cov_factor', &
               'Illegal value of "select_localization" in cov_cutoff_mod namelist', &
                source, revision, revdate )
      end select
   endif

endif
!---------------------------------------------------------

if(present(localization_override)) then
   localization_type = localization_override
else
   localization_type = select_localization
endif

z = abs(z_in)

!----------------------------------------------------------

if(localization_type == 1) then ! Standard Gaspari Cohn localization

   if( z >= c*2.0_r8 ) then

      comp_cov_factor = 0.0_r8

   else if( z <= c ) then
      r = z / c
      comp_cov_factor = &
           ( ( ( -0.25_r8*r +0.5_r8 )*r +0.625_r8 )*r -5.0_r8/3.0_r8 )*r**2 + 1.0_r8
!!$           r**5 * (-0.25_r8 ) + &
!!$           r**4 / 2.0_r8 +              &
!!$           r**3 * 5.0_r8/8.0_r8 -       &
!!$           r**2 * 5.0_r8/3.0_r8 + 1.0_r8
   else

      r = z / c
      comp_cov_factor = &
           ( ( ( ( r/12.0_r8 -0.5_r8 )*r +0.625_r8 )*r +5.0_r8/3.0_r8 )*r -5.0_r8 )*r &
!!$           r**5 / 12.0_r8  -  &
!!$           r**4 / 2.0_r8   +  &
!!$           r**3 * 5.0_r8 / 8.0_r8 + &
!!$           r**2 * 5.0_r8 / 3.0_r8 - 5.0_r8*r &
           + 4.0_r8 - 2.0_r8 / (3.0_r8 * r)
   endif

else if(localization_type == 2) then ! BOXCAR localization

   if(z < 2.0_r8 * c) then
      comp_cov_factor = 1.0_r8
   else
      comp_cov_factor = 0.0_r8
   endif

else if(localization_type == 3) then ! Ramped localization

   if(z >= 2.0_r8 * c) then
      comp_cov_factor = 0.0_r8
   else if(z >= c .and. z < 2.0_r8 * c) then
      comp_cov_factor = (2.0_r8 * c - z) / c
   else
      comp_cov_factor = 1.0_r8
   endif

else ! Otherwise namelist parameter is illegal; this is an error

     call error_handler(E_ERR,'comp_cov_factor', &
              'Illegal value of "localization" in cov_cutoff_mod namelist', &
               source, revision, revdate )

endif

end function comp_cov_factor

end module cov_cutoff_mod

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
! ===========================================================================================================
module xx_integral_hh
  use xx_kinds
  use xx_public_variables
  use xx_pass_integ
  implicit none
  private :: COUPLING_HH
  public :: IFCT_HH_X12
contains
! ------------------------------
  function IFCT_HH_X12(dum) result(dsig)
    real(kind=double), dimension(dim(ii)), intent(in) :: dum ! vegas integration variable
    real(kind=double), dimension(dim(ii))             :: var ! internal integration variable
    real(kind=double)                  :: dsig
    integer                            :: inlo, iq
    real(kind=double), dimension(1:30) :: massin,massix,massin_1,massin_2
    real(kind=double), dimension(1:20) :: C
    real(kind=double)  :: m1,m2,s,qf,qr,beta,del_s4,gamt,beta1,beta2
    real(kind=double)  :: sw,alpha,alpha_s,mu_tgb,nlo,ALPHAS
    real(kind=double)  :: x1m,x1p,x1,x1_jac,x2m,x2p,x2,x2_jac,t2m,t2p,t2,t2_jac,s4m,s4p,s4,s4_jac
    real(kind=double)  :: z4m,z4p,z4,s4s4,theta_s4,prop_s4,betas4,t2s4m,t2s4p,t2s4,t2s4_jac
    real(kind=double)  :: z3m,z3p,z3,s3s3,theta_s3,prop_s3,beta1s3,beta2s3,t2s3m,t2s3p,t2s3,t2s3_jac
    real(kind=double)  :: s3m,s3p,s3,s3_jac,betax,t2xm,t2xp,t2x,t2x_jac,s4s3m,s4s3p,s4s3,s4s3_jac
    real(kind=double)  :: s4xm,s4xp,s4x,s4x_jac
    real(kind=double)  :: s4_1m,s4_1p,s4_1_jac,s4_1
    real(kind=double)  :: t2_1m,t2_1p,t2_1_jac,t2_1
    real(kind=double)  :: s3_1m,s3_1p,s3_1_jac,s3_1
    real(kind=double)  :: z3_1m,z3_1p,z3_1
    real(kind=double)  :: beta1_1,beta2_1
    real(kind=double)  :: s4_2m,s4_2p,s4_2_jac,s4_2
    real(kind=double)  :: t2_2m,t2_2p,t2_2_jac,t2_2
    real(kind=double)  :: s3_2
    real(kind=double)  :: beta_1,beta1_2,beta2_2
    real(kind=double)  :: LUMI
    real(kind=double)  :: HH_QBB,HH_QBV,HH_QBSY,HH_QBS,HH_QBD,HH_QBH,HH_QGH,HH_QGOS,HH_GBH,HH_GBOS

    if (ii>8) then                                             ! finish early
       dsig = 0.0
       return
    end if

    var(1:dim(ii)) = dum(1:dim(ii)) * ( 1.0 - 2.0*cut ) + cut  ! cut off the integration in general

    iq = -1                                                    ! no impact for the bottom density

    massin(1:30)   = 0.0                                       ! initialize the massin arrays
    massix(1:30)   = 0.0
    massin_1(1:30) = 0.0
    massin_2(1:30) = 0.0

    m1 = mch                                                   ! assign the final state masses
    m2 = mch
    mu_tgb = mu_susy * tan_b                                   ! for the delta mb corrections

    sw    = sqrt( 1.0 - mw**2/mz**2 )                          ! weak parameters in the on-shell scheme
    alpha = sqrt(2.0) * mw**2 * sw**2 /pi * gf

    x1m    = (m1+m2)**2 /sc                                    ! x1-x2 integration, map x->log(x)
    x1p    = 1.0
    x1     = x1m * (x1p/x1m)**var(1)
    x1_jac = x1 * log(x1p/x1m)

    x2m    = (m1+m2)**2 /sc /x1
    x2p    = 1.0
    x2     = x2m * (x2p/x2m)**var(2)
    x2_jac = x2 * log(x2p/x2m)

    s = x1 * x2 * sc                                           ! partonic cm energy

    if (iscaling==1) s = (m1+m2)**2 * (eta+1)                  ! overwrite integration for scaling fct

    if (isca==0) then                                          ! renormalization/factorization scale, factorization scale low
       qf = scafac  * (m1+m2)/4.0
       qr = scafac  * (m1+m2)/2.0
    else if (isca==1) then
       qf = scafac * sqrt(s)/4.0
       qr = scafac * sqrt(s)
    end if

    if (qf < 30.0) qf = 30.0                                   ! too small scales are not good
    if (qr < 30.0) qr = 30.0

    if (iscaling==0) then                                      ! nlo factor [always nlo alpha_s]
       alpha_s = ALPHAS(qr,2)
       nlo     = 4.0 * pi * alpha_s
    else if (iscaling==1) then
       nlo = 1.0
    end if

    if (ii<=0) then                                            ! coupling factor alpha_s not always nlo
       inlo = 0
    else if (ii>0) then
       inlo = 1
    end if

    theta_s4 = 0.0                                             ! the os subtraction theta function
    theta_s3 = 0.0
    gamt     = ewi * mt                                        ! always a small width needed for arcustan

    if ((s>(mt+m2)**2).and.(mt>m1)) theta_s4 = 1.0             ! intermediate top
    if ((s>(mt+m1)**2).and.(mt>m2)) theta_s3 = 1.0             ! intermediate top

    del_s4 = eps_sli * min(m1**2,m2**2)                      & ! unit is m^2
                   * (1.0-(m1+m2)**2/s)*(1.0-(m1-m2)**2/s)     ! rescale the s4 cutoff
    if (ii/=4) del_s4 = 0.0

    s4p    = 0.0                                                               ! only for real corrections

    select case (ii)                                                           ! all the phase spaces
    case(-1,0,1,2,3)                                                           ! born_lo, born_nlo, virt

       beta   = sqrt(1.0-(m1+m2)**2/s)                       &                 ! t2 integration
               *sqrt(1.0-(m1-m2)**2/s)
       t2m    = -1.0/2.0 * ( s + m2**2 - m1**2 + s*beta )
       t2p    = -1.0/2.0 * ( s + m2**2 - m1**2 - s*beta )
       t2     = var(3) * (t2p-t2m) + t2m
       t2_jac = t2p-t2m

    case(4)                                                                    ! real(qb)

       beta   = ((1.0-del_s4/s)**2-(m1+m2)**2/s)              &                ! t2 integration shifted
               *((1.0-del_s4/s)**2-(m1-m2)**2/s)
       if (beta>=0.0) then
          beta = sqrt(beta)
       else
          print*, " IFCT_HH_X12: sqrt(beta) not defined "
          beta = 0.0
       end if
       t2m    = -1.0/2.0 * ( s - del_s4 + m2**2 - m1**2 + s*beta )
       t2p    = -1.0/2.0 * ( s - del_s4 + m2**2 - m1**2 - s*beta )
       t2     = var(3) * (t2p-t2m) + t2m
       t2_jac = t2p-t2m

       s4m    = del_s4                                                         ! s4 mapped to log
       s4p    = s + t2 + m2**2 - m1**2 + s*m2**2/t2                            !  -> only a slight improvement
!       s4     = s4m * (s4p/s4m)**var(4)
!       s4_jac = s4 * log(s4p/s4m)
       if (s4p<s4m) then
          n_faulty = n_faulty+1
          dsig = 0.0
          return
       end if
       s4     = var(4) * (s4p-s4m) + s4m
       s4_jac = s4p-s4m

    case(5)                                                                    ! crossed(qg), no os subtraction

       beta   = ((1.0-del_s4/s)**2-(m1+m2)**2/s)              &                ! t2 integration shifted
               *((1.0-del_s4/s)**2-(m1-m2)**2/s)
       if (beta>=0.0) then
          beta = sqrt(beta)
       else
          print*, " IFCT_HH_X12: sqrt(beta) not defined "
          beta = 0.0
       end if
       t2m    = -1.0/2.0 * ( s - del_s4 + m2**2 - m1**2 + s*beta )
       t2p    = -1.0/2.0 * ( s - del_s4 + m2**2 - m1**2 - s*beta )
       t2     = var(3) * (t2p-t2m) + t2m
       t2_jac = t2p-t2m

       s4m    = del_s4                                                         ! s4 mapped to log
       s4p    = s + t2 + m2**2 - m1**2 + s*m2**2/t2                            !  -> only a slight improvement
       s4     = var(4) * (s4p-s4m) + s4m
       s4_jac = s4p-s4m

       s3_1m    = 0.0                                                          ! for the s3 integration, complete phase space
       s3_1p    = s + m1**2 - m2**2 - 2.0 * sqrt(s*m1**2)
       s3_1     = var(5) * (s3_1p-s3_1m) + s3_1m
       s3_1_jac = s3_1p-s3_1m

       beta1_1    = (s-s3_1-m2**2+m1**2)**2 - 4.0*m1**2*s
       if (beta1_1<0.0) then
          n_faulty = n_faulty+1
          dsig = 0.0
          return
       else
          beta1_1 = sqrt( beta1_1 )
       end if
       s4_1m    = s3_1/2.0/(s3_1+m2**2) * (s - s3_1 - m1**2 - m2**2 - beta1_1)
       s4_1p    = s3_1/2.0/(s3_1+m2**2) * (s - s3_1 - m1**2 - m2**2 + beta1_1)
       s4_1     = var(6) * (s4_1p-s4_1m) + s4_1m
       s4_1_jac = s4_1p-s4_1m

       beta2_1    = (s-s4_1-m1**2+m2**2)**2 - 4.0*s*m2**2
       if (beta2_1<0.0) then
          n_faulty = n_faulty+1
          dsig = 0.0
          return
       else
          beta2_1 = sqrt( beta2_1 )
       end if
       t2_1m    = -1.0/2.0 * ( s - s4_1 - m1**2 + m2**2 + beta2_1 )
       t2_1p    = -1.0/2.0 * ( s - s4_1 - m1**2 + m2**2 - beta2_1 )
       t2_1     = var(7) * (t2_1p-t2_1m) + t2_1m
       t2_1_jac = t2_1p-t2_1m

       t2_1_jac = t2_1_jac * 2.0/s4_1*(s4_1+m1**2) / beta2_1                   ! including the over-all jacobian

    case(6)                                                                    ! on-shell(qg)

       if (theta_s3==0.0) then
          dsig = 0.0
          return
       end if

       s3_1m    = 0.0                                                          ! for the s3 integration, complete phase space
       s3_1p    = s + m1**2 - m2**2 - 2.0 * sqrt(s*m1**2)
       z3_1m    = atan( (s3_1m + m2**2 - mt**2)/mt/gamt )
       z3_1p    = atan( (s3_1p + m2**2 - mt**2)/mt/gamt )
       z3_1     = var(3) * (z3_1p-z3_1m) + z3_1m

       s3_1     = mt*gamt*tan(z3_1) + mt**2 - m2**2
       s3_1_jac = ((s3_1+m2**2-mt**2)**2/mt/gamt+mt*gamt)*(z3_1p-z3_1m)

       beta1_1    = (s-s3_1-m2**2+m1**2)**2 - 4.0*m1**2*s
       if (beta1_1<0.0) then
          n_faulty = n_faulty+1
          dsig = 0.0
          return
       else
          beta1_1 = sqrt( beta1_1 )
       end if
       s4_1m    = s3_1/2.0/(s3_1+m2**2) * (s - s3_1 - m1**2 - m2**2 - beta1_1)
       s4_1p    = s3_1/2.0/(s3_1+m2**2) * (s - s3_1 - m1**2 - m2**2 + beta1_1)
       s4_1     = var(4) * (s4_1p-s4_1m) + s4_1m
       s4_1_jac = s4_1p-s4_1m

       beta2_1    = (s-s4_1-m1**2+m2**2)**2 - 4.0*s*m2**2
       if (beta2_1<0.0) then
          n_faulty = n_faulty+1
          dsig = 0.0
          return
       else
          beta2_1 = sqrt( beta2_1 )
       end if
       t2_1m    = -1.0/2.0 * ( s - s4_1 - m1**2 + m2**2 + beta2_1 )
       t2_1p    = -1.0/2.0 * ( s - s4_1 - m1**2 + m2**2 - beta2_1 )
       t2_1     = var(5) * (t2_1p-t2_1m) + t2_1m
       t2_1_jac = t2_1p-t2_1m

       t2_1_jac = t2_1_jac * 2.0/s4_1*(s4_1+m1**2) / beta2_1                   ! including the over-all jacobian

       s3_2 = mt**2 - m2**2                                                    ! go to constrained phase space

       beta1_2    = (s-s3_2-m2**2+m1**2)**2 - 4.0*m1**2*s
       if (beta1_2<0.0) then
          n_faulty = n_faulty+1
          dsig = 0.0
          return
       else
          beta1_2 = sqrt( beta1_2 )
       end if
       s4_2m    = s3_2/2.0/(s3_2+m2**2) * (s - s3_2 - m1**2 - m2**2 - beta1_2)
       s4_2p    = s3_2/2.0/(s3_2+m2**2) * (s - s3_2 - m1**2 - m2**2 + beta1_2)
       s4_2     = var(4) * (s4_2p-s4_2m) + s4_2m
       s4_2_jac = s4_2p-s4_2m

       beta2_2    = (s-s4_2-m1**2+m2**2)**2 - 4.0*s*m2**2
       if (beta2_2<0.0) then
          n_faulty = n_faulty+1
          dsig = 0.0
          return
       else
          beta2_2 = sqrt( beta2_2 )
       end if
       t2_2m    = -1.0/2.0 * ( s - s4_2 - m1**2 + m2**2 + beta2_2 )
       t2_2p    = -1.0/2.0 * ( s - s4_2 - m1**2 + m2**2 - beta2_2 )
       t2_2     = var(5) * (t2_2p-t2_2m) + t2_2m
       t2_2_jac = t2_2p-t2_2m

       t2_2_jac = t2_2_jac * 2.0/s4_2*(s4_2+m1**2) / beta2_2                   ! including the over-all jacobian

       prop_s3 = mt**2*gamt**2 / ((s3_1+m2**2-mt**2)**2+mt**2*gamt**2)         ! compensate for the wrong breit-wigner

    case(7,8)                                                                  ! crossed(gb) including on-shell

       if ( (ii==8).and.(theta_s4==0.0) ) then
          dsig = 0.0
          return
       end if

       s4m    = 0.0
       s4p    = s + m2**2 - m1**2 - 2.0 * sqrt(s*m2**2)

       z4m    = atan( (s4m + m1**2 - mt**2)/mt/gamt )                          ! s4 mapped to atan
       z4p    = atan( (s4p + m1**2 - mt**2)/mt/gamt )                          ! since s4t=s4+m1^2-mt^2
       z4     = var(3) * (z4p-z4m) + z4m

       s4     = mt*gamt*tan(z4) + mt**2 - m1**2
       s4_jac = ((s4+m1**2-mt**2)**2/mt/gamt+mt*gamt)*(z4p-z4m)

       beta   = sqrt( dabs( (s-s4-m1**2+m2**2)**2 - 4.0*s*m2**2 ) )            ! absolute value because of numerics
       t2m    = -1.0/2.0 * ( s - s4 - m1**2 + m2**2 + beta )                   ! t2 integration as usual
       t2p    = -1.0/2.0 * ( s - s4 - m1**2 + m2**2 - beta )
       t2     = var(4) * (t2p-t2m) + t2m
       t2_jac = t2p-t2m

       if (theta_s4==1.0) then                                                 ! restricted phase space for theta_s4=1
          s4_1 = mt**2 - m1**2

          beta_1   = sqrt( (s-s4_1-m1**2+m2**2)**2 - 4.0*s*m2**2 )
          t2_1m    = -1.0/2.0 * ( s - s4_1 - m1**2 + m2**2 + beta_1 )
          t2_1p    = -1.0/2.0 * ( s - s4_1 - m1**2 + m2**2 - beta_1 )
          t2_1     = var(4) * (t2_1p-t2_1m) + t2_1m
          t2_1_jac = t2_1p-t2_1m

          prop_s4 = mt**2*gamt**2 / ((s4+m1**2-mt**2)**2+mt**2*gamt**2)        ! compensate for the wrong breit-wigner
       end if

    end select

    massin(1)  = s                                             ! assign the mass arrays
    massin(2)  = t2
    massin(3)  = s4                                            ! s4 only for 2->3 processes
    massin(6)  = m1
    massin(7)  = mt
    massin(8)  = mz
    massin(9)  = mh1
    massin(10) = mh2
    massin(12) = qr                                            ! renormalization scale
    massin(13) = qf                                            ! factorization scale
    massin(15) = mg
    massin(16) = msq(-5)
    massin(17) = msq( 5)
    massin(18) = msq(-6)
    massin(19) = msq( 6)
    massin(20) = del_s4                                        ! slicing parameter only fort ii=4
    massin(21) = s4p                                           ! s4^max for the log(delta) terms
    massin(25) = gamt                                          ! width of the top
    massin(26) = 0.0                                           ! gamma/m of s channel particles in breit wigner

    if ( (ii==5).or.(ii==6).or.(ii==8) ) then                  ! s3 subtraction phase space
       massin_1(1:30) = massin(1:30)
       massin_1(2)    = t2_1
       massin_1(3)    = s4_1
       massin_1(4)    = s3_1
    end if

    if (ii==6) then                                            ! restricted phace space for s3 subtraction
       massin_2(1:30) = massin(1:30)
       massin_2(2)    = t2_2
       massin_2(3)    = s4_2
       massin_2(4)    = s3_2
    end if

    call COUPLING_HH(qr,C)

    if ( (m1+m2) < mz )  C(9)  = 0.0                           ! switch off s-channel couplings
    if ( (m1+m2) < mh1 ) C(10) = 0.0
    if ( (m1+m2) < mh2 ) C(11) = 0.0

    dsig = 0.0

    select case (ii)
    case(-1,0,1)                                                                             ! born
       dsig = dsig +       LUMI(inlo,60,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *   HH_QBB(massin,C)     * t2_jac
    case(2)                                                                                ! virtual
       dsig = dsig + nlo * LUMI(inlo,60,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *   HH_QBV(massin,C)     * t2_jac
    case(3)                                                                                ! virtual-susy
       dsig = dsig +       LUMI(inlo,60,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *   HH_QBSY(massin,C,mu_tgb,alpha_s) * t2_jac
    case(4)                                                                                ! real(qb)
       dsig = dsig + nlo * LUMI(inlo,60,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *(  HH_QBH(massin,C)                                                 &
                          +HH_QBD(massin,C) )   * t2_jac   * s4_jac
    case(5)                                                                                ! crossed(qg)
       dsig = dsig + nlo * LUMI(inlo,70,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *   HH_QGH(massin,C)     * t2_jac   * s4_jac
       if (theta_s3==0.0)                                                                   &
       dsig = dsig + nlo * LUMI(inlo,70,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *   HH_QGOS(massin_1,C)  * t2_1_jac * s4_1_jac * s3_1_jac
    case(6)                                                                                ! on-shell(qg)
       dsig = dsig + nlo * LUMI(inlo,70,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *(  HH_QGOS(massin_1,C)  * t2_1_jac * s4_1_jac * s3_1_jac            &
                         - HH_QGOS(massin_2,C)  * t2_2_jac * s4_2_jac * s3_1_jac * prop_s3 )
    case(7)                                                                                ! crossed(gb)
       dsig = dsig + nlo * LUMI(inlo,80,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *   HH_GBH(massin,C)     * t2_jac   * s4_jac
       if (theta_s4==0.0)                                                                   &
       dsig = dsig + nlo * LUMI(inlo,80,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *   HH_GBOS(massin,C)    * t2_jac   * s4_jac
    case(8)                                                                                ! on-shell(qg)
       dsig = dsig + nlo * LUMI(inlo,80,icoll,idub,iq,x1,x2,qf)/s**2                        &
                       *(  HH_GBOS(massin,C)    * t2_jac   * s4_jac                         &
                         - HH_GBOS(massin_1,C)  * t2_1_jac * s4_jac * prop_s4 )
    case default
       dsig = 0.0
    end select

    if (iscaling==0)                                         & ! some jacobians
         dsig = dsig * x1_jac * x2_jac

    if (iscaling==0)                                         &  ! what is missing from the matrix elements
        dsig = dsig * 4.0/(m1+m2)**2 * (4.D0*pi*alpha)**2

    if (iscaling==0)                                         &  ! result in pb
         dsig = dsig * gevpb

    ii_done(ii) = 1

  end function IFCT_HH_X12
! ------------------------------
  subroutine COUPLING_HH(qr,C)
    real(kind=double),                  intent(in)   :: qr
    real(kind=double), dimension(1:20), intent(out)  :: C

    integer              :: ic
    real(kind=double)    :: e,g
    real(kind=double)    :: t3,qq,sw2,cw2,sw,cw,c2w,sqrt2
    real(kind=double)    :: yb,yt
    real(kind=double)    :: sa,ca,sb,cb,sbma,cbma,sbpa,cbpa,c2b
    real(kind=double)    :: RUNM_EXT

    sqrt2 = sqrt(2.0)

    t3 = -0.5                                                  ! bottom quark quantum numbers
    qq = -1.0/3.0

    sw = sqrt( 1.0 - mw**2/mz**2 )                             ! weak mixing angle
    sw2 = sw**2
    cw2 = 1.0 - sw2
    cw  = sqrt(cw2)
    c2w = cw**2-sw**2

    sa = sin_a                                                 ! scalar higgs mixing angle
    ca = cos_a

    cb   = 1.D0/dsqrt(1.D0+tan_b**2)                           ! beta as in tan(beta)
    sb   = tan_b*cb

    sbma = sb*ca - cb*sa
    cbma = cb*ca + sb*sa
    sbpa = sb*ca + cb*sa
    cbpa = cb*ca - sb*sa
    c2b  = cb**2-sb**2

    e = 1.0                                                    ! over-all factor
    g = e/sw

    if (ii<=0) then
       yb = RUNM_EXT(qr,5,1)
       yt = RUNM_EXT(qr,6,1)
    else
       yb = RUNM_EXT(qr,5,2)
       yt = RUNM_EXT(qr,6,2)
    end if

    C(1:20) = 0.0                                              ! initialize arrays

    C(1)  = e/sw/cw * ( t3 - qq*sw2 )                          ! lq
    C(2)  = e/sw/cw * (    - qq*sw2 )                          ! rq
    C(3)  = e*qq                                               ! pq

    C(4)  =  g * yb*tan_b/mw                                   ! hl
    C(5)  =  g * yt/tan_b/mw                                   ! hr
    C(6)  =  g/sqrt2 * yb/mw * sa/cb                           ! h1
    C(7)  = -g/sqrt2 * yb/mw * ca/cb                           ! h2

    C(8)  = -e                                                 ! ssp
    C(9)  = -g*c2w/2.0/cw                                      ! ssz
    C(10) = -g*( mw*sbma + mz/2.0/cw*c2b*sbpa )                ! lambda1 (checked with ME)
    C(11) = -g*( mw*cbpa - mz/2.0/cw*c2b*cbma )                ! lambda2 (checked with ME)

  end subroutine COUPLING_HH

end module xx_integral_hh

max_stars_repo_name: sliem/docker-hep

max_issues_repo_name: sliem/docker-hep

max_forks_repo_name: sliem/docker-hep


! { dg-do run }
!
! PR 42647: Missed initialization/dealloc of allocatable scalar DT with allocatable component
!
! Contributed by Tobias Burnus <burnus@gcc.gnu.org>

type t
 integer, allocatable :: p
end type t
type(t), allocatable :: a

deallocate(a,stat=istat)
if (istat == 0) STOP 1
end


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

!>    @brief writes gradients in tecplot-format
!>    @param[in] ident index number for file name
!>    @param[in] ismpl local sample index
      SUBROUTINE write_gtecplot(ident,ismpl)
        use arrays
#ifdef AD
        use g_arrays
#endif
#ifdef AD_RM
        use arrays_ad
#endif
        use mod_genrl
        use mod_genrlc
        use mod_data
        use mod_time
        use mod_flow
        use mod_temp
        use mod_inverse
        use mod_linfos
        use mod_OMP_TOOLS
        IMPLICIT NONE
        integer :: ismpl
        integer :: i, j, k
        INCLUDE 'OMP_TOOLS.inc'

        DOUBLE PRECISION tu
        INTEGER i1, i2, i3, i4, i1s, i2s, ident, lblank, lout, locstr, &
          clast
        EXTERNAL lblank, locstr, clast

        CHARACTER filename*256, snumber*8, sline*1024

        IF ( .NOT. tec_out) RETURN
#ifdef NOPLT
        RETURN
#endif

!     get his own file discriptor index
        CALL omp_new_file_handler(lout,16)

        tu = tunit

        CALL chln(project,i1,i2)
        CALL chln(project_sfx(ismpl),i1s,i2s)
        IF (ident>=0) THEN
          WRITE(snumber,'(1I7)') ident
        ELSE IF (ident==-1) THEN
          WRITE(snumber,'(A8)') 'final'
        ELSE IF (ident==-2) THEN
          WRITE(snumber,'(A8)') 'debug'
        ELSE IF (ident==-3) THEN
          WRITE(snumber,'(A8)') 'ens_mean'
        ELSE IF (ident==-4) THEN
          WRITE(snumber,'(A8)') 'mean'
        ELSE IF (ident==-5) THEN
          WRITE(snumber,'(A8)') 'ens_mean'
        END IF
        CALL chln(snumber,i3,i4)

        IF (i1s==0) THEN
          filename = project(i1:i2) // '_G_p' // snumber(i3:i4) // '.plt'
          OPEN(lout,file=filename,status='unknown',blank='null')
          WRITE(lout,'(3A)') 'title = "', project(i1:i2) // '_' // &
            snumber(i3:i4), '"'
        ELSE
          filename = project(i1:i2) // project_sfx(ismpl) (i1s:i2s) // &
            '_G_p' // snumber(i3:i4) // '.plt'
          OPEN(lout,file=filename,status='unknown',blank='null')
          WRITE(lout,'(3A)') 'title = "', project(i1:i2) // &
            project_sfx(ismpl) (i1s:i2s) // '_' // snumber(i3:i4), '"'
        END IF

        IF (linfos(3)>=1) THEN
          WRITE(*,'(3A)') '  [W] : gradients Tecplot to "', &
            filename(1:lblank(filename)), '"'
        END IF

        WRITE(lout,'(13A)') 'variables = "x",            "y",            "z", "uindex",', &
          '  "i",  "j",  "k",', &
#ifndef AD_RM
          '         "g_head",         "g_temp",','         "g_conc"','         "g_pres"'
#else
          '         "head_ad",         "temp_ad",','         "conc_ad"','         "pres_ad"'

#endif
!     use sline for ZONE string
        i1 = locstr(project_sfx(ismpl),'_time_')
        i2 = locstr(project_sfx(ismpl),'_out')
        IF (i1>=1 .AND. i2>=1) THEN
          WRITE(sline,'(1A,1I7,3A,1I7)') 'zone T="', ident, &
            '.time step", SolutionTime=', project_sfx(ismpl) (i1+6:i2- &
            1), ', StrandID=', ident
          WRITE(lout,'(1A,3(A,I5),1A)') sline(1:clast(sline)), &
            ', i=', i0, ', j=', j0, ', k=', k0, ', f=point'
        ELSE
          WRITE(lout,'(3(A,I5),1A)') 'zone i=', i0, ', j=', j0, &
            ', k=', k0, ', f=point'
        END IF

!       output of box-centred data:
        DO k = 1, k0
          DO j = 1, j0
            DO i = 1, i0
              WRITE(lout, '(3(g15.6,1x),i10,3i6,4(e17.8,1x))') &
                delxa(i), delya(j), delza(k), uindex(i,j,k), i, j, k, &
!#ifndef AD_RM
                g_head(i,j,k,ismpl), g_temp(i,j,k,ismpl), g_conc(i,j,k,1,ismpl), g_pres(i,j,k,ismpl)
!#else
!                head_ad(i,j,k,ismpl), temp_ad(i,j,k,ismpl), conc_ad(i,j,k,1,ismpl), pres_ad(i,j,k,ismpl)
!#endif

            END DO
          END DO
        END DO

        CLOSE(lout)
        CALL omp_del_file_handler(lout)
        CALL compress_file(compress_out,filename)

        RETURN
      END


text: #redirect The Posh Bagel

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C    %W% %G%
C****************************************************************
C
C   File: p_report.f
C   Purpose: IPF shell program to process /REPORTS commands
C
C   Author: Walt Powell
C
C****************************************************************
C
	integer function p_report (in_buffer, out_buffer)
        character in_buffer *(*), out_buffer *(*)

        include 'ipfinc/parametr.inc'
        include 'ipfinc/blank.inc'
        include 'ipfinc/jobctl.inc'
        include 'ipfinc/lfiles.inc'
        include 'ipfinc/filnam.inc'
        include 'ipfinc/errorsw.inc'
        include 'ipfinc/pfstates.inc'
        include 'ipfinc/prt.inc'
        include 'ipfinc/errmsg.inc'
        include 'ipfinc/errorx.inc'
        include 'ipfinc/update.inc'
        include 'ipfinc/usranl.inc'

        character null * 1, linefeed * 1, capital * 40,
     &            text * 80, word(50) * 60, out_filename * 60
        integer apdoutbuf, o1, o2, status, vltdifrpt, lfodifrpt,
     &          areaintrpt, chglisrpt, inttierpt, phshftrpt,
     &          findstr, scrfil, open_file, open_apnd, capdifrpt,
     &          usranlrpt, busdifrpt, aredifrpt, lindifrpt, ownintchg,
     &          zonintchg, arsummrpt, genrpt, regxrpt, areatotrpt,
     &          qutdifrpt
        logical apndfil, opened

        p_report = 0     ! default return SUCCESS state
        numerr = 0       ! reinitialize error count
        null = char(0)
        linefeed = char(10)
        max_buf = len( out_buffer )
        out_buffer(1:1) = null
        apndfil = .true.
        opened = .false.
        out_filename = ' '

        i1 = 1
        i2 = index (in_buffer, null)
        num_word = 0

c*************** Now done by command parser ****************************
c***c
c***c       Replace "*[EOM]" with "(END)"
c***c
c***        last = index (in_buffer, null)
c***        ix = findstr (in_buffer(1:last), '*[EOM]')
c***        if (ix .gt. 0) then
c***           in_buffer(ix:ix+5) = '(END)' // null
c***        endif
c***
c**********************************************************************

        next = nxt_term(in_buffer)
        if (in_buffer(next:next) .eq. linefeed) then
           next = nxt_term(in_buffer(next+1:)) + next
        endif
        inrcd = in_buffer(1:next-1)
        call uscan(in_buffer(2:next-1), word, nwrd, '=',' ,'//linefeed)
c
c       Capitalize all word() except <filename>
c
        i = 0
        do while ( i .lt. nwrd )
           i = i + 1
           if ( word(i) .eq. '=' ) i = i + 1
           word(i) = capital(word(i))
           if ( word(i) .eq. 'OUTPUT' ) word(i) = 'APPEND'
           if ( word(i) .eq. 'OVERWRITE'  .or.
     &          word(i) .eq. 'APPEND' ) then
              i = i + 1
              if ( word(i) .eq. '=' ) i = i + 1
           endif
        enddo

        out_filename = ' '
        scrfil = 0
        isave_lprt = 0
        iwrd = 1
        do while (iwrd .le. nwrd)
           if ( word(iwrd) .eq. 'OVERWRITE'  .or.
     &          word(iwrd) .eq. 'APPEND'           ) then
              if ( word(iwrd) .eq. 'OVERWRITE' ) apndfil = .false.
              iwrd = iwrd + 1
              if (word(iwrd) .eq. '=') iwrd = iwrd + 1
              out_filename = word(iwrd)
              if (out_filename .ne. ' ') then

c                temporarily borrow "lprt" so reports can use "prtpkg"

                 if (lprt .eq. wscfil) then
                    write (errbuf(1),11)
  11                format(' programmer error in P_REPORT -- ',
     &                 'logical units are improperly assigned')
                    call prterx ('F',1)
                    goto 300
                 endif
                 isave_lprt = lprt
                 lprt = wscfil
                 scrfil = wscfil
                 call close_file (scrfil)
                 if ( apndfil ) then
                    status = open_apnd( scrfil, out_filename )
                 else
                    status = open_file(scrfil,out_filename, 'F', 'W',
     &                                 iost)
                 endif
                 if (status .ne. 0) then
                    scrfil = 0
                    lc = lastch( out_filename )
                    if ( lc .eq. 0 ) lc = 1
                    write( errbuf(1), 111 ) out_filename(1:lc)
  111               format(' Error opening file: ', a )
                    call prterx ('E', 1)
                    write( out_buffer, 113 ) out_filename(1:lc), null
  113               format(' !!! ERROR OPENING FILE :  ',a,' !!!',a)
                    p_report = 1
                    goto 300
                 else
                    opened = .true.
                 endif
              endif
           endif
           iwrd = iwrd + 1
        enddo

	iwrd = 1
        do while (iwrd .le. nwrd)
           if (word(iwrd)(1:6) .eq. 'SELECT') then
              iwrd = iwrd + 1

              if (word(iwrd)(1:9) .eq. 'BUS_INPUT') then
                 itext = findstr (in_buffer, 'BUS_INPUT')
     &                   + lastch (word(iwrd))
                 call businrpt(in_buffer(itext:), out_buffer, scrfil)
                 iwrd = nwrd + 1

              else if (word(iwrd)(1:12) .eq. 'BUS_BR_INPUT') then
                 itext = findstr (in_buffer, 'BUS_BR_INPUT')
     &                   + lastch (word(iwrd))
                 call busbrinrpt(in_buffer(itext:), out_buffer, scrfil)
                 iwrd = nwrd + 1

              else if (word(iwrd)(1:13) .eq. 'BUS_BR_OUTPUT') then
                 itext = findstr (in_buffer, 'BUS_BR_OUTPUT')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
  144               format(' Request rejected - unsolved base data in',
     &                     ' residence')
                    call prterx ('W', 1)
                 else
                    call busbrotrpt(in_buffer(itext:), out_buffer,
     &                              scrfil)
                 endif

              else if (word(iwrd)(1:8) .eq. 'LINEFLOW') then
                 itext = findstr (in_buffer, 'LINEFLOW')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call lineflowxl(in_buffer(itext:), out_buffer,
     &                              scrfil)
                 endif

              else if (word(iwrd)(1:11) .eq. 'AREA_TOTALS') then
                 itext = findstr (in_buffer, 'AREA_TOTALS')
     1                 + lastch (word(iwrd))
                 iwrd = nwrd + 1
                 if (ostates .le. 4) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call updzon
                    status = areatotrpt(in_buffer(itext:), out_buffer,
     &                                  scrfil)
                 endif

              else if (word(iwrd)(1:16) .eq. 'OVERLOADED_LINES') then
                 itext = findstr (in_buffer, 'OVERLOADED_LINES')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call ovldlnsrpt(in_buffer(itext:), out_buffer,
     &                              scrfil)
                 endif

              else if (word(iwrd)(1:14) .eq. 'OVERLOADED_TXS') then
                 itext = findstr (in_buffer, 'OVERLOADED_TXS')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call ovldtxsrpt(in_buffer(itext:), out_buffer,
     &                              scrfil)
                 endif

              else if (word(iwrd)(1:8) .eq. 'BUS_UVOV') then
                 itext = findstr (in_buffer, 'BUS_UVOV')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call busuvovrpt(in_buffer(itext:), out_buffer,
     &                              scrfil)
                 endif

              else if (word(iwrd)(1:9) .eq. 'VOLTAGE_C') then
                 itext = findstr (in_buffer, 'VOLTAGE_C')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = vltdifrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:9) .eq. 'CAPACITOR_C') then
                 itext = findstr (in_buffer, 'CAPACITOR_C')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = capdifrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:6) .eq. 'LINE_C') then
                 itext = findstr (in_buffer, 'LINE_C')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = lfodifrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:13) .eq. 'FILT_OUTPUT') then
                 itext = findstr (in_buffer, 'FILT_OUTPUT')
     1                 + lastch (word(iwrd))
                 iwrd = nwrd + 1
                 if (ostates .le. 4) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call filtoutrpt(in_buffer(itext:), out_buffer,
     &                              scrfil)
                 endif

              else if (word(iwrd)(1:10) .eq. 'AI_SUMMARY') then
                 itext = findstr (in_buffer, 'AI_SUMMARY')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call updzon
                    status = areaintrpt(in_buffer(itext:), out_buffer,
     &                                  scrfil)
                 endif

              else if (word(iwrd)(1:9) .eq. 'AREA_SUMM') then
                 itext = findstr (in_buffer, 'AREA_SUMM')
     1                 + lastch (word(iwrd))
                 iwrd = nwrd + 1
                 if (ostates .le. 4) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call updzon
                    status = arsummrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:8) .eq. 'TIE_LINE') then
                 itext = findstr (in_buffer, 'TIE_LINE')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call updzon
                    status = inttierpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:9) .eq. 'AREA_SUMM') then
                 itext = findstr (in_buffer, 'AREA_SUMM')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call updzon
                    status = arsummrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:11) .eq. 'NETWORK_CHA') then
                 itext = findstr (in_buffer, 'NETWORK_CHA')
     &                   + lastch (word(iwrd))
                 status = chglisrpt(in_buffer(itext:), out_buffer,
     &                              scrfil)
                 iwrd = nwrd + 1

              else if (word(iwrd)(1:6) .eq. 'BUS_OV') then
                 itext = findstr (in_buffer, 'BUS_OV')
     1                 + lastch (word(iwrd))
                 iwrd = nwrd + 1
                 if (ostates .le. 4) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call busovrpt(in_buffer(itext:), out_buffer,
     &                              scrfil)
                 endif

              else if (word(iwrd)(1:6) .eq. 'BUS_UV') then
                 itext = findstr (in_buffer, 'BUS_UV')
     1                 + lastch (word(iwrd))
                 iwrd = nwrd + 1
                 if (ostates .le. 4) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    call busuvrpt(in_buffer(itext:), out_buffer,
     &                              scrfil)
                 endif

              else if (word(iwrd)(1:8) .eq. 'GEN_SUMM') then
                 itext = findstr (in_buffer, 'GEN_SUMM')
     1                 + lastch (word(iwrd))
                 iwrd = nwrd + 1
                 if (ostates .le. 4) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = genrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:11) .eq. 'PHASE_SHIFT') then
                 itext = findstr (in_buffer, 'PHASE_SHIFT')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = phshftrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:8) .eq. 'REG_XFMR') then
                 itext = findstr (in_buffer, 'REG_XFMR')
     1                 + lastch (word(iwrd))
                 iwrd = nwrd + 1
                 if (ostates .le. 4) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = regxrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:13) .eq. 'USER_ANALYSIS') then
                 itext = findstr (in_buffer, 'USER_ANALYSIS')
     &                   + lastch (word(iwrd))
                 iwrd = nwrd + 1
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
c
c                   Any prior opened USER_ANALYSIS output must be
C                   closed
c
                    if (opened) then
                      call close_file (26)
                    endif
                    status = usranlrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                    if (opened) then
                      usrfil(numusr) = out_filename
                      call close_file (scrfil)
                    endif
                 endif

              else if (word(iwrd)(1:7) .eq. 'AR_COMP') then
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = aredifrpt(in_buffer, out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:8) .eq. 'BUS_COMP') then
                 itext = findstr (in_buffer, 'BUS_COMP')
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = busdifrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:7) .eq. 'BR_COMP') then
                 itext = findstr (in_buffer, 'BR_COMP')
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = lindifrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:9) .eq. 'OWNER_INT') then
                 itext = findstr (in_buffer, 'OWNER_INT')
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = ownintchg(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:8) .eq. 'ZONE_INT') then
                 itext = findstr (in_buffer, 'ZONE_INT')
                 iwrd = nwrd + 1
c                 if (ostates .le. 4) then
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = zonintchg(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else if (word(iwrd)(1:13) .eq. 'REACTIVE_UTIL') then
                 itext = findstr (in_buffer, 'REACTIVE_UTIL')
                 iwrd = nwrd + 1
                 if ( lskp .eq. 0 ) then
                    p_report = 1
                    write (errbuf(1), 144)
                    call prterx ('W', 1)
                 else
                    status = qutdifrpt(in_buffer(itext:), out_buffer,
     &                                 scrfil)
                 endif

              else
                 itext = findstr (in_buffer, 'SELECT')
     &                   + len ('SELECT')
                 call outreport(in_buffer(itext:), out_buffer, scrfil)
                 iwrd = nwrd + 1
              endif
           else
              iwrd = iwrd + 1
           endif
        enddo

        call prtime('REPORTS')

  300   continue

        if (isave_lprt .ne. 0) lprt = isave_lprt

        if (scrfil .gt. 0) then
           call close_file (scrfil)
           scrfil = 0
        endif

        inrcd = buf
c
c       Append error messages to buffer
c
        j = 1
        length = 1
        o2 = index (out_buffer,null)
        do while (j .le. numerr .and. length .gt. 0)
           length = apdoutbuf(o2, errm(j), out_buffer(o2:))
           o2 = o2 + length
           j = j + 1
        enddo
c
c 	Append summary
c
        if ( o2 .gt. max_buf - 50 ) then
           o2 = max_buf - 50
           do while ( o2 .gt. 1  .and.
     &                out_buffer(o2:o2) .ne. linefeed )
              o2 = o2 - 1
           enddo
           out_buffer(o2:o2) = null
        endif
        write (text, 340) 'p_report.f', p_report, ostates
  340   format ('/', a, ' return status: ', i2, ' IPF state: ', i2)
        length = apdoutbuf(o2, text, out_buffer(o2:))
        o2 = o2 + length
c
c       Reset error flag
c
        call setercnt (0, ' ')
c
c       Debug printout invoked with IPF command /TRACE, CHANGE = ON
c
        if (kase1(27) .ne. 0) then
           o1 = 1
           do while (o1 .lt. max_buf .and.
     &              (out_buffer(o1:o1) .ne. null))
              next = nxt_term (out_buffer(o1+1:)) + o1
              o2 = min (o1+131, next)
              write (dbug, 830) out_buffer(o1:o2)
  830         format (1x, a)
              o1 = next
              if (out_buffer(o1:o1) .eq. linefeed) o1 = o1 + 1
           enddo
        endif
        return
	end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      SUBROUTINE chemstep(I,J,change)
!@sum chemstep Calculate new concentrations after photolysis & chemistry
!@auth Drew Shindell (modelEifications by Greg Faluvegi)
!@ver  1.0 (based on ds3ch4_chem_calc_jun1202_M23)
!@calls rates,chem1,chem1prn
c
C**** GLOBAL parameters and variables:
C
      USE RESOLUTION, only    : im,jm,lm,ls1
      USE DYNAMICS, only      : am, byam
      USE GEOM, only          : BYDXYP,dxyp
      USE TRACER_COM, only: n_CH4,n_CH3OOH,n_Paraffin,n_PAN,n_Isoprene,
     &                   n_AlkylNit,n_Alkenes,n_N2O5,n_NOx,n_HO2NO2,
     &                   n_Ox,n_HNO3,n_H2O2,n_CO,n_HCHO,trm,ntm
      USE TRCHEM_Shindell_COM, only: chemrate,photrate,mass2vol,
     &                   yCH3O2,yC2O3,yXO2,yXO2N,yRXPAR,yAldehyde,
     &                   yROR,nCH3O2,nC2O3,nXO2,nXO2N,nRXPAR,
     &                   nAldehyde,nROR,nr,nn,dt2,nss,ks,dest,prod,
     &                   ny,nhet,rr,nO1D,nOH,nNO,nHO2,ta,nM,ss,
     &                   nO3,nNO2,nNO3,prnrts,jprn,iprn,lprn,ay,
     &                   prnchg,y,bymass2vol,kss,nps,kps,nds,kds,
     &                   npnr,nnr,ndnr,kpnr,kdnr,nH2O,changeAldehyde
c
      IMPLICIT NONE
c
C**** Local parameters and variables and arguments:
!@var change change due to chemistry in mass/time
!@var I,J passed horizontal spatial indicies
!@var L,iter dummy loop variable
!@var maxl highest level with chemistry
!@var qqqCH3O2,CH3O2loss,XO2_NO,XO2N_HO2,RXPAR_PAR,ROR_CH2,C2O3prod,
!@+   C2O3dest,XO2prod,XO2dest,XO2_XO2,XO2Nprod,XO2Ndest,RXPARprod,
!@+   RXPARdest,Aldehydeprod,Aldehydedest,RORprod,RORdest,total,
!@+   rnewval,sumN,dNOx,ratio,sumD,newD,ratioD,newP,ratioP,changeA
!@+   sumP dummy temp variables
!@var tempiter,tempiter2 temp vars for equilibrium calcs iterations
      REAL*8, DIMENSION(IM,JM,LM,ntm)     :: change
      INTEGER, INTENT(IN) :: I,J
      INTEGER L,iter,maxl,igas
      REAL*8 qqqCH3O2,CH3O2loss,XO2_NO,XO2N_HO2,RXPAR_PAR,ROR_CH2,
     & C2O3prod,C2O3dest,XO2prod,XO2dest,XO2_XO2,XO2Nprod,XO2Ndest,
     & RXPARprod,RXPARdest,Aldehydeprod,Aldehydedest,RORprod,RORdest,
     & total,rnewval,sumN,dNOx,ratio,sumD,newD,ratioD,newP,ratioP,
     & changeA,sumP,tempiter,tempiter2
C
C     TROPOSPHERIC CHEMISTRY ONLY:
      maxl=LS1-1

      do L=1,maxl
       y(nCH3O2,L)   =yCH3O2(I,J,L)
       y(nC2O3,L)    =yC2O3(I,J,L)
       y(nXO2,L)     =yXO2(I,J,L)
       y(nXO2N,L)    =yXO2N(I,J,L)
       y(nRXPAR,L)   =yRXPAR(I,J,L)
       y(nAldehyde,L)=yAldehyde(I,J,L)
       y(nROR,L)     =yROR(I,J,L)
      enddo
C
cc    calculate reaction rates with present concentrations
      call rates(maxl,I,J)
c
c     chem1 call sample:
c     (klist,l,numel,nlist,ndlist,rate,change,multip)
c              numel=number of elements in reaction list nlist (1 or 2)
c              change=dest or prod array
c              multip=1(prod) or -1(dest)
c     chemical destruction:
      call chem1(kdnr,maxl,2,nn,ndnr,chemrate,dest,-1)
c     chemical production:
      call chem1(kpnr,maxl,2,nnr,npnr,chemrate,prod,1)
c     photolytic destruction:
      call chem1(kds,maxl,1,ks,nds,photrate,dest,-1)
c     photolytic production:
      call chem1(kps,maxl,2,kss,nps,photrate,prod,1)
C
c     Oxidation of Isoprene and Alkenes produces less than one
c      HCHO, Alkenes, and CO per rxn, correct here following Houweling
      do L=1,maxl
       prod(n_CO,L)=prod(n_CO,L)-0.63*rr(35,I,J,L)*y(n_Alkenes,L)*
     &  y(nO3,L)*dt2
       prod(n_HCHO,L)=prod(n_HCHO,L)-0.36*rr(35,I,J,L)*y(n_Alkenes,L)*
     &  y(nO3,L)*dt2
       prod(n_HCHO,L)=prod(n_HCHO,L)-0.39*rr(30,I,J,L)*y(n_Isoprene,L)*
     &  y(nOH,L)*dt2
       prod(n_Alkenes,L)=prod(n_Alkenes,L)-0.42*rr(30,I,J,L)*
     &  y(n_Isoprene,L)*y(nOH,L)*dt2
       prod(n_HCHO,L)=prod(n_HCHO,L)-0.10*rr(31,I,J,L)*y(n_Isoprene,L)*
     &  y(nO3,L)*dt2
       prod(n_Alkenes,L)=prod(n_Alkenes,L)-0.45*rr(31,I,J,L)*
     &  y(n_Isoprene,L)*y(nO3,L)*dt2
      enddo
c
c     set CH3O2 values (concentration = production/specific loss)
      do L=1,maxl
        iter=1
        qqqCH3O2=(rr(11,I,J,L)*y(nO1D,L)+rr(12,I,J,L)*y(nOH,L))
     &  *y(n_CH4,L)+rr(23,I,J,L)*y(n_CH3OOH,L)*y(nOH,L)

       tempiter=rr(20,I,J,L)*y(nNO,L)+rr(22,I,J,L)*y(nHO2,L)
 10    CH3O2loss=tempiter+rr(27,I,J,L)*yCH3O2(I,J,L)
       if(CH3O2loss.gt.1E-7)then
          y(nCH3O2,L)=qqqCH3O2/CH3O2loss
        else
          y(nCH3O2,L)=1.0
        endif
        yCH3O2(I,J,L)=y(nCH3O2,L)
        iter=iter+1
        if(iter.le.7)goto10 ! replace with while loop?
      enddo
c
      do L=1,maxl
c
c       Set C2O3, XO2, XO2N, RXPAR, Aldehyde & ROR values
c
c        First set various specific loss rates
         XO2_NO=y(nNO,L)*4.2E-12*exp(180./ta(L))
         XO2N_HO2=y(nHO2,L)*y(nNO,L)*rr(44,I,J,L)*
     &   rr(43,I,J,L)/XO2_NO
         RXPAR_PAR=y(n_Paraffin,L)*8.E-11
         ROR_CH2=1.6E3
c
c       Set value for C2O3
        iter=1
        C2O3prod=rr(38,I,J,L)*yAldehyde(I,J,L)*y(nOH,L)+
     &  (rr(29,I,J,L)*y(nM,L)+ss(15,I,J,L))*y(n_PAN,L)+0.15*
     &  rr(31,I,J,L)*y(nO3,L)*y(n_Isoprene,L)
        tempiter=rr(39,I,J,L)*y(nNO,L)+rr(54,I,J,L)*y(nNO2,L)+
     &   rr(41,I,J,L)*y(nHO2,L)
 20     C2O3dest=tempiter+rr(40,I,J,L)*yC2O3(I,J,L)
        if(C2O3dest.gt.1E-7)then
          y(nC2O3,L)=(C2O3prod/C2O3dest)
        else
          y(nC2O3,L)=1.0
        endif
        yC2O3(I,J,L)=y(nC2O3,L)
        iter=iter+1
        if(iter.le.7)goto20     ! replace with while loop?
c
c       Set value for XO2
        iter=1
        XO2prod=ss(16,I,J,L)*yAldehyde(I,J,L)+
     &  y(nC2O3,L)*(rr(39,I,J,L)*y(nNO2,L)+rr(40,I,J,L)*
     &  y(nC2O3,L)*2.+rr(41,I,J,L)*y(nHO2,L))
     &  +rr(42,I,J,L)*yROR(I,J,L)*0.96
     &  +y(nOH,L)*(rr(37,I,J,L)*y(n_Paraffin,L)*0.87+rr(34,I,J,L)*
     &  y(n_Alkenes,L)+rr(30,I,J,L)*y(n_Isoprene,L)*0.85+
     &  rr(33,I,J,L)*y(n_AlkylNit,L))+
     &  y(nO3,L)*(rr(35,I,J,L)*y(n_Alkenes,L)*0.29+
     &  rr(31,I,J,L)*y(n_Isoprene,L)*0.18)
        tempiter=XO2_NO+rr(43,I,J,L)*y(nHO2,L)
        tempiter2=1.7E-14*exp(1300./ta(L))
 30     XO2_XO2=yXO2(I,J,L)*tempiter2
        XO2dest=tempiter+XO2_XO2
        if(XO2dest.gt.1E-7.and.ss(16,I,J,L).gt.5E-5)then
          y(nXO2,L)=(XO2prod/XO2dest)
        else
          y(nXO2,L)=1.0
        endif
        yXO2(I,J,L)=y(nXO2,L)
        iter=iter+1
        if(iter.le.7)goto30   ! replace with while loop?
c
c       Set value for XO2N
        XO2Nprod=rr(37,I,J,L)*y(n_Paraffin,L)*y(nOH,L)*0.13+
     &  rr(42,I,J,L)*yROR(I,J,L)*0.04+rr(30,I,J,L)*y(n_Isoprene,L)*
     &  y(nOH,L)*0.15
        XO2Ndest=XO2N_HO2+rr(44,I,J,L)*y(nNO,L)
        if(XO2Ndest.gt.1E-7)then
          y(nXO2N,L)=(XO2Nprod/XO2Ndest)
        else
          y(nXO2N,L)=1.0
        endif
        yXO2N(I,J,L)=y(nXO2N,L)
c
c       Set value for RXPAR
        RXPARprod=rr(37,I,J,L)*y(n_Paraffin,L)*y(nOH,L)*0.11+
     &  rr(34,I,J,L)*yROR(I,J,L)*2.1+rr(35,I,J,L)*y(n_Alkenes,L)*
     &  y(nO3,L)*0.9
        RXPARdest=RXPAR_PAR
        if(RXPARdest.gt.0.)then
          y(nRXPAR,L)=(RXPARprod/RXPARdest)
        else
          y(nRXPAR,L)=1.0
        endif
        yRXPAR(I,J,L)=y(nRXPAR,L)
c
c       Set value for Aldehyde
        Aldehydeprod=rr(37,I,J,L)*y(n_Paraffin,L)*y(nOH,L)*0.11+
     &  rr(34,I,J,L)*y(n_Alkenes,L)*y(nOH,L)+
     &  rr(42,I,J,L)*yROR(I,J,L)*1.1+rr(35,I,J,L)*y(n_Alkenes,L)*
     &  y(nO3,L)*0.44
        Aldehydedest=rr(38,I,J,L)*y(nOH,L)+ss(16,I,J,L)

c       Check for equilibrium
        if(Aldehydedest*y(nAldehyde,L)*dt2.lt.y(nAldehyde,L))then
         changeAldehyde=
     &   (Aldehydeprod-y(nAldehyde,L)*Aldehydedest)*dt2
         if(changeAldehyde.gt.y(nAldehyde,L))
     &   changeAldehyde=y(nAldehyde,L)
         y(nAldehyde,L)=y(nAldehyde,L)+changeAldehyde
         if(y(nAldehyde,L).lt.0.)y(nAldehyde,L)=1.0
        else
          y(nAldehyde,L)=(Aldehydeprod/(Aldehydedest+0.5E-5))
        endif
        yAldehyde(I,J,L)=y(nAldehyde,L)
c
c       Set value for ROR
        RORprod=rr(37,I,J,L)*y(n_Paraffin,L)*y(nOH,L)*0.76+
     &  rr(42,I,J,L)*yROR(I,J,L)*0.02
        RORdest=rr(42,I,J,L)+ROR_CH2
        if(RORdest.gt.0.)then
          y(nROR,L)=(RORprod/RORdest)
        else
          y(nROR,L)=1.0
        endif
        yROR(I,J,L)=y(nROR,L)
c
c       Add in parrafin loss term via rxpar reaction
        dest(n_Paraffin,L)=dest(n_Paraffin,L)-y(nRXPAR,L)*RXPAR_PAR*dt2

c       Add in CH3OOH production via XO2N + HO2
        prod(n_CH3OOH,L)=prod(n_CH3OOH,L)+XO2N_HO2*y(nXO2N,L)*dt2
c
      end do  ! L
c
c     If NOx in equil with N2O5, HO2NO2 or PAN, remove from changes
      do L=1,maxl
       if(-dest(n_N2O5,L).ge.y(n_N2O5,L))then
        dest(n_NOx,L)=dest(n_NOx,L)+2.*chemrate(52,L)
        prod(n_NOx,L)=prod(n_NOx,L)-2.*(chemrate(46,L)+photrate(7,L))
       endif
       if(-dest(n_HO2NO2,L).ge.y(n_HO2NO2,L))then
        dest(n_NOx,L)=dest(n_NOx,L)+chemrate(51,L)
        prod(n_NOx,L)=prod(n_NOx,L)-(chemrate(18,L)+chemrate(45,L)+
     &   photrate(10,L)+photrate(11,L))
       endif
       if(-dest(n_PAN,L).ge.y(n_PAN,L))then
        dest(n_NOx,L)=dest(n_NOx,L)+chemrate(54,L)
        prod(n_NOx,L)=prod(n_NOx,L)-(chemrate(29,L)+photrate(15,L))
       endif
      enddo
c
c     print rxtn rates if desired (chem1prn : argument before multip is
c     index = number of call):
c
       if(prnrts.and.J.eq.jprn.and.I.eq.iprn)then
        do igas=1,ntm
         total=0.
         write(6,108) ' Species: ',ay(1,igas),ay(2,igas)
c
         call chem1prn
     *   (kdnr,2,nn,ndnr,chemrate,1,-1,igas,total,maxl,I,J)
c
         if(igas.eq.n_Ox)write(6,109)
     *   ' losses from 10 and 11 reset by fraction of O1D via rxn,',
     *   'dy = ?'
         if(igas.eq.n_NOx.and.-dest(n_HO2NO2,lprn).ge.y(n_HO2NO2,lprn))
     *    write(6,110)'loss by reaction 51 (HO2NO2 formation) removed',
     *    rr(51,I,J,lprn)*y(nHO2,lprn)*y(nNO2,lprn)*dt2
         if(igas.eq.n_NOx.and.-dest(n_N2O5,lprn).ge.y(n_N2O5,lprn))
     *    write(6,110)'losses by reaction 52 (N2O5 formation) removed',
     *    2*rr(52,I,J,lprn)*y(nNO3,lprn)*y(nNO2,lprn)*dt2
         if(igas.eq.n_NOx.and.-dest(n_PAN,lprn).ge.y(n_PAN,lprn))
     *    write(6,110)'losses by reaction 54 (PAN formation) removed',
     *    chemrate(54,lprn)
c
         call chem1prn
     &   (kpnr,2,nnr,npnr,chemrate,2,1,igas,total,maxl,I,J)
c
         if(igas.eq.n_NOx.and.-dest(n_HO2NO2,lprn).ge.y(n_HO2NO2,lprn))
     *    write(6,110)'gain by reactions 18 & 45 (from HO2NO2) rmoved',
     *    (rr(18,I,J,lprn)*y(nOH,L)+rr(45,I,J,lprn)*y(nM,lprn)+
     *    ss(10,I,J,lprn)+ss(11,I,J,lprn))*y(n_HO2NO2,lprn)*dt2
         if(igas.eq.n_NOx.and.-dest(n_N2O5,lprn).ge.y(n_N2O5,lprn))
     *        write(6,110)
     *        'gains by reaction 46 (N2O5 decomposition) removed',
     *        2.*(rr(46,I,J,lprn)*y(nM,lprn))*y(n_N2O5,lprn)*dt2
         if(igas.eq.n_NOx.and.-dest(n_PAN,lprn).ge.y(n_PAN,lprn))
     *    write(6,110)'gain by reaction 29 (from PAN) removed',
     *    rr(29,I,J,lprn)*y(nM,lprn)*y(n_PAN,lprn)*dt2
c
         call chem1prn(kds,1,ks,nds,photrate,3,-1,igas,total,maxl,I,J)
c
         call chem1prn(kps,2,kss,nps,photrate,4,1,igas,total,maxl,I,J)
c
         if(igas.eq.n_NOx.and.-dest(n_N2O5,lprn).ge.y(n_N2O5,lprn))
     *    write(6,110)'gains by reaction 7 (N2O5 photolysis) removed',
     *    ss(7,I,J,lprn)*y(n_N2O5,lprn)*dt2
         if(igas.eq.n_NOx.and.-dest(n_N2O5,lprn).ge.y(n_N2O5,lprn))
     *        write(6,110)
     *        'net change due to N2O5 is ',
     *        2.*(y(n_N2O5,lprn)-(rr(52,I,J,lprn)*y(nNO3,lprn)*
     *        y(nNO2,lprn))/
     *        (rr(46,I,J,lprn)*y(nM,lprn)+ss(7,I,J,lprn)))
         if(igas.eq.n_NOx.and.-dest(n_HO2NO2,lprn).ge.y(n_HO2NO2,lprn))
     *    write(6,110)'gain by rxns 10 & 11 (HO2NO2 photolysis) rmoved'
     *    ,(ss(10,I,J,lprn)+ss(11,I,J,lprn))*y(n_HO2NO2,lprn)*dt2
         if(igas.eq.n_NOx.and.-dest(n_HO2NO2,lprn).ge.y(n_HO2NO2,lprn))
     *        write(6,110)'net change due to HO2NO2 is ',
     *        y(n_HO2NO2,lprn)-((rr(51,I,J,lprn)*y(nHO2,lprn)*
     *       y(nNO2,lprn))/(rr(18,I,J,lprn)*
     *     y(nOH,lprn)+rr(45,I,J,lprn)*y(nM,lprn)+ss(10,I,J,lprn)
     *     +ss(11,I,J,lprn)))
         if(igas.eq.n_NOx.and.-dest(n_PAN,lprn).ge.y(n_PAN,lprn))
     *        write(6,110)'net change due to PAN is ',
     *        y(n_PAN,lprn)-((rr(54,I,J,lprn)*y(nC2O3,lprn)*
     *        y(nNO2,lprn))/
     *        (rr(29,I,J,lprn)*y(nM,lprn)+ss(15,I,J,lprn)))
         if(igas.eq.n_Ox.or.igas.eq.n_NOx)total=
     *    100.*(dest(igas,lprn)+prod(igas,lprn))/y(igas,lprn)
         if(igas.eq.n_CH3OOH) write(6,'(a48,a6,e10.3)')
     *    'production from XO2N + HO2 ','dy = ',
     *    y(nHO2,lprn)*y(nNO,lprn)*rr(44,I,J,lprn)*rr(43,I,J,lprn)/
     *    (y(nNO,lprn)*4.2E-12*exp(180./ta(lprn)))
     *    *y(nXO2N,lprn)*dt2
         if(igas.eq.n_Paraffin) write(6,'(a48,a6,e10.3)')
     *    'destruction from RXPAR ',
     *    'dy = ',-y(nRXPAR,lprn)*y(n_Paraffin,lprn)*8.E-11*dt2
         write(6,118) ' Total change in ',ay(1,igas),ay(2,igas),
     *  ' is ',total,' percent; dy= ',dest(igas,lprn)+prod(igas,lprn)
         write(6,*)
        enddo ! igas
       endif  ! chem diags
 108  format(a10,2x,2a4)
 109  format(a55,7x,a6,e10.3)
 110  format(a68,e10.3)
 118  format(a17,2a4,a4,f10.0,a14,e12.3)
 119  format(a45,7x,f10.0,a5,e12.3)
c
      if(prnchg.and.J.eq.jprn.and.I.eq.iprn)write(6,'(a35,3(2x,i2))')
     * ' Total change by species at I, J, L',i,j,lprn
c
cc    Loop to calculate tracer changes
      do igas=1,ntm
       do L=1,maxl
c
         change(I,J,L,igas)=
     &   (dest(igas,L)+prod(igas,L))*dxyp(J)*AM(L,I,J)*
     &   bymass2vol(igas)/y(nM,L)
c
c        Set N2O5 to equilibrium when necessary (near ground,
c        N2O5 is thermally unstable, has a very short lifetime)
         if(igas.eq.n_N2O5.and.-dest(igas,L).ge.y(n_N2O5,L))then
           rnewval=(rr(52,I,J,L)*y(nNO3,L)*y(nNO2,L))/
     &     (rr(46,I,J,L)*y(nM,L)+ss(7,I,J,L))
           if(rnewval.lt.1.)rnewval=1.
           change(I,J,L,igas)=(rnewval-y(n_N2O5,L))*dxyp(J)*AM(L,I,J)*
     &     bymass2vol(igas)/y(nM,L)
c          endif
         endif
c
c        Conserve NOx with respect to N2O5
         if(igas.eq.n_NOx.and.-dest(n_N2O5,L).ge.y(n_N2O5,L))then
          rnewval=(rr(52,I,J,L)*y(nNO3,L)*y(nNO2,L))/
     &    (rr(46,I,J,L)*y(nM,L)+ss(7,I,J,L))
          change(I,J,L,igas)=
     &    change(I,J,L,igas)+2.*(y(n_N2O5,L)-rnewval)
     &    *dxyp(J)*AM(L,I,J)*bymass2vol(igas)/y(nM,L)
         endif
c
c        Set HO2NO2 to equil when necessary
         if(igas.eq.n_HO2NO2.and.-dest(igas,L).ge.y(n_HO2NO2,L))then
          rnewval=(rr(51,I,J,L)*y(nHO2,L)*y(nNO2,L))/
     *    (rr(18,I,J,L)*y(nOH,L)+rr(45,I,J,L)*y(nM,L)+ss(10,I,J,L)
     *    +ss(11,I,J,L))
          if(rnewval.lt.1.)rnewval=1.
          change(I,J,L,igas)=(rnewval-y(n_HO2NO2,L))*dxyp(J)*AM(L,I,J)*
     &    bymass2vol(igas)/y(nM,L)
         endif
c
c        Conserve NOx with respect to HO2NO2
         if(igas.eq.n_NOx.and.-dest(n_HO2NO2,L).ge.y(n_HO2NO2,L))then
          rnewval=(rr(51,I,J,L)*y(nHO2,L)*y(nNO2,L))/
     *    (rr(18,I,J,L)*y(nOH,L)+rr(45,I,J,L)*y(nM,L)+ss(10,I,J,L)
     *    +ss(11,I,J,L))
          change(I,J,L,igas)=
     &    change(I,J,L,igas)+(y(n_HO2NO2,L)-rnewval)*
     &    dxyp(J)*AM(L,I,J)*bymass2vol(igas)/y(nM,L)
         endif
c
c        Set PAN to equilibrium when necessary (near ground,
c         PAN is thermally unstable, has a very short lifetime)
         if(igas.eq.n_PAN.and.-dest(igas,L).ge.y(n_PAN,L))then
          rnewval=(rr(54,I,J,L)*y(nC2O3,L)*y(nNO2,L))/
     *    (rr(29,I,J,L)*y(nM,L)+ss(15,I,J,L))
          if(rnewval.lt.1.)rnewval=1.
          change(I,J,L,igas)=(rnewval-y(n_PAN,L))*dxyp(J)*AM(L,I,J)
     &    *bymass2vol(igas)/y(nM,L)
         endif
c
c        Conserve NOx with respect to PAN
         if(igas.eq.n_NOx.and.-dest(n_PAN,L).ge.y(n_PAN,L))then
          rnewval=(rr(54,I,J,L)*y(nC2O3,L)*y(nNO2,L))/
     *    (rr(29,I,J,L)*y(nM,L)+ss(15,I,J,L))
          if(rnewval.lt.1.)rnewval=1.
          change(I,J,L,igas)=change(I,J,L,igas)+(y(n_PAN,L)-rnewval)*
     &    dxyp(J)*AM(L,I,J)*bymass2vol(igas)/y(nM,L)
         endif
c
       end do ! L
      end do  ! igas
c
cc    Insure nitrogen conservation
      if(prnchg.and.J.eq.jprn.and.I.eq.iprn)then
       write(*,*) 'changes before nitrogen conservation routine'
       write(*,*) 'NOx, N2O5, HO2NO2, HNO3, PAN, AlkylNit'
       write(*,*) change(I,J,lprn,n_NOx),change(I,J,lprn,n_N2O5),
     & change(I,J,lprn,n_HO2NO2),change(I,J,lprn,n_HNO3),
     & change(I,J,lprn,n_PAN),change(I,J,lprn,n_AlkylNit)
      endif
c
      do L=1,maxl
c
c       First check for nitrogen loss > 100%
        if(-change(I,J,L,n_NOx).gt.trm(I,J,L,n_NOx))
     &  change(I,J,L,n_NOx)=1.-trm(I,J,L,n_NOx)
        if(-change(I,J,L,n_N2O5).gt.trm(I,J,L,n_N2O5))
     &  change(I,J,L,n_N2O5)=1.-trm(I,J,L,n_N2O5)
        if(-change(I,J,L,n_HO2NO2).gt.trm(I,J,L,n_HO2NO2))
     &  change(I,J,L,n_HO2NO2)=1.-trm(I,J,L,n_HO2NO2)
        if(-change(I,J,L,n_HNO3).gt.trm(I,J,L,n_HNO3))
     &  change(I,J,L,n_HNO3)=1.-trm(I,J,L,n_HNO3)
        if(-change(I,J,L,n_PAN).gt.trm(I,J,L,n_PAN))
     &  change(I,J,L,n_PAN)=1.-trm(I,J,L,n_PAN)
        if(-change(I,J,L,n_AlkylNit).gt.trm(I,J,L,n_AlkylNit))
     &  change(I,J,L,n_AlkylNit)=1.-trm(I,J,L,n_AlkylNit)
c
c       Next insure balance between dNOx and sum of dOthers
        sumN=(2.*change(I,J,L,n_N2O5))*mass2vol(n_N2O5)+
     *  (change(I,J,L,n_HNO3))*mass2vol(n_HNO3)+
     *  (change(I,J,L,n_HO2NO2))*mass2vol(n_HO2NO2)+
     *  (change(I,J,L,n_PAN))*mass2vol(n_PAN)+
     *  (change(I,J,L,n_AlkylNit))*mass2vol(n_AlkylNit)
        dNOx=change(I,J,L,n_NOx)*mass2vol(n_NOx)
        if(prnchg.and.J.eq.jprn.and.I.eq.iprn.and.L.eq.lprn)
     &  write(*,*) 'sumN, dNOx = ',sumN,dNOx
        ratio=-sumN/dNOx
        IF(ratio.le.0.999 .OR. ratio.ge.1.001) THEN
         if(dNOx.gt.0.)then
c         NOx being produced (net positive change)
          if (ratio.gt.1.)then
           sumD=0.
c          reduce N destruction to match NOx prodcution
           if(change(I,J,L,n_N2O5).lt.0.)   sumD=sumD+
     *     2.*change(I,J,L,n_N2O5)*mass2vol(n_N2O5)
           if(change(I,J,L,n_HO2NO2).lt.0.) sumD=sumD+
     *     change(I,J,L,n_HO2NO2)*mass2vol(n_HO2NO2)
           if(change(I,J,L,n_HNO3).lt.0.)   sumD=sumD+
     *     change(I,J,L,n_HNO3)*mass2vol(n_HNO3)
           if(change(I,J,L,n_PAN).lt.0.)    sumD=sumD+
     *     change(I,J,L,n_PAN)*mass2vol(n_PAN)
           if(change(I,J,L,n_AlkylNit).lt.0.)sumD=sumD+
     *     change(I,J,L,n_AlkylNit)*mass2vol(n_AlkylNit)
           newD=(sumN/ratio)+sumD-sumN
           ratioD=newD/sumD
           if(change(I,J,L,n_N2O5).lt.0.)    change(I,J,L,n_N2O5)=
     *     change(I,J,L,n_N2O5)    *ratioD
           if(change(I,J,L,n_HO2NO2).lt.0.)  change(I,J,L,n_HO2NO2)=
     *     change(I,J,L,n_HO2NO2)  *ratioD
           if(change(I,J,L,n_HNO3).lt.0.)    change(I,J,L,n_HNO3)=
     *     change(I,J,L,n_HNO3)    *ratioD
           if(change(I,J,L,n_PAN).lt.0.)     change(I,J,L,n_PAN)=
     *     change(I,J,L,n_PAN)     *ratioD
           if(change(I,J,L,n_AlkylNit).lt.0.)change(I,J,L,n_AlkylNit)=
     *     change(I,J,L,n_AlkylNit)*ratioD
          endif
c
          if (ratio.le.1..and.ratio.gt.0.)then
c          reduce NOx production to match N loss
           change(I,J,L,n_NOx)=change(I,J,L,n_NOx)*ratio
          endif
         else
c         NOx being destroyed (net change is negative)
          if (ratio.gt.1.)then
           sumP=0
c          reduce N production to match NOx loss
           if(change(I,J,L,n_N2O5).gt.0.)    sumP=sumP+
     *     2.*change(I,J,L,n_N2O5)*mass2vol(n_N2O5)
           if(change(I,J,L,n_HO2NO2).gt.0.)  sumP=sumP+
     *     change(I,J,L,n_HO2NO2)*mass2vol(n_HO2NO2)
           if(change(I,J,L,n_HNO3).gt.0.)    sumP=sumP+
     *     change(I,J,L,n_HNO3)*mass2vol(n_HNO3)
           if(change(I,J,L,n_PAN).gt.0.)     sumP=sumP+
     *     change(I,J,L,n_PAN)*mass2vol(n_PAN)
           if(change(I,J,L,n_AlkylNit).gt.0.)sumP=sumP+
     *     change(I,J,L,n_AlkylNit)*mass2vol(n_AlkylNit)
           newP=(sumN/ratio)+sumP-sumN
           ratioP=newP/sumP
           if(change(I,J,L,n_N2O5).gt.0)     change(I,J,L,n_N2O5)=
     *        change(I,J,L,n_N2O5)*ratioP
           if(change(I,J,L,n_HO2NO2).gt.0.)  change(I,J,L,n_HO2NO2)=
     *        change(I,J,L,n_HO2NO2)*ratioP
           if(change(I,J,L,n_HNO3).gt.0.)    change(I,J,L,n_HNO3)=
     *        change(I,J,L,n_HNO3)*ratioP
           if(change(I,J,L,n_PAN).gt.0.)     change(I,J,L,n_PAN)=
     *        change(I,J,L,n_PAN)*ratioP
           if(change(I,J,L,n_AlkylNit).gt.0.)change(I,J,L,n_AlkylNit)=
     *        change(I,J,L,n_AlkylNit)*ratioP
          endif
          if (ratio.le.1.and.ratio.gt.0)then
c          increase N production to match NOx loss (12/4/2001)
           sumP=0.
           if(change(I,J,L,n_N2O5).gt.0.)    sumP=sumP+
     *     2.*change(I,J,L,n_N2O5)*mass2vol(n_N2O5)
           if(change(I,J,L,n_HO2NO2).gt.0.)  sumP=sumP+
     *     change(I,J,L,n_HO2NO2)*mass2vol(n_HO2NO2)
           if(change(I,J,L,n_HNO3).gt.0.)    sumP=sumP+
     *     change(I,J,L,n_HNO3)*mass2vol(n_HNO3)
           if(change(I,J,L,n_PAN).gt.0.)     sumP=sumP+
     *     change(I,J,L,n_PAN)*mass2vol(n_PAN)
           if(change(I,J,L,n_AlkylNit).gt.0.)sumP=sumP+
     *     change(I,J,L,n_AlkylNit)*mass2vol(n_AlkylNit)
           newP=(sumN/ratio)+sumP-sumN
           ratioP=newP/sumP
           if(change(I,J,L,n_N2O5).gt.0.)    change(I,J,L,n_N2O5)=
     *     change(I,J,L,n_N2O5)*ratioP
           if(change(I,J,L,n_HO2NO2).gt.0.)  change(I,J,L,n_HO2NO2)=
     *     change(I,J,L,n_HO2NO2)*ratioP
           if(change(I,J,L,n_HNO3).gt.0.)    change(I,J,L,n_HNO3)=
     *     change(I,J,L,n_HNO3)*ratioP
           if(change(I,J,L,n_PAN).gt.0.)     change(I,J,L,n_PAN)=
     *     change(I,J,L,n_PAN)*ratioP
           if(change(I,J,L,n_AlkylNit).gt.0.)change(I,J,L,n_AlkylNit)=
     *     change(I,J,L,n_AlkylNit)*ratioP
          endif
         endif

        END IF ! skipped section above if ratio very close to one
C
        if(prnchg.and.J.eq.jprn.and.I.eq.iprn.and.L.eq.lprn)
     &  write(*,*) 'ratio for conservation =',ratio
C
      end do ! L
C
cc    Print chemical changes in a particular grid box if desired
      if(prnchg.and.J.eq.jprn.and.I.eq.iprn)then
        do igas=1,ntm
         changeA=change(I,J,Lprn,igas)*y(nM,lprn)*mass2vol(igas)*
     *   bydxyp(J)*byam(lprn,I,J)
         if(y(igas,lprn).eq.0.)then
            write(6,156) ay(1,igas),ay(2,igas),': ',changeA
     *            ,' molecules;  y=0'
         else
            write(6,155) ay(1,igas),ay(2,igas),': ',changeA
     *           ,' molecules produced; ',
     *   (100.*changeA)/y(igas,lprn),' percent of'
     *   ,y(igas,lprn),'(',1.E9*y(igas,lprn)/y(nM,lprn),' ppbv)'
         endif
c
         if(igas.eq.ntm)then
          write(6,'(a10,58x,e13.3,6x,f10.3,a5)')
     *    ' H2O     :',y(nH2O,LPRN),(y(nH2O,LPRN)/
     *    y(nM,LPRN))*1.E6,' ppmv'
          write(6,'(a10,58x,e13.3,6x,f10.3,a5)')
     *    ' CH3O2   :',yCH3O2(I,J,LPRN),(yCH3O2(I,J,LPRN)/
     *    y(nM,LPRN))*1.E9,' ppbv'
          write(6,'(a10,58x,e13.3,6x,f10.3,a5)')
     *    ' C2O3    :',y(nC2O3,LPRN),(y(nC2O3,LPRN)/y(nM,LPRN))*1.E9,
     *    ' ppbv'
          write(6,'(a10,58x,e13.3,6x,f10.3,a5)')
     *    ' XO2     :',y(nXO2,LPRN),(y(nXO2,LPRN)/y(nM,LPRN))*1.E9,
     *    ' ppbv'
          write(6,'(a10,58x,e13.3,6x,f10.3,a5)')
     *    ' XO2N    :',y(nXO2N,LPRN),(y(nXO2N,LPRN)/
     *    y(nM,LPRN))*1.E9,' ppbv'
          write(6,'(a10,58x,e13.3,6x,f10.3,a5)')
     *    ' RXPAR   :',y(nRXPAR,LPRN),(y(nRXPAR,LPRN)/
     *    y(nM,LPRN))*1.E9,' ppbv'
          write(6,'(a10,58x,e13.3,6x,f10.3,a5)')
     *    ' Aldehyde:',y(nAldehyde,LPRN),(y(nAldehyde,LPRN)/
     *    y(nM,LPRN))*1.E9,' ppbv'
          write(6,'(a10,58x,e13.3,6x,f10.3,a5)')
     *    ' ROR     :',y(nROR,LPRN),(y(nROR,LPRN)/
     *    y(nM,LPRN))*1.E9,' ppbv'
         endif
        enddo
      endif  !end of prnchg loop
c
c*** tracer masses & slopes are now updated in apply_tracer_3Dsource ***
      do igas=1,ntm
       do L=1,maxl
        if(change(I,J,L,igas).gt.1.E20 .OR.
     &    -change(I,J,L,igas).gt.trm(I,J,L,igas)) THEN
           WRITE(6,*)'>>change set to 0 in chemstep: I,J,L,igas,change'
     &    ,I,J,L,igas,change(I,J,L,igas)
           change(I,J,L,igas) = 0.
        endif
c
c Leaving this here to reinstate it later in some form. GSF 2/14/02 :
c         if(change(I,J,L,igas).le.1.E20) THEN
c           if(L.eq.1.and.igas.eq.1)then
c             changeA=(change(I,J,1,igas)*y(nM,1)*mass2vol(igas))*
c     *       bydxyp(J)*byam(1,I,J)
c             TAIJS(I,J,1)=TAIJS(I,J,1)+
c     *       1.E9*(changeA)/y(nM,L) !O3 change in ppbv
c           endif
c           if(igas.eq.ntm)then
c            if(y(nOH,L).gt.0..and.y(nOH,L).lt.1.E20)tempJLS(J,L,nOH)=
c     *       tempJLS(J,L,nOH)+y(nOH,L)
c            tempJLS(J,L,nH2O)=tempJLS(J,L,nH2O)+(y(nH2O,L)/y(nM,L))
c           endif
c         endif
c
       end do    ! L
      end do     ! igas

 155  format(1x,2a4,a2,e13.3,a21,f10.0,a11,2x,e13.3,3x,a1,f12.5,a6)
 156  format(1x,2a4,a2,e13.3,a16)
 157  format(1x,2a4,a32,e13.3,a11)
c
      return
      end SUBROUTINE chemstep
cc    __________________________________________________________________
      SUBROUTINE rates(maxl,I,J)
!@sum rates calculate reaction rates with present concentrations
!@auth Drew Shindell (modelEifications by Greg Faluvegi)
!@ver  1.0 (based on ds3ch4_chem_calc_jun1202_M23)
c
C**** GLOBAL parameters and variables:
C
      USE TRCHEM_Shindell_COM, only: nr,chemrate,photrate,rr,y,nn,dt2,
     &                          ss,ks,ny,dest,prod,JPPJ
c
      IMPLICIT NONE
c
C**** Local parameters and variables and arguments:
!@var kalt local dummy L-loop variable
!@var maxl passed highest chemistry level
!@var ireac,igas dummy loop variables
!@var I,J passed horixontal spatial indicies
      INTEGER kalt, ireac, igas, maxl
      INTEGER, INTENT(IN) :: I,J
c
      do kalt=1,maxl
        do ireac=1,nr
          chemrate(ireac,kalt)=rr(ireac,I,J,kalt)*y(nn(1,ireac),kalt)*
     &    y(nn(2,ireac),kalt)*dt2
        enddo
        do ireac=1,JPPJ
          photrate(ireac,kalt)=ss(ireac,I,J,kalt)*y(ks(ireac),kalt)*dt2
        enddo
c       Initialize change arrays
        do igas=1,ny
          dest(igas,kalt)=0.0
          prod(igas,kalt)=0.0
        enddo
      enddo
      return
      end SUBROUTINE rates

cc    __________________________________________________________________
      SUBROUTINE chem1(kdnr,maxl,numel,nn,ndnr,chemrate,dest,multip)
!@sum chem1 calculate chemical destruction/production
!@auth Drew Shindell (modelEifications by Greg Faluvegi)
!@ver  1.0 (based on ds3ch4_chem_calc_jun1202_M23)
c
C**** GLOBAL parameters and variables:
C
      USE TRCHEM_Shindell_COM, only: p_2, p_3, p_4, n_igas, numfam,nfam
c
      IMPLICIT NONE
c
C**** Local parameters and variables and arguments:
!@var maxl passed highest chemistry level
!@var numel first index of nn array (shouldn't code like this anymore?)
!@var kdnr kdnr,kpnr,kds, or kps    passed from chemstep
!@var nn nn,nnr,ks, or kss          passed from chemstep
!@var ndnr ndnr,npnr,nds, or nps    passed from chemstep
!@var chemrate chemrate or photrate passed from chemstep
!@var dest dest or prod             passed from chemstep
!@var multip -1 for destruction, +1 for production
!@var i,ireac,igas,ial,nbeg,nend dummy loop variables
!@var dk dummy variable
      INTEGER ireac,igas,ial,i,dk,nbeg,nend
      INTEGER, INTENT(IN)            :: maxl, numel,multip
      INTEGER, DIMENSION(p_4)        :: kdnr
      INTEGER, DIMENSION(numel,p_2)  :: nn       ! can't do anymore?
      INTEGER, DIMENSION(p_3)        :: ndnr
      REAL*8,  DIMENSION(p_2,maxl)   :: chemrate ! can't do anymore?
      REAL*8,  DIMENSION(n_igas,maxl):: dest     ! can't do anymore?
C
      ireac=0
c     Reactive families
      do igas=1,numfam
        dk=kdnr(igas+1)-kdnr(igas)
        if(dk.ge.1) then
          do i=1,dk
           ireac=ireac+1
           do ial=1,maxl
            if(nn(1,ndnr(ireac)).ge.nfam(igas).and.nn(1,ndnr(ireac)).lt.
     *       nfam(igas+1))then
             dest(igas,ial)=dest(igas,ial)+multip
     &       *chemrate(ndnr(ireac),ial)
c            Save change array for individual family elements
             dest(nn(1,ndnr(ireac)),ial)=dest(nn(1,ndnr(ireac)),ial)+
     *          multip*chemrate(ndnr(ireac),ial)
            endif
            if(numel.eq.2)then
             if(nn(2,ndnr(ireac)).ge.nfam(igas).and.nn(2,ndnr(ireac))
     *        .lt.nfam(igas+1))then
              dest(igas,ial)=dest(igas,ial)+
     *        multip*chemrate(ndnr(ireac),ial)
              dest(nn(2,ndnr(ireac)),ial)=dest(nn(2,ndnr(ireac)),ial)+
     *        multip*chemrate(ndnr(ireac),ial)
             endif
            endif
           end do ! ial
          end do  ! i
        end if
      end do      ! igas
C
c     INDIVIDUAL SPECIES
      nbeg=numfam+1
      nend=nfam(1)-1
      do igas=nbeg,nend
        dk=kdnr(igas+1)-kdnr(igas)
        if(dk.ge.1) then
          do i=1,dk
           ireac=ireac+1
           do ial=1,maxl
            dest(igas,ial)=dest(igas,ial)+multip*
     &      chemrate(ndnr(ireac),ial)
           end do
          end do
        end if
      end do
      return
      end SUBROUTINE chem1

cc    __________________________________________________________________
      SUBROUTINE chem1prn(kdnr,numel,nn,ndnr,chemrate,
     &                    index,multip,igas,total,maxl,I,J)
!@sum chem1prn for printing out the chemical reactions
!@auth Drew Shindell (modelEifications by Greg Faluvegi)
!@ver  1.0 (based on ds3ch4_chem_calc_jun1202_M23)
c
C**** GLOBAL parameters and variables:
C
      USE TRCHEM_Shindell_COM, only: ay, lprn, nfam, p_4, numfam, y,
     &                              p_2, p_3
c
      IMPLICIT NONE
c
C**** Local parameters and variables and arguments:
!@var kdnr kdnr,kpnr,kds, or kps from chemstep
!@var numel first index of nn array (shouldn't code like this anymore?)
!@var nn nn,nnr,ks, or kss from chemstep
!@var ndnr ndnr,npnr,nds, or nps from chemstep
!@var chemrate chemrate or photrate from chemstep
!@var index passed index to know which call this is... {1,2,3,4}
!@var multip 1 for production, -1 for destruction
!@var igas passed index for gas number
!@var total dummy summation
!@var maxl highest chemistry level
!@var I,J passed horizontal spatial indicies
!@var label character string for printing
!@var irec dummy loop variables
!@var per dummy temp variable
      INTEGER, INTENT(IN) :: igas,I,J,maxl,multip,index,numel
      INTEGER ireac
      REAL*8 total,per
      REAL*8,  DIMENSION(p_2,maxl)   :: chemrate !can't do this anymore?
      INTEGER, DIMENSION(p_3)        :: ndnr
      INTEGER, DIMENSION(numel,p_2)  :: nn       !can't do this anymore?
      INTEGER, DIMENSION(p_4)        :: kdnr
      character*17 label

c     skip this section during (within) family chemisty calls
      if(igas.le.numfam) then
c
c     FAMILIES
      if(kdnr(igas+1)-kdnr(igas).lt.1)goto40
      do ireac=kdnr(igas),kdnr(igas+1)-1
         if(index.le.2)label=' chem reaction # '
         if(index.gt.2)label=' phot reaction # '
         if(nn(1,ndnr(ireac)).ge.nfam(igas).and.nn(1,ndnr(ireac)).lt.
     *     nfam(igas+1))then
           per=0.
           if(y(igas,lprn).ne.0.)per=multip*100.*
     *     chemrate(ndnr(ireac),lprn)/y(igas,lprn)
           write(6,177) label,ndnr(ireac),' percent change from ',
     *     ay(1,nn(1,ndnr(ireac))),ay(2,nn(1,ndnr(ireac))),' = ',per,
     *     ' dy=',multip*chemrate(ndnr(ireac),lprn)
           total=total+per
         endif
         if(numel.eq.2)then
          if(nn(2,ndnr(ireac)).ge.nfam(igas).and.nn(2,ndnr(ireac)).lt.
     *     nfam(igas+1))then
           per=0.
           if(y(igas,lprn).ne.0.)per=multip*100.*
     *     chemrate(ndnr(ireac),lprn)/y(igas,lprn)
           write(6,177) label,ndnr(ireac),' percent change from ',
     *     ay(1,nn(2,ndnr(ireac))),ay(2,nn(2,ndnr(ireac))),' = ',per,
     *     ' dy=',multip*chemrate(ndnr(ireac),lprn)
           total=total+per
          endif
         endif  ! end of numel=2 loop
      end do  !end of ireac loop
      goto40

      end if ! non family chem only loop

c     INDIVIDUAL SPECIES
      if(kdnr(igas+1)-kdnr(igas).lt.1)goto40
      do ireac=kdnr(igas),kdnr(igas+1)-1
        if(index.le.2)label=' chem reaction # '
        if(index.gt.2)label=' phot reaction # '
c       skip same reaction if written twice
        if(ireac.eq.1 . or. ndnr(ireac).ne.ndnr(ireac-1) .or.
     &  total.eq.0.) then
         if(nn(1,ndnr(ireac)).eq.igas)then
           per=0.
           if(y(igas,lprn).ne.0.)per=100.*multip*
     *     chemrate(ndnr(ireac),lprn)/y(igas,lprn)
           write(6,106) label,ndnr(ireac),' percent change = ',per,
     *     ' dy=',multip*chemrate(ndnr(ireac),lprn)
           total=total+per
         endif
         if(numel.eq.2)then
          if(nn(2,ndnr(ireac)).eq.igas)then
           per=0.
           if(y(igas,lprn).ne.0.)per=100.*multip*
     *     chemrate(ndnr(ireac),lprn)/y(igas,lprn)
           write(6,106) label,ndnr(ireac),' percent change = ',per,
     *     ' dy=',multip*chemrate(ndnr(ireac),lprn)
           total=total+per
          endif
         endif  !end numel=2 loop
 106     format(a17,i3,a18,f10.0,a4,e12.3)
 177     format(a17,i3,a21,2a4,a3,f10.0,a4,e12.3)
  20    end if
  30  end do ! ireac
  40  continue
      return
      end SUBROUTINE chem1prn


text:
*
*
      subroutine denbin(nup,nb,denb)
*
*     calculate binary number density in the case of primordial binares
*     --------------------------------------------------------------------
*     density is calculated for binaries which position is inside r(lmax).
*     --------------------------------------------------------------------
*     lmax = nzon(nup)
*     ----------------
*
*
      include 'common.h'
*
cChanged following ftncheck
c      real*8 rb(nmax),denb(nmax),co
      real*8 rb(nmax),denb(nmax),co
*
      integer nup,lmax,i,im1,ib,nspan,nspan2,nmin1,nmax1,nb(nmax)
*
      nspan = 2
      nspan2 = 2*nspan
      co = 0.75d0*(nspan2 + 1)/pi
      lmax = nzst(nup)
      ib = 0
*
      do 10 i = 1,lmax
*
         nb(i) = 0
         im1 = iname(i)
         if(ikind(im1).eq.2) then
           ib = ib + 1
           rb(ib) = r(i)
           nb(i) = ib
           denb(ib) = 0.d0
         endif
*
 10   continue
*
      do 20 i = 1,ib
*
         if(i.le.nspan) then
           nmin1 = 1
           nmax1 = nspan2 + 1
           denb(i) = co/(rb(nmax1)**3 - rb(nmin1)**3)
*
         else
*
           if(i+nspan.gt.ib) then
c             nmin1 = 2*i - ib - nspan
c             denb(i) = co/(rb(ib)**3 - rb(nmin1)**3)
             nmin1 = ib - 2*nspan
             if(nmin1.le.0) then
               print*,' wrong-nmin0 i,ib,nmin,nspan =',i,ib,nmin,nspan
               call flush(6)
               nmin1 = 1
             endif
             nmax1 = ib
             denb(i) = co/(rb(nmax1)**3 - rb(nmin1)**3)
*
           else
*
             nmin1 = i - nspan
             if(nmin1.le.0) then
               print*,' wrong-nmin1 i,ib,nmin,nspan =',i,ib,nmin,nspan
               call flush(6)
               nmin1 = 1
             endif
             nmax1 = i + nspan
             denb(i) = co/(rb(nmax1)**3 - rb(nmin1)**3)
*
           endif
*
         endif
*
 20      continue
*
*       maybe is good to smooth the binary density
*
         return
*
         end
*
*

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse