text:
      Complex*16 FUNCTION FPIBEL(W,flpar)
C****************************************************************************
C                       PION FORM FACTOR FROM BELLE
C
C            Rho(770)+rho(1450)+rho(1700) model is used as in:
C             M.Fujikawa et al., Phys.Rev.D 78 (2008) 072006
C                          -------------------
C       * flpar: = 0 - (all free) fit result for par(1...11)
C                = 1 - (par(1)=F_pi(0)=1-fixed) fit result for par(2...11)
C
C       * parameter : par(1)=overall norm. factor (=F_pi(0))
C                   : par(2)=rho(770) mass
C                   : par(3)=rho(770) width
C                   : par(4)=rho(1450) mass
C                   : par(5)=rho(1450) width
C                   : par(6)=rho(1450) admixture abs. value (|beta|)
C                   : par(7)=rho(1450) admixture phase (phi=arg(beta))
C                   : par(8)=rho(1700) mass
C                   : par(9)=rho(1700) width
C                   : par(10)=rho(1700) admixture abs. value (|gamma|)
C                   : par(11)=rho(1700) admixture phase (phi3=arg(gamma))
C
C       * x:=s=(M_pipi0)^2=W^2
C
C       * Notice: the following code was extracted (and checked) directly
C                 from the Belle fit function with minor cosmetic changes
C                 (this is internal comment of Belle)
C
C       * Notice: the following code was extracted (and checked) directly
C                 from  tauola_with_belle_fpi.F sent by D. Epifanov and
C                 it is an identical copy of
C                 Complex*16 FUNCTION FPIBEL(W,flpar)
C                 used in tauola_with_belle_fpi.F
C
C       * Called: CURR_PIPI0  if (FF2PIRHO = 2  or FF2PIRHO =3)
C                 in ../RChL-currents/value_parameter.f
C                 FF2PIRHO = 2 is for flpar = 0
C                 FF2PIRHO = 3 is for flpar = 1
C
C****************************************************************************
      implicit none
      integer flpar
      real W
      real*8 x,par(11),pi,pimas,pi0mas,taumas
      real*8 beta,berho1,berho2,berho3
      real*8 ps,prho1,prho2,prho3
      real*8 gamma1,gamma2,gamma3
      real*8 hs,h1,h2,h3,dhds1,dhds2,dhds3
      real*8 d1,d2,d3,fs1,fs2,fs3
      real*8 a1,b1,c1,bw_re1,bw_im1
      real*8 a2,b2,c2,bw_re2,bw_im2
      real*8 a3,b3,c3,bw_re3,bw_im3
      real*8 sinphi,cosphi,sinphi3,cosphi3
      complex*16 mbw1,mbw2,mbw3,mff,mephi,mephi3
      parameter(pi=3.141592653589)

      x=dble(W**2)
C------------- PARAMETERS --------------------
C     Particle mass (unit: GeV)

      pimas  = 0.1395702
      pi0mas = 0.13498
      taumas = 1.77699

C---------------------------------------------
C Opt. par. from Table VII of PRD78 (2008) 072006

      if(flpar.eq.0) then
         par(1) = 1.02
         par(2) = 0.7749
         par(3) = 0.1486
         par(4) = 1.428
         par(5) = 0.413
         par(6) = 0.13
         par(7) = 197
         par(8) = 1.694
         par(9) = 0.135
         par(10)= 0.028
         par(11)= -3
      else
         par(1) = 1.0
         par(2) = 0.7746
         par(3) = 0.1481
         par(4) = 1.446
         par(5) = 0.434
         par(6) = 0.15
         par(7) = 202
         par(8) = 1.728
         par(9) = 0.164
         par(10)= 0.037
         par(11)= 24
      endif

C========= Beta (s, rho(770), rho(1450), rho(1700)) ==================

       beta  =sqrt( (1.0 - (pimas - pi0mas)**2/x)
     +             *(1.0 - (pimas + pi0mas)**2/x) )

c      ----- rho(770) -----
       berho1=sqrt( ( 1.0 - (pimas - pi0mas)**2/(par(2)*par(2)) )
     +             *( 1.0 - (pimas + pi0mas)**2/(par(2)*par(2)) ) )

c      ----- rho(1450) -----
       berho2=sqrt( ( 1.0 - (pimas - pi0mas)**2/(par(4)*par(4)) )
     +             *( 1.0 - (pimas + pi0mas)**2/(par(4)*par(4)) ) )

c      ----- rho(1700) -----
       berho3=sqrt( ( 1.0 - (pimas - pi0mas)**2/(par(8)*par(8)) )
     +             *( 1.0 - (pimas + pi0mas)**2/(par(8)*par(8)) ) )

C========= Momentum (s, rho(770), rho(1450), rho(1700)) ==============
C
C       ps    : momentum of pi at s
C       prho1 : momentum of pi at s=rho(770)**2
C       prho2 : momentum of pi at s=rho(1450)**2
C       prho3 : momentum of pi at s=rho(1700)**2

       ps  = 0.5 * sqrt(x) * beta
c      ----- rho(770) -----
       prho1 = 0.5* par(2) * berho1
c      ----- rho(1450) -----
       prho2 = 0.5* par(4) * berho2
c      ----- rho(1700) -----
       prho3 = 0.5* par(8) * berho3

C========= Width (rho(770), rho(1450), rho(1700)) ====================

C      ----- rho(770) -----
       gamma1 = par(3) * ( par(2)*par(2)/x ) * (ps/prho1)**3
C      ----- rho(1450) -----
       gamma2 = par(5) * ( par(4)*par(4)/x ) * (ps/prho2)**3
C      ----- rho(1700) -----
       gamma3 = par(9) * ( par(8)*par(8)/x ) * (ps/prho3)**3

C============== h(rho(770), rho(1450), rho(1700)) ====================
C
C       hs : h(s)...s(s=x)
C       h1 : h(s) for s=rho(770)**2
C       h2 : h(s) for s=rho(1450)**2
C       h3 : h(s) for s=rho(1700)**2

        hs = (2.0/pi)*(ps/sqrt(x))*
     +       log( (sqrt(x) + 2.0*ps)/(2.0*pimas) )
C      ----- rho(770) -----
        h1 = (2.0/pi)*(prho1/par(2))*
     +       log( (par(2) + 2.0*prho1)/(2.0*pimas) )
C      ----- rho(1450) -----
        h2 = (2.0/pi)*(prho2/par(4))*
     +       log( (par(4) + 2.0*prho2)/(2.0*pimas) )
C      ----- rho(1700) -----
        h3 = (2.0/pi)*(prho3/par(8))*
     +       log( (par(8) + 2.0*prho3)/(2.0*pimas) )

