text:
C    @(#)kmpymx.f	20.3 2/13/96
        function  kmpymx (m,n)

        include 'ipfinc/parametr.inc'

        include 'ipfinc/red7.inc'

        kmpymx = ymtrx(1,m) - ymtrx(1,n)
        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:
 Adam Dead Baby Cow Veal
Poet, Artist, Chef, Lover

Adam used to be the morning cook at Delta of Venus. He made the best food. He is the creator of the Adams Special. I dont think he eats veal, because that would be wrong.

Adam is a cool guy, one time in line at Davis Food CoOp the Coop the credit/debit machines went down and I was graced with his presence. It was nice to have a goodlooking guy ask me what my sign was, instead of some older grump complaining about the machines. He also sometimes likes pie.  Users/MichelleAccurso

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c
c -----------------------------------------------------
c
      subroutine valout (mgrid, lst, lend, time, nvar)
c
      implicit double precision (a-h,o-z)
      logical            graf
      parameter  (maxgr = 192, maxlv=12)
      common /nodal/ rnode(12,maxgr),node(17,maxgr),lstart(maxlv),
     *  newstl(maxlv),
     *  listsp(maxlv),intrat(maxlv),tol,bzonex,bzoney,mstart,ndfree,
     *  lfine,iorder,mxnest,kcheck,lwidth,
     *  maxcl, graf, lhead
      include "calloc.i"
      common /userdt/ cfl,gamma,gamma1,xprob,yprob,alpha,Re,iprob,
     .                ismp,gradThreshold
      logical         flag
c
c valout = graphics output of soln values for contour or surface plots.
c    if mgrid <> 0,  output only that grid, (should have
c                          lst = lend here)
c
      call basic (time, lst, lend, 3 )
c
      write(3,100) mgrid,  nvar
100   format(10h*VALS     ,2i10)
c
      flag = .false.
      if (lst .eq. lend .and. mgrid .ne. 0) flag = .true.
c
      level = lst
10    if (level .gt. lend) go to 99
          mptr = lstart(level)
          if (flag) mptr = mgrid
20        if (mptr .eq. 0) go to 50
              maxi   = node(5,mptr)
              maxj   = node(6,mptr)
	      xlow = rnode(1,mptr)
	      ylow = rnode(2,mptr)
              maxip1 = maxi + 1
              maxjp1 = maxj + 1
              loc = node(7,mptr)
	      locirr = node(14,mptr)
	      mitot = maxi-1+2*lwidth
	      mjtot = maxj-1+2*lwidth
	      locsm = igetsp(mitot*mjtot*nvar)
              call outvar(alloc(loc),maxip1,maxjp1,nvar,mptr,
     1                    alloc(locirr),mitot,mjtot,lwidth,
     2                    node(17,mptr),rnode(9,mptr),
     3                    rnode(10,mptr),xlow,ylow,alloc(locsm))
	      call reclam(locsm,mitot*mjtot*nvar)
              mptr = node(10,mptr)
              if (flag) go to 50
          go to 20
50        continue
          level = level + 1
      go to 10
c
 99   return
      end


!-----------------------------------------------------------------------
      program test
      implicit none
      include 'ceedf.h'

      integer ceed,err
      integer x,i,n
      real*8 a(10)
      real*8 norm,diff
      integer*8 aoffset
      character arg*32

      call getarg(1,arg)

      call ceedinit(trim(arg)//char(0),ceed,err)

      n=10

      call ceedvectorcreate(ceed,n,x,err)

      do i=1,10
        if (mod(i,2) == 0) then
          a(i)=i-1
        else
          a(i)=-(i-1)
        endif
      enddo
      aoffset=0
      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,a,aoffset,err)

      call ceedvectornorm(x,ceed_norm_1,norm,err)
      diff = norm - 45.
      if (abs(diff)>1.0D-14) then
! LCOV_EXCL_START
        write(*,*) 'Error L1 norm ',norm,' != 45.'
! LCOV_EXCL_STOP
      endif

      call ceedvectornorm(x,ceed_norm_2,norm,err)
      diff = norm - sqrt(285.)
      if (abs(diff)>1.0D-6) then
! LCOV_EXCL_START
        write(*,*) 'Error L2 norm ',norm,' != sqrt(285.)'
! LCOV_EXCL_STOP
      endif

      call ceedvectornorm(x,ceed_norm_max,norm,err)
      diff = norm - 9.
      if (abs(diff)>1.0D-14) then
! LCOV_EXCL_START
        write(*,*) 'Error Max norm ',norm,' != 9.'
! LCOV_EXCL_STOP
      endif

      call ceedvectordestroy(x,err)
      call ceeddestroy(ceed,err)

      end
!-----------------------------------------------------------------------


text:
subroutine readpartele_solid(partele)
  use solid_variables, only: ne_solid, ne_solid_1, n_solid
  implicit none

  integer :: partele(ne_solid)
  integer :: file,i
  integer j

  file=21
  open(file, FILE="partele_solid.in", STATUS="old")
  do i=1,ne_solid_1
     read(file,*) partele(i)
  enddo

  do j=2,n_solid
  	do i=1,ne_solid_1
	partele((j-1)*ne_solid_1 + i) = partele(i)
	end do
  end do

  close(file)

  return
end subroutine readpartele_solid

subroutine readpartnode_solid(partnode)
  use solid_variables, only: nn_solid, nn_solid_1, n_solid
  implicit none

  integer :: partnode(nn_solid)
  integer :: file,i
  integer j
  file=21
  open(file, FILE="partnode_solid.in", STATUS="old")
  do i=1,nn_solid_1
     read(file,*) partnode(i)
  enddo

  do j=2,n_solid
  	do i=1,nn_solid_1
	partnode((j-1)*nn_solid_1 + i) = partnode(i)
	end do
  end do

  close(file)

  return
end subroutine readpartnode_solid

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


      program demo_show
      use M_kracken, only : kracken, show
      implicit none

      call kracken('demo', ' default keyword -i 10 -j 20.20 -value my default string')
      call show('demo',.false.,0)

      end program demo_show


text:
C     LAST UPDATE 28/03/96
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE SETUP
      IMPLICIT NONE
C
C Purpose: Puts angles in radians and works out detector orientation
C
C Calls   2: ARSET , AISET
C Called by: FTOREC
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C
C Common arrays:
C
      REAL RDAT(2500*2500),RECSPA(1250*1250),STDEV(1250*1250),PHI(3,3)
      INTEGER NPOINT(1250*1250)
C
C Common scalars:
C
      REAL DIST,XC,YC,ROTX,ROTY,ROTZ,TILT,WAVE,XYRAT
      REAL DMIN,DMAX,SIGMIN,SIGMAX,RMIN,RMAX,ZMIN,ZMAX,RBAK,ZBAK
      REAL DELR,DELZ,DELSIG,DELD
      REAL ABSBAS,ABSEML,ABSPAP,C0,PAPN,FILN
      INTEGER NFBEFR,NPBEFR
      INTEGER NXMIN,NXMAX,NYMIN,NYMAX,NRD,NZSIG
C
C Local variables:
C
      REAL CRX,SRX,CRY,SRY,CRZ,SRZ
C
C Common blocks:
C
      COMMON /ARRAYS/ RDAT,RECSPA,STDEV,NPOINT
      COMMON /ORIENT/ PHI
      COMMON /TRANSF/ DIST,XC,YC,ROTX,ROTY,ROTZ,TILT,WAVE,XYRAT
      COMMON /RLIMIT/ DMIN,DMAX,SIGMIN,SIGMAX,RMIN,RMAX,ZMIN,ZMAX,
     &                RBAK,ZBAK
      COMMON /DELTAS/ DELR,DELZ,DELSIG,DELD
      COMMON /FILMAB/ ABSBAS,ABSEML,ABSPAP,NFBEFR,NPBEFR,C0,PAPN,FILN
      COMMON /ILIMIT/ NXMIN,NXMAX,NYMIN,NYMAX,NRD,NZSIG
C
C-----------------------------------------------------------------------
C
C========Make tilt have same sign as for GUCKMAL.
C
      TILT = -TILT/57.296
C
      ROTX = ROTX/57.296
      ROTY = ROTY/57.296
      ROTZ = ROTZ/57.296
      SIGMIN = SIGMIN/57.296
      SIGMAX = SIGMAX/57.296
      DELSIG = DELSIG/57.296
      RMAX = FLOAT(NRD-1)*DELR + RMIN
      ZMAX = FLOAT(NZSIG-1)*DELZ + ZMIN
C
      NYMIN = NYMIN + 1
      NXMIN = NXMIN + 1
      NYMAX = NYMAX + 1
      NXMAX = NXMAX + 1
      CALL ARSET(RECSPA,1250*1250,0.0)
      CALL ARSET(STDEV,1250*1250,0.0)
      CALL AISET(NPOINT,1250*1250,0)
C
C========Set up film absorption coefficients
C
      C0 = (1.0-EXP(-ABSEML))*(1.0+EXP(-ABSEML-ABSBAS))
      PAPN = FLOAT(NPBEFR)*ABSPAP
      FILN = FLOAT(NFBEFR)*(2.0*ABSEML+ABSBAS)
C
C========Form sines and cosines of detector rotation angles
C
      CRX = COS(ROTX)
      SRX = SIN(ROTX)
      CRY = COS(ROTY)
      SRY = SIN(ROTY)
      CRZ = COS(ROTZ)
      SRZ = SIN(ROTZ)
C
C========Form detector rotation matrix RzRyRx
C
      PHI(1,1) = CRZ*CRY
      PHI(1,2) = CRZ*SRY*SRX + SRZ*CRX
      PHI(1,3) = -CRZ*SRY*CRX + SRZ*SRX
      PHI(2,1) = -SRZ*CRY
      PHI(2,2) = -SRZ*SRY*SRX + CRZ*CRX
      PHI(2,3) = SRZ*SRY*CRX + CRZ*SRX
      PHI(3,1) = SRY
      PHI(3,2) = -CRY*SRX
      PHI(3,3) = CRY*CRX
C
      RETURN
      END


text: #redirect Traffic Circle


text:
A number of products have their origins in Davis, including:

  Afghan Pepper Company Afghan Pepper Company hot sauce  not in production yet
  BRÜBAR Energy Bars
  http://diamondtrustgame.com/buy.php Diamond Trust of London  Nintendo DS game (http://www.kickstarter.com/projects/1443658586/diamondtrustoflondon More Info)
  Moller International Flying cars (https://www.youtube.com/watch?vJ7dJvFoTx3Q sort of...)
  Schilling Robotics Robotic arms and Schilling Robotics Remotely operated submersible vehicles
  Serenade agricultural fungicide, created by AgraQuest
  Embedded JK microsystems Single Board Computers
  Mori Seiki http://www.dmgmoriseikiusa.com/nhxseries/nhx4000 NHX4000 horizontal machining centers
  The Philosophers Stoneground  sprouted almond butter
  Sudwerk Sudwerk beer, including Aggie Lager
  UC Davis Olive Oil
  Users/BillKendrick Tux Paint

Former Davis products:

  FlavrSavr tomato, created by CalGene LLC
  Square Tomato
  YummyDummy Chocolate Company YummyDummy Chocolate
  Canned Fruit Cocktail was invented by Dr. Cruess, for whom Cruess Hall is named


subroutine refit_1fc(spec, freqs, fc, omega, n, oo, ff, nf, res)
    implicit none
    ! Argument declarations
    integer :: ff, oo, nf

    double precision, intent(in) :: omega(oo), fc(ff), n(ff), spec(nf), freqs(nf)
    double precision, intent(out) :: res(3)

    ! Variable declarations
    double precision :: synth(nf), best_ssq, ssq, weights(nf)
    integer :: i, k
    best_ssq = HUGE(ssq)
    weights = 1!/freqs
    do i = 1, ff
        do k = 1, oo
            synth = omega(k) / (1 + (freqs/fc(i))**n(i))
            ssq = sum(((synth-spec)*weights)**2)
            if (ssq < best_ssq) then
                best_ssq = ssq
                res(1) = fc(i)
                res(2) = omega(k)
                res(3) = n(i)
            end if
        end do
    end do
end subroutine refit_1fc

subroutine refit_2fc(spec, freqs, fc1, fc2, omega, n, oo, ff, gg, nf, res)
    implicit none
    ! Argument declarations
    integer :: ff, oo, nf, gg
    double precision, intent(in) :: omega(oo), n(ff), spec(nf), freqs(nf)
    double precision, intent(in) :: fc1(ff), fc2(gg)
    double precision, intent(out) :: res(4)

    ! Variable declarations
    double precision :: synth(nf), best_ssq, ssq, weights(nf)
    integer :: i, k, l
    best_ssq = HUGE(ssq)
    do i = 1, ff
        do k = 1, oo
            do l = 1, gg
                !synth = log10(omega(k)) + log10(1 + (freqs/fc2(l))**n(i))
                !synth = synth - log10(1 + (freqs/fc1(i))**n(i))
                synth = (1 + (freqs/fc2(l))**n(i))
                synth = synth / (1 + (freqs/fc1(i))**n(i))
                synth = synth * omega(k)
                ssq = sum((synth-spec)**2)
                if (ssq < best_ssq) then
                    best_ssq = ssq
                    res(1) = fc1(i)
                    res(2) = omega(k)
                    res(3) = n(i)
                    res(4) = fc2(l)
                end if
            end do
        end do
    end do

end subroutine refit_2fc

subroutine fit_1fc(spec, freqs, fc, omega, n, oo, nn, ff, nf, res)
    implicit none
    ! Argument declarations
    integer :: ff, nn, oo, nf
    double precision, intent(in) :: omega(oo), fc(ff), n(nn), spec(nf), freqs(nf)
    double precision, intent(out) :: res(3)

    ! Variable declarations
    double precision :: synth(nf), best_ssq, ssq, weights(nf)
    integer :: i, j, k
    best_ssq = HUGE(ssq)
    weights = 1/freqs
    do i = 1, ff
        do j = 1, nn
            do k = 1, oo
                synth = omega(k) / (1 + (freqs/fc(i))**n(j))
                ssq = sum(((synth-spec))**2)
                if (ssq < best_ssq) then
                    best_ssq = ssq
                    res(1) = fc(i)
                    res(2) = omega(k)
                    res(3) = n(j)
                end if
            end do
        end do
    end do

end subroutine fit_1fc

subroutine fit_2fc(spec, freqs, fc1, fc2, omega, n, oo, nn, ff, nf, res, gg, Y)
    implicit none
    ! Argument declarations
    integer :: ff, nn, oo, nf, gg
    double precision, intent(in) :: omega(oo), n(nn), spec(nf), freqs(nf)
    double precision, intent(in) :: fc1(ff), fc2(gg), Y
    double precision, intent(out) :: res(4)

    ! Variable declarations
    double precision :: synth(nf), best_ssq, ssq, weights(nf)
    integer :: i, j, k, l
    best_ssq = HUGE(ssq)
    do i = 1, ff
        do j = 1, nn
            do k = 1, oo
                do l = 1, gg
                    synth = log10(omega(k)) + (1./Y)*log10(1 + &
                    (freqs/fc2(l))**(n(j)*Y))
                    synth = synth - (1./Y)*log10(1 + (freqs/fc1(i))**(Y*n(j)))
                   ! synth = (1 + (freqs/fc2(l))**n(j))
                   ! synth = synth / (1 + (freqs/fc1(i))**n(j))
                   ! synth = synth * omega(k)
                    ssq = sum((synth-spec)**2)
                    if (ssq < best_ssq) then
                        best_ssq = ssq
                        res(1) = fc1(i)
                        res(2) = omega(k)
                        res(3) = n(j)
                        res(4) = fc2(l)
                    end if
                end do
            end do
        end do
    end do

end subroutine fit_2fc


text:
! { dg-do compile { target { { i?86-*-* x86_64-*-* } && { ! { ia32 } } } } }

! { dg-options "-Ofast -mavx512f -ffixed-xmm1 -ffixed-xmm2 -ffixed-xmm3 -ffixed-xmm4 -ffixed-xmm5 -ffixed-xmm6 -ffixed-xmm7 -ffixed-xmm8 -ffixed-xmm9 -ffixed-xmm10 -ffixed-xmm11 -ffixed-xmm12 -ffixed-xmm13 -ffixed-xmm14 -ffixed-xmm15" }

      IMPLICIT REAL*8(A-H,O-Z)
      ALLOCATABLE DD1(:), DD2(:), WY(:,:)
      ALLOCATE( DD1(MAX), DD2(MAX), WY(MAX,MAX))
         DO J = J1,J2
            DO I = I1, I2
               DD1(I) = D1 * (WY(I-2,J) - WY(I+2,J) +
     >              (WY(I+1,J) - WY(I-1,J)))
            END DO
            DO I = I1, INT(D2 * D3(I))
            END DO
         END DO
      END

! { dg-final { scan-assembler-not "vbroadcastsd\[ \\t\]+%xmm\[0-9\]+, %ymm\[0-9\]+" } }


text:
C*************************************************************************
C  DRY MATTER PARTITIONING
       SUBROUTINE SC_PARTIT (WaterBal, Part, Climate, CaneCrop,
     -                    Growth, Out, Control, ISWITCH,
     -                    STGDOY, DSUCyest, CELLSE_DM, LEAFTT)
c     Use the modules:
      USE CNG_MODULEDEFS
      USE MODULEDEFS

      IMPLICIT NONE
      SAVE

c     Declare modules in use:
c     :::::::::::::::::::::::
      TYPE(WATERTYPE)    WaterBal
      TYPE(PARTTYPE)     Part
      TYPE(CLIMATETYPE)  Climate
      TYPE(CANECROPTYPE) CaneCrop
      TYPE(OUTTYPE)      Out
      TYPE(GROTHTYPE)    Growth
      TYPE(SwitchType) ISWITCH ! MJ: no longer used

c     Dssat control variable:
      TYPE(ControlType) Control

c     Parameter variables:
c     ::::::::::::::::::::
c      INTEGER DOY
c      REAL    IRR
      INTEGER STGDOY(20)

c     Common 'OUTBLOCK' variables in use:
c     :::::::::::::::::::::::::::::::::::
c     Thermal time for leaf development --> start of stalk elongation
      ! REAL    CHUPI, LEAFTT
      REAL    LEAFTT
      REAL    ROOTDM
      REAL    TRASDM

c     Common PART variables in use:
c     :::::::::::::::::::::::::::::
      REAL    AERDWDT
      REAL    AERLDM
      REAL    APFMX
      REAL    CHUPIBASE
      REAL    DELTAOUT
      REAL    DELTTMAX
      REAL    FDWDT
      REAL    FTCON
      REAL    NSUCMAS
      REAL    PCB
      REAL    STKDM
      REAL    STKPFMAX
      REAL    SUCA
      REAL    SUCMAS
      REAL    SURCON
      REAL    SWDF2AMP
c      REAL    TBDELTT
      REAL    TBFT
      REAL    TDMAS
      REAL    TOPDM

c     FIBRE PARTITIONING
c     MJ, June 2015
c     Note: Additional state variables not in the PART block
c     Cellulosic DM (t/ha)
      REAL, INTENT(OUT) ::  CELLSE_DM
c     Genetic parameters related to fibre partitioning
c     Fraction of stalk fibre that is cellulose
      REAL FF_CELLSE

c     Common WATER variables in use:
c     ::::::::::::::::::::::::::::::
      REAL    SWDF2

c     Common CLIMT variables in use:
c     ::::::::::::::::::::::::::::::
      REAL    TEMPMN
      REAL    TEMPMX

c     Common CANCRP variables in use:
c     :::::::::::::::::::::::::::::::
      REAL    TMEANDEDLF
      REAL    TMEANLF

c     Local variables:
c     ::::::::::::::::
      REAL    DELTA
      REAL    DELTT
      REAL    DNSTKDT
      REAL    DNSUCDT
      REAL    DNSUCTEMP
      REAL    DSTKDT
      REAL    DSUCDT
      REAL    DSUCTEMP
      REAL    FTDELTA
      REAL    FWDELTA
      REAL    MEANTEMP
      REAL    NNSTK
      REAL    NSUCCAP
      REAL    ROOTDWDT
      REAL    ROOTPF
      REAL    STKPF
      REAL    SUCCAP
      REAL    SUCEQ
      REAL    SUCPF

c     Maximum root partitioning fraction
      REAL    MAX_ROOTPF

c     For reading coeffs:
      LOGICAL CF_ERR, SPC_ERR

c     Yesterday's value of delta sucrose mass:
c     t/ha/d
      REAL, INTENT(OUT) :: DSUCyest

c     *******************************************************************
c     Init:
      IF (CONTROL%DYNAMIC .EQ. RUNINIT) THEN
c       MJ, Mar 2010: zero all of the partitioning variables in the composite variable:
         Part%AERDWDT  = 0.0
         Part%AERLDM  = 0.0
         Part%AMPLIF  = 0.0
         Part%APFMX  = 0.0
         Part%BRXPDM  = 0.0
         Part%CHUPIBASE  = 0.0
         Part%DELTAOUT  = 0.0
         Part%DELTTMAX  = 0.0
         Part%DMMAS(150)  = 0.0
         Part%FDWDT  = 0.0
         Part%FLODGE  = 0.0
         Part%FTCON  = 0.0
         Part%NSUCMAS  = 0.0
         Part%PARCE  = 0.0
         Part%PCB  = 0.0
         Part%PURITY  = 0.0
         Part%RESPCF  = 0.0
         Part%RFA  = 0.0
         Part%RFB  = 0.0
         Part%RFC  = 0.0
         Part%ROOTF  = 0.0
         Part%STKDM  = 0.0
         Part%STKPFMAX  = 0.0
         Part%STKWM  = 0.0
         Part%SUCA  = 0.0
         Part%SUCB  = 0.0
         Part%SUCMAS  = 0.0
         Part%SUCNEWMAS  = 0.0
         Part%SURCON  = 0.0
         Part%SWDF2AMP  = 0.0
         Part%TBASEFNS  = 0.0
         Part%TBASETOT  = 0.0
         Part%TBDELTT  = 0.0
         Part%TBFT  = 0.0
         Part%TDMAS  = 0.0
         Part%TOPDM  = 0.0
         Part%TOPRATE  = 0.0
         Part%TOPTFNS  = 0.0
         Part%TOPTTOT  = 0.0
         Part%TOTBASE  = 0.0

c      Now zero local variables
         ! CHUPI = 0.
         ROOTDM = 0.
         TRASDM = 0.
         AERDWDT  = 0.0
         AERLDM  = 0.0
         APFMX  = 0.0
         CHUPIBASE  = 0.0
         DELTAOUT  = 0.0
         DELTTMAX  = 0.0
         FDWDT  = 0.0
         FTCON  = 0.0
         NSUCMAS  = 0.0
         PCB  = 0.0
         STKDM  = 0.0
         STKPFMAX  = 0.0
         SUCA  = 0.0
         SUCMAS  = 0.0
         SURCON  = 0.0
         SWDF2AMP  = 0.0
         TBFT  = 0.0
         TDMAS  = 0.0
         TOPDM  = 0.0
         SWDF2  = 0.0
         TEMPMN  = 0.0
         TEMPMX  = 0.0
         TMEANDEDLF  = 0.0
         TMEANLF  = 0.0
         DELTA  = 0.0
         DELTT  = 0.0
         DNSTKDT  = 0.0
         DNSUCDT  = 0.0
         DNSUCTEMP  = 0.0
         DSTKDT  = 0.0
         DSUCDT  = 0.0
         DSUCTEMP  = 0.0
         FTDELTA  = 0.0
         FWDELTA  = 0.0
         MEANTEMP  = 0.0
         NNSTK  = 0.0
         NSUCCAP  = 0.0
         ROOTDWDT  = 0.0
         ROOTPF  = 0.0
         STKPF  = 0.0
         SUCCAP  = 0.0
         SUCEQ  = 0.0
         SUCPF  = 0.0
         MAX_ROOTPF  = 0.0

      ENDIF

      IF (CONTROL%DYNAMIC .EQ. SEASINIT) THEN
c         These values should come from cultivar / species file

c        Species params:
c         Part%TBDELTT   = 28.
         Part%SWDF2AMP  = 0.5

c        Get these parameters from file:
c         CALL GET_SPECIES_COEFF(Part%TBDELTT,  'TBDELTT',
c     &                          Control, SPC_ERR)

         CALL GET_CULTIVAR_COEFF(Part%SWDF2AMP, 'SWDF2AMP',
     &                           Control, SPC_ERR)

c         Set default values:
          Part%APFMX     = 0.88
          Part%CHUPIBASE = 1050.
          Part%DELTTMAX  = 0.07
          Part%FTCON     = 0.32
          Part%PCB       = 0.6
          Part%STKPFMAX  = 0.65
          Part%SUCA      = 0.58
          Part%SURCON    = 0.99
          Part%TBFT      = 25.
          MAX_ROOTPF     = 0.95

c         Read from file:
c         :::::::::::::::
c             AFPMX:
              CALL GET_CULTIVAR_COEFF(Part%APFMX, 'APFMX',
     -                             CONTROL,  CF_ERR)
c             CHUPIBASE:
              CALL GET_CULTIVAR_COEFF(Part%CHUPIBASE, 'CHUPIBASE',
     -                             CONTROL,  CF_ERR)
c             DELTTMAX:
              CALL GET_CULTIVAR_COEFF(Part%DELTTMAX, 'DELTTMAX',
     -                             CONTROL,  CF_ERR)
c             FTCON:
c             :::::::::
c             MJ, Mar 2008: moved FTCON to species file:
c              CALL GET_CULTIVAR_COEFF(Part%FTCON, 'FTCON',
c     -                             CONTROL,  CF_ERR)
              CALL GET_SPECIES_COEFF(Part%FTCON, 'FTCON',
     -                             CONTROL,  SPC_ERR)
c             :::::::::
c             PCB:
c             :::::::::
c             MJ, Mar 2008: moved PCB to species file
c              CALL GET_CULTIVAR_COEFF(Part%PCB, 'PCB',
c     -                             CONTROL,  CF_ERR)

              CALL GET_SPECIES_COEFF(Part%PCB, 'PCB',
     -                             CONTROL,  SPC_ERR)
c             :::::::::

c             STKPFMAX:
              CALL GET_CULTIVAR_COEFF(Part%STKPFMAX, 'STKPFMAX',
     -                             CONTROL,  CF_ERR)
c             SUCA:
              CALL GET_CULTIVAR_COEFF(Part%SUCA, 'SUCA',
     -                             CONTROL,  CF_ERR)
c             SURCON:
c             :::::::::
c             MJ, Mar 2008: moved SURCON to species file
c              CALL GET_CULTIVAR_COEFF(Part%SURCON, 'SURCON',
c     -                             CONTROL,  CF_ERR)
              CALL GET_SPECIES_COEFF(Part%SURCON, 'SURCON',
     -                             CONTROL,  SPC_ERR)
c             :::::::::

c             TBFT:
              CALL GET_CULTIVAR_COEFF(Part%TBFT, 'TBFT',
     -                             CONTROL,  CF_ERR)

c             MAX_ROOTPF
              CALL GET_SPECIES_COEFF(MAX_ROOTPF, 'MAX_ROOTPF',
     &                               Control, SPC_ERR)

c      Fibre partitioning
c      MJ, June 2015
       FF_CELLSE = 0.50
!       CALL GET_CULTIVAR_COEFF(FF_CELLSE, 'FF_CELLSE', CONTROL, CF_ERR)
       CELLSE_DM = 0.0

c         :::::::::::::::::::::::::::::::::::::::::::::::::::::::

c         These are variables that need to be zeroed for each run/treatment:
          Part%TDMAS     = 0.
          Part%TOPDM     = 0.
          Part%SUCMAS    = 0.
          Part%STKDM     = 0.
          Part%NSUCMAS   = 0.
          Part%FDWDT     = 0.
          Part%DELTAOUT  = 0.
          Part%AERDWDT   = 0.
          Part%AERLDM    = 0.
          Out%TRASDM     = 0.
          Part%TOPDM     = 0.
          Out%ROOTDM     = 0.
          Part%ROOTF     = 0.5

c         No sucrose yesterday...
          DSUCyest = 0.0

      ELSEIF (CONTROL%DYNAMIC .EQ. RATE) THEN

c     MJ: copy values from common blocks to local variables:
c     ::::::::::::::::::::::::::::::::::::::::::::::::::::::
c     WATER:
c     ::::::
      SWDF2 = WaterBal%SWDF2

c     PART:
c     :::::
      AERDWDT   = Part%AERDWDT
      AERLDM    = Part%AERLDM
      APFMX     = Part%APFMX
      CHUPIBASE = Part%CHUPIBASE
      DELTAOUT  = Part%DELTAOUT
      DELTTMAX  = Part%DELTTMAX
      FDWDT     = Part%FDWDT
      FTCON     = Part%FTCON
      NSUCMAS   = Part%NSUCMAS
      PCB       = Part%PCB
      STKDM     = Part%STKDM
      STKPFMAX  = Part%STKPFMAX
      SUCA      = Part%SUCA
      SUCMAS    = Part%SUCMAS
      SURCON    = Part%SURCON
      SWDF2AMP  = Part%SWDF2AMP
c      TBDELTT   = Part%TBDELTT
      TBFT      = Part%TBFT
      TDMAS     = Part%TDMAS
      TOPDM     = Part%TOPDM

c     CLIMT:
c     ::::::
      TEMPMN    = Climate%TEMPMN
      TEMPMX    = Climate%TEMPMX

c     Canecrop:
c     :::::::::
      TMEANDEDLF = CaneCrop%TMEANDEDLF
      TMEANLF    = CaneCrop%TMEANLF

c     Outblock:
c     :::::::::
c     Note: rootdm and trasdm should be in different block...
      ! CHUPI  = Out%CHUPI + Out%HUPI
      ROOTDM = Out%ROOTDM
      TRASDM = Out%TRASDM

c     *******************************************************************

c         Local variables:
          DELTA     = 0.
          DELTT     = 0.
          DNSTKDT   = 0.
          DNSUCDT   = 0.
          DNSUCTEMP = 0.
          DSTKDT    = 0.
          DSUCDT    = 0.
          DSUCTEMP  = 0.
          FTDELTA   = 0.
          FWDELTA   = 0.
          MEANTEMP  = 0.
          NNSTK     = 0.
          NSUCCAP   = 0.
          ROOTDWDT  = 0.
          ROOTPF    = 0.
          STKPF     = 0.
          SUCCAP    = 0.
          SUCEQ     = 0.
          SUCPF     = 0.

 !CNB JUL 2000 {NEWPART} AERDWDT IS THE AERL DAILY ALLOCATION

c     Calculate partitioning of biomass increment to roots:
c     :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c         MJ: sample from 376 variety file:
c     ----------------------------------------------------------------
c     MAX FRAC DM TO AERLDM              APFMX       |0.880   |t/t   7
c     PARTITION COEFF.                   PCB         |0.600   |      8
c     ----------------------------------------------------------------
          ROOTDM=TDMAS-AERLDM

c         MJ, 2007-04-03: replaced 0.95 with 'MAX_ROOTPF' in species file
c	    ROOTPF=AMIN1(0.95, 1.0 - APFMX * (1.0-EXP(-PCB*TDMAS)))
	    ROOTPF=AMIN1(MAX_ROOTPF, 1.0 - APFMX * (1.0-EXP(-PCB*TDMAS)))
	    ROOTDWDT=FDWDT * ROOTPF
          Growth%GRORT   = ROOTDWDT * 100.

c     MJ: 2006/09/12:
c     changed this:
c	ROOTDM=ROOTDM * ROOTDWDT
c     to this:
          ROOTDWDT = AMAX1(ROOTDWDT, 0.000001)
          ROOTDM=ROOTDM + ROOTDWDT
          ROOTDM = AMAX1(ROOTDM, 0.000001)

c         Calculate AERDWDT, ensuring that it is at least 0.
          AERDWDT=AMAX1(FDWDT-ROOTDWDT,0.000001)

c         Update Aerial dry mass
          AERLDM=AERLDM+AERDWDT

c         Set stalk partitioning factor to maximum from cultivar file
c          STKPF=STKPFMAX

C     sdb added 16/10/2003
c     MJ, Feb 2010: replaced hard-coded value with STKPFMAX, as was meant to be:
c     MJ, Jun 2017: partitioning fraction exceeds 1 between 9 and 10 deg C.
c       I don't see anywhere that stkpf is restricted to 0.0-1.0!
	    stkpf = STKPFMAX + exp(-0.5 * (((tempmx+tempmn)/2.0)-7.5))

c         If the stalk has not yet emerged, 0 biomass is allocated to the
c         stalk:
c         ::::::
c         CNB/AS-Apr2002-partitioning to stalk now triggered by CHUPI
	    IF (LEAFTT .LT. Chupibase) THEN
              STKPF=0.000
          ELSE
c             Stalk emergence phenological phase
              IF (CaneCrop%GROPHASE .LT. 3) THEN
                  CaneCrop%GROPHASE = 3
                  STGDOY(4) = Control%YRDOY
              ENDIF
          ENDIF

c         Calculate daily increment in stalk mass
	    DSTKDT=AERDWDT*STKPF

c         Update stalk drymass
	    STKDM=STKDM+DSTKDT

c         MJ: changed 1 to 1. in following eq.
!DAILY INCREMENT IN NON-STALK (TOPS AND TRASH)
	    DNSTKDT=AERDWDT*(1.-STKPF)

c     NONSTALK IS PROPORTIONALLY DIVIDED BETWEEN TRASH AND TOPS
c     CALCULATE TOPDM AND TRASHDM IN THE FOLLOWING BLOCK
c     :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c         Update non-stalk aerial dry mass
	    NNSTK = AMAX1(DNSTKDT + TRASDM + TOPDM,0.0)

c         Calculate trash amount
!          TRASDM = 0.52 * NNSTK *  (TMEANDEDLF / (TMEANLF+0.00001))
!          MJ, June 2015: changed 0.52 to 0.70 based on 08RE14 experimental data
	    TRASDM = 0.70 * NNSTK *  (TMEANDEDLF / (TMEANLF+0.00001))
	    IF (TRASDM.GT.NNSTK) TRASDM=NNSTK
	    IF (NNSTK.LT.5.0) TRASDM=0.0

c         Set mass of tops to non-stalk matter minus trash
c         i.e. TOPDM is above ground non-stalk living biomass...
c         which is green leaves and meristems
	    TOPDM = NNSTK - TRASDM
c     :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

c     Calculate sucrose partitioning
c     :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
	    MEANTEMP=(TEMPMN+TEMPMX)/2.0

c         Following line no longer active
c  	    DELTT = DELTTMAX*(1-EXP(0.2*(MEANTEMP-TBDELTT)))

c     LOGISTIC TEMP FUNCTION
	    FTDELTA=1.0/(1.0+EXP(ftcon*(MEANTEMP-TBFT)))

c     FWDELTA=(1-SWDF9) Another option that works quite well (EP1=3 instead of EP1=2)
      	FWDELTA=(1-SWDF2)**SWDF2AMP

      	DELTA=DELTTMAX*(FWDELTA+FTDELTA-FWDELTA*FTDELTA)
	    IF (DELTA.LT.0.01) DELTA=0.01

c     NOTE: STKDM *MUST* be updated to today's integrated value before running this code.
	    IF (STKDM.LT.(SUCA/DELTA)) THEN
c             SUCPF IS THE PARTITIONING FRACTION TO SUCROSE IN THE STALK
c             ...and is only used as an output, if at all.
              SUCPF=DELTA*STKDM
              SUCEQ=(DELTA/2)*STKDM**2
	    ELSE
              SUCPF=SUCA
              SUCEQ=SUCA*(STKDM-0.5*SUCA/DELTA)
	    ENDIF

          SUCCAP = (SUCEQ-SUCMAS)
          NSUCCAP = (STKDM-SUCEQ)-NSUCMAS
          DSUCTEMP = SURCON*SUCCAP
          DNSUCTEMP=SURCON*NSUCCAP

          IF (DNSUCTEMP .GT. DSTKDT) DNSUCTEMP = DSTKDT
          IF (DNSUCTEMP .LT. 0.) DNSUCTEMP = 0.00

          DNSUCDT = DNSUCTEMP
          DSUCDT  = DSTKDT - DNSUCDT
          DSUCyest = DSUCDT
          SUCMAS  = SUCMAS + DSUCDT
          NSUCMAS = NSUCMAS + DNSUCDT

c     MJ, June 2015:
c     Partitioning of fibre into cellulosic products
c     Currently a static partitioning.
c     Non-sucrose mass: NSUCMAS
c     Leaf mass: TOPDM
c     Trash mass: TRASDM
      CELLSE_DM = (NSUCMAS + TOPDM + TRASDM) * FF_CELLSE

c     :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

! CNB 9/2002 TO OUTPUT DELTA WE NAME IT DELTAOUT BECAUSE THE ARE OTHER DELTAs
	deltaout=delta

 122   FORMAT(F5.0,10F8.3)

c     *******************************************************************

c     MJ: copy values from local variables to common blocks:
c     ::::::::::::::::::::::::::::::::::::::::::::::::::::::
c     WATER:
c     ::::::
       WaterBal%SWDF2 = SWDF2

c     PART:
c     :::::
      Part%AERDWDT   = AERDWDT
      Part%AERLDM    = AERLDM
      Part%APFMX     = APFMX
      Part%CHUPIBASE = CHUPIBASE
      Part%DELTAOUT  = DELTAOUT
      Part%DELTTMAX  = DELTTMAX
      Part%FDWDT     = FDWDT
      Part%FTCON     = FTCON
      Part%NSUCMAS   = NSUCMAS
      Part%PCB       = PCB
      Part%STKDM     = STKDM
      Part%STKPFMAX  = STKPFMAX
      Part%SUCA      = SUCA
      Part%SUCMAS    = SUCMAS
      Part%SURCON    = SURCON
      Part%SWDF2AMP  = SWDF2AMP
c      Part%TBDELTT   = TBDELTT
      Part%TBFT      = TBFT
      Part%TDMAS     = TDMAS
      Part%TOPDM     = TOPDM

c     CLIMT:
c     ::::::
      Climate%TEMPMN = TEMPMN
      Climate%TEMPMX = TEMPMX

c     Canecrop:
c     :::::::::
      CaneCrop%TMEANDEDLF = TMEANDEDLF
      CaneCrop%TMEANLF    = TMEANLF

c     Outblock:
c     :::::::::
c     Note: rootdm and trasdm should be in different block...
c      Out%CHUPI  = CHUPI
      Out%ROOTDM = ROOTDM
      Out%TRASDM = TRASDM

c     End of Dyanamic = RATE
      ENDIF
c     *******************************************************************

      RETURN
      END

ext: for


text:
!==============================================================================!
  subroutine Calculate_Sgs_Hybrid(grid)
!------------------------------------------------------------------------------!
!   Calculates SGS stresses and turbulent viscosity for 'LES'.                 !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Const_Mod
  use Flow_Mod
  use Comm_Mod
  use Les_Mod
  use Rans_Mod
  use Grid_Mod
!------------------------------------------------------------------------------!
!   near(c) is the number of corresponding cell on the nearest wall.
!   In case that, in parallel executions, the subdomain does not have
!   any nearwall cells, the near(c) is zero.
!   near(c) is calculated in NearWallCells.f90, only ones in the beginig
!   of a simulation.
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
!-----------------------------------[Locals]-----------------------------------!
  integer :: c
  real    :: lf
!==============================================================================!

  do c = 1, grid % n_cells
    lf = grid % vol(c) ** ONE_THIRD
    vis_t_sgs(c) = density                &
                 * (lf*lf)                &          ! delta^2
                 * c_dyn(c)               &          ! c_dynamic
                 * shear(c)
  end do

  call Comm_Mod_Exchange_Real(grid, vis_t_sgs)

  end subroutine

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


text:
        !COMPILER-GENERATED INTERFACE MODULE: Thu Feb  5 15:03:42 2015
        MODULE MC23AD__genmod
          INTERFACE
            SUBROUTINE MC23AD(N,ICN,A,LICN,LENR,IDISP,IP,IQ,LENOFF,IW,  &
     &IW1)
              INTEGER(KIND=4) :: LICN
              INTEGER(KIND=4) :: N
              INTEGER(KIND=4) :: ICN(LICN)
              REAL(KIND=8) :: A(LICN)
              INTEGER(KIND=4) :: LENR(N)
              INTEGER(KIND=4) :: IDISP(2)
              INTEGER(KIND=4) :: IP(N)
              INTEGER(KIND=4) :: IQ(N)
              INTEGER(KIND=4) :: LENOFF(N)
              INTEGER(KIND=4) :: IW(N,5)
              INTEGER(KIND=4) :: IW1(N,2)
            END SUBROUTINE MC23AD
          END INTERFACE
        END MODULE MC23AD__genmod

max_stars_repo_name: mtsafarzadeh/FLASHverES

max_issues_repo_name: mtsafarzadeh/FLASHverES

max_forks_repo_name: mtsafarzadeh/FLASHverES


module mesh_info_M

    USE global_parameters_M

    implicit none

!---Mesh information
    integer  :: nodes_per_elem ! Nodes per element
    integer  :: num_elem ! Number of elements in the mesh
    integer  :: max_num_nodes
    integer  :: ndf ! ndf
    integer  :: matrix_length
    integer  :: num_elem_external
    real(dp) :: elem_size
    real(dp) :: mass_elem ! mass per element
    real(dp) :: Area_Core, Area_Pipe, Area_Plenum

!---Mesh arrays
    integer, allocatable  :: conn_matrix(:,:)
    real(dp), allocatable :: elem_lengths(:) ! length of elements
    real(dp), allocatable :: elem_node_lengths(:)
    real(dp), allocatable :: global_coord(:,:)
    real(dp), allocatable :: power_soln_test(:)
    real(dp), allocatable :: area_variation(:,:)

    integer :: Fuel_Core_End
    integer :: Fuel_Core_Start
    integer :: Fuel_Outlet_End
    integer :: Fuel_Inlet_Start
    !---Heat exchanger
    real(dp) :: Area_Heat_Exchanger
    integer :: Heat_Exchanger_Start
    integer :: Heat_Exchanger_End

end module mesh_info_M


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

      MODULE UITZICHT_DATA
!***********************************************************************
!     SPECTRALE GEGEVENS VAN 400 TOT 700 NM MET STAPPEN VAN 5 NM
!***********************************************************************
!
      REAL(4), DIMENSION(61) :: AWATER, BWATER, CHLSPE, PLANCK
!
!     SPECIFIEKE ABSORPTIE ALGEN
!
      DATA CHLSPE / 0.685, 0.781, 0.828, 0.883, 0.913,
     1              0.939, 0.973, 1.001, 1.000, 0.971,
     2              0.944, 0.928, 0.917, 0.902, 0.870,
     3              0.839, 0.798, 0.773, 0.750, 0.717,
     4              0.688, 0.645, 0.618, 0.582, 0.528,
     5              0.504, 0.474, 0.444, 0.416, 0.384,
     6              0.357, 0.321, 0.294, 0.273, 0.276,
     7              0.278, 0.281, 0.279, 0.270, 0.252,
     8              0.256, 0.262, 0.261, 0.268, 0.281,
     9              0.299, 0.316, 0.328, 0.329, 0.337,
     A              0.361, 0.397, 0.457, 0.529, 0.556,
     B              0.534, 0.485, 0.411, 0.334, 0.270,
     C              0.215/
!
!     SPECTRALE VERDELING INVALLEND LICHT
      DATA PLANCK / 0.18380, 0.18836, 0.19280, 0.19710, 0.20134,
     1              0.20542, 0.20937, 0.21320, 0.21689, 0.22046,
     2              0.22389, 0.22720, 0.23037, 0.23342, 0.23630,
     3              0.23912, 0.24177, 0.24430, 0.24671, 0.24898,
     4              0.25115, 0.25320, 0.25511, 0.25691, 0.25860,
     5              0.26018, 0.26166, 0.26300, 0.26427, 0.26540,
     6              0.26650, 0.26745, 0.26830, 0.26910, 0.26978,
     7              0.27037, 0.27089, 0.27132, 0.27167, 0.27194,
     8              0.27214, 0.27226, 0.27232, 0.27231, 0.27223,
     9              0.27209, 0.27189, 0.27163, 0.27130, 0.27094,
     A              0.27050, 0.27004, 0.26950, 0.26894, 0.26833,
     B              0.26767, 0.26697, 0.26620, 0.26545, 0.26464,
     C              0.26379/
!
!     VERSTROOIING WATER
!
      DATA BWATER / 0.005290, 0.005025, 0.004776, 0.004543, 0.004323,
     1              0.004117, 0.003922, 0.003739, 0.003567, 0.003404,
     2              0.003250, 0.003105, 0.002968, 0.002838, 0.002715,
     3              0.002599, 0.002488, 0.002384, 0.002285, 0.002191,
     4              0.002102, 0.002017, 0.001936, 0.001860, 0.001787,
     5              0.001717, 0.001651, 0.001588, 0.001528, 0.001471,
     6              0.001416, 0.001364, 0.001314, 0.001267, 0.001221,
     7              0.001178, 0.001136, 0.001097, 0.001059, 0.001022,
     8              0.000987, 0.000954, 0.000922, 0.000891, 0.000862,
     9              0.000833, 0.000806, 0.000780, 0.000755, 0.000731,
     A              0.000708, 0.000686, 0.000664, 0.000644, 0.000624,
     B              0.000605, 0.000587, 0.000569, 0.000552, 0.000536,
     C              0.00052/
!
!     ABSORPTIE WATER
!
      DATA AWATER / 0.00576, 0.00617, 0.00669, 0.00727, 0.00790,
     1              0.00854, 0.00918, 0.00980, 0.01039, 0.01093,
     2              0.01144, 0.01193, 0.01241, 0.01293, 0.01353,
     3              0.01426, 0.01520, 0.01645, 0.01810, 0.02032,
     4              0.02380, 0.02818, 0.03294, 0.03768, 0.04212,
     5              0.04604, 0.04937, 0.05211, 0.05438, 0.05639,
     6              0.05849, 0.06108, 0.06470, 0.06999, 0.07768,
     7              0.08886, 0.10519, 0.12723, 0.15488, 0.18723,
     8              0.22238, 0.25724, 0.28740, 0.30690, 0.31434,
     9              0.31761, 0.31940, 0.32232, 0.32862, 0.33985,
     A              0.35650, 0.37763, 0.40052, 0.42033, 0.43334,
     B              0.44359, 0.45513, 0.47316, 0.50342, 0.55154,
     C              0.62200/
      END MODULE UITZICHT_DATA

      SUBROUTINE UIT_ZI ( DIEP1 , DIEP2 , ANGLE , C_GL1 , C_GL2 ,
     &                    C_DET , HELHUM, TAU   , CORCHL, CHLORO,
     &                    DETRIT, GLOEIR, AH_380, SECCHI, D_1   ,
     &                    EXTPAR, EXTP_D, DOSECC)
!>\file
!>       Transparency due to Chlorophyll, detritus, inorganics and humic accids

!***********************************************************************
!     SUBROUTINE UIT_ZI, GEBASEERD OP HET MODEL UITZICHT VAN
!     H. BUITEVELD, RIZA, POSTBUS 17, 8200 AA LELYSTAD (TEL 03200-70737)
!
!     BEREKEND DOORZICHT EN EXTINKTIE OP BASIS VAN CHLOROFYL, DETRITUS
!     GLOEIREST, ABSORPTIE HUMUSZUREN BIJ 380 NM EN SPECTRA OM DE 5 NM
!     VAN DE ABSORPTIE EN VERSTROOIING VAN WATER, DE VERDELING DE
!     INVALLEND LICHT EN DE SPECIFIEKE SPECIFIEKE ABSORPTIE ALGEN
!
!     CHANGES IN THE MODULE:
!     DATE   AUTHOR          DESCRIPTION
!     ------ --------------- -------------------------------------------
!     971119 Jan van Beek    Extended Ascii characters verwijderd
!     970127 Rik Sonneveldt  In subroutine BEP_D beveiliging tegen CHLORO
!                            < 0 ingebouwd (op verzoek van Maarten Ouboter).
!     911125 Andre Hendriks  Code beter leesbaar gemaakt, en variabele-
!                            namen langer dan 6 letters vervangen door korte
!                            namen.
!     910926 WOLF MOOIJ      IMPLEMENTATION IN DELWAQ-BLOOM
!                            BELANGRIJKSTE VERANDERING: INPLAATS VAN
!                            TOTAAL ZWEVEND STOF WORDT NU DIRECT
!                            DE DETRITUS CONCENTRATIE INGELEZEN
!***********************************************************************
!***********************************************************************
!     ARGUMENTS:
!     NAME   TYPE SIZE   I/O DESCRIPTION
!     ------ ---- ------ --- -------------------------------------------
!     ANGLE   R*4        IN  FUNCTION, DEFAULT constant = 30 x
!     AH_380  R*4        IN  PARAMETER: EXTINCTIE HUMUSZUREN (1/m)
!     CHLORO  R*4        IN  PARAMETER: CHLOROPHYL (mg/m3)
!     CORCHL  R*4        IN  CONSTANT, DEFAULT = 2.5
!     C_DET   R*4        IN  CONSTANT, DEFAULT = 0.026
!     C_GL1   R*4        IN  CONSTANT, DEFAULT = 0.73
!     C_GL2   R*4        IN  CONSTANT, DEFAULT = 1.0
!     DETRIT  R*4        IN  PARAMETER: GESUSPENDEERD DETRITUS (gDW/m3)
!     DIEP1   R*4        IN  CONSTANT, DEFAULT = 1.0 (m)
!     DIEP2   R*4        IN  CONSTANT, DEFAULT = 1.2 (m)
!     D_1     R*4        OUT DUMMY PARAMETER D 10% TRANSMIS 560 NM (m)
!     EXTPAR  R*4        OUT PARAMETER: EXTINCTIE OP 1M (1/m)
!     EXTP_D  R*4        OUT DUMMY PARAMETER EXTINCTIE OP D_1 (1/m)
!     GLOEIR  R*4        IN  PARAMETER: ANORGANISCH ZWEVEND STOF (gDW/m3)
!     HELHUM  R*4        IN  CONSTANT, DEFAULT = 0.014
!     SECCHI  R*4        OUT PARAMETER: DOORZICHT (m)
!     TAU     R*4        IN  CONSTANT, DEFAULT = 7.8
!***********************************************************************
!***********************************************************************
!     COMMON VARIABLES:
!     NAME   TYPE SIZE   DESCRIPTION
!     ------ ---- ------ -----------------------------------------------
!     AWATER  R*4 (61)   ABSORPTIE WATER
!     BWATER  R*4 (61)   VERSTROOING WATER
!     CHLSPE  R*4 (61)   SPECIFIEKE ABSORPTIE ALGEN
!     PLANCK  R*4 (61)   VERDELING INVALLEND LICHT
!***********************************************************************
!***********************************************************************
!     LOCAL VARIABLES:
!     NAME   TYPE SIZE   DESCRIPTION
!     ------ ---- ------ -----------------------------------------------
!     A       R*4
!     A_CHL   R*4
!     A_DET   R*4
!     A_HUM   R*4
!     B       R*4
!     B_CHL   R*4
!     B_GL    R*4
!     C_CHL   R*4
!     C_GL    R*4
!     C_MU    R*4
!     D_2     R*4
!     EXT_KI  R*4
!     I_550   I*4
!     LAMBDA  I*4
!     SOM_C   R*4
!     SOM_D1  R*4
!     SOM_D2  R*4
!     SOM_H   R*4
!     S_D1    R*4
!     S_D2    R*4
!     TELLER  I*4
!     ZW_STF  R*4
!***********************************************************************
!
!     include '..\inc\ioblck.inc'
      USE UITZICHT_DATA
!
      REAL    A     , A_CHL , A_DET , A_HUM , AH_380,
     1        ANGLE , B     , B_CHL , B_GL  ,
     2        C_CHL , C_DET , C_GL  , C_GL1 , C_GL2 ,
     3        C_MU  , CHLORO, CORCHL, D_1   , D_2   ,
     4        DETRIT, DIEP1 , DIEP2 , EXT_KI, EXTP_D,
     5        EXTPAR, GLOEIR, HELHUM, S_D1  , S_D2  ,
     6        SECCHI, SOM_C , SOM_D1, SOM_D2, SOM_H ,
     7        TAU   , ZW_STF
      INTEGER I_550 , LAMBDA, TELLER
      LOGICAL DOSECC
!
      IF ( CHLORO .GE. 0.0 .OR. DETRIT .GE. 0.0 .OR.
     1     GLOEIR .GE. 0.0) THEN
         C_MU = COS ( ANGLE * 0.0174533)
         IF ( DOSECC )
     J   CALL BEP_D  ( C_GL1 , C_GL2 , C_DET , HELHUM, CORCHL,
     1                 C_MU  , CHLORO, DETRIT, GLOEIR, AH_380,
     2                 D_1   , D_2   )
         I_550 = ((550-400)/5)+1
!
!        CHLOROFYL BUNDELVERZWAKKING
!
         C_CHL  = ( 0.058 + 0.018 * CHLORO)* CHLSPE (I_550)
         C_CHL  = ( C_CHL + 0.12 * ( CHLORO**0.63)) * CORCHL
         SOM_D1 = 0.0
         SOM_D2 = 0.0
         S_D1   = 0.0
         S_D2   = 0.0
         SOM_C  = 0.0
         SOM_H  = 0.0
         DO 100 TELLER  =  1, 61
            LAMBDA = 400+(TELLER-1)*5
!
!           HUMUSZUREN ABOSORPTIE
!
            A_HUM = AH_380 * EXP ( -HELHUM * (LAMBDA-380.0))
!
!           ALGEN ABSOROPTIE EN VERSTROOIING
!
            IF ( CHLORO .LT. 0.000001) THEN
               A_CHL = 0.0
               B_CHL = 0.0
            ELSE
               A_CHL = ( 0.058 + 0.018 * CHLORO) * CHLSPE ( TELLER) *
     1                 CORCHL
               B_CHL = C_CHL - A_CHL
            ENDIF
!
!           GLOEIREST EN DETRITUS
!
            ZW_STF = GLOEIR + DETRIT
            C_GL   = C_GL1 * (( ZW_STF**C_GL2)) * (400.0 / LAMBDA)
            A_DET  = C_DET * DETRIT * ( 400.0/ LAMBDA)
            B_GL   = C_GL - A_DET
!
!           TOTAAL ABSORPTIE EN VERSTROOIING BIJ LAMBDA
!
            A = AWATER ( TELLER) + A_HUM + A_DET + A_CHL
            B = BWATER ( TELLER) + B_GL + B_CHL
!
!           EXTINKTIE BIJ LAMBDA
!
            EXT_KI = 1 / C_MU *
     1               SQRT ( A * A + (0.425 * C_MU - 0.19) * A * B)
            SOM_D1 = SOM_D1 + PLANCK ( TELLER) * EXP ( -EXT_KI * DIEP1)
            SOM_D2 = SOM_D2 + PLANCK ( TELLER) * EXP ( -EXT_KI * DIEP2)
            IF ( DOSECC ) THEN
            S_D1   = S_D1   + PLANCK ( TELLER) * EXP ( -EXT_KI * D_1)
            S_D2   = S_D2   + PLANCK ( TELLER) * EXP ( -EXT_KI * D_2)
            SOM_H  = SOM_H  + PLANCK ( TELLER) * EXP ( -(A+B)  * D_1)
            SOM_C  = SOM_C  + PLANCK ( TELLER) * EXP ( -(A+B)  * D_2)
            ENDIF
  100    CONTINUE
         EXTPAR = ( 1.0 / ( DIEP1 - DIEP2) * LOG ( SOM_D2 / SOM_D1))
         IF ( DOSECC ) THEN
         EXTP_D = ( 1.0 / ( D_1   - D_2  ) * LOG ( S_D2   / S_D1  ))
         SOM_C  = ( 1.0 / ( D_1   - D_2  ) * LOG ( SOM_C  / SOM_H ))
         SECCHI = TAU / ( EXTP_D + SOM_C)
         ELSE
         SECCHI = -1.0
         ENDIF
      ENDIF
      RETURN
      END

      SUBROUTINE BEP_D ( C_GL1 , C_GL2 , C_DET , HELHUM, CORCHL,
     1                   C_MU  , CHLORO, DETRIT, GLOEIR, AH_380,
     2                   D_1   , D_2   )
!***********************************************************************
!     BEPAALD DIEPTE WAAR 10 % VAN HET LICHT OVER IS BIJ 550  NM
!***********************************************************************
!***********************************************************************
!     ARGUMENTS:
!     NAME   TYPE SIZE   I/O DESCRIPTION
!     ------ ---- ------ --- -------------------------------------------
!     AH_380  R*4        IN  PARAMETER: EXTINCTIE HUMUSZUREN (1/M)
!     CHLORO  R*4        IN  PARAMETER: CHLOROPHYL (fG/L)
!     CORCHL  R*4        IN  CONSTANT,  DEFAULT  =  2.5
!     C_DET   R*4        IN  CONSTANT,  DEFAULT  =  0.026
!     C_GL1   R*4        IN  CONSTANT,  DEFAULT  =  0.73
!     C_GL2   R*4        IN  CONSTANT,  DEFAULT  =  1.0
!     C_MU    R*4        IN  CONSTANT,  COSINUS VAN ANGLE (DEFAULT  =  30)
!     DETRIT  R*4        IN  PARAMETER: GESUSPENDEERD DETRITUS (MG/L)
!     D_1     R*4        OUT DUMMY PARAMETER D 10% TRANSMIS 560 NM (M)
!     D_2     R*4        OUT DUMMY PARAMETER D_1 + 0.1 (M)
!     GLOEIR  R*4        IN  PARAMETER: ANORGANISCH ZWEVEND STOF (MG/L)
!     HELHUM R*4        IN  CONSTANT,  DEFAULT  =  0.014
!***********************************************************************
!***********************************************************************
!     COMMON VARIABLES:
!     NAME   TYPE SIZE   DESCRIPTION
!     ------ ---- ------ -----------------------------------------------
!     AWATER  R*4 (61)   ABSORPTIE WATER
!     BWATER  R*4 (61)   VERSTROOING WATER
!     CHLSPE  R*4 (61)   SPECIFIEKE ABSORPTIE ALGEN
!     PLANCK  R*4 (61)   VERDELING INVALLEND LICHT
!***********************************************************************
!***********************************************************************
!     LOCAL VARIABLES:
!     NAME   TYPE SIZE   DESCRIPTION
!     ------ ---- ------ -----------------------------------------------
!     A      R*4
!     A_CHL  R*4
!     A_DET  R*4
!     A_HUM  R*4
!     B      R*4
!     B_CHL  R*4
!     B_GL   R*4
!     C_CHL  R*4
!     C_GL   R*4
!     EXT_KI R*4
!     I_550  I*4
!     LAMBDA I*4
!     TELLER I*4
!     ZW_STF R*4
!***********************************************************************

      USE UITZICHT_DATA
!
      REAL    A     , A_CHL , A_DET , A_HUM , AH_380,
     1        B     , B_CHL , B_GL  , C_CHL , C_DET ,
     2        C_GL  , C_GL1 , C_GL2 , C_MU  , CHLORO,
     3        CORCHL, D_1   , D_2   , DETRIT, EXT_KI,
     4        GLOEIR, HELHUM, ZW_STF
      INTEGER I_550 , LAMBDA, TELLER
!
      I_550 = ((550-400)/5)+1
!
!     Beveiliging tegen negatieve waarde CHLORO (RS27jan97 voor Maarten O.)
!
      CHLORO = MAX(0.0,CHLORO)
!
!     CHLOROFYL BUNDEL VERZWAKKING
!
      C_CHL = ( 0.058 + 0.018 * CHLORO) * CHLSPE ( I_550)
      C_CHL = ( C_CHL + 0.12 * ( CHLORO**0.63)) * CORCHL
      TELLER = (( 560 - 400) / 5) + 1
      LAMBDA = 400 + ( TELLER - 1) * 5
!
!     ABSORPTIE HUMUSZUREN
!
      A_HUM = AH_380 * EXP ( -HELHUM * ( LAMBDA - 380))
!
!     CHLOROFYL
!
      IF ( CHLORO .LT. 0.000001) THEN
         A_CHL = 0.0
         B_CHL = 0.0
      ELSE
         A_CHL = ( 0.058 + 0.018 * CHLORO) * CHLSPE ( TELLER) * CORCHL
         B_CHL = C_CHL - A_CHL
      ENDIF
!
!     GLOEIREST EN DETRITUS
!
      ZW_STF = GLOEIR + DETRIT
      C_GL   = C_GL1 * ( ZW_STF**C_GL2) * ( 400.0 / LAMBDA)
      A_DET  = C_DET * DETRIT * ( 400.0 / LAMBDA)
      B_GL   = C_GL - A_DET
      A      = AWATER ( TELLER) + A_HUM + A_DET + A_CHL
      B      = BWATER ( TELLER) + B_GL  + B_CHL
      EXT_KI = 1 / C_MU *
     1         SQRT (( A * A + ( 0.425 * C_MU - 0.19) * A * B))
!
!     DIEPTE 10 % TRANSMISSIE 560 NM
!
      D_1 = 2.3 / EXT_KI
      D_2 = D_1 + 0.1
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      program aa04

      real x(10), y(10), z(10), u(10,10)

      do i = 1, 10
         x(11-i) = 0.
      enddo

      do i = 1, 10
         x(11-i) = real(i)
      enddo

      end


text:
submodule (stats) stats_sm

contains

    module real(8) function get_bin_size(this)
        class(histogram),  intent(in) :: this
        get_bin_size = (this%bin_max-this%bin_min)/this%nbins
    end function get_bin_size

    ! Create a histogram given the number of bins
    module subroutine create_hist1(this, col, mask, nbins)
        class(histogram), intent(out) :: this
        real(8), intent(in)           :: col(:)
        logical, intent(in), optional :: mask(:)
        integer, intent(in), optional :: nbins
        logical, allocatable :: hist_mask(:)      ! Mask selecting those entries in col that are in the current interval
        real(8)              :: min_col, max_col, bin_size
        real(8), allocatable :: col_data(:)
        integer              :: i

        if (present(nbins)) then
            this%nbins = nbins
        end if

        ! Create local contiguous copy of the data being processed
        if (present(mask)) then
            col_data = pack(col,mask)
        else
            col_data = col
        end if

        ! Initialize the frequency
        allocate(this%frequency(0:this%nbins-1))

        ! Calculate the range of data in 'col' for each bin
        max_col = maxval(col_data)
        min_col = minval(col_data)
        if (max_col == min_col) then
            max_col = max_col + 0.5d0
            min_col = min_col - 0.5d0
        end if
        bin_size = (max_col-min_col)/this%nbins
        do i=0, this%nbins-1
            hist_mask  = (min_col + i*bin_size) < col_data .and. col_data <= (min_col + (i+1)*bin_size)
            this%frequency(i) = count(hist_mask)/bin_size
        end do
        this%bin_max = max_col
        this%bin_min = min_col
    end subroutine create_hist1

    ! Create a histogram given the bin size and range
    module subroutine create_hist2(this, col, bin_size, bin_range, mask)
        class(histogram), intent(out) :: this
        real(8), intent(in)           :: col(:)
        real(8), intent(in)           :: bin_size
        real(8), intent(in)           :: bin_range(2)
        logical, intent(in), optional :: mask(:)
        real(8), allocatable :: col_data(:)
        logical, allocatable :: hist_mask(:)      ! Mask selecting those entries in col that are in the current interval
        integer              :: i

        ! Create local contiguous copy of the data being processed
        if (present(mask)) then
            col_data = pack(col,mask)
        else
            col_data = col
        end if

        ! Initialize the frequency
        this%bin_max = bin_range(2)
        this%bin_min = bin_range(1)
        this%nbins = nint((this%bin_max-this%bin_min)/bin_size)
        allocate(this%frequency(0:this%nbins-1))

        ! Calculate the range of data in 'col' for each bin
        do i=0, this%nbins-1
            hist_mask  = (this%bin_min + i*bin_size) < col_data .and. &
                          col_data <= (this%bin_min + (i+1)*bin_size)
            this%frequency(i) = count(hist_mask)/bin_size
        end do
    end subroutine create_hist2

    ! Create a histogram given the bin size and range
    module subroutine create_hist3(this, col, bin_size, bin_range, mask)
        class(histogram), intent(out) :: this
        real(8), intent(in)           :: col(:)
        integer, intent(in)           :: bin_size
        integer, intent(in)           :: bin_range(2)
        logical, intent(in), optional :: mask(:)
        real(8)          :: r_bin_range(2)

        ! Check the range
        if (mod(bin_range(2)-bin_range(1),bin_size) /= 0) then
            stop '***Error[create_hist3]: inconsistent bin definition'
        end if

        r_bin_range = real(bin_range,8)
        call this%create(col, real(bin_size,8), r_bin_range, mask)
        this%int_bin_size = .true.
    end subroutine create_hist3

    ! Write histogram
    module subroutine write_histogram(this, u)
        class(histogram), intent(in) :: this
        integer, intent(in)          :: u
        integer :: i
        real(8) :: bin_size

        bin_size = (this%bin_max - this%bin_min)/this%nbins
        associate (bin_start=>this%bin_min, freq=>this%frequency)
            if (this%int_bin_size) then
                do i=0, this%nbins-1
                    write(u,'(2(a,i0),a,f0.2,a,f0.2)') '(',nint(bin_start+i*bin_size), &
                    '->',nint(bin_start+(i+1)*bin_size),'],',bin_start+i*bin_size+bin_size/2,',',freq(i)
                end do
            else
                do i=0, this%nbins-1
                    write(u,'(2(a,es9.2),a,es9.2,a,es9.2)') '(',bin_start+i*bin_size,'->', &
                    bin_start+(i+1)*bin_size,'],',bin_start+i*bin_size+bin_size/2,',',freq(i)
                end do
            end if
        end associate
    end subroutine write_histogram

    ! Calculate mean and standard deviation
    module subroutine calculate_stats(col, mean, mask, standard_deviation, skewness)
        real(8), intent(in)            :: col(:)
        real(8), intent(out)           :: mean
        logical, intent(in), optional  :: mask(:)
        real(8), intent(out), optional :: standard_deviation
        real(8), intent(out), optional :: skewness
        real(8) :: mv(3), ssdev
        integer :: i, n

        mv = 0
        if (present(mask)) then
            n = count(mask)
            if (n == 0) return

            !$omp parallel do shared(col, mask) reduction(+:mv)
            do i=1,size(col)
                if (mask(i)) then
                    mv(1) = mv(1) + col(i)
                    if (present(standard_deviation) .or. present(skewness)) then
                        mv(2) = mv(2) + col(i)**2
                        if (present(skewness)) then
                            mv(3)= mv(3)+ col(i)**3
                        end if
                    end if
                end if
            end do
            !$omp end parallel do

        else
            n = size(col)
            !$omp parallel do shared(col) reduction(+:mv)
            do i=1,n
                mv(1) = mv(1) + col(i)
                if (present(standard_deviation) .or. present(skewness)) then
                    mv(2) = mv(2) + col(i)**2
                    if (present(skewness)) then
                        mv(3)= mv(3)+ col(i)**3
                    end if
                end if
            end do
            !$omp end parallel do

        end if

        ! Mean and standard deviation
        mv = mv/n
        mean = mv(1)
        if (present(standard_deviation) .or. present(skewness)) then
            ssdev = sqrt(abs(mv(2) - mv(1)**2))  ! Rounding could make it < 0 so use abs to be safe
            if (present(skewness)) then
                skewness = (mv(3) - mean*(3*ssdev**2 + mean**2))/ssdev**3
            end if
            if (present(standard_deviation)) then
                standard_deviation = ssdev
            end if
        end if

    end subroutine calculate_stats

    module function random() result(ran)
!=======================================================================
!
!    This is an adapted version of subroutine RANECU written by
!    F. James (Comput. Phys. Commun. (1990) 60, 329-344), which
!    has been modified to give a single random number at each call.
!    The seeds iseed1 and iseed2 must be initialized in the main
!    program and transferred through the named common block /REED/.
!
!=======================================================================

        real(8) :: ran

        ! ---   Declaracion de variables   ---
        integer :: i1, i2, iz

        ! Saved data. Note that inline initialization implies 'save'
        integer :: iseed1 = 860934
        integer :: iseed2 = 542039

        ! ---   Parametros   ---
        real(8), parameter :: uscale=1.0d0/2.0d0**31

!     ---   Calculo del numero aleatorio   ---

        i1 = iseed1/53668
        iseed1 = 40014*(iseed1-i1*53668)-i1*12211

        if (iseed1 < 0) then
            iseed2 = iseed1+2147483563
        end if

        i2 = iseed2/52774
        iseed2 = 40692*(iseed2-i2*52774)-i2*3791

        if (iseed2 < 0) then
            iseed2 = iseed2+2147483399
        end if

        iz = iseed1-iseed2

        if (iz < 1) then
            iz = iz+2147483562
        end if

        ran = iz*uscale
        return

    end function random

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

    module function gasdev()
!=======================================================================
!
!     Returns a normally distributed deviate with zero mean and unit
!     variance.
!
!=======================================================================

        real(8) :: gasdev

        ! ---   Variable declaration   ---
        real(8)       :: fac, r
        real(8)       :: v1, v2

        ! ---   Use static memory for persistence   ---
        real(8), save :: gset = 0
        logical, save :: have_extra_deviate = .false.

!       ---   We don't have an extra deviate handy, so ...   ---
        if (.not. have_extra_deviate) then

            uniform_square: do

                ! ---   pick two uniform numbers in the square (-1,+1)

                v1 = 2*random()-1
                v2 = 2*random()-1

                ! ---   See if they are in the unit circle   ---

                r = v1*v1+v2*v2

                ! ---   and if they are not, try again   ---

                if (r < 1 .and. r > 0) exit uniform_square

            end do uniform_square

!           ---  Make Box-Muller polar transformation to get 2 normal deviates   ---

            fac = sqrt(-2*log(r)/r)

!           ---   return one and save the other for the next time   ---

            gset = v1*fac
            gasdev = v2*fac
            have_extra_deviate = .true.

        else

!           ---   We have an extra deviate so return it and unset flag   ---

            gasdev = gset
            have_extra_deviate = .false.

        end if

        return

    end function gasdev

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

end submodule stats_sm

max_stars_repo_name: sgeard/code

max_issues_repo_name: sgeard/code

max_forks_repo_name: sgeard/code


text:
      SUBROUTINE REP_ENTRY( S, REPORT )

!***********************************************************************
!* Writes Standard UDS Entry to Report File
!*
!* Language: Fortran
!*
!* Author: Stuart G. Mentzer
!*
!* Date: 1999/08/20
!***********************************************************************

      ! Headers
      INCLUDE 'uds_str.fi'
      INCLUDE 'uds_fxn.fi'
      INCLUDE 'file.fi'

      ! Arguments ______________________________________________________

      RECORD /UDS_SPEC/  S ! UDS_SPEC object

      RECORD /FILE/  REPORT ! Report file object

      ! Variables ______________________________________________________

      INTEGER IOS

      CHARACTER VEHICLE*30, CUTOFF*20

      ! Functions ______________________________________________________

      CHARACTER RJUST*12

      EXTERNAL RJUST

      ! Return if report file not open
      IF ( .NOT. REPORT.OPEN ) RETURN

      ! Set up entries
      VEHICLE = ' '
      WRITE( VEHICLE(1:4), '(I4)', IOSTAT=IOS ) S.YEAR
      VEHICLE(6:) = S.MAKE(:L_TRIM(S.MAKE))//' '//S.MODEL
      CUTOFF = ' '
      IF ( S.FCUT .GT. 0. ) THEN
        WRITE( CUTOFF, '(G20.6)', IOSTAT=IOS ) S.FCUT
      ELSE IF ( S.PREF .GT. 0. ) THEN
        WRITE( CUTOFF, '(G20.6)', IOSTAT=IOS ) S.PREF
      END IF

      ! Write the UDS entry
      WRITE( REPORT.UNIT,
     & '(A12,2X,5X,I1,5X,2X,A30,3X,A2,3X,A20,2X,A20,2X,A20)',
     & IOSTAT=IOS )
     & RJUST( S.TSTNAM ), S.VEHNO, VEHICLE, S.AXIS,
     & S.YTYP, S.XTYP, CUTOFF
      WRITE( REPORT.UNIT,
     & '(A12,2X,4X,A2,5X,2X,A30,3X,2X,3X,A20,2X,A20,2X,A20)',
     & IOSTAT=IOS )
     & RJUST( S.CURNAM ), S.SENLOC, S.SENATT,
     & S.YUNITS, S.XUNITS, S.STATUS

      RETURN
      END

max_stars_repo_name: DeadParrot/NHTSA-Tools

max_issues_repo_name: DeadParrot/NHTSA-Tools

max_forks_repo_name: DeadParrot/NHTSA-Tools


text:
      program t
      implicit none
! io-control-spec read stmt end label
      character*4::out
      integer::ierr
      open (95, file='infile', status='old')
      read (95, '(a)') out
      read (95, '(a)',end=100) out
      print *, out
      print *,'i wish this'

 100  print *,'test was successful'
      endprogram t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
      SUBROUTINE BOBYQB (N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,
     1  MAXFUN,XBASE,XPT,FVAL,XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,
     2  SL,SU,XNEW,XALT,D,VLAG,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(*),XL(*),XU(*),XBASE(*),XPT(NPT,*),FVAL(*),
     1  XOPT(*),GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),
     2  SL(*),SU(*),XNEW(*),XALT(*),D(*),VLAG(*),W(*)
C
C     The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, IPRINT and MAXFUN
C       are identical to the corresponding arguments in SUBROUTINE BOBYQA.
C     XBASE holds a shift of origin that should reduce the contributions
C       from rounding errors to values of the model and Lagrange functions.
C     XPT is a two-dimensional array that holds the coordinates of the
C       interpolation points relative to XBASE.
C     FVAL holds the values of F at the interpolation points.
C     XOPT is set to the displacement from XBASE of the trust region centre.
C     GOPT holds the gradient of the quadratic model at XBASE+XOPT.
C     HQ holds the explicit second derivatives of the quadratic model.
C     PQ contains the parameters of the implicit second derivatives of the
C       quadratic model.
C     BMAT holds the last N columns of H.
C     ZMAT holds the factorization of the leading NPT by NPT submatrix of H,
C       this factorization being ZMAT times ZMAT^T, which provides both the
C       correct rank and positive semi-definiteness.
C     NDIM is the first dimension of BMAT and has the value NPT+N.
C     SL and SU hold the differences XL-XBASE and XU-XBASE, respectively.
C       All the components of every XOPT are going to satisfy the bounds
C       SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when
C       XOPT is on a constraint boundary.
C     XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the
C       vector of variables for the next call of CALFUN. XNEW also satisfies
C       the SL and SU constraints in the way that has just been mentioned.
C     XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW
C       in order to increase the denominator in the updating of UPDATE.
C     D is reserved for a trial step from XOPT, which is usually XNEW-XOPT.
C     VLAG contains the values of the Lagrange functions at a new point X.
C       They are part of a product that requires VLAG to be of length NDIM.
C     W is a one-dimensional array that is used for working space. Its length
C       must be at least 3*NDIM = 3*(NPT+N).
C
C     Set some constants.
C
      HALF=0.5D0
      ONE=1.0D0
      TEN=10.0D0
      TENTH=0.1D0
      TWO=2.0D0
      ZERO=0.0D0
      NP=N+1
      NPTM=NPT-NP
      NH=(N*NP)/2
C
C     The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
C     BMAT and ZMAT for the first iteration, with the corresponding values of
C     of NF and KOPT, which are the number of calls of CALFUN so far and the
C     index of the interpolation point at the trust region centre. Then the
C     initial XOPT is set too. The branch to label 720 occurs if MAXFUN is
C     less than NPT. GOPT will be updated if KOPT is different from KBASE.
C
      CALL PRELIM (N,NPT,X,XL,XU,RHOBEG,IPRINT,MAXFUN,XBASE,XPT,
     1  FVAL,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,KOPT)
      XOPTSQ=ZERO
      DO 10 I=1,N
      XOPT(I)=XPT(KOPT,I)
   10 XOPTSQ=XOPTSQ+XOPT(I)**2
      FSAVE=FVAL(1)
      IF (NF .LT. NPT) THEN
          IF (IPRINT .GT. 0) PRINT 390
          GOTO 720
      END IF
      KBASE=1
C
C     Complete the settings that are required for the iterative procedure.
C
      RHO=RHOBEG
      DELTA=RHO
      NRESC=NF
      NTRITS=0
      DIFFA=ZERO
      DIFFB=ZERO
      ITEST=0
      NFSAV=NF
C
C     Update GOPT if necessary before the first iteration and after each
C     call of RESCUE that makes a call of CALFUN.
C
   20 IF (KOPT .NE. KBASE) THEN
          IH=0
          DO 30 J=1,N
          DO 30 I=1,J
          IH=IH+1
          IF (I .LT. J) GOPT(J)=GOPT(J)+HQ(IH)*XOPT(I)
   30     GOPT(I)=GOPT(I)+HQ(IH)*XOPT(J)
          IF (NF .GT. NPT) THEN
              DO 50 K=1,NPT
              TEMP=ZERO
              DO 40 J=1,N
   40         TEMP=TEMP+XPT(K,J)*XOPT(J)
              TEMP=PQ(K)*TEMP
              DO 50 I=1,N
   50         GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
          END IF
      END IF
C
C     Generate the next point in the trust region that provides a small value
C     of the quadratic model subject to the constraints on the variables.
C     The integer NTRITS is set to the number "trust region" iterations that
C     have occurred since the last "alternative" iteration. If the length
C     of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to
C     label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW.
C
   60 CALL TRSBOX (N,NPT,XPT,XOPT,GOPT,HQ,PQ,SL,SU,DELTA,XNEW,D,
     1  W,W(NP),W(NP+N),W(NP+2*N),W(NP+3*N),DSQ,CRVMIN)
      DNORM=DMIN1(DELTA,DSQRT(DSQ))
      IF (DNORM .LT. HALF*RHO) THEN
          NTRITS=-1
          DISTSQ=(TEN*RHO)**2
          IF (NF .LE. NFSAV+2) GOTO 650
C
C     The following choice between labels 650 and 680 depends on whether or
C     not our work with the current RHO seems to be complete. Either RHO is
C     decreased or termination occurs if the errors in the quadratic model at
C     the last three interpolation points compare favourably with predictions
C     of likely improvements to the model within distance HALF*RHO of XOPT.
C
          ERRBIG=DMAX1(DIFFA,DIFFB,DIFFC)
          FRHOSQ=0.125D0*RHO*RHO
          IF (CRVMIN .GT. ZERO .AND. ERRBIG .GT. FRHOSQ*CRVMIN)
     1       GOTO 650
          BDTOL=ERRBIG/RHO
          DO 80 J=1,N
          BDTEST=BDTOL
          IF (XNEW(J) .EQ. SL(J)) BDTEST=W(J)
          IF (XNEW(J) .EQ. SU(J)) BDTEST=-W(J)
          IF (BDTEST .LT. BDTOL) THEN
              CURV=HQ((J+J*J)/2)
              DO 70 K=1,NPT
   70         CURV=CURV+PQ(K)*XPT(K,J)**2
              BDTEST=BDTEST+HALF*CURV*RHO
              IF (BDTEST .LT. BDTOL) GOTO 650
          END IF
   80     CONTINUE
          GOTO 680
      END IF
      NTRITS=NTRITS+1
C
C     Severe cancellation is likely to occur if XOPT is too far from XBASE.
C     If the following test holds, then XBASE is shifted so that XOPT becomes
C     zero. The appropriate changes are made to BMAT and to the second
C     derivatives of the current model, beginning with the changes to BMAT
C     that do not depend on ZMAT. VLAG is used temporarily for working space.
C
   90 IF (DSQ .LE. 1.0D-3*XOPTSQ) THEN
          FRACSQ=0.25D0*XOPTSQ
          SUMPQ=ZERO
          DO 110 K=1,NPT
          SUMPQ=SUMPQ+PQ(K)
          SUM=-HALF*XOPTSQ
          DO 100 I=1,N
  100     SUM=SUM+XPT(K,I)*XOPT(I)
          W(NPT+K)=SUM
          TEMP=FRACSQ-HALF*SUM
          DO 110 I=1,N
          W(I)=BMAT(K,I)
          VLAG(I)=SUM*XPT(K,I)+TEMP*XOPT(I)
          IP=NPT+I
          DO 110 J=1,I
  110     BMAT(IP,J)=BMAT(IP,J)+W(I)*VLAG(J)+VLAG(I)*W(J)
C
C     Then the revisions of BMAT that depend on ZMAT are calculated.
C
          DO 150 JJ=1,NPTM
          SUMZ=ZERO
          SUMW=ZERO
          DO 120 K=1,NPT
          SUMZ=SUMZ+ZMAT(K,JJ)
          VLAG(K)=W(NPT+K)*ZMAT(K,JJ)
  120     SUMW=SUMW+VLAG(K)
          DO 140 J=1,N
          SUM=(FRACSQ*SUMZ-HALF*SUMW)*XOPT(J)
          DO 130 K=1,NPT
  130     SUM=SUM+VLAG(K)*XPT(K,J)
          W(J)=SUM
          DO 140 K=1,NPT
  140     BMAT(K,J)=BMAT(K,J)+SUM*ZMAT(K,JJ)
          DO 150 I=1,N
          IP=I+NPT
          TEMP=W(I)
          DO 150 J=1,I
  150     BMAT(IP,J)=BMAT(IP,J)+TEMP*W(J)
C
C     The following instructions complete the shift, including the changes
C     to the second derivative parameters of the quadratic model.
C
          IH=0
          DO 170 J=1,N
          W(J)=-HALF*SUMPQ*XOPT(J)
          DO 160 K=1,NPT
          W(J)=W(J)+PQ(K)*XPT(K,J)
  160     XPT(K,J)=XPT(K,J)-XOPT(J)
          DO 170 I=1,J
          IH=IH+1
          HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J)
  170     BMAT(NPT+I,J)=BMAT(NPT+J,I)
          DO 180 I=1,N
          XBASE(I)=XBASE(I)+XOPT(I)
          XNEW(I)=XNEW(I)-XOPT(I)
          SL(I)=SL(I)-XOPT(I)
          SU(I)=SU(I)-XOPT(I)
  180     XOPT(I)=ZERO
          XOPTSQ=ZERO
      END IF
      IF (NTRITS .EQ. 0) GOTO 210
      GOTO 230
C
C     XBASE is also moved to XOPT by a call of RESCUE. This calculation is
C     more expensive than the previous shift, because new matrices BMAT and
C     ZMAT are generated from scratch, which may include the replacement of
C     interpolation points whose positions seem to be causing near linear
C     dependence in the interpolation conditions. Therefore RESCUE is called
C     only if rounding errors have reduced by at least a factor of two the
C     denominator of the formula for updating the H matrix. It provides a
C     useful safeguard, but is not invoked in most applications of BOBYQA.
C
  190 NFSAV=NF
      KBASE=KOPT
      CALL RESCUE (N,NPT,XL,XU,IPRINT,MAXFUN,XBASE,XPT,FVAL,
     1  XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,DELTA,KOPT,
     2  VLAG,W,W(N+NP),W(NDIM+NP))
C
C     XOPT is updated now in case the branch below to label 720 is taken.
C     Any updating of GOPT occurs after the branch below to label 20, which
C     leads to a trust region iteration as does the branch to label 60.
C
      XOPTSQ=ZERO
      IF (KOPT .NE. KBASE) THEN
          DO 200 I=1,N
          XOPT(I)=XPT(KOPT,I)
  200     XOPTSQ=XOPTSQ+XOPT(I)**2
      END IF
      IF (NF .LT. 0) THEN
          NF=MAXFUN
          IF (IPRINT .GT. 0) PRINT 390
          GOTO 720
      END IF
      NRESC=NF
      IF (NFSAV .LT. NF) THEN
          NFSAV=NF
          GOTO 20
      END IF
      IF (NTRITS .GT. 0) GOTO 60
C
C     Pick two alternative vectors of variables, relative to XBASE, that
C     are suitable as new positions of the KNEW-th interpolation point.
C     Firstly, XNEW is set to the point on a line through XOPT and another
C     interpolation point that minimizes the predicted value of the next
C     denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL
C     and SU bounds. Secondly, XALT is set to the best feasible point on
C     a constrained version of the Cauchy step of the KNEW-th Lagrange
C     function, the corresponding value of the square of this function
C     being returned in CAUCHY. The choice between these alternatives is
C     going to be made when the denominator is calculated.
C
  210 CALL ALTMOV (N,NPT,XPT,XOPT,BMAT,ZMAT,NDIM,SL,SU,KOPT,
     1  KNEW,ADELT,XNEW,XALT,ALPHA,CAUCHY,W,W(NP),W(NDIM+1))
      DO 220 I=1,N
  220 D(I)=XNEW(I)-XOPT(I)
C
C     Calculate VLAG and BETA for the current choice of D. The scalar
C     product of D with XPT(K,.) is going to be held in W(NPT+K) for
C     use when VQUAD is calculated.
C
  230 DO 250 K=1,NPT
      SUMA=ZERO
      SUMB=ZERO
      SUM=ZERO
      DO 240 J=1,N
      SUMA=SUMA+XPT(K,J)*D(J)
      SUMB=SUMB+XPT(K,J)*XOPT(J)
  240 SUM=SUM+BMAT(K,J)*D(J)
      W(K)=SUMA*(HALF*SUMA+SUMB)
      VLAG(K)=SUM
  250 W(NPT+K)=SUMA
      BETA=ZERO
      DO 270 JJ=1,NPTM
      SUM=ZERO
      DO 260 K=1,NPT
  260 SUM=SUM+ZMAT(K,JJ)*W(K)
      BETA=BETA-SUM*SUM
      DO 270 K=1,NPT
  270 VLAG(K)=VLAG(K)+SUM*ZMAT(K,JJ)
      DSQ=ZERO
      BSUM=ZERO
      DX=ZERO
      DO 300 J=1,N
      DSQ=DSQ+D(J)**2
      SUM=ZERO
      DO 280 K=1,NPT
  280 SUM=SUM+W(K)*BMAT(K,J)
      BSUM=BSUM+SUM*D(J)
      JP=NPT+J
      DO 290 I=1,N
  290 SUM=SUM+BMAT(JP,I)*D(I)
      VLAG(JP)=SUM
      BSUM=BSUM+SUM*D(J)
  300 DX=DX+D(J)*XOPT(J)
      BETA=DX*DX+DSQ*(XOPTSQ+DX+DX+HALF*DSQ)+BETA-BSUM
      VLAG(KOPT)=VLAG(KOPT)+ONE
C
C     If NTRITS is zero, the denominator may be increased by replacing
C     the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if
C     rounding errors have damaged the chosen denominator.
C
      IF (NTRITS .EQ. 0) THEN
          DENOM=VLAG(KNEW)**2+ALPHA*BETA
          IF (DENOM .LT. CAUCHY .AND. CAUCHY .GT. ZERO) THEN
              DO 310 I=1,N
              XNEW(I)=XALT(I)
  310         D(I)=XNEW(I)-XOPT(I)
              CAUCHY=ZERO
              GO TO 230
          END IF
          IF (DENOM .LE. HALF*VLAG(KNEW)**2) THEN
              IF (NF .GT. NRESC) GOTO 190
              IF (IPRINT .GT. 0) PRINT 320
  320         FORMAT (/5X,'Return from BOBYQA because of much',
     1          ' cancellation in a denominator.')
              GOTO 720
          END IF
C
C     Alternatively, if NTRITS is positive, then set KNEW to the index of
C     the next interpolation point to be deleted to make room for a trust
C     region step. Again RESCUE may be called if rounding errors have damaged
C     the chosen denominator, which is the reason for attempting to select
C     KNEW before calculating the next value of the objective function.
C
      ELSE
          DELSQ=DELTA*DELTA
          SCADEN=ZERO
          BIGLSQ=ZERO
          KNEW=0
          DO 350 K=1,NPT
          IF (K .EQ. KOPT) GOTO 350
          HDIAG=ZERO
          DO 330 JJ=1,NPTM
  330     HDIAG=HDIAG+ZMAT(K,JJ)**2
          DEN=BETA*HDIAG+VLAG(K)**2
          DISTSQ=ZERO
          DO 340 J=1,N
  340     DISTSQ=DISTSQ+(XPT(K,J)-XOPT(J))**2
          TEMP=DMAX1(ONE,(DISTSQ/DELSQ)**2)
          IF (TEMP*DEN .GT. SCADEN) THEN
              SCADEN=TEMP*DEN
              KNEW=K
              DENOM=DEN
          END IF
          BIGLSQ=DMAX1(BIGLSQ,TEMP*VLAG(K)**2)
  350     CONTINUE
          IF (SCADEN .LE. HALF*BIGLSQ) THEN
              IF (NF .GT. NRESC) GOTO 190
              IF (IPRINT .GT. 0) PRINT 320
              GOTO 720
          END IF
      END IF
C
C     Put the variables for the next calculation of the objective function
C       in XNEW, with any adjustments for the bounds.
C
C
C     Calculate the value of the objective function at XBASE+XNEW, unless
C       the limit on the number of calculations of F has been reached.
C
  360 DO 380 I=1,N
      X(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XNEW(I)),XU(I))
      IF (XNEW(I) .EQ. SL(I)) X(I)=XL(I)
      IF (XNEW(I) .EQ. SU(I)) X(I)=XU(I)
  380 CONTINUE
      IF (NF .GE. MAXFUN) THEN
          IF (IPRINT .GT. 0) PRINT 390
  390     FORMAT (/4X,'Return from BOBYQA because CALFUN has been',
     1      ' called MAXFUN times.')
          GOTO 720
      END IF
      NF=NF+1
      CALL CALFUN (N,X,F)
      IF (IPRINT .EQ. 3) THEN
          PRINT 400, NF,F,(X(I),I=1,N)
  400      FORMAT (/4X,'Function number',I6,'    F =',1PD18.10,
     1       '    The corresponding X is:'/(2X,5D15.6))
      END IF
      IF (NTRITS .EQ. -1) THEN
          FSAVE=F
          GOTO 720
      END IF
C
C     Use the quadratic model to predict the change in F due to the step D,
C       and set DIFF to the error of this prediction.
C
      FOPT=FVAL(KOPT)
      VQUAD=ZERO
      IH=0
      DO 410 J=1,N
      VQUAD=VQUAD+D(J)*GOPT(J)
      DO 410 I=1,J
      IH=IH+1
      TEMP=D(I)*D(J)
      IF (I .EQ. J) TEMP=HALF*TEMP
  410 VQUAD=VQUAD+HQ(IH)*TEMP
      DO 420 K=1,NPT
  420 VQUAD=VQUAD+HALF*PQ(K)*W(NPT+K)**2
      DIFF=F-FOPT-VQUAD
      DIFFC=DIFFB
      DIFFB=DIFFA
      DIFFA=DABS(DIFF)
      IF (DNORM .GT. RHO) NFSAV=NF
C
C     Pick the next value of DELTA after a trust region step.
C
      IF (NTRITS .GT. 0) THEN
          IF (VQUAD .GE. ZERO) THEN
              IF (IPRINT .GT. 0) PRINT 430
  430         FORMAT (/4X,'Return from BOBYQA because a trust',
     1          ' region step has failed to reduce Q.')
              GOTO 720
          END IF
          RATIO=(F-FOPT)/VQUAD
          IF (RATIO .LE. TENTH) THEN
              DELTA=DMIN1(HALF*DELTA,DNORM)
          ELSE IF (RATIO .LE. 0.7D0) THEN
              DELTA=DMAX1(HALF*DELTA,DNORM)
          ELSE
              DELTA=DMAX1(HALF*DELTA,DNORM+DNORM)
          END IF
          IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