C============ dhds(rho(770), rho(1450), rho(1700)) ===================
C
C       dhds = dh/ds = h(M_rho)[1/(8*p(M_rho)^2) - 1/(2M_rho^2)
C                                              + 1/(2pi*M_rho^2)
C       dhds1 : dh/ds for M_rho=rho(770)
C       dhds2 : dh/ds for M_rho=rho(1450)
C       dhds3 : dh/ds for M_rho=rho(1700)

C      ----- rho(770) -----
        dhds1 = h1*( 1.0/(8.0*prho1**2) - 1.0/(2.0*par(2)**2) )
     +                       +  1.0/(2.0*pi*par(2)**2)
C      ----- rho(1450) -----
        dhds2 = h2*( 1.0/(8.0*prho2**2) - 1.0/(2.0*par(4)**2) )
     +                       +  1.0/(2.0*pi*par(4)**2)
C      ----- rho(1700) -----
        dhds3 = h3*( 1.0/(8.0*prho3**2) - 1.0/(2.0*par(8)**2) )
     +                       +  1.0/(2.0*pi*par(8)**2)

C============ d(rho(770), rho(1450), rho(1700)) ======================
C
C       d = 3/pi * (m_pi^2/p(M_rho)^2) *
C            ln {(M_rho+2*p(M_rho))/2M_pi} + M_rho/(2pi*p(M_rho))
C                               - (m_pi^2 * M_rho)/(pi*p(M_rho)^3)
C       d1 : d for M_rho=rho(770)
C       d2 : d for M_rho=rho(1450)
C       d2 : d for M_rho=rho(1700)

C      ----- rho(770) -----
        d1 = (3.0/pi)*(pimas**2/prho1**2)
     +         * log((par(2) + 2.0*prho1)/(2.0*pimas))
     +         + par(2)/(2.0*pi*prho1)
     +         - ((pimas**2)*par(2))/(pi*prho1**3)
C      ----- rho(1450) -----
        d2 = (3.0/pi)*(pimas**2/prho2**2)
     +         * log((par(4) + 2.0*prho2)/(2.0*pimas))
     +         + par(4)/(2.0*pi*prho2)
     +         - ((pimas**2)*par(4))/(pi*prho2**3)
C      ----- rho(1700) -----
        d3 = (3.0/pi)*(pimas**2/prho3**2)
     +         * log((par(8) + 2.0*prho3)/(2.0*pimas))
     +         + par(8)/(2.0*pi*prho3)
     +         - ((pimas**2)*par(8))/(pi*prho3**3)

C================ f(s) (rho(770), rho(1450), rho(1700) ) ==========
C
C       f(s) = gamma-rho * M_rhp**2/p(M_rho)[ p(s)^2(h(s)-h(M_rho)
C                                 + (M_rho^2 -s)p(M_rho)^2 * dh/ds]
C       fs1 : f(s) for M_rho=rho(770)
C       fs2 : f(s) for M_rho=rho(1450)
C       fs3 : f(s) for M_rho=rho(1700)
C
C      ----- rho(770) -----
        fs1 = par(3) * (par(2)**2/prho1**3) *
     +        ( ps**2 *(hs - h1) + (par(2)**2 -x)*prho1**2 * dhds1 )
C      ----- rho(1450) -----
        fs2 = par(5) * (par(4)**2/prho2**3) *
     +        ( ps**2 *(hs - h2) + (par(4)**2 -x)*prho2**2 * dhds2 )
C      ----- rho(1700) -----
        fs3 = par(9) * (par(8)**2/prho3**3) *
     +        ( ps**2 *(hs - h3) + (par(8)**2 -x)*prho3**2 * dhds3 )

C====== BW form G&S model (rho(770) + rho(1450) + rho(1700)) ======
C    bw_re - real part ; bw_im - imaginary part

C      ----- rho(770) -----
       a1=par(2)*par(2) - x + fs1
       b1=sqrt(x)*gamma1
       c1=par(2)*par(2) + d1*par(2)*par(3)

       bw_re1= (a1*c1)/(a1**2 + b1**2)
       bw_im1= (b1*c1)/(a1**2 + b1**2)
       mbw1=dcmplx(bw_re1,bw_im1)

c      ----- rho(1450) -----
       a2=par(4)*par(4) - x + fs2
       b2=sqrt(x)*gamma2
       c2=par(4)*par(4)+ d2*par(4)*par(5)

       bw_re2= (a2*c2)/(a2**2 + b2**2)
       bw_im2= (b2*c2)/(a2**2 + b2**2)
       mbw2=dcmplx(bw_re2,bw_im2)

c      ----- rho(1700) -----
       a3=par(8)*par(8) - x + fs3
       b3=sqrt(x)*gamma3
       c3=par(8)*par(8)+ d3*par(8)*par(9)

       bw_re3= (a3*c3)/(a3**2 + b3**2)
       bw_im3= (b3*c3)/(a3**2 + b3**2)
       mbw3=dcmplx(bw_re3,bw_im3)

C====== admixtures of rho(1450) and rho(1700) ======

C --------- exp(i*arg(beta)) ----------
       sinphi=sin(par(7)*(pi/180.0))  ! use degree
       cosphi=cos(par(7)*(pi/180.0))  ! use degree
       mephi=dcmplx(cosphi,sinphi)

C --------- exp(i*arg(gamma)) ---------
       sinphi3=sin(par(11)*(pi/180.0))  ! use degree
       cosphi3=cos(par(11)*(pi/180.0))  ! use degree
       mephi3=dcmplx(cosphi3,sinphi3)

C================= Form factor =====================

       mff = 1.0/(1.0+par(6)*mephi+par(10)*mephi3)
       FPIBEL=par(1)*mff*(mbw1+par(6)*mephi*mbw2+par(10)*mephi3*mbw3)
      RETURN
      END

max_stars_repo_name: AllaMaevskaya/AliRoot

max_issues_repo_name: AllaMaevskaya/AliRoot

max_forks_repo_name: AllaMaevskaya/AliRoot


text:
lvalue/test.cpp
lvalue/tb.cpp
lvalue/monitor.cpp
lvalue/main.cpp

ext: f


!-------------------------------------------------------------------------------
!>
!! Operator module
!!
!! @par Description
!!         This module contains the subroutines for differential oeprators.
!!
!! @author  H.Tomita
!!
!! @par History
!! @li      2004-02-17 (H.Tomita)    Imported from igdc-4.33
!! @li      2006-04-17 (H.Tomita)    Add sub[OPRT_divergence]
!! @li      2006-08-11 (H.Tomita)    Implementatio of miura scheme(2004) with Thurbuen(1996)'s limiter.
!!                                   Add sub[OPRT_divergence2]
!!                                       sub[OPRT_divergence2_prep]
!!                                       sub[OPRT_divergence2_all].
!! @li      2006-08-22 (Y.Niwa)      divide the rows for the calc. of clap and clap_pl due to the rule of XLF.
!! @li      2007-01-26 (H.Tomita)    Optimization of sub[oprt_diffusion].
!! @li      2007-11-28 (T.Mitsui)    bugfix in oprt_divergence2, _all
!! @li      2008-01-24 (Y.Niwa)      add OPRT_divergence2{_prep,,_all}_rev
!! @li      2008-04-28 (T.Mitsui)    bug fix in OPRT_divergence2{_all}_rev
!! @li      2009-09-04 (H.Taniguchi) bug fix in OPRT_divergence2{_all}_rev Zero clear of wrk[_pl] is needed.
!! @li      2010-06-08 (S.Iga)       new grid is implemented (see, string XTMS)
!! @li      2011-09-27 (T.Seiki)     merge optimization by RIST and M.Terai
!! @li      2012-01-17 (M.Terai)     update optimization(case6) in div2rev
!! @li      2012-05-01 (T.Yamaura)   bug fix in div2rev
!! @li      2012-06-28 (M.Terai)     Removed wrapper subroutine to invoked the operators directly macro var.
!!
!<
module mod_oprt
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use mod_precision
  use mod_misc
!ESC!  use mod_stdio
!ESC!  use mod_prof
!ESC!  use mod_adm, only: &
!ESC!     ADM_nxyz,           &
!ESC!     TI    => ADM_TI,    &
!ESC!     TJ    => ADM_TJ,    &
!ESC!     AI    => ADM_AI,    &
!ESC!     AIJ   => ADM_AIJ,   &
!ESC!     AJ    => ADM_AJ,    &
!ESC!     K0    => ADM_KNONE, &
!ESC!     vlink => ADM_vlink, &
!ESC!     ADM_lall,           &
!ESC!     ADM_lall_pl,        &
!ESC!     ADM_kall,           &
!ESC!     ADM_jall,           &
!ESC!     ADM_iall,           &
!ESC!     ADM_gall_pl
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: OPRT_diffusion

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine OPRT_diffusion( &
       dscl,      dscl_pl,      &
       scl,       scl_pl,       &
       kh,        kh_pl,        &
       coef_intp, coef_intp_pl, &
       coef_diff, coef_diff_pl  )
!ESC!    use mod_adm, only: &
!ESC!       ADM_have_pl,  &
!ESC!       ADM_have_sgp, &
!ESC!       ADM_imin,     &
!ESC!       ADM_imax,     &
!ESC!       ADM_jmin,     &
!ESC!       ADM_jmax,     &
!ESC!       ADM_gslf_pl,  &
!ESC!       ADM_gmin_pl,  &
!ESC!       ADM_gmax_pl
!ESC!    use mod_grd, only: &
!ESC!       XDIR => GRD_XDIR, &
!ESC!       YDIR => GRD_YDIR, &
!ESC!       ZDIR => GRD_ZDIR
    implicit none

    real(RP), intent(out) :: dscl        (ADM_iall,ADM_jall,ADM_kall,ADM_lall   )
    real(RP), intent(out) :: dscl_pl     (ADM_gall_pl      ,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: scl         (ADM_iall,ADM_jall,ADM_kall,ADM_lall   )
    real(RP), intent(in)  :: scl_pl      (ADM_gall_pl      ,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: kh          (ADM_iall,ADM_jall,ADM_kall,ADM_lall   )
    real(RP), intent(in)  :: kh_pl       (ADM_gall_pl      ,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: coef_intp   (ADM_iall,ADM_jall,1:3,    ADM_nxyz,TI:TJ,ADM_lall   )
    real(RP), intent(in)  :: coef_intp_pl(ADM_gall_pl      ,1:3,    ADM_nxyz,      ADM_lall_pl)
    real(RP), intent(in)  :: coef_diff   (ADM_iall,ADM_jall,1:6,    ADM_nxyz,      ADM_lall   )
    real(RP), intent(in)  :: coef_diff_pl(              1:ADM_vlink,ADM_nxyz,      ADM_lall_pl)

    real(RP) :: vt   (ADM_iall,ADM_jall,ADM_nxyz,TI:TJ)
    real(RP) :: vt_pl(ADM_gall_pl      ,ADM_nxyz)

    integer  :: imin, imax, jmin, jmax, kall, lall
    integer  :: ij, ijp1, ijm1

    integer  :: i, j, k, l, n, v, d
    !---------------------------------------------------------------------------

    call DEBUG_rapstart('OPRT_diffusion')

    imin = ADM_imin
    imax = ADM_imax
    jmin = ADM_jmin
    jmax = ADM_jmax
    kall = ADM_kall
    lall = ADM_lall

    !$omp parallel default(none),private(i,j,k,l,d), &
    !$omp shared(imin,imax,jmin,jmax,kall,lall,ADM_have_sgp,dscl,scl,kh,vt,coef_intp,coef_diff)
    do l = 1, lall
    do k = 1, kall
       !$omp do schedule(static) collapse(2)
       do d = XDIR, ZDIR
       do j = jmin-1, jmax
       do i = imin-1, imax
          vt(i,j,d,TI) = ( ( + 2.0_RP * coef_intp(i,j,1,d,TI,l) &
                             - 1.0_RP * coef_intp(i,j,2,d,TI,l) &
                             - 1.0_RP * coef_intp(i,j,3,d,TI,l) ) * scl(i  ,j  ,k,l) &
                         + ( - 1.0_RP * coef_intp(i,j,1,d,TI,l) &
                             + 2.0_RP * coef_intp(i,j,2,d,TI,l) &
                             - 1.0_RP * coef_intp(i,j,3,d,TI,l) ) * scl(i+1,j  ,k,l) &
                         + ( - 1.0_RP * coef_intp(i,j,1,d,TI,l) &
                             - 1.0_RP * coef_intp(i,j,2,d,TI,l) &
                             + 2.0_RP * coef_intp(i,j,3,d,TI,l) ) * scl(i+1,j+1,k,l) &
                         ) / 3.0_RP
       enddo
       enddo
       enddo
       !$omp end do nowait

       !$omp do schedule(static) collapse(2)
       do d = XDIR, ZDIR
       do j = jmin-1, jmax
       do i = imin-1, imax
          vt(i,j,d,TJ) = ( ( + 2.0_RP * coef_intp(i,j,1,d,TJ,l) &
                             - 1.0_RP * coef_intp(i,j,2,d,TJ,l) &
                             - 1.0_RP * coef_intp(i,j,3,d,TJ,l) ) * scl(i  ,j  ,k,l) &
                         + ( - 1.0_RP * coef_intp(i,j,1,d,TJ,l) &
                             + 2.0_RP * coef_intp(i,j,2,d,TJ,l) &
                             - 1.0_RP * coef_intp(i,j,3,d,TJ,l) ) * scl(i+1,j+1,k,l) &
                         + ( - 1.0_RP * coef_intp(i,j,1,d,TJ,l) &
                             - 1.0_RP * coef_intp(i,j,2,d,TJ,l) &
                             + 2.0_RP * coef_intp(i,j,3,d,TJ,l) ) * scl(i  ,j+1,k,l) &
                         ) / 3.0_RP
       enddo
       enddo
       enddo
       !$omp end do

       if ( ADM_have_sgp(l) ) then ! pentagon
          !$omp master
          vt(imin-1,jmin-1,:,TI) = vt(imin,jmin-1,:,TJ)
          !$omp end master
       endif

       !$omp do schedule(static)
       do j = jmin, jmax
       do i = imin, imax
          dscl(i,j,k,l) = ( coef_diff(i,j,1,XDIR,l) * ( vt(i  ,j  ,XDIR,TI) + vt(i  ,j  ,XDIR,TJ) ) &
                          + coef_diff(i,j,1,YDIR,l) * ( vt(i  ,j  ,YDIR,TI) + vt(i  ,j  ,YDIR,TJ) ) &
                          + coef_diff(i,j,1,ZDIR,l) * ( vt(i  ,j  ,ZDIR,TI) + vt(i  ,j  ,ZDIR,TJ) ) &
                          ) * 0.5_RP * ( kh(i  ,j  ,k,l) + kh(i+1,j+1,k,l) )
       enddo
       enddo
       !$omp end do

       !$omp do schedule(static)
       do j = jmin, jmax
       do i = imin, imax
          dscl(i,j,k,l) = dscl(i,j,k,l) &
                        + ( coef_diff(i,j,2,XDIR,l) * ( vt(i  ,j  ,XDIR,TJ) + vt(i-1,j  ,XDIR,TI) ) &
                          + coef_diff(i,j,2,YDIR,l) * ( vt(i  ,j  ,YDIR,TJ) + vt(i-1,j  ,YDIR,TI) ) &
                          + coef_diff(i,j,2,ZDIR,l) * ( vt(i  ,j  ,ZDIR,TJ) + vt(i-1,j  ,ZDIR,TI) ) &
                          ) * 0.5_RP * ( kh(i  ,j  ,k,l) + kh(i  ,j+1,k,l) )
       enddo
       enddo
       !$omp end do

       !$omp do schedule(static)
       do j = jmin, jmax
       do i = imin, imax
          dscl(i,j,k,l) = dscl(i,j,k,l) &
                        + ( coef_diff(i,j,3,XDIR,l) * ( vt(i-1,j  ,XDIR,TI) + vt(i-1,j-1,XDIR,TJ) ) &
                          + coef_diff(i,j,3,YDIR,l) * ( vt(i-1,j  ,YDIR,TI) + vt(i-1,j-1,YDIR,TJ) ) &
                          + coef_diff(i,j,3,ZDIR,l) * ( vt(i-1,j  ,ZDIR,TI) + vt(i-1,j-1,ZDIR,TJ) ) &
                          ) * 0.5_RP * ( kh(i-1,j  ,k,l) + kh(i  ,j  ,k,l) )
       enddo
       enddo
       !$omp end do

       !$omp do schedule(static)
       do j = jmin, jmax
       do i = imin, imax
          dscl(i,j,k,l) = dscl(i,j,k,l) &
                        + ( coef_diff(i,j,4,XDIR,l) * ( vt(i-1,j-1,XDIR,TJ) + vt(i-1,j-1,XDIR,TI) ) &
                          + coef_diff(i,j,4,YDIR,l) * ( vt(i-1,j-1,YDIR,TJ) + vt(i-1,j-1,YDIR,TI) ) &
                          + coef_diff(i,j,4,ZDIR,l) * ( vt(i-1,j-1,ZDIR,TJ) + vt(i-1,j-1,ZDIR,TI) ) &
                          ) * 0.5_RP * ( kh(i-1,j-1,k,l) + kh(i  ,j  ,k,l) )
       enddo
       enddo
       !$omp end do

       !$omp do schedule(static)
       do j = jmin, jmax
       do i = imin, imax
          dscl(i,j,k,l) = dscl(i,j,k,l) &
                        + ( coef_diff(i,j,5,XDIR,l) * ( vt(i-1,j-1,XDIR,TI) + vt(i  ,j-1,XDIR,TJ) ) &
                          + coef_diff(i,j,5,YDIR,l) * ( vt(i-1,j-1,YDIR,TI) + vt(i  ,j-1,YDIR,TJ) ) &
                          + coef_diff(i,j,5,ZDIR,l) * ( vt(i-1,j-1,ZDIR,TI) + vt(i  ,j-1,ZDIR,TJ) ) &
                          ) * 0.5_RP * ( kh(i  ,j-1,k,l) + kh(i  ,j  ,k,l) )
       enddo
       enddo
       !$omp end do

       !$omp do schedule(static)
       do j = jmin, jmax
       do i = imin, imax
          dscl(i,j,k,l) = dscl(i,j,k,l) &
                        + ( coef_diff(i,j,6,XDIR,l) * ( vt(i  ,j-1,XDIR,TJ) + vt(i  ,j  ,XDIR,TI) ) &
                          + coef_diff(i,j,6,YDIR,l) * ( vt(i  ,j-1,YDIR,TJ) + vt(i  ,j  ,YDIR,TI) ) &
                          + coef_diff(i,j,6,ZDIR,l) * ( vt(i  ,j-1,ZDIR,TJ) + vt(i  ,j  ,ZDIR,TI) ) &
                          ) * 0.5_RP * ( kh(i  ,j  ,k,l) + kh(i+1,j  ,k,l) )
       enddo
       enddo
       !$omp end do

       !$omp workshare
       dscl(:,jmin-1,k,l) = 0.0_RP
       dscl(:,jmax+1,k,l) = 0.0_RP
       dscl(imin-1,:,k,l) = 0.0_RP
       dscl(imax+1,:,k,l) = 0.0_RP
       !$omp end workshare
    enddo ! k loop
    enddo ! l loop
    !$omp end parallel

    if ( ADM_have_pl ) then
       n = ADM_gslf_pl

       do l = 1, ADM_lall_pl
       do k = 1, ADM_kall

          do v = ADM_gmin_pl, ADM_gmax_pl
             ij   = v
             ijp1 = v + 1
             if( ijp1 == ADM_gmax_pl+1 ) ijp1 = ADM_gmin_pl

             do d = 1, ADM_nxyz
                vt_pl(ij,d) = ( ( + 2.0_RP * coef_intp_pl(v,1,d,l) &
                                  - 1.0_RP * coef_intp_pl(v,2,d,l) &
                                  - 1.0_RP * coef_intp_pl(v,3,d,l) ) * scl_pl(n   ,k,l) &
                              + ( - 1.0_RP * coef_intp_pl(v,1,d,l) &
                                  + 2.0_RP * coef_intp_pl(v,2,d,l) &
                                  - 1.0_RP * coef_intp_pl(v,3,d,l) ) * scl_pl(ij  ,k,l) &
                              + ( - 1.0_RP * coef_intp_pl(v,1,d,l) &
                                  - 1.0_RP * coef_intp_pl(v,2,d,l) &
                                  + 2.0_RP * coef_intp_pl(v,3,d,l) ) * scl_pl(ijp1,k,l) &
                              ) / 3.0_RP
             enddo
          enddo

          dscl_pl(:,k,l) = 0.0_RP

          do v = ADM_gmin_pl, ADM_gmax_pl
             ij   = v
             ijm1 = v - 1
             if( ijm1 == ADM_gmin_pl-1 ) ijm1 = ADM_gmax_pl ! cyclic condition

             dscl_pl(n,k,l) = dscl_pl(n,k,l) &
                            + ( coef_diff_pl(v-1,XDIR,l) * ( vt_pl(ijm1,XDIR) + vt_pl(ij,XDIR) ) &
                              + coef_diff_pl(v-1,YDIR,l) * ( vt_pl(ijm1,YDIR) + vt_pl(ij,YDIR) ) &
                              + coef_diff_pl(v-1,ZDIR,l) * ( vt_pl(ijm1,ZDIR) + vt_pl(ij,ZDIR) ) &
                              ) * 0.5_RP * ( kh_pl(n,k,l) + kh_pl(ij,k,l) )
          enddo

       enddo
       enddo
    else
       dscl_pl(:,:,:) = 0.0_RP
    endif

    call DEBUG_rapend('OPRT_diffusion')

    return
  end subroutine OPRT_diffusion

end module mod_oprt
!-------------------------------------------------------------------------------


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

C   --*** PREB1 *** (BLOT) Display database element block information
C   --   Written by Amy Gilkey - revised 01/05/88
C   --
C   --PREB1 displays the element block internal information (either the
C   --connectivity or the attributes.
C   --
C   --Parameters:
C   --   OPTION - IN - '*' to print all, else print options:
C   --      'C' to print connectivity
C   --      'A' to print attributes
C   --   NOUT - IN - the output file, <=0 for standard
C   --   IEL0 - IN - the element offset of the elements in the block
C   --   NLISEL - IN - the number of selected elements
C   --   LISEL - IN - the indices of the selected elements
C   --   NLINK - IN - the number of nodes per element
C   --   NATR - IN - the number of attributes per element
C   --   LINK - IN - the connectivity array for this block
C   --   ATRIB - IN - the attribute array for this block

      CHARACTER*(*) OPTION
      INTEGER LISEL(*)
      INTEGER LINK(NLINK,*)
      REAL ATRIB(NATR,*)
      INTEGER MAPEL(*), MAPND(*)
      LOGICAL ISABRT
      LOGICAL DOCONN, DOATR

      IF (NLISEL .LE. 0) RETURN

      DOCONN = ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'C') .GT. 0))
     &   .AND. (NLINK .GT. 0)
      DOATR  = ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'A') .GT. 0))
     &   .AND. (NATR .GT. 0)
      IF (.NOT. (DOCONN .OR. DOATR)) RETURN

      IF (NOUT .GT. 0) THEN
         IF (DOCONN .AND. DOATR) THEN
            WRITE (NOUT, 10000) 'Connectivity and Attributes'
         ELSE IF (DOCONN) THEN
            WRITE (NOUT, 10000) 'Connectivity'
         ELSE IF (DOATR) THEN
            WRITE (NOUT, 10000) 'Attributes'
         END IF
      ELSE
         IF (DOCONN .AND. DOATR) THEN
            WRITE (*, 10000) 'Connectivity and Attributes'
         ELSE IF (DOCONN) THEN
            WRITE (*, 10000) 'Connectivity'
         ELSE IF (DOATR) THEN
            WRITE (*, 10000) 'Attributes'
         END IF
      END IF

      DO 100 IX = 1, NLISEL
         IEL = LISEL(IX)
         NE = IEL - IEL0
         IF (ISABRT ()) RETURN
         IF (DOCONN) THEN
            IF (NOUT .GT. 0) THEN
               WRITE (NOUT, 10010, IOSTAT=IDUM)
     &          NE, MAPEL(IEL), (MAPND(LINK(I,NE)), I=1,NLINK)
            ELSE
               WRITE (*, 10010, IOSTAT=IDUM)
     &            NE, MAPEL(IEL), (MAPND(LINK(I,NE)), I=1,NLINK)
            END IF
         END IF
         IF (DOATR) THEN
            IF (NOUT .GT. 0) THEN
               WRITE (NOUT, 10020, IOSTAT=IDUM)
     &            NE, MAPEL(IEL), (ATRIB(I,NE), I=1,NATR)
            ELSE
               WRITE (*, 10020, IOSTAT=IDUM)
     &            NE, MAPEL(IEL), (ATRIB(I,NE), I=1,NATR)
            END IF
         END IF
  100 CONTINUE

      RETURN

10000  FORMAT (1X, '   Local ID  Global ID       ', A)
10010  FORMAT (1X, I11, I11, 5X, 8I11, :, /,
     &   (18X, 8I11))
10020  FORMAT (1X, I11, I11, 3X, 4 (2X, 1PE11.4), :, /,
     &   (16X, 4 (2X, 1PE11.4)))
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


program devnull
use, intrinsic:: iso_fortran_env, only: int64
implicit none (type, external)

integer(int64) :: tic,toc,rate
character(*), parameter :: nulunix='/dev/null', nulwin='NUL',fout='out.txt'
integer, parameter :: N=1000
integer :: ios,u,i
real :: tnul, tscr

! --- benchmark NUL
call system_clock(tic,count_rate=rate)

open(newunit=u,file=nulunix,status='replace',iostat=ios, action='write')
if (ios /= 0) open(newunit=u,file=nulwin,status='replace',iostat=ios, action='write')
if (ios /= 0) error stop 'could not open a NULL file handle'

do i = 1,N
    write(u,*) 'blah blah blah'
    flush(u)
enddo
close(u)
call system_clock(toc)

tnul = (toc-tic)/real(rate)

print *,tnul,' seconds to write to NUL'

!---- benchmark scratch
call system_clock(tic)
open(newunit=u, status='scratch')
do i = 1,N
    write(u,*) 'blah blah blah'
    flush(u)