C
C     Recalculate KNEW and DENOM if the new F is less than FOPT.
C
          IF (F .LT. FOPT) THEN
              KSAV=KNEW
              DENSAV=DENOM
              DELSQ=DELTA*DELTA
              SCADEN=ZERO
              BIGLSQ=ZERO
              KNEW=0
              DO 460 K=1,NPT
              HDIAG=ZERO
              DO 440 JJ=1,NPTM
  440         HDIAG=HDIAG+ZMAT(K,JJ)**2
              DEN=BETA*HDIAG+VLAG(K)**2
              DISTSQ=ZERO
              DO 450 J=1,N
  450         DISTSQ=DISTSQ+(XPT(K,J)-XNEW(J))**2
              TEMP=DMAX1(ONE,(DISTSQ/DELSQ)**2)
              IF (TEMP*DEN .GT. SCADEN) THEN
                  SCADEN=TEMP*DEN
                  KNEW=K
                  DENOM=DEN
              END IF
  460         BIGLSQ=DMAX1(BIGLSQ,TEMP*VLAG(K)**2)
              IF (SCADEN .LE. HALF*BIGLSQ) THEN
                  KNEW=KSAV
                  DENOM=DENSAV
              END IF
          END IF
      END IF
C
C     Update BMAT and ZMAT, so that the KNEW-th interpolation point can be
C     moved. Also update the second derivative terms of the model.
C
      CALL UPDATE (N,NPT,BMAT,ZMAT,NDIM,VLAG,BETA,DENOM,KNEW,W)
      IH=0
      PQOLD=PQ(KNEW)
      PQ(KNEW)=ZERO
      DO 470 I=1,N
      TEMP=PQOLD*XPT(KNEW,I)
      DO 470 J=1,I
      IH=IH+1
  470 HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J)
      DO 480 JJ=1,NPTM
      TEMP=DIFF*ZMAT(KNEW,JJ)
      DO 480 K=1,NPT
  480 PQ(K)=PQ(K)+TEMP*ZMAT(K,JJ)