enddo
close(u)
call system_clock(toc)

tscr = (toc-tic)/real(rate)

print *,tscr,' seconds to write to scratch file'

print '(A,F7.3,A)','NUL is ',tscr/tnul,' times faster than scratch.'

end program


text:
The Internal Affairs Commission (IAC) is a subordinate body of the ASUCD ASUCD Senate Senate. IAC reviews legislation regarding matters of association policy and procedure before they go on to the Senate of the undergraduate student government.  As such, many of ASUCDs policy wonks get their start in IAC. Meetings are largely spent reviewing exciting stuff such as bylaws, LongRange Plans, and other internal matters.  The commission is made of nine members and a few alternate members.

IAC has open meetings on Mondays at 5:10 PM in the MU Mee Room.  The meetings are typically longer than that of any other ASUCD Commission, as they have so many bills to review.  IAC is also mostly responsible for writing the LongRange Plans for all ASUCD units.

The commission has been described as Infernal Affairs Commission Infernal Affairs, as they get all of the hard legislation and none of the glory that comes from running programs. However, IAC does hold programs like ASUCD legislationwriting clinics. The most recent one was successful but poorly attended, as only two of the six newly elected senators attended.

 Members

Commission Chair: Users/SpencerMcManus Spencer McManus

Commissioners

 History

 Past Chairs
  Users/SpencerMcManus Spencer McManus (Fall 2013Spring 2014)
  Amy Martin (Winter 2011)
  Users/SergioCano Sergio Cano (Winter 2010Fall 2010) and (Winter 2011Spring 2013)
  Users/MatthewShannon Matt Shannon (Spring 2009)
  Amy Hartstein (Fall 2008Winter 2009) and (Fall 2009Winter 2010)
  Users/JillWeinstein Jill Weinstein (Spring 2007)
  Mike Rivera (Fall 2006Winter 2007) and (Fall 2007Spring 2008)
  Paul Harms (Spring 2006)
  Kai SavareeRuess (Spring 2005Winter 2006)
  Marvin Zamora (Spring 2005)
  Kahliah Laney (Spring 2004Winter 2005)
  Users/BobGill Bob Gill (Fall 2003Spring 2004)
  Jakob Suski (Winter 2003Spring 2003)
  Pelin Erdal (Fall 2002)
  Erik Canales (Fall 2001Spring 2002)
  Allen Shiu (Fall 2000Spring 2001)
  Anthony Volkar (Fall 1999Spring 2000)
  Melanie Sengupta (Winter 1998Spring 1999)
  Scott Reed (Winter 1997Winter 1998)
  Matt Malone (Fall 1996Fall 1997)
  Maya Kwiat (Winter 1995Fall 1995)
  Peter Nguyen (Fall 1993Winter 1995)
  Larissa Adams (Fall 1992Spring 1993)
  Donny Anderson (Fall 1991Spring 1992)

 IAC: Birthplace of Power