C
C     Include the new interpolation point, and make the changes to GOPT at
C     the old XOPT that are caused by the updating of the quadratic model.
C
      FVAL(KNEW)=F
      DO 490 I=1,N
      XPT(KNEW,I)=XNEW(I)
  490 W(I)=BMAT(KNEW,I)
      DO 520 K=1,NPT
      SUMA=ZERO
      DO 500 JJ=1,NPTM
  500 SUMA=SUMA+ZMAT(KNEW,JJ)*ZMAT(K,JJ)
      SUMB=ZERO
      DO 510 J=1,N
  510 SUMB=SUMB+XPT(K,J)*XOPT(J)
      TEMP=SUMA*SUMB
      DO 520 I=1,N
  520 W(I)=W(I)+TEMP*XPT(K,I)
      DO 530 I=1,N
  530 GOPT(I)=GOPT(I)+DIFF*W(I)
C
C     Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT.
C
      IF (F .LT. FOPT) THEN
          KOPT=KNEW
          XOPTSQ=ZERO
          IH=0
          DO 540 J=1,N
          XOPT(J)=XNEW(J)
          XOPTSQ=XOPTSQ+XOPT(J)**2
          DO 540 I=1,J
          IH=IH+1
          IF (I .LT. J) GOPT(J)=GOPT(J)+HQ(IH)*D(I)
  540     GOPT(I)=GOPT(I)+HQ(IH)*D(J)
          DO 560 K=1,NPT
          TEMP=ZERO
          DO 550 J=1,N
  550     TEMP=TEMP+XPT(K,J)*D(J)
          TEMP=PQ(K)*TEMP
          DO 560 I=1,N
  560     GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
      END IF