There is a high incidence of powerful people coming from the ASUCD Internal Affairs Commission.

A few examples:

President
  Peter Nguyen, former ASUCD President 199596 and coauthor of the Democracy at Davis Constitutional Amendment
  Scott Reed, former ASUCD President 199899
  Sara Henry, former ASUCD President 20032004
  Darnell Holloway, former ASUCD President 20062007

 Vice President
  Aaron Bloom, former ASUCD Vice President 199596 and coauthor of the Democracy at Davis Constitutional Amendment
  Theron Spurgeon, former ASUCD Vice President 19992000
  Darnell Holloway, former ASUCD Vice President 20052006
  Molly Fluet, former ASUCD Vice President 20082009
  Previn Witana, former ASUCD Vice President 20102011
  Maxwell Kappes, ASUCD Vice President 20142015

 Senator
  Brett Leggett, former ASUCD Senator 199596
  Maya Kwiat, former ASUCD Senator 199596
  Erica Rios, former ASUCD Senator 199596
  Amos Tubb, former ASUCD Senator 199596
  Jason Capitan, former ASUCD Senator 199596
  Users/BobGill Robert Gill, former ASUCD Senator
  Janine Fiel, former ASUCD Senator 20042005
  Darnell Holloway, former ASUCD Senator 20042005
  Marvin Zamora, former ASUCD Senator 20052006
  Michael Lay, former ASUCD Senator 20062007
  Molly Fluet, former ASUCD Senator 20062007
  Justin Patrizio, former ASUCD Senator 20082009
  Erin Lebe, former ASUCD Senator 20082009
  Kevin Massoudi, former ASUCD Senator 20082009
  Previn Witana, former ASUCD Senator 20082009
  Don Ho, former ASUCD Senator 20092010
  Amy Martin, former ASUCD Senator 20102011
  Kabir Kapur, former ASUCD Senator 20112012
  Maxwell Kappes, former ASUCD Senator 20122013
  Gareth Smythe, ASUCD Senator 20132014