C
C     Calculate the parameters of the least Frobenius norm interpolant to
C     the current data, the gradient of this interpolant at XOPT being put
C     into VLAG(NPT+I), I=1,2,...,N.
C
      IF (NTRITS .GT. 0) THEN
          DO 570 K=1,NPT
          VLAG(K)=FVAL(K)-FVAL(KOPT)
  570     W(K)=ZERO
          DO 590 J=1,NPTM
          SUM=ZERO
          DO 580 K=1,NPT
  580     SUM=SUM+ZMAT(K,J)*VLAG(K)
          DO 590 K=1,NPT
  590     W(K)=W(K)+SUM*ZMAT(K,J)
          DO 610 K=1,NPT
          SUM=ZERO
          DO 600 J=1,N
  600     SUM=SUM+XPT(K,J)*XOPT(J)
          W(K+NPT)=W(K)
  610     W(K)=SUM*W(K)
          GQSQ=ZERO
          GISQ=ZERO
          DO 630 I=1,N
          SUM=ZERO
          DO 620 K=1,NPT
  620     SUM=SUM+BMAT(K,I)*VLAG(K)+XPT(K,I)*W(K)
          IF (XOPT(I) .EQ. SL(I)) THEN
              GQSQ=GQSQ+DMIN1(ZERO,GOPT(I))**2
              GISQ=GISQ+DMIN1(ZERO,SUM)**2
          ELSE IF (XOPT(I) .EQ. SU(I)) THEN
              GQSQ=GQSQ+DMAX1(ZERO,GOPT(I))**2
              GISQ=GISQ+DMAX1(ZERO,SUM)**2
          ELSE
              GQSQ=GQSQ+GOPT(I)**2
              GISQ=GISQ+SUM*SUM
          END IF
  630     VLAG(NPT+I)=SUM
C
C     Test whether to replace the new quadratic model by the least Frobenius
C     norm interpolant, making the replacement if the test is satisfied.
C
          ITEST=ITEST+1
          IF (GQSQ .LT. TEN*GISQ) ITEST=0
          IF (ITEST .GE. 3) THEN
              DO 640 I=1,MAX0(NPT,NH)
              IF (I .LE. N) GOPT(I)=VLAG(NPT+I)
              IF (I .LE. NPT) PQ(I)=W(NPT+I)
              IF (I .LE. NH) HQ(I)=ZERO
              ITEST=0
  640         CONTINUE
          END IF
      END IF
C
C     If a trust region step has provided a sufficient decrease in F, then
C     branch for another trust region calculation. The case NTRITS=0 occurs
C     when the new interpolation point was reached by an alternative step.
C
      IF (NTRITS .EQ. 0) GOTO 60
      IF (F .LE. FOPT+TENTH*VQUAD) GOTO 60
C
C     Alternatively, find out if the interpolation points are close enough
C       to the best point so far.
C
      DISTSQ=DMAX1((TWO*DELTA)**2,(TEN*RHO)**2)
  650 KNEW=0
      DO 670 K=1,NPT
      SUM=ZERO
      DO 660 J=1,N
  660 SUM=SUM+(XPT(K,J)-XOPT(J))**2
      IF (SUM .GT. DISTSQ) THEN
          KNEW=K
          DISTSQ=SUM
      END IF
  670 CONTINUE