Controller
  Kai SavareeRuess, former ASUCD Controller 20062007
  Users/PaulHarms Paul Harms, former ASUCD Controller 20072008 and 20082009
  Joey Chen, former ASUCD Controller 20102011
  Don Ho, former ASUCD Controller 20112012
  Melanie Maemura, former ASUCD Controller 20122013

Unit Directors
  Brett Leggett, Director of the ASUCD Academic Affairs Office 199697
  Dorothy Yee, former Director of the ASUCD City/County Affairs Office 199596 and former Assistant Director of the ASUCD Academic Affairs Office 199697
  Scott Reed, former Director of the Book Exchange 19971998
  Derick Lennox, former Director of ASUCD Lobby Corps 20062008
  Matt Shannon, former Director of ASUCD University Affairs 20072009
  Kabir Kapur, Director of Office of Advocacy and Student Representation 20142015

Elections Committee Chairs
  Sergio Cano Sergio Cano, former ASUCD Elections Committee Chairperson Winter 2011
  Aaron Hsu, former ASUCD Elections Committee Chairperson 20122013
  Eric Renslo, ASUCD Elections Committee Chairperson 20132014

Chief of Staff
  Users/KevinPowers Kevin Powers, ASUCD Chief of Staff 20072008
  Amy Hartstein, ASUCD Chief of Staff 20082009
  Francisco Lara, ASUCD Chief of Staff 20142015

Other
  Tom Aguilar, GSA Vice Chair 20072008, GSA Chair 20062007, GSA External Chair 20052006 (was actually the Internal Affairs Representative to the Legislative Assembly, back in the long, long ago)
  ASUCD Business and Finance Commission, Greg Webb Gregory D. Webb
  Users/MissyWhitney Missy Lyla Whitney, ASUCD Court Chief Justice 20082010

IAC was very influential under the leadership of one of its longestserving chairs, Peter Nguyen.  Serving first as member and Vice Chair of the IAC from 199293, Nguyen went on to chair the commission from Fall 1993 through Winter of 1995 when he was elected ASUCD President.  During this era, a large number of ASUCD leaders emerged from IAC (see above),  and IAC was the primary progressive resistance to the Davis College Republicansdominated Executive Council.  In fact, so many Council bills died in IAC that, for purely political reasons, the Executive Council called a closed session to discuss the job performance of Nguyen as chair in Fall 1994.  Nguyen opened up the closed session to the public, exposed the political gambit to the press, forced the Council majority to back down from removing him as Chair, and caused a backlash against the Davis College Republicans which led to their fall from power in the Winter 1995 election.

 Dramatic incidents

  Rather dramatically, Users/RevChad was forced to resign from IAC when he joined the 2005/2006 Elections Committee as you can only hold one Senateconfirmed position at any given time.
  Pledge of Allegiance Controversies (20062007): Steve Ostrowski has been a frequent cause of drama, as has the commissioners overreactions to Ostrowski.  Steve introduced a bill to add the U.S. pledge of allegiance to the Senate agenda three times, it failed three times.  Due to the surprise of Users/PaulHarms various Users/AndrewPeake ASUCD officials at being videotaped, the meeting was shortly adjourned. The tape was turned into a http://www.youtube.com/watch?vEicfEjI4YA seven minute Youtube video with music, stock images, and written commentary in the form of title cards.  The cameraman was funded by the conservative http://www.leadershipinstitute.org/ Leadership Institute.  There was also a letter written by former DPD Officer and former Controller  Cameron Menezes that Infernal Affairs Commission attacked the commissioners position on the Pledge.

20070207 11:53:39 nbsp Where can I find the bills that have been proposed before ASUCD regarding the Pledge of Allegiance? Users/CarlosOverstreet
 http://asucd.ucdavis.edu/government/govdocuments/senatebills/bills.html Bill #12 is one of the previous bills to have the pledge mandatory in senate meetings.  ~Matt
 http://asucd.ucdavis.edu/government/govdocuments/senatebills/20052006 The 0506 year, bill #84.
   I see how its only a minor change in the agenda of the meeting, however, I would like to know why it was oppossed with such vehemence, besides the fact Steve Ostrowski was the author?  Users/CarlosOverstreet
    There were a variety of reasons, including that it went against the uc davis principles of community. Im not in a good state to answer your question.
    I read through the principles of community and I dont see how the pledge of allegiance (without a real loose interpretation of the principles) contradicts this.  Users/CarlosOverstreet
     Its a violation of the Principles of Community because the Pledge of Allegiance is, in part, a form of blasphemy in my religion.  As a Discordian, Im sure that Under God really means Under Goddess, so Im cool with that.  But it says To the republic, for which it stands.  This is clearly false, as wiki:wikipedia:Emperor_Norton Emperor Norton, our saint, has clearly declared the United States an empire.  To say anything different would sully the name of Norton I, Emperor of the United States and Protector of Mexico, and people should be ashamed to recite this desecrating nonsense.   Users/BrentLaabs
     Do you refrain from saying the pledge of allegiance at all?  Users/CarlosOverstreet
      In the last several Senate meetings, Andrew Peake Peake or Michael A. Molnar Molnar have been leading the Pledge of Allegiance, during public discussion, which is when I normally read off the names. There have been varying degrees of participation, and varying degrees of Under God, but in public discussion it is completely voluntary, and people are allowed to say and think what they want. Users/MaxMikalonis
     Emperor Norton did tolerate the continued presence of Congress, after they did not follow his orders and break up. That can be viewed as him tolerating the duel existance of the Republic. Users/JosephBleckman
    This bill parallels the placement of the flags at Senate. It has the possibility of alienating those that 1) do not identify themselves as American or 2) do not want to or feel it is necessary to pledge to any flag or country. The pledge can also be a very intimidating event for anybody that is attending Senate who is not comfortable with reciting or even hearing the pledge. The flags are a passive display, whereas reciting the pledge is an invasive event to some. Granted, recital of the pledge is voluntary, but listening to the pledge is not. Taken in account with bylaws that require certain people to be present at Senate, the passing of the bill can create a very unfair environment to some. So in the end, this bill can serve to alienate those that are mandated to attend Senate meetings and may even discourage student participation in ASUCD.  Users/EricWu
     The pledge of allegiance is intimidating? Oh, please, Eric. Pathetic, and indicative of how very little you understand or appreciate the United States of America. Users/CameronMenezes

20070207 12:21:27 nbsp The IAC shouldve known better, Mason Harrison anyone?  I dont think its a big deal but then again I wasnt there. Users/GregWebb

As such I want to prevent any speculation or claims regarding this incident, I dislike drama and would like to not see it on this wiki. ~Users/DavidPoole

Linking the video is rather pointless due to how the video itself is rather misleading. The only substantial part of the entire clip is when someone who isnt even on the commission tells the camera man to turn it off. Most people who watch it have absolutely no idea what is really going on, and Id rather not have misleading clip posted here to confuse people. Aside from that, if people feel that the event is worth a paragraph on this page so be it. ~Matt Shannon
  Im inclined to second Commissioner Shannons motion. The clip implies that both formerChairman Harms and Senator Peake are IAC Commissioners, which is an unfortunate implication to draw in this case. Unless this footage actually does make some semimajor media/news providers, the misrepresentation of the Commission is rather shady to have on this page. Now, if the footage does make Fox News, or something like that, then itd have a place here, due to the noteworthiness of it. Users/JosephBleckman Commissioner Bleckman
   Is there a law in ASCUD prohibiting the video recording of ASUCD Commission meetings? Why was steve recording in the first place?  Users/CarlosOverstreet
   I agree with all of you, I linked it to mitigate an edit war, I guess I rather not have pointless reverts etc. I think the issue has passed though.

20070228 16:53:32 nbsp Wow. That video is ridiculously stupid. Users/TusharRawat

20120514 00:06:35 nbsp Its Joshua Herskovitz not Herskoivtz Users/JoshuaHerskovitz

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine chdxtz(delxi,dlxmx0,dzvc0,iodb,iopt,jordlm,kdim,
     $ kmax,km1,kord,kxt,nodbmx,noptmx,nordmx,nordzs,noutpt,
     $ nrd1mx,nsscmx,scalim,scfczs,sscrew,qmin,smp100,uzvec1,
     $ zklogu,zvec0,zvclg0)