C
C     If KNEW is positive, then ALTMOV finds alternative new positions for
C     the KNEW-th interpolation point within distance ADELT of XOPT. It is
C     reached via label 90. Otherwise, there is a branch to label 60 for
C     another trust region iteration, unless the calculations with the
C     current RHO are complete.
C
      IF (KNEW .GT. 0) THEN
          DIST=DSQRT(DISTSQ)
          IF (NTRITS .EQ. -1) THEN
              DELTA=DMIN1(TENTH*DELTA,HALF*DIST)
              IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
          END IF
          NTRITS=0
          ADELT=DMAX1(DMIN1(TENTH*DIST,DELTA),RHO)
          DSQ=ADELT*ADELT
          GOTO 90
      END IF
      IF (NTRITS .EQ. -1) GOTO 680
      IF (RATIO .GT. ZERO) GOTO 60
      IF (DMAX1(DELTA,DNORM) .GT. RHO) GOTO 60
C
C     The calculations with the current value of RHO are complete. Pick the
C       next values of RHO and DELTA.
C
  680 IF (RHO .GT. RHOEND) THEN
          DELTA=HALF*RHO
          RATIO=RHO/RHOEND
          IF (RATIO .LE. 16.0D0) THEN
              RHO=RHOEND
          ELSE IF (RATIO .LE. 250.0D0) THEN
              RHO=DSQRT(RATIO)*RHOEND
          ELSE
              RHO=TENTH*RHO
          END IF
          DELTA=DMAX1(DELTA,RHO)
          IF (IPRINT .GE. 2) THEN
              IF (IPRINT .GE. 3) PRINT 690
  690         FORMAT (5X)
              PRINT 700, RHO,NF
  700         FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of',
     1          ' function values =',I6)
              PRINT 710, FVAL(KOPT),(XBASE(I)+XOPT(I),I=1,N)
  710         FORMAT (4X,'Least value of F =',1PD23.15,9X,
     1          'The corresponding X is:'/(2X,5D15.6))
          END IF
          NTRITS=0
          NFSAV=NF
          GOTO 60
      END IF