c
c     This subroutine chooses a step size and order according to the
c     Taylor's series accuracy criterion, examining the z vector and
c     its associated smoothed or average derivatives. Subroutine
c     chdxtz.f performs the same function for the r vector and its
c     associated average derivatives.
c
c     This subroutine is called by:
c
c       EQ6/path.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c
c     Principal output:
c
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer kmax,nodbmx,noptmx,nordmx,nrd1mx,nsscmx
c
      integer noutpt
c
      integer iodb(nodbmx),iopt(noptmx)
c
      integer jordlm,kdim,km1,kord,kxt,nordzs
c
      logical qmin
c
      character(len=48) uzvec1(kmax)
c
      real(8) dzvc0(nrd1mx,kmax),sscrew(nsscmx),zvclg0(kmax),zvec0(kmax)
c
      real(8) delxi,dlxmx0,scalim,scfczs,smp100,zklogu
c
c-----------------------------------------------------------------------
c
c     Local variable declarations with global dimensions.
c
      integer isv_nordmx
c
      SAVE isv_nordmx
c
      integer, dimension(:), allocatable :: kkk
c
      SAVE kkk
c
      character(len=24), dimension(:), allocatable :: uscal
c
      SAVE uscal
c
      real(8), dimension(:), allocatable :: scmax,scraw
c
      SAVE scmax,scraw
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,ip1,j2,j3,k,kcol,kcoli
c
      integer ilnobl
c
      character(len=24) ubas,uzt,uzn
c
      real(8) adx,delchg,dlxii,fctr,fctrmn,lxx0,mxx0,relchg,xip1
c
      real(8) fctrl
c
c-----------------------------------------------------------------------
c
      data ubas /'Basis variable          '/
c
c-----------------------------------------------------------------------
c
c     Allocate or reallocate local work arrays as needed.
c
      if (.not.ALLOCATED(kkk)) then
c
c       Local work arrays are not allocated. Zero the saved
c       array size variables. Note that only one array is tested
c       to see if it is allocated. It is assumed that all local
c       work arrays are either allocated or not.
c
        isv_nordmx = 0
      else
c
c       Local work arrays are allocated. Check to see if any of the
c       array size variables have changed. If so, deallocate
c       the corresponding local work arrays and zero the corresponding
c       saved size variables.
c
        if (nordmx .ne. isv_nordmx) then
          DEALLOCATE(kkk,uscal,scmax,scraw)
          isv_nordmx = 0
        endif
      endif
c
c     At this point, the saved array size values are zero if the
c     corresponding arrays need to be allocated.
c
      if (isv_nordmx .eq. 0) then
        ALLOCATE(kkk(0:nordmx))
        ALLOCATE(uscal(0:nordmx))
        ALLOCATE(scmax(0:nordmx),scraw(0:nordmx))
        isv_nordmx = nordmx
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Zero the contents of the local work arrays.
c
      do k = 0,nordmx
        kkk(k) = 0
      enddo
c
      do k = 0,nordmx
        uscal(k) = ' '
      enddo
c
      do k = 0,nordmx
        scmax(k) = 0.
        scraw(k) = 0.
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (kord .le. 0) then
        nordzs = 0
        scfczs = dlxmx0/delxi
        go to 999
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Set up the regular choice of orders. Here i is the order, which
c     technically runs from 0 to jordlm. Here jordlm + 1 is the hidden
c     extra order, which is used to constrain the error in using order
c     jordlm.
c
      do i = 0,jordlm
        scmax(i) = 0.
        scraw(i) = 0.
      enddo
c
c     For i = 0, take delxi = dlxmx0. This keeps delxi greater
c     than or equal to dlxmx0, except within searches, as for phase
c     boundaries.
c
      scmax(0) = dlxmx0/delxi
      scraw(0) = scmax(0)
      kkk(0) = 0
      uscal(0) = 'None'
c
c     Process the higher orders.
c
      do i = 1,jordlm
        ip1 = i + 1
        dlxii = (delxi**(ip1))/fctrl(ip1)
        fctrmn = 1.e+38
        kkk(i) = 0
        uscal(i) = 'None'
c
c       Constrain the step size to satisfy the tolerance on the
c       fractional error of matrix variables.
c
        do kcol = 1,kdim
          qmin = kcol.ge.km1 .and. kcol.le.kxt
          lxx0 = zvclg0(kcol)
c
c         Skip the constraint for a mineral with a very small mass.
c
          if (.not.qmin .or. lxx0.gt.zklogu) then
c
            adx = abs(dzvc0(ip1,kcol))
            delchg = adx*dlxii
            if (delchg .ne. 0.) then
c
c             Taylor's series in a linear quantity (e.g., mass).
c
              mxx0 = zvec0(kcol)
              if (mxx0 .ge. smp100) then
                relchg = delchg/mxx0
                fctr = sscrew(1)/relchg
              else
                relchg = 1.e+38
                fctr = 1.e+38
              endif
c
              if (fctr .lt. fctrmn) then
                fctrmn = fctr
                uscal(i)(1:24) = ubas(1:24)
                kkk(i) = kcol
              endif
            endif
          endif
        enddo
c
        xip1 = ip1
        scmax(i) = fctrmn**(1./xip1)
        scraw(i) = scmax(i)
        if (scmax(i) .gt. scalim) scmax(i) = scalim
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (iodb(5) .ge. 2) then
        write (noutpt,1000)
 1000   format(/' --- Order and Scale Factor Data (chdxtz) ---',/)
        do i = 0,jordlm
          uzt = uscal(i)
          kcoli = kkk(i)
          if (uzt(1:24) .eq. ubas(1:24)) then
            uzn = uzvec1(kcoli)(1:24)
          else
            uzn = ' '
          endif
c
          write (noutpt,1010) i,scmax(i),scraw(i)
 1010     format('   Order= ',i2,', scale= ',f10.5,', raw scale= ',
     $    f10.5)
          j2 = ilnobl(uzt)
          j3 = ilnobl(uzn)
          if (j3 .gt. 0) then
            write (noutpt,1020) uzt(1:j2),uzn(1:j3)
 1020       format(5x,'Constrained by ',a,1x,a)
          else
            write (noutpt,1030) uzt(1:j2)
 1030       format(5x,'Constrained by ',a)
          endif
        enddo
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Choose the order that gives the highest scale factor.
c
      scfczs = scmax(0)
      nordzs = 0
      do i = 1,jordlm
        if (scmax(i) .gt. scfczs) then
          scfczs = scmax(i)
          nordzs = i
        endif
      enddo
c
  999 continue
      end


text:
Describe Users/SASSAFRASBEAUTY here.

20120827 15:38:27 nbsp Welcome to the Wiki/Business Owner Welcome to the Wiki, Sassafrasbeauty! I reworked your Sassafras Beauty page to clean it up a bit. I labelled the pics in a cleaner manner as well. Hope you like it Users/PeteB