C
C     Return from the calculation, after another Newton-Raphson step, if
C       it is too short to have been tried before.
C
      IF (NTRITS .EQ. -1) GOTO 360
  720 IF (FVAL(KOPT) .LE. FSAVE) THEN
          DO 730 I=1,N
          X(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XOPT(I)),XU(I))
          IF (XOPT(I) .EQ. SL(I)) X(I)=XL(I)
          IF (XOPT(I) .EQ. SU(I)) X(I)=XU(I)
  730     CONTINUE
          F=FVAL(KOPT)
      END IF
      IF (IPRINT .GE. 1) THEN
          PRINT 740, NF
  740     FORMAT (/4X,'At the return from BOBYQA',5X,
     1      'Number of function values =',I6)
          PRINT 710, F,(X(I),I=1,N)
      END IF
      RETURN
      END


text:
! 8.24 в учебнике

program ex_8_24
    implicit none
    integer, parameter      :: R_ = 4
    integer                 :: In = 0, Out = 0, N = 0, i = 0
    real(R_)                :: x, a, b, h, sum, integral, t
    character(*), parameter :: output_file = "output.txt", &
                               input_file = "../data/input.txt"

    open (file=input_file, newunit=In)
        read(In, *) a, b
        read(In, *) N
        read(In, *) x
    close (In)

    if (mod(N,2) == 1 .or. N <= 0) then
        write(*,*) 'Число N должно быть положительным и четным. Отмена...'
        stop
    endif

    h = (b - a) / N
    t = a
    sum = F(t)

    i = 1
    do
        t = t + h
        sum = sum + 4 * F(t)

        i = i + 2

        if (i >= N) exit

        t = t + h
        sum = sum + 2 * F(t)
    enddo

    sum = sum + F(b)
    integral = h/3 * sum

    open (file=output_file, newunit=Out)
        write(Out,*) 'Ответ ', integral
    close (Out)

contains

real function F(x)
    implicit none
    real, intent(in) :: x

    F = (q(x) ** 2 + 1.57) / sqrt(q(x) + 2.15)

end function F

pure real function q(x)
    implicit none
    real, intent(in) :: x
    real             :: sum, old_E, new_E
    integer          :: i, j

    sum   = 1
    new_E = 1

    i   = 1
    do
        old_E = new_E
        new_E = old_E * x**i

        do j = i*2+1, 2, -1
            new_E = new_E / j
        enddo

        sum = sum + new_E
        if (new_E == old_E) exit
    enddo

    q = sum

end function q

end program ex_8_24

max_stars_repo_name: taxnuke/algorithm-exercises

max_issues_repo_name: taxnuke/algorithm-exercises

max_forks_repo_name: taxnuke/algorithm-exercises


text:
subroutine setup_stressperm
!***********************************************************************
! Copyright 2011 Los Alamos National Security, LLC  All rights reserved
! Unless otherwise indicated,  this information has been authored by an
! employee or employees of the Los Alamos National Security, LLC (LANS),
! operator of the  Los  Alamos National  Laboratory  under Contract  No.
! DE-AC52-06NA25396  with  the U. S. Department  of  Energy.  The  U. S.
! Government   has   rights  to  use,  reproduce,  and  distribute  this
! information.  The  public may copy  and  use this  information without
! charge, provided that this  Notice and any statement of authorship are
! reproduced on all copies.  Neither  the  Government nor LANS makes any
! warranty,   express   or   implied,   or   assumes  any  liability  or
! responsibility for the use of this information.
!***********************************************************************

  use comai
  use combi
  use comdi
  use comdti
  use comsi
  implicit none
  integer jj1,i,iispmd,i1,i2,jj, kb
  integer kbxmin,kbxmax,kbymin,kbymax,kbzmin,kbzmax
  real*8 xdmin,xdmax,ydmin,ydmax,zdmin,zdmax,xi,yi,zi
  real*8 xdmint,xdmaxt,ydmint,ydmaxt,zdmint,zdmaxt
  real*8 xkb,ykb,zkb,xd1,yd1,zd1,dis,dist,disx,disy,disz
  real*8 dis_tol

  parameter (dis_tol = 1.d-8)

  jj1 =0
  do i = 1,n0
     iispmd = ispm(i)
     ispmd = ispmt(iispmd)
     if(ispmd.eq.3.or.ispmd.eq.5.or.ispmd.eq.7 &
          .or.ispmd.eq.8) then
        xi = cord(i,1)
        yi = cord(i,2)
        zi = cord(i,3)
        xdmax = -1.d30
        xdmin = -1.d30
        ydmax = -1.d30
        ydmin = -1.d30
        zdmax = -1.d30
        zdmin = -1.d30
        xdmaxt = -1.d30
        xdmint = 1.d30
        ydmaxt = -1.d30
        ydmint = 1.d30
        zdmaxt = -1.d30
        zdmint = 1.d30
        kbxmax = i
        kbxmin = i
        kbymax = i
        kbymin = i
        kbzmax = i
        kbzmin = i
        i1 = nelm(i)+1
        i2 = nelm(i+1)
        do jj = i1,i2
           kb = nelm(jj)
           xkb = cord(kb,1)
           ykb = cord(kb,2)
           zkb = cord(kb,3)
           if(xkb-xi.gt.0) then
              dis = abs(xkb-xi)
              dist = abs(xkb-xi) + abs(ykb-yi) + abs(zkb-zi)
              if(dis.gt.xdmax) then
                 kbxmax = kb
                 xdmax = dis
                 xdmaxt = dist
              else if(abs(dis-xdmax).lt.dis_tol) then
                 if(dist.lt.xdmaxt) then
                    xdmaxt = dist
                    kbxmax = kb
                 endif
              endif
           else if(xkb-xi.lt.0) then
              dist = abs(xkb-xi) + abs(ykb-yi) + abs(zkb-zi)
              dis = abs(xkb-xi)
              if(dis.gt.xdmin) then
                 kbxmin = kb
                 xdmin = dis
                 xdmint = dist
              else if(abs(dis-xdmin).lt.dis_tol) then
                 if(dist.lt.xdmint) then
                    xdmint = dist
                    kbxmin = kb
                 endif
              endif
           endif
           if(ykb-yi.gt.0) then
              dist = abs(xkb-xi) + abs(ykb-yi) + abs(zkb-zi)
              dis = abs(ykb-yi)
              if(dis.gt.ydmax) then
                 kbymax = kb
                 ydmax = dis
                 ydmaxt = dist
              else if(abs(dis-ydmax).lt.dis_tol) then
                 if(dist.lt.ydmaxt) then
                    ydmaxt = dist
                    kbymax = kb
                 endif
              endif
           else if(ykb-yi.lt.0) then
              dist = abs(xkb-xi) + abs(ykb-yi) + abs(zkb-zi)
              dis = abs(ykb-yi)
              if(dis.gt.ydmin) then
                 kbymin = kb
                 ydmin = dis
                 ydmint = dist
              else if(abs(dis-ydmin).lt.dis_tol) then
                 if(dist.lt.ydmint) then
                    ydmint = dist
                    kbymin = kb
                 endif
              endif
           endif
           if(zkb-zi.gt.0) then
              dist = abs(xkb-xi) + abs(ykb-yi) + abs(zkb-zi)
              dis = abs(zkb-zi)
              if(dis.gt.zdmax) then
                 kbzmax = kb
                 zdmax = dis
                 zdmaxt = dist
              else if(abs(dis-zdmax).lt.dis_tol) then
                 if(dist.lt.zdmaxt) then
                    zdmaxt = dist
                    kbzmax = kb
                 endif
              endif
           else if(zkb-zi.lt.0) then
              dist = abs(xkb-xi) + abs(ykb-yi) + abs(zkb-zi)
              dis = abs(zkb-zi)
              if(dis.gt.zdmin) then
                 kbzmin = kb
                 zdmin = dis
                 zdmint = dist
              else if(abs(dis-zdmin).lt.dis_tol) then
                 if(dist.lt.zdmint) then
                    zdmint = dist
                    kbzmin = kb
                 endif
              endif
           endif
        enddo
        !
        !     store max and min of each direction (may create (slightly)larger stencil)
        !
        ipermx(i,1) = kbxmin
        ipermx(i,2) = kbxmax
        ipermy(i,1) = kbymin
        ipermy(i,2) = kbymax
        ipermz(i,1) = kbzmin
        ipermz(i,2) = kbzmax
        if(cord(kbxmax,1)-cord(kbxmin,1).le.0.0) then
           jj1 =1
           write(ierr,*) &
                'dis(x) failed, node ',i,' model 3 or 5 sub stres_perm'
        endif
        if(cord(kbymax,2)-cord(kbymin,2).le.0.0) then
           jj1 =1
           write(ierr,*) &
                'dis(y) failed, node ',i,' model 3 or 5 sub stres_perm'
        endif
        if(cord(kbzmax,3)-cord(kbzmin,3).le.0.0) then
           jj1 =1
           write(ierr,*) &
                'dis(z) failed, node ',i,' model 3 or 5 sub stres_perm'
        endif
        !     stop for zero distances
        if(jj1.ne.0) stop
     endif
  enddo
  go to 2001
  !
  !     test code
  !
  open(99,file='stress_perm',status = 'unknown')
  do i = 1,n0
     iispmd = ispm(i)
     ispmd = ispmt(iispmd)
     if(ispmd.eq.3.or.ispmd.eq.5) then
        write(99,*) 'node = ',i
        write(99,*)  &
             ' x node pair', ipermx(i,2), ipermx(i,1),  &
             ' dis ', cord(ipermx(i,2),1)-cord(ipermx(i,1),1)
        if(cord(ipermx(i,2),1)-cord(ipermx(i,1),1).le.0.001)  &
             write (99,*) '>>> zero distance x <<<<'
        write(99,*)  &
             ' y node pair', ipermy(i,2), ipermy(i,1),  &
             ' dis ', cord(ipermy(i,2),2)-cord(ipermy(i,1),2)
        if(cord(ipermy(i,2),2)-cord(ipermy(i,1),2).le.0.001)  &
             write (99,*) '>>> zero distance y <<<<'
        write(99,*)  &
             ' z node pair', ipermz(i,2), ipermz(i,1),  &
             ' dis ', cord(ipermz(i,2),3)-cord(ipermz(i,1),3)
        if(cord(ipermz(i,2),3)-cord(ipermz(i,1),3).le.0.001)  &
             write (99,*) '>>> zero distance z <<<<'
     endif
  enddo
  close(99)

2001 continue

  return

end subroutine setup_stressperm
!......................................................................

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
+incdir+$SVUNIT_INSTALL/examples/modules/clk_and_reset/someLibDir
+incdir+$SVUNIT_INSTALL/svunit_base/util

$SVUNIT_INSTALL/examples/modules/clk_and_reset/my_and_module.v

max_stars_repo_name: umarcor/svunit

max_issues_repo_name: umarcor/svunit

max_forks_repo_name: umarcor/svunit