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

! Fortran implementation of ACORN random number generator of order less
! than or equal to 12 (higher orders can be obtained by increasing the
! parameter value MAXORD).

! NOTES: 1. The variable idum is a dummy variable. The common block
!           IACO is used to transfer data into the function.

!        2. Before the first call to ACORN the common block IACO must
!           be initialised by the user, as follows. The values of
!           variables in the common block must not subsequently be
!           changed by the user.

!             KORDEI - order of generator required ( must be =< MAXORD)

!             MAXINT - modulus for generator, must be chosen small
!                      enough that 2*MAXINT does not overflow

!             ixv(1) - seed for random number generator
!                      require 0 < ixv(1) < MAXINT

!             (ixv(I+1),I=1,KORDEI)
!                    - KORDEI initial values for generator
!                      require 0 =< ixv(I+1) < MAXINT

!        3. After initialisation, each call to ACORN generates a single
!           random number between 0 and 1.

!        4. An example of suitable values for parameters is

!             KORDEI   = 10
!             MAXINT   = 2**30
!             ixv(1)   = an odd integer in the (approximate) range
!                        (0.001 * MAXINT) to (0.999 * MAXINT)
!             ixv(I+1) = 0, I=1,KORDEI

! Author: R.S.Wikramaratna,                           Date: October 1990
!-----------------------------------------------------------------------
double precision function acorni(idum)
    implicit double precision (a-h,o-z)
    parameter (KORDEI=12,MAXOP1=KORDEI+1,MAXINT=2**30)
    common/iaco/ ixv(MAXOP1)
    do i=1,KORDEI
        ixv(i+1)=(ixv(i+1)+ixv(i))
        if(ixv(i+1) >= MAXINT) ixv(i+1)=ixv(i+1)-MAXINT
    end do
    acorni=dble(ixv(KORDEI+1))/MAXINT
    return
end function


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

!>    @brief (weighted) jacoby-matrix output for each physical value
!>    @param[in] ismpl local sample index
      SUBROUTINE write_jacw(ismpl)
        use arrays
#ifndef JACOBI_FREE
#ifndef AD_RM
        use g_arrays
#else
        use arrays_ad
#endif
        use mod_genrl
        use mod_genrlc
        use mod_inverse
        use mod_data
        use mod_linfos
        IMPLICIT NONE
        integer :: ismpl
        integer :: i, j, k, l

        DOUBLE PRECISION dx, dy, dz, errd
        DOUBLE PRECISION, ALLOCATABLE :: jac_row(:)
        INTEGER i1, i2, i3, ii, kk, type, ll_pv, lout
        INTEGER i1_, i2_, i3_, i4_

        CHARACTER filename*80, snumber*8
        INTEGER lblank
        EXTERNAL lblank

        IF (write_disable) RETURN
#ifdef NOWJAC
        RETURN
#endif

        IF (mpara>=9999) THEN
!       needs modifying on the write-format
          WRITE(*,*) 'error : in "write_jacw" need "mpara"<9999 !'
          STOP
        END IF

        CALL omp_new_file_handler(lout,16)

        ALLOCATE(jac_row(mpara))
        CALL chln(project,i1_,i2_)
        WRITE(snumber,'(1I7)') iter_inv
        CALL chln(snumber,i3_,i4_)

        DO ll_pv = 1, npv
          kk = 0
          IF (ll_pv==pv_head) kk = ndata_h
          IF (ll_pv==pv_temp) kk = ndata_t
          IF (ll_pv==pv_conc) kk = ndata_c
          IF (ll_pv==pv_pres) kk = ndata_p

          IF (kk>0) THEN
            filename = project(i1_:i2_) // '_JW' // pv_name(ll_pv) // &
              '_' // snumber(i3_:i4_) // '.plt'
            OPEN(lout,file=filename,status='unknown',blank='null')

            IF (linfos(2)>=0) THEN
              WRITE(*,'(5A)') '  [W] : ', pv_name(ll_pv), &
                '-JacW to "', filename(1:lblank(filename)), '"'
            END IF

! 'x y z uid pv-value por1 akx1 aky1 kz1 alx1 aly1 lz1 q1 por2 akx2 aky2 kz2 alx2 aly2 lz2 q2'
            CALL head_tecpl(lout,kk,1,1,pv_name(ll_pv),ismpl)

            DO l = 1, ndata
              type = idata(l,cid_pv)
              IF (type==ll_pv) THEN
                i = idata(l,cid_i)
                j = idata(l,cid_j)
                k = idata(l,cid_k)
                errd = 1.D0/ddata(l,cdd_w)
                dx = delxa(i)
                dy = delya(j)
                dz = delza(k)
                DO ii = 1, mpara
                  jac_row(ii) = jac(l,ii)
                END DO
                WRITE(lout,'(3(e16.7),1I8,9999(e16.7,1X))') dx, dy, &
                  dz, uindex(i,j,k), sdata(l,ismpl), &
                  (errd*jac_row(ii),ii=1,mpara)
              END IF
            END DO
            CLOSE(lout)
          END IF
        END DO
        DEALLOCATE(jac_row)

        CALL omp_del_file_handler(lout)
#endif
        RETURN
      END


    Program ztzrzf_example

!     ZTZRZF Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use blas_interfaces, Only: dznrm2, ztrsm
      Use lapack_example_aux, Only: nagf_file_print_matrix_complex_gen_comp
      Use lapack_interfaces, Only: zgeqp3, ztzrzf, zunmqr, zunmrz
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Complex (Kind=dp), Parameter :: one = (1.0_dp, 0.0_dp)
      Complex (Kind=dp), Parameter :: zero = (0.0_dp, 0.0_dp)
      Integer, Parameter :: inc1 = 1, nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: tol
      Integer :: i, ifail, info, j, k, lda, ldb, lwork, m, n, nrhs
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:, :), tau(:), work(:)
      Real (Kind=dp), Allocatable :: rnorm(:), rwork(:)
      Integer, Allocatable :: jpvt(:)
      Character (1) :: clabs(1), rlabs(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs
!     .. Executable Statements ..
      Write (nout, *) 'ZTZRZF Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n, nrhs
      lda = m
      ldb = m
      lwork = (n+1)*nb
      Allocate (a(lda,n), b(ldb,nrhs), tau(n), work(lwork), rnorm(n), &
        rwork(2*n), jpvt(n))

!     Read A and B from data file

      Read (nin, *)(a(i,1:n), i=1, m)
      Read (nin, *)(b(i,1:nrhs), i=1, m)

!     Initialize JPVT to be zero so that all columns are free

      jpvt(1:n) = 0

!     Compute the QR factorization of A with column pivoting as
!     A = Q*(R11 R12)*(P**T)
!            ( 0  R22)
      Call zgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)

!     Compute C = (C1) = (Q**H)*B, storing the result in B
!                  (C2)
      Call zunmqr('Left', 'Conjugate transpose', m, nrhs, n, a, lda, tau, b, &
        ldb, work, lwork, info)

!     Choose TOL to reflect the relative accuracy of the input data

      tol = 0.01_dp

!     Determine and print the rank, K, of R relative to TOL

loop: Do k = 1, n
        If (abs(a(k,k))<=tol*abs(a(1,1))) Then
          Exit loop
        End If
      End Do loop
      k = k - 1

      Write (nout, *) 'Tolerance used to estimate the rank of A'
      Write (nout, 100) tol
      Write (nout, *) 'Estimated rank of A'
      Write (nout, 110) k
      Write (nout, *)
      Flush (nout)

!     Compute the RZ factorization of the K by K part of R as
!     (R1 R2) = (T 0)*Z
      Call ztzrzf(k, n, a, lda, tau, work, lwork, info)

!     Compute least squares solutions of triangular problems by
!     back substitution in T*Y1 = C1, storing the result in B
      Call ztrsm('Left', 'Upper', 'No transpose', 'Non-Unit', k, nrhs, one, a, &
        lda, b, ldb)

!     Compute estimates of the square roots of the residual sums of
!     squares (2-norm of each of the columns of C2)

      Do j = 1, nrhs
        rnorm(j) = dznrm2(m-k, b(k+1,j), inc1)
      End Do

!     Set the remaining elements of the solutions to zero (to give
!     the minimum-norm solutions), Y2 = 0

      b(k+1:n, 1:nrhs) = zero

!     Form W = (Z**H)*Y
      Call zunmrz('Left', 'Conjugate transpose', n, nrhs, k, n-k, a, lda, tau, &
        b, ldb, work, lwork, info)

!     Permute the least squares solutions stored in B to give X = P*W

      Do j = 1, nrhs
        Do i = 1, n
          work(jpvt(i)) = b(i, j)
        End Do
        b(1:n, j) = work(1:n)
      End Do

!     Print least squares solutions

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call nagf_file_print_matrix_complex_gen_comp('General', ' ', n, nrhs, b, &
        ldb, 'Bracketed', 'F7.4', 'Least squares solution(s)', 'Integer', &
        rlabs, 'Integer', clabs, 80, 0, ifail)

!     Print the square roots of the residual sums of squares

      Write (nout, *)
      Write (nout, *) 'Square root(s) of the residual sum(s) of squares'
      Write (nout, 100) rnorm(1:nrhs)

100   Format (3X, 1P, 7E11.2)
110   Format (1X, I6)
    End Program


text:
      subroutine kscomp (work)
      implicit double precision (a-h,o-z)
      common /kscomm/ rdf   ,adf   ,fdl   ,fdm   ,rho   ,drho  ,rhomax,
     1                fun0  ,slope ,delx  ,alpha ,alpmax,a1    ,a2    ,
     2                a3    ,a4    ,f1    ,f2    ,f3    ,f4    ,alim  ,
     3                atest ,ftest ,ifscl ,ifoff ,isx   ,isx0  ,isxlb ,
     4                isxub ,iscl  ,ig0   ,idf   ,islp  ,iobj0 ,iy    ,
     5                ip    ,ih    ,ihess ,iside ,isact ,idobj ,idg   ,
     6                itmp1 ,itmp2 ,inext ,jnext ,jsel  ,itcnt ,icntr ,
     7                icnta ,isdflg,isdrst,ifncl ,nunit ,ndv   ,ncon  ,
     8                nobj  ,nside ,nscale,iprnt ,itmax ,igrad ,limit
      dimension work(*)
c
c          copy variables from common block kscomm to work array
c
c          author   - Gregory A. Wrenn
c          location - Lockheed Engineering and Sciences Co.
c                     144 Research Drive
c                     Hampton, Va. 23666
c
c          last modification - 19 July 1996
c
      work(1)  = rdf
      work(2)  = adf
      work(3)  = fdl
      work(4)  = fdm
      work(5)  = rho
      work(6)  = drho
      work(7)  = rhomax
      work(8)  = fun0
      work(9)  = slope
      work(10) = delx
      work(11) = alpha
      work(12) = alpmax
      work(13) = a1
      work(14) = a2
      work(15) = a3
      work(16) = a4
      work(17) = f1
      work(18) = f2
      work(19) = f3
      work(20) = f4
      work(21) = alim
      work(22) = atest
      work(23) = ftest
      work(24) = float(ifscl)
      work(25) = float(ifoff)
      work(26) = float(isx)
      work(27) = float(isx0)
      work(28) = float(isxlb)
      work(29) = float(isxub)
      work(30) = float(iscl)
      work(31) = float(ig0)
      work(32) = float(idf)
      work(33) = float(islp)
      work(34) = float(iobj0)
      work(35) = float(iy)
      work(36) = float(ip)
      work(37) = float(ih)
      work(38) = float(ihess)
      work(39) = float(iside)
      work(40) = float(isact)
      work(41) = float(idobj)
      work(42) = float(idg)
      work(43) = float(itmp1)
      work(44) = float(itmp2)
      work(45) = float(inext)
      work(46) = float(jnext)
      work(47) = float(jsel)
      work(48) = float(itcnt)
      work(49) = float(icntr)
      work(50) = float(icnta)
      work(51) = float(isdflg)
      work(52) = float(isdrst)
      work(53) = float(ifncl)
      work(54) = float(nunit)
      work(55) = float(ndv)
      work(56) = float(ncon)
      work(57) = float(nobj)
      work(58) = float(nside)
      work(59) = float(nscale)
      work(60) = float(iprnt)
      work(61) = float(itmax)
      work(62) = float(igrad)
      work(63) = float(limit)
c
      return
      end

ext: f

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


        SUBROUTINE modelCumulTTFrom(calendarMomentst1,calendarCumulst1,cumulTT,cumulTTFromZC65, cumulTTFromZC39, cumulTTFromZC91)
                        IMPLICIT NONE
                        CHARACTER(65), ALLOCATABLE , DIMENSION(:), INTENT(IN) ::
                                calendarMomentst1
                        REAL, ALLOCATABLE , DIMENSION(:), INTENT(IN) :: calendarCumulst1
                        REAL, INTENT(IN) :: cumulTT
                        REAL, INTENT(OUT) :: cumulTTFromZC65
                        REAL, INTENT(OUT) :: cumulTTFromZC39
                        REAL, INTENT(OUT) :: cumulTTFromZC91
                        cumulTTFromZC65 = 0.0
                        cumulTTFromZC39 = 0.0
                        cumulTTFromZC91 = 0.0
                        IF(ANY(calendarMomentst1 .EQ. 'Anthesis')) THEN
                            cumulTTFromZC65 = cumulTT -
                            calendarCumulst1(indice(calendarMomentst1, 'Anthesis'))
                        END IF
                        IF(ANY(calendarMomentst1 .EQ. 'FlagLeafLiguleJustVisible')) THEN
                            cumulTTFromZC39 = cumulTT -
                                    calendarCumulst1(indice(calendarMomentst1,
                                    'FlagLeafLiguleJustVisible'))
                        END IF
                        IF(ANY(calendarMomentst1 .EQ. 'EndGrainFilling')) THEN
                            cumulTTFromZC91 = cumulTT -
                                    calendarCumulst1(indice(calendarMomentst1, 'EndGrainFilling'))
                        END IF
                    END SUBROUTINE modelCumulTTFrom

        END MODULE


text:
      double complex function Fcc_qpqmgpgm(j1,j2,j3,j4,j5,j6,za,zb)
      implicit none
      integer j1,j2,j3,j4,j5,j6
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      double complex L0,Lsm1,Lsm1_2mh,Lnrat,I3m
      double precision t

      Fcc_qpqmgpgm=
     .Lsm1(-s(j1,j2),-t(j1,j2,j4),-s(j1,j4),-t(j1,j2,j4))*
     .(-((zb(j1,j2)**2*(-(za(j1,j5)*zb(j1,j4))-za(j2,j5)*zb(j2,j4))**2)/
     .(za(j5,j6)*zb(j1,j4)*zb(j2,j4)**2*
     .(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))))+
     .(zb(j1,j4)*(-(za(j1,j5)*zb(j1,j2))+za(j4,j5)*zb(j2,j4))**2)/
     .(za(j5,j6)*zb(j2,j4)**2*(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4
     .))*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))))-
     .(2d0*L0(-t(j2,j3,j4),-s(j5,j6))*za(j1,j5)*zb(j1,j3)*
     .(za(j1,j2)*zb(j2,j3)-za(j1,j4)*zb(j3,j4))*
     .(-(za(j2,j5)*zb(j2,j3))+za(j4,j5)*zb(j3,j4)))/
     .(s(j5,j6)*za(j5,j6)*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))*
     .zb(j3,j4)*(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)))+
     .Lsm1(-s(j1,j2),-t(j1,j2,j3),-s(j2,j3),-t(j1,j2,j3))*
     .((za(j1,j2)**2*(-(za(j1,j3)*zb(j1,j6))-za(j2,j3)*zb(j2,j6))**2)/
     .(za(j1,j3)**2*za(j2,j3)*(-(za(j1,j3)*zb(j1,j4))-
     .za(j2,j3)*zb(j2,j4))*(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))*
     .zb(j5,j6))-(za(j2,j3)*
     .(za(j1,j2)*zb(j2,j6)+za(j1,j3)*zb(j3,j6))**2)/
     .(za(j1,j3)**2*(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))*zb(j5,j6)))-
     .(2d0*L0(-t(j1,j3,j4),-s(j5,j6))*za(j2,j4)*
     .(-(za(j1,j4)*zb(j1,j2))+za(j3,j4)*zb(j2,j3))*zb(j2,j6)*
     .(-(za(j1,j4)*zb(j1,j6))-za(j3,j4)*zb(j3,j6)))/
     .(s(j5,j6)*za(j3,j4)*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*zb(j5,j6))

      Fcc_qpqmgpgm=Fcc_qpqmgpgm+
     .I3m(s(j1,j4),s(j2,j3),s(j5,j6))*
     .((za(j1,j5)*za(j2,j4)*zb(j1,j3)*zb(j2,j6)+
     .((-(za(j1,j4)*zb(j1,j3))-za(j2,j4)*zb(j2,j3))*
     .(za(j1,j3)*za(j2,j5)*zb(j1,j6)*zb(j2,j4)+
     .za(j1,j3)*zb(j1,j4)*
     .(-(za(j2,j5)*zb(j2,j6))+za(j3,j5)*zb(j3,j6))))/
     .(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))+
     .((-(za(j1,j4)*zb(j1,j3))+za(j2,j4)*zb(j2,j3))*
     .(-(za(j1,j5)*zb(j1,j6))-za(j4,j5)*zb(j4,j6)))/2)/
     .((-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)))-
     .(za(j4,j5)*zb(j1,j3)*(-(za(j1,j2)*zb(j1,j6))+za(j2,j3)*zb(j3,j6)))
     ./
     .((-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*t(j1,j2,j3)))+
     .Lsm1_2mh(s(j1,j4),t(j1,j2,j3),s(j2,j3),s(j5,j6))*
     .(-(((-(za(j1,j2)*zb(j1,j4))+za(j2,j3)*zb(j3,j4))**2*
     .(za(j1,j2)*zb(j2,j6)+za(j1,j3)*zb(j3,j6))**2)/
     .(za(j2,j3)*(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))**3*zb(j5,j6)))+
     .(za(j1,j2)**2*zb(j4,j6)**2*t(j1,j2,j3)**2)/
     .(za(j2,j3)*(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))**3*zb(j5,j6)))+
     .Lsm1_2mh(s(j3,j4),t(j1,j2,j3),s(j1,j2),s(j5,j6))*
     .(((-(za(j1,j3)*zb(j1,j6))-za(j2,j3)*zb(j2,j6))**2*
     .(-(za(j1,j2)*zb(j1,j4))+za(j2,j3)*zb(j3,j4))**2)/
     .(za(j2,j3)*(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))**3*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))*zb(j5,j6))-
     .(za(j2,j3)*zb(j4,j6)**2*t(j1,j2,j3)**2)/
     .((-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))**3*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))*zb(j5,j6)))+
     .(2d0*Lnrat(-s(j2,j3),-s(j5,j6))*(-2d0*(s(j1,j2)-s(j3,j4))*za(j2,j4
     .)*zb(j2,j3)*
     .(-(za(j2,j5)*zb(j2,j6))-za(j3,j5)*zb(j3,j6))+
     .(za(j2,j5)*zb(j1,j3)*((s(j1,j2)-s(j3,j4))*
     .((-s(j1,j4)+s(j2,j3)-s(j5,j6))*za(j3,j5)*zb(j2,j3)+
     .(-s(j1,j4)-s(j2,j3)+s(j5,j6))*za(j5,j6)*zb(j2,j6))+
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*
     .(-((-s(j1,j4)+s(j2,j3)-s(j5,j6))*za(j2,j5)*zb(j2,j3))+
     .(-s(j1,j4)-s(j2,j3)+s(j5,j6))*za(j5,j6)*zb(j3,j6))))/
     .(za(j5,j6)*zb(j1,j4))-
     .(za(j2,j4)*(zb(j1,j3)*
     .(za(j1,j2)*zb(j2,j6)+za(j1,j3)*zb(j3,j6))-
     .zb(j3,j4)*(-(za(j2,j4)*zb(j2,j6))-za(j3,j4)*zb(j3,j6))-
     .zb(j2,j3)*(-(za(j1,j2)*zb(j1,j6))+za(j2,j4)*zb(j4,j6)))*
     .(za(j3,j5)*zb(j2,j3)*zb(j5,j6)+zb(j2,j6)*t(j1,j2,j3)))/zb(j5,j6)
     .))/((s(j1,j4)**2-2d0*s(j1,j4)*s(j2,j3)+s(j2,j3)**2-2d0*s(j1,j4)*s(
     .j5,j6)-
     .2d0*s(j2,j3)*s(j5,j6)+s(j5,j6)**2)*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)))

      Fcc_qpqmgpgm=Fcc_qpqmgpgm+
     .I3m(s(j2,j3),s(j1,j4),s(j5,j6))*
     .((za(j1,j5)*za(j2,j4)*zb(j1,j3)*zb(j2,j6)+
     .((za(j1,j4)*zb(j1,j3)-za(j2,j4)*zb(j2,j3))*
     .(-(za(j2,j5)*zb(j2,j6))-za(j3,j5)*zb(j3,j6)))/2+
     .((-(za(j1,j4)*zb(j1,j3))-za(j2,j4)*zb(j2,j3))*
     .(za(j1,j3)*za(j2,j5)*zb(j1,j6)*zb(j2,j4)+
     .za(j2,j3)*zb(j2,j4)*
     .(-(za(j1,j5)*zb(j1,j6))+za(j4,j5)*zb(j4,j6))))/
     .(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4)))/
     .((-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)))-
     .(za(j2,j4)*(za(j2,j5)*zb(j1,j2)+za(j4,j5)*zb(j1,j4))*zb(j3,j6))/
     .((-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*t(j1,j2,j4)))+
     .Lsm1_2mh(s(j2,j3),t(j1,j2,j4),s(j1,j4),s(j5,j6))*
     .(((za(j2,j3)*zb(j1,j2)-za(j3,j4)*zb(j1,j4))**2*
     .(-(za(j1,j5)*zb(j1,j2))+za(j4,j5)*zb(j2,j4))**2)/
     .(za(j5,j6)*zb(j1,j4)*(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))**3)-
     .(za(j3,j5)**2*zb(j1,j2)**2*t(j1,j2,j4)**2)/
     .(za(j5,j6)*zb(j1,j4)*(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))**3))+
     .Lsm1_2mh(s(j3,j4),t(j1,j2,j4),s(j1,j2),s(j5,j6))*
     .(-(((za(j2,j3)*zb(j1,j2)-za(j3,j4)*zb(j1,j4))**2*
     .(-(za(j1,j5)*zb(j1,j4))-za(j2,j5)*zb(j2,j4))**2)/
     .(za(j5,j6)*zb(j1,j4)*(-(za(j1,j3)*zb(j1,j4))-
     .za(j2,j3)*zb(j2,j4))**3*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))))+
     .(za(j3,j5)**2*zb(j1,j4)*t(j1,j2,j4)**2)/
     .(za(j5,j6)*(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))**3*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))))+
     .(2d0*Lnrat(-s(j1,j4),-s(j5,j6))*(-2d0*(s(j1,j2)-s(j3,j4))*za(j1,j4
     .)*zb(j1,j3)*
     .(-(za(j1,j5)*zb(j1,j6))-za(j4,j5)*zb(j4,j6))-
     .(za(j2,j4)*zb(j1,j6)*((s(j1,j2)-s(j3,j4))*
     .((s(j1,j4)-s(j2,j3)-s(j5,j6))*za(j1,j4)*zb(j4,j6)-
     .(-s(j1,j4)-s(j2,j3)+s(j5,j6))*za(j1,j5)*zb(j5,j6))+
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))*
     .(-((s(j1,j4)-s(j2,j3)-s(j5,j6))*za(j1,j4)*zb(j1,j6))-
     .(-s(j1,j4)-s(j2,j3)+s(j5,j6))*za(j4,j5)*zb(j5,j6))))/
     .(za(j2,j3)*zb(j5,j6))+
     .(zb(j1,j3)*(-(za(j1,j4)*
     .(za(j2,j5)*zb(j1,j2)+za(j3,j5)*zb(j1,j3)))+
     .za(j2,j4)*(-(za(j1,j5)*zb(j1,j2))+za(j4,j5)*zb(j2,j4))+
     .za(j3,j4)*(-(za(j1,j5)*zb(j1,j3))+za(j4,j5)*zb(j3,j4)))*
     .(-(za(j1,j4)*za(j5,j6)*zb(j4,j6))+za(j1,j5)*t(j1,j2,j4)))/
     .za(j5,j6)))/
     .((s(j1,j4)**2-2d0*s(j1,j4)*s(j2,j3)+s(j2,j3)**2-2d0*s(j1,j4)*s(j5,
     .j6)-
     .2d0*s(j2,j3)*s(j5,j6)+s(j5,j6)**2)*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)))

      Fcc_qpqmgpgm=Fcc_qpqmgpgm+
     .Lsm1(-s(j1,j4),-t(j1,j3,j4),-s(j3,j4),-t(j1,j3,j4))*
     .((za(j3,j4)**2*(za(j1,j3)*zb(j3,j6)+za(j1,j4)*zb(j4,j6))**2)/
     .(za(j1,j3)**3*(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*
     .zb(j5,j6)*t(j1,j3,j4))-
     .(za(j1,j4)**2*(-(za(j1,j3)*zb(j1,j6))+za(j3,j4)*zb(j4,j6))**2)/
     .(za(j1,j3)**3*(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*zb(j5,j
     .6)*
     .t(j1,j3,j4)))+(2d0*L0(-t(j1,j3,j4),-s(j3,j4))*za(j1,j4)*zb(j1,j3)*
     .(-(za(j1,j4)*zb(j1,j6))-za(j3,j4)*zb(j3,j6))*
     .(za(j1,j3)*zb(j3,j6)+za(j1,j4)*zb(j4,j6)))/
     .(s(j3,j4)*za(j1,j3)*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))*
     .zb(j5,j6)*t(j1,j3,j4))+(2d0*L0(-t(j1,j3,j4),-s(j1,j4))*za(j1,j4)*z
     .b(j1,j3)*
     .(-(za(j1,j4)*zb(j1,j6))-za(j3,j4)*zb(j3,j6))*
     .(-(za(j1,j3)*zb(j1,j6))+za(j3,j4)*zb(j4,j6)))/
     .(s(j1,j4)*za(j1,j3)*(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*
     .zb(j5,j6)*t(j1,j3,j4))+Lsm1_2mh(s(j1,j2),t(j1,j3,j4),s(j3,j4),
     .s(j5,j6))*(-((za(j2,j5)**2*zb(j1,j3)**3)/
     .(za(j5,j6)*zb(j1,j4)*(-(za(j2,j3)*zb(j1,j3))-
     .za(j2,j4)*zb(j1,j4))*zb(j3,j4)*t(j1,j3,j4)))-
     .((-(za(j1,j4)*zb(j1,j2))+za(j3,j4)*zb(j2,j3))**3*
     .(za(j1,j3)*zb(j3,j6)+za(j1,j4)*zb(j4,j6))**2)/
     .(za(j3,j4)*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))**3*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*zb(j5,j6)*t(j1,j3,j4)
     .)+(za(j1,j4)**2*(-(za(j1,j4)*zb(j1,j2))+za(j3,j4)*zb(j2,j3))*zb(j2
     .,j6)**2*
     .t(j1,j3,j4))/
     .(za(j3,j4)*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))**3*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*zb(j5,j6)))+
     .(I3m(s(j1,j4),s(j2,j3),s(j5,j6))*zb(j1,j3)*
     .(-(s(j5,j6)*za(j2,j4)*((-s(j1,j4)-s(j2,j3)+s(j5,j6))*za(j1,j5)*
     .zb(j2,j6)-2d0*za(j1,j4)*za(j3,j5)*zb(j2,j3)*zb(j4,j6)))+
     .2d0*za(j1,j4)*za(j1,j5)*za(j2,j5)*
     .((-s(j1,j4)-s(j2,j3)+s(j5,j6))*zb(j1,j2)-
     .2d0*za(j3,j4)*zb(j1,j4)*zb(j2,j3))*zb(j5,j6)-
     .2d0*za(j1,j4)*za(j4,j5)*zb(j4,j6)*
     .(s(j5,j6)*(-s(j1,j4)-s(j2,j3)+s(j5,j6))+
     .(s(j1,j4)-s(j2,j3)-s(j5,j6))*t(j1,j3,j4))))/
     .((s(j1,j4)**2-2d0*s(j1,j4)*s(j2,j3)+s(j2,j3)**2-2d0*s(j1,j4)*s(j5,
     .j6)-
     .2d0*s(j2,j3)*s(j5,j6)+s(j5,j6)**2)*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)))

      Fcc_qpqmgpgm=Fcc_qpqmgpgm+
     .Lsm1(-s(j2,j3),-t(j2,j3,j4),-s(j3,j4),-t(j2,j3,j4))*
     .(-(((za(j3,j5)*zb(j2,j3)+za(j4,j5)*zb(j2,j4))**2*zb(j3,j4)**2)/
     .(za(j5,j6)*zb(j2,j4)**3*(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))*
     .t(j2,j3,j4)))+(zb(j2,j3)**2*
     .(-(za(j2,j5)*zb(j2,j4))-za(j3,j5)*zb(j3,j4))**2)/
     .(za(j5,j6)*zb(j2,j4)**3*(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))*
     .t(j2,j3,j4)))-(2d0*L0(-t(j2,j3,j4),-s(j3,j4))*za(j2,j4)*zb(j2,j3)*
     .(za(j3,j5)*zb(j2,j3)+za(j4,j5)*zb(j2,j4))*
     .(-(za(j2,j5)*zb(j2,j3))+za(j4,j5)*zb(j3,j4)))/
     .(s(j3,j4)*za(j5,j6)*zb(j2,j4)*
     .(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))*t(j2,j3,j4))-
     .(2d0*L0(-t(j2,j3,j4),-s(j2,j3))*za(j2,j4)*zb(j2,j3)*
     .(-(za(j2,j5)*zb(j2,j4))-za(j3,j5)*zb(j3,j4))*
     .(-(za(j2,j5)*zb(j2,j3))+za(j4,j5)*zb(j3,j4)))/
     .(s(j2,j3)*za(j5,j6)*zb(j2,j4)*(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3
     .,j4))*
     .t(j2,j3,j4))+Lsm1_2mh(s(j1,j2),t(j2,j3,j4),s(j3,j4),s(j5,j6))*
     .(-(((za(j3,j5)*zb(j2,j3)+za(j4,j5)*zb(j2,j4))**2*
     .(za(j1,j2)*zb(j2,j3)-za(j1,j4)*zb(j3,j4))**3)/
     .(za(j5,j6)*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))**3*
     .zb(j3,j4)*(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))*t(j2,j3,j4))
     .)-(za(j2,j4)**3*zb(j1,j6)**2)/
     .(za(j2,j3)*za(j3,j4)*(-(za(j2,j3)*zb(j1,j3))-za(j2,j4)*zb(j1,j4))*
     .zb(j5,j6)*t(j2,j3,j4))+
     .(za(j1,j5)**2*zb(j2,j3)**2*(za(j1,j2)*zb(j2,j3)-za(j1,j4)*zb(j3,j4
     .))*
     .t(j2,j3,j4))/
     .(za(j5,j6)*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))**3*zb(j3,j
     .4)*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))))+
     .(I3m(s(j2,j3),s(j1,j4),s(j5,j6))*za(j2,j4)*
     .(-2d0*za(j5,j6)*zb(j1,j6)*zb(j2,j3)*zb(j2,j6)*
     .(-((-s(j1,j4)-s(j2,j3)+s(j5,j6))*za(j1,j2))+
     .2d0*za(j1,j4)*za(j2,j3)*zb(j3,j4))-
     .s(j5,j6)*zb(j1,j3)*((-s(j1,j4)-s(j2,j3)+s(j5,j6))*za(j1,j5)*
     .zb(j2,j6)-2d0*za(j1,j4)*za(j3,j5)*zb(j2,j3)*zb(j4,j6))-
     .2d0*za(j3,j5)*zb(j2,j3)*zb(j3,j6)*
     .(s(j5,j6)*(-s(j1,j4)-s(j2,j3)+s(j5,j6))+
     .(-s(j1,j4)+s(j2,j3)-s(j5,j6))*t(j2,j3,j4))))/
     .((s(j1,j4)**2-2d0*s(j1,j4)*s(j2,j3)+s(j2,j3)**2-2d0*s(j1,j4)*s(j5,
     .j6)-
     .2d0*s(j2,j3)*s(j5,j6)+s(j5,j6)**2)*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)))

      Fcc_qpqmgpgm=Fcc_qpqmgpgm+
     .I3m(s(j1,j2),s(j3,j4),s(j5,j6))*
     .((za(j3,j5)*zb(j2,j3)*(s(j1,j4)*za(j1,j4)*zb(j1,j6)-
     .s(j2,j4)*za(j2,j4)*zb(j2,j6))-
     .5d0*za(j1,j3)*za(j2,j4)*za(j3,j5)*zb(j1,j2)*zb(j2,j3)*zb(j3,j6)-
     .(-s(j1,j3)+s(j2,j4))*za(j3,j4)*
     .(-(za(j1,j5)*zb(j1,j2))+za(j4,j5)*zb(j2,j4))*zb(j3,j6)-
     .za(j1,j4)*zb(j1,j2)*(-(s(j3,j4)*za(j1,j5)*zb(j1,j6))-
     .2d0*(-((s(j1,j3)+s(j2,j3))*za(j2,j5)*zb(j2,j6))-
     .za(j2,j3)*za(j4,j5)*zb(j2,j6)*zb(j3,j4))-
     .s(j2,j4)*za(j3,j5)*zb(j3,j6)+
     .(2d0*za(j2,j3)*zb(j1,j2)*zb(j3,j4)*
     .(-(za(j2,j5)*zb(j2,j6))-za(j4,j5)*zb(j4,j6)))/zb(j1,j4)+
     .3d0*(za(j2,j5)*za(j3,j4)*zb(j2,j4)*zb(j3,j6)+
     .s(j1,j3)*za(j4,j5)*zb(j4,j6))))/
     .(2d0*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))*
     .(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4)))-
     .(za(j2,j4)*zb(j1,j2)*(za(j2,j5)*zb(j1,j2)+za(j4,j5)*zb(j1,j4))*
     .zb(j3,j6))/
     .(zb(j1,j4)*(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))*t(j1,j2,j4
     .))
     .-(za(j2,j5)*zb(j1,j3)**2*(-(za(j1,j4)*zb(j1,j6))-za(j3,j4)*zb(j3,j
     .6)))/
     .(zb(j1,j4)*t(j1,j3,j4)**2)+
     .(za(j2,j5)*zb(j1,j3)**2*(-(za(j2,j5)*zb(j2,j3))+
     .za(j4,j5)*zb(j3,j4)))/
     .(2d0*za(j5,j6)*zb(j1,j4)*zb(j3,j4)*t(j1,j3,j4))+
     .(za(j2,j4)*zb(j1,j3)*zb(j1,j6)*
     .(-(za(j1,j4)*zb(j1,j6))-za(j3,j4)*zb(j3,j6)))/
     .(2d0*za(j3,j4)*zb(j1,j4)*zb(j5,j6)*t(j1,j3,j4))+
     .((za(j1,j5)**2*zb(j1,j3)**3*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4)))/
     .(za(j5,j6)*zb(j1,j4)*zb(j3,j4))+
     .(zb(j1,j2)*((-s(j1,j2)+2d0*s(j1,j3)+s(j3,j4)-s(j5,j6))*
     .za(j4,j5)-2d0*za(j2,j4)*za(j5,j6)*zb(j2,j6))*zb(j3,j6))/
     .zb(j1,j4)+(zb(j1,j3)*
     .(-(za(j1,j3)*za(j1,j5)*zb(j1,j2)*zb(j1,j3))+
     .za(j1,j3)*za(j2,j5)*zb(j1,j2)*zb(j2,j3)-
     .za(j1,j5)*za(j3,j4)*zb(j1,j3)*zb(j2,j4)+
     .za(j2,j5)*za(j3,j4)*zb(j2,j3)*zb(j2,j4))*zb(j3,j6))/
     .(zb(j1,j4)*zb(j3,j4))+
     .(za(j2,j4)**2*(-((s(j1,j2)-s(j3,j4)-s(j5,j6))*zb(j1,j2))+
     .za(j1,j3)*zb(j1,j2)*zb(j1,j3)+
     .za(j3,j4)*zb(j1,j3)*zb(j2,j4))*zb(j2,j6)**2)/
     .(za(j3,j4)*zb(j1,j4)*zb(j5,j6))+
     .(4d0*za(j1,j4)*za(j2,j5)*zb(j1,j3)*
     .(-(za(j1,j4)*zb(j1,j2))+za(j3,j4)*zb(j2,j3))*zb(j2,j6))/
     .t(j1,j3,j4)+((-(za(j1,j4)*zb(j1,j2))+za(j3,j4)*zb(j2,j3))*
     .(-(za(j1,j4)*zb(j1,j6))-za(j3,j4)*zb(j3,j6))**2*
     .(2d0*s(j3,j4)*s(j5,j6)+(s(j1,j2)-s(j3,j4)-s(j5,j6))*t(j1,j3,j4))
     .)/(za(j3,j4)*zb(j5,j6)*t(j1,j3,j4)**2)+
     .(za(j4,j5)*zb(j1,j2)*(za(j1,j4)*zb(j1,j6)-za(j2,j4)*zb(j2,j6))*
     .(s(j1,j4)+t(j2,j3,j4)))/(za(j3,j4)*zb(j1,j4)))/
     .(2d0*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))*
     .(-(za(j1,j3)*zb(j1,j2))-za(j3,j4)*zb(j2,j4))))

      Fcc_qpqmgpgm=Fcc_qpqmgpgm+
     .I3m(s(j1,j2),s(j3,j4),s(j5,j6))*
     .(((s(j1,j3)-s(j2,j4))*za(j4,j5)*zb(j3,j4)*
     .(za(j1,j2)*zb(j2,j6)+za(j1,j3)*zb(j3,j6))+
     .za(j1,j4)*(-(s(j1,j3)*za(j1,j5)*zb(j1,j3))+
     .s(j2,j3)*za(j2,j5)*zb(j2,j3))*zb(j4,j6)+
     .5d0*za(j1,j2)*za(j1,j4)*za(j4,j5)*zb(j1,j3)*zb(j2,j4)*zb(j4,j6)+
     .za(j1,j2)*zb(j2,j3)*(-(s(j3,j4)*za(j2,j5)*zb(j2,j6))+
     .(2d0*za(j1,j2)*za(j3,j4)*zb(j1,j4)*
     .(-(za(j1,j5)*zb(j1,j6))-za(j3,j5)*zb(j3,j6)))/za(j2,j3)+
     .3d0*(-(za(j1,j3)*za(j4,j5)*zb(j1,j6)*zb(j3,j4))+
     .s(j2,j4)*za(j3,j5)*zb(j3,j6))-
     .2d0*(-((s(j1,j4)+s(j2,j4))*za(j1,j5)*zb(j1,j6))+
     .za(j1,j5)*za(j3,j4)*zb(j1,j4)*zb(j3,j6))-
     .s(j1,j3)*za(j4,j5)*zb(j4,j6)))/
     .(2d0*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))*
     .(-(za(j1,j3)*zb(j1,j4))-za(j2,j3)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)))+
     .(za(j1,j2)*za(j4,j5)*zb(j1,j3)*
     .(-(za(j1,j2)*zb(j1,j6))+za(j2,j3)*zb(j3,j6)))/
     .(za(j2,j3)*(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))*t(j1,j2,j3))-
     .(za(j2,j4)**2*zb(j1,j6)*(-(za(j2,j5)*zb(j2,j3))+za(j4,j5)*zb(j3,j4
     .)))/
     .(za(j2,j3)*t(j2,j3,j4)**2)+
     .(za(j2,j4)*za(j2,j5)*zb(j1,j3)*
     .(-(za(j2,j5)*zb(j2,j3))+za(j4,j5)*zb(j3,j4)))/
     .(2d0*za(j2,j3)*za(j5,j6)*zb(j3,j4)*t(j2,j3,j4))+
     .(za(j2,j4)**2*zb(j1,j6)*(-(za(j1,j4)*zb(j1,j6))-za(j3,j4)*zb(j3,j6
     .)))/
     .(2d0*za(j2,j3)*za(j3,j4)*zb(j5,j6)*t(j2,j3,j4))+
     .((za(j1,j5)**2*zb(j1,j3)**2*
     .((s(j1,j2)-s(j3,j4)-s(j5,j6))*za(j1,j2)-
     .za(j1,j2)*za(j2,j4)*zb(j2,j4)-za(j1,j3)*za(j2,j4)*zb(j3,j4)
     .))/(za(j2,j3)*za(j5,j6)*zb(j3,j4))-
     .(za(j2,j4)*za(j4,j5)*(-(za(j1,j2)*za(j1,j4)*zb(j1,j6)*
     .zb(j2,j4))+za(j1,j2)*za(j2,j4)*zb(j2,j4)*zb(j2,j6)-
     .za(j1,j3)*za(j1,j4)*zb(j1,j6)*zb(j3,j4)+
     .za(j1,j3)*za(j2,j4)*zb(j2,j6)*zb(j3,j4)))/
     .(za(j2,j3)*za(j3,j4))+
     .(za(j2,j4)**3*zb(j2,j6)**2*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)))/
     .(za(j2,j3)*za(j3,j4)*zb(j5,j6))-
     .(za(j1,j2)*za(j4,j5)*((-s(j1,j2)+2d0*s(j2,j4)+s(j3,j4)-s(j5,j6))*
     .zb(j3,j6)+2d0*za(j1,j5)*zb(j1,j3)*zb(j5,j6)))/za(j2,j3)+
     .(za(j1,j2)*(-(za(j1,j5)*zb(j1,j3))+za(j2,j5)*zb(j2,j3))*
     .zb(j3,j6)*(s(j2,j3)+t(j1,j3,j4)))/(za(j2,j3)*zb(j3,j4))+
     .(4d0*za(j1,j5)*za(j2,j4)*zb(j1,j6)*zb(j2,j3)*
     .(za(j1,j2)*zb(j2,j3)-za(j1,j4)*zb(j3,j4)))/t(j2,j3,j4)+
     .((za(j1,j2)*zb(j2,j3)-za(j1,j4)*zb(j3,j4))*
     .(-(za(j2,j5)*zb(j2,j3))+za(j4,j5)*zb(j3,j4))**2*
     .(2d0*s(j3,j4)*s(j5,j6)+(s(j1,j2)-s(j3,j4)-s(j5,j6))*t(j2,j3,j4)))/
     .(za(j5,j6)*zb(j3,j4)*t(j2,j3,j4)**2))/
     .(2d0*(-(za(j1,j3)*zb(j2,j3))-za(j1,j4)*zb(j2,j4))*
     .(za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4))))

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


    !use crop2mlModules
        REAL::oldcanopyShootNumber, calc
        INTEGER:: emergedLeaves, i
        oldcanopyShootNumber = canopyShootNumber
        emergedLeaves = INT(MAX(1, CEILING(leafNumber - 1)))
        shoots = fibonacci(emergedLeaves)
        canopyShootNumber = MIN(shoots * sowingDensity, INT(targetFertileShoot))
        averageShootNumberPerPlant = canopyShootNumber / sowingDensity
        IF (canopyShootNumber /= oldcanopyShootNumber) THEN
            CALL AddToList(tilleringProfile,canopyShootNumber - oldcanopyShootNumber)
        END IF
         tillerNumber = SIZE(tilleringProfile)
        DO i=SIZE(leafTillerNumberArray), CEILING(leafNumber)-1
            CALL AddToListInt(leafTillerNumberArray,tillerNumber)
        END DO


text:
	real*4 function zx(tx, h, q, lamb)

c	Returns the value of Z = ZX (redshift) corresponding to a given
c	light travel time TX (measured in Gyr).

c	H = Ho in km/sec/Mpc
c	Q = qo (if problems with qo = 0, try 0.0001)

	implicit none
	real*4 ltt, z(47), age, tx, h, q, hh0, dz, t, omega0, omegainv, lamb, zi
	integer j
	external t
	data z /0.,.001,.002,.004,.006,.008,.01,.02,.04,.06,.08,.1,.2,.3,
     &      .4,.5,.545,.6,.7,.8,.9,.945,1.,1.2,1.4,1.6,1.8,2.,3.,4.,5.,6.,
     &      7.,8.,9.,10.,12.,14.,16.,18.,20.,30.,40.,60.,80.,100.,1000./

c     data z / 0., .001, .002, .004, .006, .008, .01, .02, .04, .06, .08
c    &, .1, .2, .3, .4, .5, .545, .6, .7, .8, .9, .945, 1., 1.2, 1.4,
c    &1.6, 1.8, 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 14., 16., 18.
c    &, 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80.
c    &, 85., 90., 95., 100., 500., 1000. /

c	Define light travel time
	ltt(q,zi) = age - t(h,q,zi,lamb)

c	Compute omega
	if (lamb .ne. 0.) then
		omega0 = (2. * (q + 1.)) / 3.
		omegainv = 1. / omega0
	endif

c	Check for zero value
	zx = 0.
	if (tx .eq. 0.) return

c	Check for maximum age of universe
	age = t(h,q,0.,lamb)
	zx = -2.
	if (tx .ge. age) return

c	Express Ho in billion years ** (-1)
	hh0 = h * 0.001022

c	Check for q=0.5 case
	if (q .eq. 0.5) then
		zx = ((1. - (((3. * hh0) * tx) / 2.)) ** (- (2. / 3.))) - 1.
		return
	endif

c	General case
	do j = 1, 47
	if (tx .le. ltt(q,z(j))) goto 2
	enddo
2	if (j > 1) then
		zx = z(j - 1)
	else
		zx = z(1)
	endif
	zi=zx

c	Iterate to find best value
	do j = 1, 100000
	if (q .ge. 0.) then
c		zero cosmological constant
		dz = ((hh0 * ((tx - t(h,q,0.,lamb)) + t(h,q,zx,lamb))) *
     *		     ((1. + zx) ** 2)) * ((1. + ((2. * q) * zx)) ** 0.5)
	elseif (q .lt. 0.) then
c		non-zero cosmological constant
		dz = ((hh0 * ((tx - t(h,q,0.,lamb)) + t(h,q,zx,lamb))) *
     *		     (1. + zx)) * sqrt(omega0)
		dz = dz * sqrt((((1. + zx) ** 3.) + omegainv) - 1.)
	endif
c	if (abs(dz / zx) .lt. 0.00001) return
c	if (abs(dz / zx) .lt. 0.00005) return
c	if (abs(dz / zx) .lt. 0.00010)  return
	if (abs(dz / zx) .lt. 0.00003) return
	zx = zx + dz
	enddo
	return
	end

max_stars_repo_name: jchavesmontero/galaxpy

max_issues_repo_name: jchavesmontero/galaxpy

max_forks_repo_name: jchavesmontero/galaxpy


text:
SUBROUTINE get_proj_shadows(hx, hy, ns, nl, &
    htol, phi_sun, sun_zen, zmax, zmin, a, mask, h_offset, &
    n_inc, m_inc, aoff_x, aoff_y, nsA, nlA, k_setting, &
    dem_nr, dem_nc, nlA_ori, nsA_ori, ierr)

    implicit none

    integer ns, nl, dem_nr, dem_nc, nlA_ori, nsA_ori, k_setting
    integer k_max, n_add, m_add, ierr, az_case, i, j
!   NOTE: n_inc and m_inc are Floating Point arrays
    real n_inc(k_setting), m_inc(k_setting)
    real h_offset(k_setting)
    real phi_sun,zmax,zmin,sun_zen,htol
    real*8 hx, hy
    real d, d0, a(dem_nr, dem_nc)
    integer ncho, aoff_x, aoff_y, nsA, nlA
    integer t_aoff_x,t_aoff_y,t_nsA,t_nlA
    integer tval(8)
    logical rstatus
    integer*2 mask(nlA_ori, nsA_ori)
    integer*2 rmax, rmin

    real pi,r2d,d2r,eps
    common/base/pi,r2d,d2r,eps
!
!   calculate the border info for the sun position
!   In Australia and in the south in particular case=1
!   for Landsat since Landsat overpass is around 10:am
!   local time.
!   That is the sun azimuth is between East and North

    call set_borderf(rstatus,phi_sun,zmax,zmin,sun_zen,hx,hy,az_case, &
        d,d0,k_max,h_offset,n_inc,m_inc,n_add,m_add,k_setting,k_setting,ierr)

    if (.not.rstatus) then
        goto 99
    endif

!   define the maximum sized subset that can be processed
    if (az_case.eq.1) then
        t_aoff_x=0
        t_aoff_y=m_add
        t_nsA=ns-n_add
        t_nlA=nl-m_add
    else if (az_case.eq.2) then
        t_aoff_x=n_add
        t_aoff_y=m_add
        t_nsA=ns-n_add
        t_nlA=nl-m_add
    else if (az_case.eq.3) then
        t_aoff_x=n_add
        t_aoff_y=0
        t_nsA=ns-n_add
        t_nlA=nl-m_add
    else if (az_case.eq.4) then
        t_aoff_x=0
        t_aoff_y=0
        t_nsA=ns-n_add
        t_nlA=nl-m_add
    endif

!   Set the sub-matrix A where shade is to be found
!   aoff_x is the offset in samples
!   aoff_y is the offset in lines
!   pos in line in image for A(i,j) is aoff_x+j
!   pos in lines in image for A(i,j) is aoff_y+i

!   check the submatrix A is valid in various ways

!   first that the individual components are valid
    if (((aoff_x.lt.0) .or. (aoff_x.ge.ns)) .or. ((aoff_y.lt.0) .or. &
        (aoff_y.ge.nl)) .or. &
        ((nsA.lt.1) .or. (nsA.gt.ns)) .or. &
        ((nlA.lt.1) .or. (nlA.gt.nl))) then
            ierr = 71
            goto 99
    endif

!   Check A is embedded in the whole image
    if((aoff_x+nsA.gt.ns) .or. (aoff_y+nlA.gt.nl)) then
        ierr = 72
        goto 99
    endif

!   check the sub-image A plus the border area
!   needed to test A still fits inside the main image
!   with the buffer available
!   NOTE: treat the four cases in pairs

    if(az_case.eq.1.or.az_case.eq.2) then
        if (aoff_y.lt.m_add) then
            ierr = 73
            goto 99
        endif
    else if (az_case.eq.3.or.az_case.eq.4) then
        if (aoff_y+nlA+m_add.gt.nl) then
            ierr = 73
            goto 99
        endif
    endif

    if(az_case.eq.2.or.az_case.eq.3) then
        if (aoff_x.lt.n_add) then
            ierr = 74
            goto 99
        endif
    else if (az_case.eq.1.or.az_case.eq.4) then
        if (aoff_x+nsA+n_add.gt.ns) then
            ierr = 74
            goto 99
        endif
    endif
!    now set up the mask image to record shade pixels in
!    A NOTE: mask has the dimensions of A and not the
!    DEM - set as a 1-D array so that it can be indexed
!    in the subroutine without worrying about set bounds

!     Mask is a long integer (integer*4) due to recl issue
!     First set to 1 so zero will represent deep shadow

    do i=1,nlA
        do j=1,nsA
            mask(i,j)=1
        enddo
    enddo

    rmax=maxval(mask(1:nlA,1:nsA))
    rmin=minval(mask(1:nlA,1:nsA))

!   proj_terrain does the job of checking for occlusion
!   along the vector from a pixel to the sun
!   if any terrain obstructs the path the pixel is set to
!   zero in Mask

    call proj_terrain(ns, nl, nsA, nlA, a, mask, Aoff_x, Aoff_y, &
    k_max, n_inc, m_inc, h_offset, zmax, htol, dem_nr, dem_nc, nlA_ori, nsA_ori)

!   Description of proj_terrain
!   ===========================
!   call proj_terrain(n_max,m_max,n,m,z,mask,n_off,m_off,k_max,
!     . n_inc,m_inc,h_offset,zmax,htol,k_setting)
!
!   subroutine to construct mask of shade pixels
!
!   z(m_max,n_max) is the main array of heights
!   A(m,n) is the (sub-)matrix of target heights in:
!      z(m_off+1,n_off+1) to z(m_off+m,n_off+n)
!      mask(m,n) is the output mask
!      on input assumed to be 1 where valid data exist 0 else
!      k_max is the number of lags in the projection
!      n_inc(k_max) real increments column number for projection
!      m_inc(k_max)  real increments row number for the projection
!      h_offset(k_max) is the height of the projection hor each lag
!      zmax is the maximum altitude in the whole array
!      htol is a tolerance (m) for the test for a hit (RMS error in z)

99  continue
    return
END SUBROUTINE get_proj_shadows

max_stars_repo_name: Oceancolour-RG/wagl

max_issues_repo_name: Oceancolour-RG/wagl

max_forks_repo_name: Oceancolour-RG/wagl


text:
! S08C-cant-assert-keyword-nonnull.f90
! In 'c_action_actual_arg_spec', assertion 'keyword' != NULL is incorrect.
! CASE C: rule for 'actual_arg_spec'.
!
! An alternate return specifier following a keyword argument in an
! actual argument list causes the front end to fail an assertion.

! The subroutine definition for reference below is commented out because
! it causes front end to fail another assertion before showing this bug.
subroutine g(k, *, *)
   integer :: k
   return 2
end subroutine

program p
    ! The following call is incorrect because F90 requires that
    ! positional arguments come before keyword arguments.
    ! However, OFP accepts it and ROSE fails an assertion on it.
    call g(k=0, *100, *200)  ! assertion failure: 'keyword' is null
100 continue
200 continue
end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
program t
INTERFACE SWITCH
SUBROUTINE INT_SWITCH (X, Y)
INTEGER, INTENT (INOUT) :: X, Y
END SUBROUTINE INT_SWITCH
SUBROUTINE REAL_SWITCH (X, Y)
REAL, INTENT (INOUT) :: X, Y
END SUBROUTINE REAL_SWITCH
SUBROUTINE COMPLEX_SWITCH (X, Y)
COMPLEX, INTENT (INOUT) :: X, Y
END SUBROUTINE COMPLEX_SWITCH
END INTERFACE SWITCH
PRINT *,'ok'
end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
      subroutine heat
!     dummy to be modified for special calculations of heat transfer
      use httr, only:qfun,qfunw
      implicit none
      qfunw=qfun
      end subroutine heat

max_stars_repo_name: aldorona/contur

max_issues_repo_name: aldorona/contur

max_forks_repo_name: aldorona/contur


text:
      subroutine steady(iflg,inflow,inflowe)
!***********************************************************************
!  Copyright, 1993, 2004,  The  Regents of the University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C***********************************************************************
CD1
CD1 PURPOSE
CD1
CD1 Determime if steady state solution achieved
CD1
C***********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 18-OCT-93    Z. Dash        22      Add prolog.
CD2              G. Zyvoloski           Initial implementation.
CD2
CD2 $Log:   /pvcs.config/fehm90/src/steady.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:06   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:17:24   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:14:46   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:27:02   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:10:34   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:46:00 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.2   Mon Feb 05 11:30:58 1996   hend
CD2 Updated Requirements Traceability
CD2
CD2    Rev 1.1   03/18/94 15:55:52   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.0   01/20/94 10:28:18   pvcs
CD2 original version in process of being certified
CD2
C***********************************************************************
CD3
CD3 INTERFACES
CD3
CD3 Formal Calling Parameters
CD3
CD3   None
CD3
CD3 Interface Tables
CD3
CD3   None
CD3
CD3 Files
CD3
CD3   Name                     Use  Description
CD3
CD3   isave                    I/O  Unit number for restart file (to write)
CD3
C***********************************************************************
CD4
CD4 GLOBAL OBJECTS
CD4
CD4 Global Constants
CD4
CD4   None
CD4
CD4 Global Types
CD4
CD4   None
CD4
CD4 Global Variables
CD4
CD4                            COMMON
CD4   Identifier      Type     Block  Description
CD4
CD4   am0             REAL*8   faar   Initial mass in problem
CD4   ame             REAL*8   faar   ?
CD4   astmo           REAL*8   faar   Initial steam mass in problem
CD4   cord            REAL*8   fbs    Contains the coordinates of each node
CD4   deneh           REAL*8   fdd    Last time step energy accumulation term
CD4                                     at each node
CD4   denei           REAL*8   fcc    Energy accumulation term
CD4   denej           REAL*8   fdd    Last time step energy accumulation time
CD4                                     derivative at each node
CD4   denh            REAL*8   fdd    Last time step mass accumulation term at
CD4                                     each node
CD4   deni            REAL*8   fcc    Mass accumulation term
CD4   denj            REAL*8   fdd    Last time step mass accumulation time
CD4                                     derivative at each node
CD4   dstm            REAL*8   fcc    Steam mass
CD4   dtot            REAL*8   faar   Current time step size in seconds
CD4   eflow           REAL*8   fdd    Energy flow at each source node
CD4   grav            REAL*8   faar   Value of gravity
CD4   icnl            INT      faai   Problem dimension
CD4   isave           INT      faai   Unit number for restart file (to write)
CD4   ka              INT      fbb    Contains boundary type information for
CD4                                     each node
CD4   n               INT      faai   Total number of nodes
CD4   pflow           REAL*8   fdd    Flowing pressure at each source node
CD4   phi             REAL*8   fdd    Pressure at each node
CD4   phini           REAL*8   fdd2   Initial pressure at each node
CD4   pho             REAL*8   fdd    Last time step pressure at each node
CD4   pnx             REAL*8   fdd    ?Permeability in the x-direction, liquid
CD4                                     velocity in the x-direction, vapor
CD4                                     velocity in the x-direction
CD4   pny             REAL*8   fdd    ?Permeability in the y-direction liquid
CD4                                     velocity in the y direction, vapor
CD4                                     velocity in the y-direction
CD4   pnz             REAL*8   fdd    ?Permeability in the z-direction, liquid
CD4                                     velocity in the z-direction, vapor
CD4                                     velocity in the z-direction
CD4   ps              REAL*8   fdd    Porosity at each node
CD4   psini           REAL*8   fdd2   Initial porosity at each node
CD4   qflux           REAL*8   fdd    Heat flux at each node
CD4   qflxm           REAL*8   fdd    Heat flux impedance at each node
CD4   qh              REAL*8   fdd    Energy source term at each node
CD4   s               REAL*8   fdd    Liquid saturation at each node
CD4   sk              REAL*8   fdd    Source strength of each node
CD4   so              REAL*8   fdd    Last time step saturation at each node
CD4   t               REAL*8   fdd    Temperature at each node
CD4   thx             REAL*8   fdd    Thermal conductivity x-direction
CD4   thy             REAL*8   fdd    Thermal conductivity y-direction
CD4   thz             REAL*8   fdd    Thermal conductivity z-direction
CD4   tini            REAL*8   fdd2   Initial temperature at each node
CD4   to              REAL*8   fdd    Last time step temperature at each node
CD4   volume          REAL*8   fdd    Volume associated at each node
CD4   wellim          REAL*8   fdd    Well impedance at each source node
CD4
CD4 Global Subprograms
CD4
CD4   Identifier      Type     Description
CD4
CD4   bnswer                   Initiate thermodynamics and solution routines
CD4   enthp           REAL*8   Calculate enthalpy as a function of pressure
CD4                              and temperature
CD4
C***********************************************************************
CD5
CD5 LOCAL IDENTIFIERS
CD5
CD5 Local Constants
CD5
CD5   None
CD5
CD5 Local Types
CD5
CD5   None
CD5
CD5 Local variables
CD5
CD5   Identifier      Type     Description
CD5
CD5   cztop           REAL*8   Maximum coordinate
CD5   etoll           REAL*8   ?Coordinate tolerance
CD5   i               INT      Loop index
CD5   ir              INT      Problem dimension
CD5
CD5 Local Subprograms
CD5
CD5   None
CD5
C***********************************************************************
CD6
CD6 FUNCTIONAL DESCRIPTION
CD6
C***********************************************************************
CD7
CD7 ASSUMPTIONS AND LIMITATIONS
CD7
CD7 None
CD7
C***********************************************************************
CD8
CD8 SPECIAL COMMENTS
CD8
CD8 This routine will be used in the future but is not currently
CD8 used by the code.
CD8
C***********************************************************************
CD9
CD9 REQUIREMENTS TRACEABILITY
CD9
CD9 Not Applicable.  See Special Comments.
CD9
C***********************************************************************

      use comai
      use combi
      use comci
      use comdi
      use comdti
      use comfi
      use comsplitts
      use comsteady
      use davidi
      implicit none

      integer i, iflg, itt
      real*8 divisor, divisore, inflow, inflowe, btol_save, tolde_save
      real*8 flow_max, ratio2, balfac1, balfac2, balfac3, balfac12
      real*8 balfac13, tol_divisor, tmch1_min, tmch1_max
      character*20 dummy
      logical null1
      parameter (tol_divisor= 1.d-8)

      if(isteady.eq.0) return
      if(iflg .eq. 0) then
c     Reading steady input parameters
         toldh = tolerance
         toldp = tolerance
         tolds = tolerance
         toldt = tolerance
         toldc = tolerance
         tolde = tolerance
         tacc = tolerance
         tol_str = tolerance
         spercent = .false.
         divisor = 1.
         divisore = 1.
         smult = 1.
         sday = 0.
         shtl = 0.0
         sdmx = 0.0
         stmch = 0.0
         balance_tol = tolerance
         tmch_old = 1.d20
         svar_flag = .false.
         sflux_flag = .false.
         info_string = 'absolute difference'
         time_ss = 1.d20
         sminstep = 2
         if (ifree .ne. 0)  then
            snstep = 150
            smult = 2.0
         else if (ico2 .lt. 0 ) then
            snstep = 500
            smult = 2.0
         else if (ico2 .eq. 0) then
            snstep = 200
            smult = 2.0
         else if (ico2 .gt. 0) then
            snstep = 200
            smult = 1.5
         end if

         do
            read (inpt, '(a80)') wdd1
            if ( null1(wdd1) .or. wdd1(1:3) .eq. 'end' ) exit
            select case (wdd1(1:4))
            case ('shea')
               read (wdd1, *) dummy, toldh
               toldp = toldh * 997. * 9.8d-6
               svar_flag = .true.
            case('spre')
               read (wdd1, *) dummy, toldp
               svar_flag = .true.
            case ('stem')
               read (wdd1, *) dummy, toldt
               svar_flag = .true.
            case('ssat')
               read (wdd1, *) dummy, tolds
               svar_flag = .true.
            case('sair')
               read (wdd1, *) dummy, toldc
               svar_flag = .true.
            case('sflu')
               read (wdd1, *) dummy, balance_tol
               sflux_flag = .true.
            case('sent')
               read (wdd1, *) dummy, tolde
               sflux_flag = .true.
            case ('stim')
               read (wdd1, *) dummy, time_ss
            case ('smul')
               read (wdd1, *) dummy, smult
            case('sday')
               read (wdd1, *) dummy, sday
            case('sdmx')
               read (wdd1, *) dummy, sdmx
            case('snst')
               read (wdd1, *) dummy, snstep
            case('smst')
               read (wdd1, *) dummy, sminstep
            case('shtl')
               read (wdd1, *) dummy, shtl
            case('sacc')
               read (wdd1, *) dummy, tacc
            case('stmc')
               read (wdd1, *) dummy, stmch
            case('smas')
               read (wdd1, *) dummy, smass
            case('tstr')
               read (wdd1, *) dummy, tol_str
            case('stss')
               read (wdd1, *) dummy, stsstr
            case('sper')
               spercent = .true.
               info_string = 'fractional change'
            end select
         end do

      else if(iflg.eq.-1) then
c 'steady macro overrides usual stopping time tims (gaz 032405)
       tims_trans = tims
       daymax_save = daymax
       tims = time_ss
       if(sdmx.gt.0.0) then
        daymax = sdmx
       else
        daymax = tims
       endif
	 daycs_save=daycs
	 daycf_save=daycf
	 dayhf_save=dayhf
	 dayhs_save=dayhs
       dayhf = 1.d50
	 dayhs = 1.d50
       daycs = 1.d50
	 day_save_ss = day
         if (sday .ne. 0) day = sday
         if (shtl .ne. 0.0) shtl0 = head_tol
         if (stmch .ne. 0.0) stmch0 = tmch
         nstep_save = nstep
         nstep = snstep
         aiaa_save = aiaa
         aiaa = smult

      else if(iflg.eq.1) then
        divisor = 1.
        divisore = 1.
        i_pdiff = 0
        pmax_i = 0.
         if (svar_flag) then
            isteady=-1
         else
            isteady = 1
            return
         endif
         if(irdof.eq.13) then
            pdifmax = 0.0
            accmax = 0.0
            accdif_i = 0
            do i = 1, n
               if (spercent) then
                  if (pho(i) .ne. 0.) then
                     divisor = pho(i)
                  else if (phi(i) .ne. 0.) then
                     divisor = phi(i)
                  else
                     divisor = 1.
                  end if
               end if
               pdifmax = max(abs((phi(i)-pho(i))/divisor),pdifmax)
c 082319 gaz
               accdif_i = abs((deni(i)*sx1(i))*dtot)
               if(accdif_i.gt.accmax) then
                i_accdif = i
                accmax = accdif_i
               endif
            end do
            if(pdifmax.gt.toldp) go to 100
            if(accmax.gt.tacc) go to 100
         else if(ico2.lt.0) then
            pdifmax = 0.0
            sdifmax = 0.0
            accmax = 0.0
            amass = 0.0
            i_pdiff = 0
            i_sdiff = 0
            i_accdif = 0
            do i = 1, n
               if (spercent) then
                  if (pho(i) .ne. 0) then
                     divisor = pho(i)
                  else if (phi(i) .ne. 0) then
                     divisor = phi(i)
                  else
                     divisor = 1.
                  end if
               end if
               pdiff_i = abs((phi(i)-pho(i))/divisor)
               if(pdiff_i.gt.pdifmax) then
                i_pdiff = i
                pmax_i = phi(i)
                pmax_io = pho(i)
                pdifmax = pdiff_i
               endif
               if (spercent) then
                     divisor = 1.
               end if
               sdiff_i = abs((s(i)-so(i))/1.)
               if(sdiff_i.gt.sdifmax) then
                i_sdiff = i
                smax_i = s(i)
                smax_io = so(i)
                sdifmax = sdiff_i
               endif
c 082319 gaz
               accdif_i = abs((deni(i)*sx1(i))*dtot)
               if(accdif_i.gt.accmax) then
                i_accdif = i
                accmax = accdif_i
               endif
            amass = amass + denh(i) * volume(i)
            end do
            if(l.eq.1) then
             amass0 = am0
             amass_ch =(amass-amass0)
             amass0 = amass
            else
             amass_ch =(amass-amass0)
             amass0 = amass
            endif
            if(pdifmax.gt.toldp) go to 100
            if(sdifmax.gt.tolds) go to 100
            if(accmax.gt.tacc) go to 100
         else if(ico2.eq.0) then
            pdifmax = 0.0
            sdifmax = 0.0
            tdifmax = 0.0
            i_accdif = 0
            accmax = 0.0
            amass = 0
            do i = 1, n
               if (spercent) then
                  if (pho(i) .ne. 0.) then
                     divisor = pho(i)
                  else if (phi(i) .ne. 0.) then
                     divisor = phi(i)
                  else
                     divisor = 1.
                  end if
               end if
               pdifmax = max(abs((phi(i)-pho(i))/divisor),pdifmax)
               if (spercent) then
                     divisor = 1.
               end if
               sdifmax = max(abs((s(i)-so(i))/1.),sdifmax)
               if (spercent) then
                  if (to(i) .ne. 0.) then
                     divisor = to(i)
                  else if (t(i) .ne. 0.) then
                     divisor = t(i)
                  else
                     divisor = 1.
                  end if
               end if
               tdifmax = max(abs((t(i)-to(i))/divisor),tdifmax)
c 082319 gaz
               accdif_i = abs((deni(i)*sx1(i))*dtot)
               if(accdif_i.gt.accmax) then
                i_accdif = i
                accmax = accdif_i
               endif
              amass = amass + denh(i) * volume(i)
            end do
            if(l.eq.1) then
             amass0 = am0
             amass_ch =(amass-amass0)
             amass0 = amass
            else
             amass_ch =(amass-amass0)
             amass0 = amass
            endif
            if(pdifmax.gt.toldp) go to 100
            if(sdifmax.gt.tolds) go to 100
            if(tdifmax.gt.toldt) go to 100
            if(accmax.gt.tacc) go to 100
         else if(ico2.gt.0) then
            pdifmax = 0.0
            sdifmax = 0.0
            tdifmax = 0.0
            pcidifmax = 0.0
            accmax = 0.0
            do i = 1, n
               if (spercent) then
                  if (pho(i) .ne. 0.) then
                     divisor = pho(i)
                  else if (phi(i) .ne. 0.) then
                     divisor = phi(i)
                  else
                     divisor = 1.
                  end if
               end if
               pdifmax = max(abs((phi(i)-pho(i))/divisor),pdifmax)
               if (spercent) then
                     divisor = 1.
               end if
               sdifmax = max(abs((s(i)-so(i))/1.),sdifmax)
               if (spercent) then
                  if (to(i) .ne. 0.) then
                     divisor = to(i)
                  else if (t(i) .ne. 0.) then
                     divisor = t(i)
                  else
                     divisor = 1.
                  end if
               end if
               tdifmax = max(abs((t(i)-to(i))/divisor),tdifmax)
               if (spercent) then
                  if (pcio(i) .ne. 0.) then
                     divisor = pcio(i)
                  else if (pci(i) .ne. 0.) then
                     divisor = pci(i)
                  else
                     divisor = 1.
                  end if
               end if
               pcidifmax = max(abs((pci(i)-pcio(i))/divisor),pcidifmax)
               accmax = max(abs(deni(i)*sx1(i)),accmax)
            end do
            if(pdifmax.gt.toldp) go to 100
            if(sdifmax.gt.tolds) go to 100
            if(tdifmax.gt.toldt) go to 100
            if(pcidifmax.gt.toldc) go to 100
            if(accmax.gt.tacc) go to 100
         endif
         isteady=1
 100     continue
      else if(iflg.eq.2) then
         if(time_ss.gt.0.0) then
	    itt = 1
         else
	    itt = 0
         endif
         if(days.ge.time_ss.and.itt.eq.1) then
            isteady = 2
            btol_save = balance_tol
            balance_tol = tolerance
            tolde_save = tolde
            tolde = tolerance
            go to 20
         endif
         if(.not. sflux_flag) then
            if (isteady .eq. 1) then
               isteady = 2
               itt = 0
            end if
            go to 10
         endif

 20      continue
         if(ico2.lt.0) then
            if (spercent) then
               if (qtoti .ne. 0) then
                  divisor = abs(qtoti)
               else if (inflow .ne. 0.) then
                  divisor = abs(inflow)
               else
                  divisor = 1.
               end if
            else
                divisor = day*8.64e4
            end if
c gaz 082819 always go through flowrate adjustment of tmch
            flow_rate = (abs(inflow) - abs(qtoti))/divisor
            if(abs(flow_rate) .le. balance_tol) then
               if(isteady.eq.-1)then
                  isteady = 1
               else if(isteady.eq.1) then
                  isteady = 2
                  itt = 0
               endif
            endif

               if (shtl.ne.0) then
c gaz 082619
c                  ratio = max(1.0d0,abs(flow_rate/balance_tol))
c                  head_tol = min(ratio*shtl0,shtl)
                   ratio = 1.
                   head_tol = min(ratio*shtl0,shtl)
               endif
               if (stmch.ne.0) then
                  flow_max = abs(flow_rate)
c                  ratio = flow_max/balance_tol
                  balfac13 = 1.5*balance_tol
                  balfac12 = 2.*balance_tol
                  balfac1 = 4*balance_tol
                  balfac2 = 100.*balance_tol
                  balfac3 = 1000.*balance_tol
                  if(flow_max.gt.balfac3.or
     &                   .days.lt.tol_str) then
                   tmch1 = 100.*stmch0
          write(ierr,*) 'balfac3 flow_max ', flow_max,' tmch1 ',tmch1
                  else if(flow_max.gt.balfac2) then
                   tmch1_min = 50.*stmch0
                   tmch1_max = 100.*stmch0
                   tmch1 = ((tmch1_max-tmch1_min)/(balfac3 -balfac2))*
     &               flow_max + tmch1_min
          write(ierr,*) 'balfac2 flow_max ', flow_max,' tmch1 ',tmch1
                  else if(flow_max.gt.balfac1) then
                   tmch1_min = 10.*stmch0
                   tmch1_max = 50.*stmch0
                   tmch1 = ((tmch1_max-tmch1_min)/(balfac2 -balfac1))*
     &               flow_max + tmch1_min
          write(ierr,*) 'balfac2 flow_max ', flow_max,' tmch1 ',tmch1
                  else if(flow_max.gt.balfac12) then
                   tmch1_min = 5.*stmch0
                   tmch1_max = 10.*stmch0
                   tmch1 = ((tmch1_max-tmch1_min)/(balfac1 -balfac12))*
     &               flow_max + tmch1_min
          write(ierr,*) 'balfac12 flow_max ', flow_max,' tmch1 ',tmch1
                  else if(flow_max.gt.balfac13) then
                   tmch1_min = 1.5*stmch0
                   tmch1_max = 5.*stmch0
                   tmch1 = ((tmch1_max-tmch1_min)/(balfac12 -balfac13))*
     &               flow_max + tmch1_min
            write(ierr,*) 'balfac13 flow_max ', flow_max,' tmch1 ',tmch1
                  else if(flow_max.gt.balance_tol) then
                   tmch1_min = stmch0
                   tmch1_max = 1.5*stmch0
                   tmch1=((tmch1_max-tmch1_min)/(balfac13-balance_tol))*
     &               flow_max + tmch1_min
            write(ierr,*) 'bal_tol flow_max ', flow_max,' tmch1 ',tmch1
                  else
                   tmch1 = stmch0
                  endif
c gaz 080519 tolerance should decrease monotonically
                  tmch = min(tmch1,tmch_old)
             write(ierr,*) ' tmch ', tmch,' tmch_old ',tmch_old
                  tmch_old = tmch
               endif
c            endif
         else if(ico2.ge.0) then
            if (spercent) then
               if (qtoti .ne. 0) then
                  divisor = abs(qtoti)
               else if (inflow .ne. 0.) then
                  divisor = abs(inflow)
               else
                  divisor = 1.
               end if
               if (qtotei .ne. 0) then
                  divisore = abs(qtotei)
               else if (inflowe .ne. 0.) then
                  divisore = abs(inflowe)
               else
                  divisore = 1.
               end if
            else
               divisor = day*8.64e4
               divisore = day*8.64e4
            end if
            flow_rate = (abs(inflow) - abs(qtoti))/divisor
            enth_rate = (abs(inflowe) - abs(qtotei))/divisore
            if(abs(flow_rate) .le. balance_tol .and.
     &           abs(enth_rate) .le. tolde) then
               if(isteady.eq.-1)then
                  isteady = 1
               else if(isteady.eq.1) then
                  isteady = 2
                  itt = 0
               endif
            else
               if (stmch.ne.0) then
                  ratio = max(1.0d0,abs(flow_rate/balance_tol),
     &                 abs(enth_rate/tolde))
                  tmch1 = min(ratio*stmch0,stmch)
                  ratio = tims/days
                  tmch2 = min(ratio*stmch0,stmch)
                  tmch = min(tmch1,tmch2)
               end if

            endif
         end if

 10      continue
         if(isteady .eq. 2) then
            if (l .lt. sminstep) then
               isteady = 1
c               return
            end if
c     If we have reached steady state or reached our maximum
c     time limit, make sure we output information
            if(itt .eq. 0) then
               if (iout .ne. 0) write(iout,50) days, trim(info_string)
               if(iptty.ne.0) write(iptty,50) days, trim(info_string)
            else
               if (iout .ne. 0) then
                  write (iout, 52) snstep, time_ss, days
                  write (iout, 53) trim(info_string)
               end if
               if(iptty.ne.0) then
                  write (iptty, 52) snstep, time_ss, days
                  write (iptty, 53) trim(info_string)
               end if
               balance_tol = btol_save
               tolde = tolde_save
            endif
         end if

         call write_steady_state

         if (isteady .eq. 2) then
            ntty = 2

            if (isty .eq. 0) then
               ifinsh = 1
               isteady = 0
               istea_pest = 1
            else
               ifinsh = 2
               isteady = 0
               tims = tims_trans
               daycs=daycs_save
               daycf=daycf_save
               dayhf=dayhf_save
               dayhs=dayhs_save
               day = day_save_ss
               daymax = daymax_save
               nstep = nstep_save
               aiaa = aiaa_save
            end if
         end if
      else if(iflg.eq.3) then
         if(isteady .ne. 2) then
            if (iout .ne. 0) write(iout,53) trim(info_string)
            if (iptty .ne. 0) write(iptty,53) trim(info_string)

            call write_steady_state

            if (iout .ne. 0) write(iout,54) snstep, time_ss, days
            if (iptty.ne.0) write(iptty,54) snstep, time_ss, days
         endif
      endif

      return

 50   format(/,'>>>>> Steady State achieved at ', g16.9,
     &     ' days (', a, ') <<<<<')
 52   format(/,'>>>>> Steady state time exceeded <<<<<',/,1p,
     &     'Max timesteps: ', i6, ' Specified Time: ',
     &     g15.6,' Simulated Time: ',g15.6)
 53   format(/,'>>>>> Steady State NOT achieved with (', a, ') <<<<<')
 54   format(/,'>>>>> Time control information <<<<<',/,1p,
     &     'Max timesteps: ', i6, ' Specified Time: ',
     &     g15.6,' Simulated Time: ',g15.6, /)
 56   format('>>>>> Total Mass ',1x,1p,g15.6
     &     'Max timesteps: ', i6, ' Specified Time: ',
     &     g15.6,' Simulated Time: ',g15.6, /)

      contains

      subroutine write_steady_state

      if (toldh .ne. tolerance)
     &     hdifmax = pdifmax /( 997. * 9.8d-6)
      if (iout .ne. 0) then
         if (isteady  .ne. 2 .and. iflg .ne. 3)
     &        write (iout, 41) trim(info_string)
         if (toldp .ne. tolerance)
     &        write(iout, 40) 'pressure', pdifmax, toldp
         if (toldp .ne. tolerance. and. i_pdiff. ne. 0)
     &        write(iout, 43) i_pdiff, pmax_i, pmax_io,izonef(i_pdiff)
         if (toldh .ne. tolerance)
     &        write(iout, 40) 'head', hdifmax, toldh
         if (tolds .ne. tolerance)
     &        write(iout, 40) 'saturation', sdifmax, tolds
         if (tolds .ne. tolerance. and. i_sdiff. ne. 0)
     &        write(iout, 44) i_sdiff, smax_i, smax_io, izonef(i_sdiff)
         if (toldt .ne. tolerance)
     &        write(iout, 40) 'temperature', tdifmax, toldt
         if (tacc .ne. tolerance)
     &        write(iout, 40) 'nodal mass', accmax, tacc
         if (tacc .ne. tolerance. and. i_accdif. ne. 0)
     &    write(iout,45)i_accdif,accmax/(sx1(i_accdif)*denh(i_accdif)),
     &            amass_ch,izonef(i_accdif)
         if (toldc .ne. tolerance)
     &        write(iout, 40) 'air pressure', pcidifmax, toldc
         if (balance_tol .ne. tolerance) write(iout, 40)
     &        'flow balance (IN-OUT)', flow_rate, balance_tol
         if (tolde .ne. tolerance) write(iout, 40)
     &        'enthalpy balance (IN-OUT)', enth_rate, tolde
            write(iout,42) tmch
         if (isty .ne. 0) write(iout,51)
      end if
      if(iptty.ne.0) then
         if (isteady .ne. 2 .and. iflg .ne. 3)
     &        write (iptty, 41) trim(info_string)
         if (toldp .ne. tolerance)
     &        write(iptty, 40) 'pressure', pdifmax, toldp
         if (toldp .ne. tolerance. and. i_pdiff. ne. 0)
     &        write(iptty, 43) i_pdiff, pmax_i, pmax_io, izonef(i_pdiff)
         if (toldh .ne. tolerance)
     &        write(iptty, 40) 'head', hdifmax, toldh
         if (tolds .ne. tolerance)
     &        write(iptty, 40) 'saturation', sdifmax, tolds
         if (tolds .ne. tolerance. and. i_sdiff. ne. 0)
     &        write(iptty, 44) i_sdiff, smax_i, smax_io, izonef(i_sdiff)
         if (toldt .ne. tolerance)
     &        write(iptty, 40) 'temperature', tdifmax, toldt
         if (tacc .ne. tolerance)
     &        write(iptty, 40) 'nodal mass', accmax, tacc
         if (tacc .ne. tolerance. and. i_accdif. ne. 0)
     &    write(iptty,45)i_accdif,accmax/(sx1(i_accdif)*denh(i_accdif)),
     &            amass_ch,izonef(i_accdif)
         if (toldc .ne. tolerance)
     &        write(iptty, 40) 'air pressure', pcidifmax, toldc
         if (balance_tol .ne. tolerance) write(iptty, 40)
     &        'flow balance (IN-OUT)', flow_rate, balance_tol
         if (tolde .ne. tolerance) write(iptty, 40)
     &        'enthalpy balance (IN-OUT)', enth_rate, tolde
         write(iptty,42) tmch
         if (isty .ne. 0) write(iptty,51)
      endif

 40   format ('>>>>> Max change in ', a, 1pg12.4, ' tolerance',
     &     1pg12.4, ' <<<<<')
 43   format ('>>>>> Gridblock Max Pres change ',i8,
     &     ' Pres(new) ',1p,g12.4,' Pres(old) ',g12.4,
     &      ' Zone ',i6,' <<<<<')
 44   format ('>>>>> Gridblock Max Sat change ',i8,
     &     ' Sat(new) ',1p,g12.4,' Sat(old) ',g12.4,
     &     ' Zone ',i6,' <<<<<')
  45  format ('>>>>> Gridblock Max Mass change ',i8,' Frac change ',1p,
     &       g12.4,' Total Mass change ',g14.6,' Zone ',i6,' <<<<<')
 41   format (/, '>>>>> Steady state check (', a, ') <<<<<')
 42   format ('>>>>> Equation tolerance (',1p,g12.4, ') <<<<<')
 51   format(/,'>>>>> Starting transient with steady solution <<<<<',/)

      end subroutine write_steady_state

      end subroutine steady

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
Image(MSC_Sign.JPG, left, thumbnail, also known as Med Sci Cafe) The Scrubs MedSci Cafe caters mostly to the med/vetmed crowd. They sell Starbucks coffee, cookies, pastries, and some breakfast foods. They have a place that makes deli sandwiches, as well as a grill serving everything from burgers (beef, chicken, garden) to cheesesteaks to even grilled cheese sandwiches, as well as hot dogs and sausages. You could also get more unique, tasty things such as the pineapple chicken teriyaki sandwich, or some kind of taco and bean salad.  The TexMex grill serves Quesadillas, tacos, and other mexican food. They do have premade sandwiches as well as a variety of salad types.  The cafe is large, clean, has plenty of seating space, a couple of tvs, and vending machines as well for those in a rush.  It also has MOObile Net wireless internet to those with a UCD login.

20070313 01:01:44 nbsp I love the staff.  They are extraordinairily friendly and helpful.  I just wish it were open on Sundays around noon so I could get a cup of coffee and feed my caffeine addiction. Users/KirstenOnell

20080602 12:45:57 nbsp I have mixed feelings about this place. I like it b/c its convienent, since theres nothing else close by. However the prices for coffee and the sandwhiches irk me. The venti iced coffee drinks are not standard Starbucks venti size the coffee is 2030 cents more than Starbucks. Users/MyaBrn

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
       subroutine autnls
c
c   this prepares for and calls the canned subroutine lmdif1
c   which will solve the unconstrained non-linear least squares
c   fitting problem. lmdif1 uses a levenberg-marquardt algorithm,
c   and requires an external subprogram to evaluate the function
c   to minimize. the subroutine autfun is used to evaluate this
c   function.
c   the function allows the ordinates of the breakpoints of
c   the spline to vary such that when the optimal b-spline (as
c   discussed in de boor) is put through these breakpoints, and
c   the chi(r) is found using this spline as the background is
c   optimized at low-r.
c
c   copyright 1992  university of washington :          matt newville
c------------------------------------------------------------------
       include 'autobk.h'
c------------------------------------------------------------------
       integer   lenwrk, loop, lminfo, isplin, iup, ilo
       integer   iwork(mtknot), ifft, nemin, nemax, nofx, ne0
       integer   i, nr, ndoc, iend, im, istrln, ntest
       parameter(lenwrk = 2*maxpts*mtknot + 20*mtknot + 2*maxpts   )
       character*128  messg , doc(maxdoc),  type*10,  file*40
       double precision work(lenwrk),  fvect(2*maxpts), varys(mtknot)
       double precision esplin(mtknot), bscoef(mtknot)
       double precision etmp(maxpts), small, anorm(3), toler, tolfac
       double precision sumsqr
       double precision e0, delq, efromq, qmn, qmx, qknot
       double precision e0tst, qmntst, qmxtst, edfmin
       parameter (small  = 1.d-10, tolfac = 1.d-4, edfmin= 100.d0)
       external autfun, nofx, istrln, sumsqr
c-----------------------------------------------------------------------
c----initialization for fitting:
c   -the starting value for the fitting tolerance (1e-5) is empirical.
c    the user will be allowed to change it, if necessary, with usertl.
c    (right now usertl is set to 1).
c   -the fitting is done in single precision. round off is a worry, but
c    i have not been able to get the variables to change enough with the
c    double precision version of lmdif1. this should still be explored.
c-----------------------------------------------------------------------
       toler  = usrtol * tolfac
       e0     = ee + e0shft
       ifft  =  1
       loop   = 0
 100   continue
       lminfo = 0
       loop   = loop + 1
       if (loop.gt.5) then
          call echo('   autobk warning: completed 5'//
     $         ' fitting loops without a stable result.')
          call echo('                   something may be wrong'//
     $         ' with the fitting,')
          call echo('                   and the background may'//
     $         ' not be reliable.')
          go to 800
       end if
c----start fitting:
       do 120 i = 1, nxmu
          spline(i) = xmudat(i)
 120   continue
       e0    = ee  + e0shft
c   move emin and emax if e0-shifted
       emin  = emin + e0shft
       emax  = emax + e0shft
       if (emax.le.emin) emax = energy(nxmu)
       nemin = nofx(emin,energy,nxmu)
       nemax = nofx(emax,energy,nxmu)
       emin  = energy(nemin)
       emax  = energy(nemax)
c  evaluate energies for the spline variables:
c  - emin and emax are not relative to the edge
c  - energy contains the input energy values, not relative to the edge.
       qmin   = sqrt(etok* abs(emin - e0 ) )
       if (e0.le.emin)    qmin   = zero
       qmax   = sqrt(etok* (emax - e0 ) )
c
c drpair is the spacing between pairs of independent points in r-space
       if (.not.gvknot)
     $      nsplin = 2 * int ( rbkg * ( qmax - qmin ) / pi )  +  1
       nsplin = min(mtknot-5, max(5,nsplin))
c
c  initialize energy values through which the first guess for the
c  spline must go (evenly spaced in q), and get the initial value
c  for the spline value at this point. we'll also get the initial
c  guesses for the variables (the b-spline coefficients) from this
c  initial spline.
c
       delq = (qmax-qmin)/(nsplin - 1)
       do 300 i = 1, nsplin
           efromq    = e0 + ( qmin + (i-1) * delq) **2 / etok
           isplin    = nofx(efromq,energy,nxmu)
           esplin(i) = energy(isplin)
           iup       = min(nxmu, isplin + 5)
           ilo       = max(1,    isplin - 5)
           bscoef(i) = (2*spline(isplin) + spline(iup) + spline(ilo))/4
300    continue
       esplin(nsplin) = one + esplin(nsplin)
c the first and last korder knots in the b-spline are nearly degenerate at
c the endpoints. spstep sets the spacing between these points.  the
c default is one -- this may help eliminate "spikes" at the endpoints.
c since each knot represents a place where a derivative can break,
c  having all four of these at one place allows a complete break from
c at this point. by moving a few of the knots just off the ends, the
c spline is a little bit stiffer at the endpoints.
       do 310 i = 1, korder
           eknot(i)        = esplin(1)      - spstep * (korder-i-1)
           eknot(nsplin+i) = esplin(nsplin) + spstep * i
310    continue
       qmn       = sqrt( etok* abs(esplin(1) - e0) )
       if (e0.lt.esplin(1))  qmn = zero
       qmx = sqrt( etok* (esplin(nsplin) - e0) )
       do 320 i = korder+1, nsplin
           qknot    = (i-korder)*(qmx - qmn)/(nsplin-korder+1)
           eknot(i) = esplin(1) + qknot**2/etok
320    continue

c
c    determine the knots for the spline:
c    knots are points at which the spline has extra freedom.
cc       ntknot = nsplin + korder
       if ( (korder.lt.3).or.(nsplin.lt.korder) ) then
          call echo('           autobk error: not enough'//
     $                      ' freedom to create spline.')
          call echo('                   change fitting range, or'//
     $                      ' order of spline')
          call fstop(' not enough data to make spline')
       end if
c  the b-spline coefficients will be the variables in the fit,
c  the above estimates for the elements of bscoef are good enough
c  especially when the spline values are smoothed for the initial
c  guesses in the fit, as below.
       varys(1) = (3*bscoef(1) + bscoef(2))/ 4
       varys(nsplin) = (3*bscoef(nsplin) + bscoef(nsplin-1))/ 4
       do 380 i = 2, nsplin-1
          varys(i) =  (bscoef(i-1) + 2*bscoef(i) + bscoef(i+1) )/4
 380   continue
       nvarys = nsplin
c set up window function
       call window(winstr,windo1,windo2,qmin,qmax,qgrid,maxpts,windo)
c  if a theory file is used, do its fft now
c  and add another variable if e0 is to be shifted
       if (theory) then
          if (eevary) then
             nvarys        = nvarys + 1
             varys(nvarys) = e0shft
          end if
          nr     = nrpts
cc          print*, ' theory : ', mftfit
          call fitfft(thiq, maxpts, mftfit, wfftc, qgrid,
     $         windo, qweigh, spline, one, ifft, zero, r1st,
     $         nrpts, thifit)

          thessq = max(small, sumsqr(thifit(nrbkg),nr1st))
          if (nr.ne.nrpts)  then
             call echo('  autnls error: fitfft is broken' )
             write(messg,'(a,2i5)') 'nr, nrpts = ', nr, nrpts
             call echo('          '//messg(:40))
             call fstop(' call matt!! autobk fitfft is broken')
          end if
       end if
c
c  initialize fvect and work arrays for lmdif1
       mfit  = nrbkg
       do 400 i =1, mfit
          fvect(i) = zero
 400   continue
       do 410 i =1, lenwrk
          work(i) = zero
 410   continue
       do 420 i =1, mtknot
          iwork(i) = 0
 420   continue
c
c  lmdif1 to do levenberg-marquardt nonlinear least squares
       if (eevary) then
          call echo('   autobk: fitting background and e0'//
     $         ' over the low-r range')
       else
          call echo('   autobk: fitting background'//
     $         ' over the low-r range')
       end if
       if (iprint.ge.2) then
          iend = 0
cc          print*, ' call autfun: ', mfit, nvarys, iend
cc          print*, varys(1), fvect(1), varys(3), fvect(3)
          call autfun(mfit, nvarys, varys, fvect, iend)
          type   = 'xmu'
          file   = 'spline0.dat'
          doc(1) = ' autobk: initial spline v. e : before fitting'
 500      format (1x,a,i4)
 501      format (1x,a)
 502      format (1x,i5.2,e14.6,e14.6,1x,i5.2,e14.6,e14.6)
 503      format (1x,i5.2,e14.6,e14.6)
 504      format (1x,a,f14.6)
          write (doc(2), 500) ' number of variables = ', nsplin
          write (doc(3), 504) ' spline_step (spstep)= ', spstep
          write (doc(4), 501) ' index, esplin, ysplin, '//
     $                        ' index, esplin, ysplin '
          ndoc   = 4
          do 550 i = 1, nsplin, 2
             if (i.lt.maxdoc)  ndoc = ndoc + 1
             if ((i.le.maxdoc).and.(i.eq.nsplin)) then
                write (doc(ndoc), 503) i , esplin(i)  , varys(i)
             elseif (i.le.maxdoc) then
                write (doc(ndoc), 502) i  , esplin(i)  , varys(i)  ,
     $                                i+1, esplin(i+1), varys(i+1)
             end if
 550      continue
          call outcol(type, file, asccmt, ndoc, mdocxx, doc, nxmu,
     $         energy, spline, xmudat, chiq, thiq)
       end if

c
c the real fit!
       call lmdif1 (autfun, mfit, nvarys, varys, fvect,
     $      toler, lminfo,  iwork,  work, lenwrk)

       call autfun(mfit, nvarys, varys, fvect, iend)
       if ( (lminfo.ge.1).and.(lminfo.le.3) ) then
          messg = 'fitting is finished.'
          im = max(1, istrln(messg))
          call echo ('                   '//messg(:im) )
       else
          messg = 'lmdif finished with an error ! '
          im = max(1, istrln(messg))
          call echo ('                   '//messg(:im) )
          write(messg, '( a, i5)' ) 'error code lminfo = ',lminfo
          im = max(1, istrln(messg))
          call echo ('                   '//messg(:im) )
          call echo ('            call matt')
       end if
c
       if (iprint.ge.2) then
          file   = 'spline.dat'
          doc(1) = ' autobk: final spline v. e : before fitting'
          ndoc   = 4
          do 580 i = 1, nsplin, 2
             if (i.lt.maxdoc)  ndoc = ndoc + 1
             if ((i.le.maxdoc).and.(i.eq.nsplin)) then
                write (doc(ndoc), 503) i , esplin(i)  , varys(i)
             elseif (i.le.maxdoc) then
                write (doc(ndoc), 502) i  , esplin(i)  , varys(i)  ,
     $                                i+1, esplin(i+1), varys(i+1)
             end if
 580      continue
          call outcol(type, file, asccmt, ndoc, mdocxx, doc, nxmu,
     $         energy, spline, xmudat, chiq, thiq)
       end if

c  even though the fit was good, if an energy shift was done, the
c  q values used to determine the number and location of the spline
c  points were wrong, and so the fit must be re-done. this will be
c  repeated up to 5 times, or until e0 is stable to within 0.5 ev.
       e0tst   = ee + e0shft
       qmntst  = sqrt(etok* abs(emin - e0tst ) )
       qmxtst  = sqrt(etok* abs(emax - e0tst ) )
       if (e0tst.le.emin)    qmntst   = zero
       ntest   = 2 * int ( rbkg * abs(qmxtst - qmntst) / pi)  + 1
       if (gvknot) ntest = nsplin
       ntest   = min(mtknot-5, max(5,ntest))
       if (eevary.and.(ntest.ne.nsplin) ) then
          call echo('   autobk warning: because the'//
     $                ' energy origin was shifted in ')
          call echo('                   in the fit, kmin'//
     $                ' and kmax have been redefined')
          call echo('                   so that the number'//
     $                ' of knots in the spline has')
          call echo('                   changed. the fit '//
     $                ' will be done again.')
          go to 100
       end if
c------------------------------------------------------------------
c  finished with fit
 800   continue
c
c  now that we have the spline as good as we can get it, let's get an
c  improved estimate of the edge-step. we should have it pretty close
c  to one, but now we do a parabolic extrapolation of the background
c  spline to the edge energy value. this should cover up most
c  difficulties in getting the edge step from a linear extrapolation
c  of the absorption data itself.
       if (stfind.and.(.not.funnrm)) then
          enor1  = e0 + enor1
          enor2  = e0 + enor2
          if (enor2.gt.energy(nxmu))   enor2 = energy(nxmu)
          if (enor1.gt.energy(nxmu))   enor1 = enor2 /2
          nnorm = 3
          if (abs(enor2 - enor1).le.edfmin) nnorm = 2
          call polyft(enor1, enor2, energy, spline,nxmu,nnorm,cnorm)
          ne0  = nofx(e0, energy, nxmu)
          step = cnorm(1) + cnorm(2)*energy(ne0)
          if (nnorm.eq.3) step = step + cnorm(3)*energy(ne0)**2
          cnorm(1) = cnorm(1) + bpre
          cnorm(2) = cnorm(2) + slopre
       end if
c
c  now redo function evaluation with the final values for everything.
       do 930 i = 1, nxmu
          spline(i) = xmudat(i)
 930   continue
       do 940 i = 1, maxpts
          chiq(i) = zero
 940   continue
       final = .true.
       call autfun(mfit, nvarys, varys, fvect, iend)
c
       return
c end subroutine autnls
       end

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
CThe subroutines in this file were taken directly from RIP code written
C by Dr. Mark Stoelinga.  They were modified by Sherrie
C Fredrick(NCAR/MMM) to work with NCL February 2015.
C NCLFORTSTART
      subroutine wrf_monotonic(out,in,lvprs,cor,idir,delta,
     &                         ew,ns,nz,icorsw)
       implicit none
       integer idir,ew,ns,nz,icorsw
       double precision delta
       double precision in(ew,ns,nz),out(ew,ns,nz)
       double precision lvprs(ew,ns,nz),cor(ew,ns)
C NCLEND

       integer i,j,k,ripk,k300

       do j=1,ns
       do i=1,ew
          if (icorsw.eq.1.and.cor(i,j).lt.0.) then
             do k=1,nz
                in(i,j,k)=-in(i,j,k)
              enddo
          endif

c
c   First find k index that is at or below (height-wise) the 300 hPa
c      level.
c
      do k = 1,nz
         ripk =  nz-k+1
         if (lvprs(i,j,k) .le. 300.d0) then
            k300=k
            goto 40
         endif
      enddo
c
 40   continue

       do k = k300, 1,-1
         if (idir.eq.1) then
            out(i,j,k)=min(in(i,j,k),in(i,j,k+1)+delta)
         elseif (idir.eq.-1) then
            out(i,j,k)=max(in(i,j,k),in(i,j,k+1)-delta)
         endif
       enddo

       do k = k300+1, nz
         if (idir.eq.1) then
            out(i,j,k)=max(in(i,j,k),in(i,j,k-1)-delta)
         elseif (idir.eq.-1) then
            out(i,j,k)=min(in(i,j,k),in(i,j,k-1)+delta)
         endif
       enddo

      end do
      end do

      return
      end

c--------------------------------------------------------------------

C NCLFORTSTART
      FUNCTION wrf_intrp_value (wvalp0,wvalp1,vlev,vcp0,vcp1,icase)
      implicit none

      integer icase
      double precision wvalp0,wvalp1,vlev,vcp0,vcp1
C NCLEND
      double precision valp0,valp1,rvalue,rgas,ussalr,sclht

      double precision wrf_intrp_value,chkdiff

      rgas    = 287.04d0     !J/K/kg
      ussalr  = 0.0065d0      ! deg C per m
      sclht   = rgas*256.d0/9.81d0

      valp0 = wvalp0
      valp1 = wvalp1
      if ( icase .eq. 2) then  !GHT
           valp0=exp(-wvalp0/sclht)
           valp1=exp(-wvalp1/sclht)
      end if

      chkdiff = vcp1 - vcp0
      if(chkdiff .eq. 0) then
         print *,"bad difference in vcp's"
         stop
      end if

      rvalue = (vlev-vcp0)*(valp1-valp0)/(vcp1-vcp0)+valp0
      if (icase .eq. 2) then  !GHT
          wrf_intrp_value = -sclht*log(rvalue)
      else
          wrf_intrp_value = rvalue
      endif

      return
      end
c------------------------------------------------------------
C NOTES:
c      vcarray is the array holding the values for the vertical
c      coordinate.
c              It will always come in with the dimensions of
c              the staggered U and V grid.
C NCLFORTSTART

      subroutine  wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,
     &                       sfp,smsfp,vcarray,interp_levels,numlevels,
     &                       icase,ew,ns,nz,extrap,vcor,logp,rmsg)

      implicit none
      integer   ew,ns,nz,icase,extrap
      integer   vcor,numlevels,logp
      double precision   datain(ew,ns,nz),pres(ew,ns,nz),tk(ew,ns,nz)
      double precision   ght(ew,ns,nz)
      double precision   terrain(ew,ns),sfp(ew,ns),smsfp(ew,ns)
      double precision   dataout(ew,ns,numlevels),qvp(ew,ns,nz)
      double precision   vcarray(ew,ns,nz)
      double precision   interp_levels(numlevels),rmsg
C NCLEND
       integer   njx,niy,nreqlvs,ripk
       integer   i,j,k,itriv,kupper
       integer   ifound,miy,mjx,isign
       double precision      rlevel,vlev,diff
       double precision      tempout(ew,ns),tmpvlev
       double precision      vcp1,vcp0,valp0,valp1
       double precision      rgas,rgasmd,sclht,ussalr,cvc,eps
       double precision      qvlhsl,ttlhsl,vclhsl,vctophsl
       double precision      wrf_intrp_value
       double precision      plhsl,zlhsl,ezlhsl,tlhsl,psurf,pratio,tlev
       double precision      ezsurf,psurfsm,zsurf,qvapor,vt
       double precision      rconst,expon,exponi
       double precision      ezlev,plev,zlev,ptarget,dpmin,dp
       double precision      pbot,zbot,tbotextrap,e
       double precision      tlclc1,tlclc2,tlclc3,tlclc4
       double precision      thtecon1,thtecon2,thtecon3
       double precision      tlcl,gamma,cp,cpmd,gammamd,gammam
       character cvcord*1

       rgas    = 287.04d0     !J/K/kg
       rgasmd  = .608d0
       ussalr  = .0065d0      ! deg C per m
       sclht   = rgas*256.d0/9.81d0
       eps     = 0.622d0
       rconst  = -9.81d0/(rgas * ussalr)
       expon   =  rgas*ussalr/9.81d0
       exponi  =  1./expon
       tlclc1   = 2840.d0
       tlclc2   = 3.5d0
       tlclc3   = 4.805d0
       tlclc4   = 55.d0
       thtecon1 = 3376.d0 ! K
       thtecon2 = 2.54d0
       thtecon3 = 0.81d0
       cp       = 1004.d0
       cpmd     = 0.887d0
       gamma    = rgas/cp
       gammamd  = rgasmd-cpmd

       if(vcor .eq. 1) then
          cvcord = 'p'
       else if((vcor .eq. 2) .or. (vcor .eq. 3)) then
          cvcord = 'z'
       else if((vcor .eq. 4) .or. (vcor .eq. 5)) then
          cvcord = 't'
       end if

       miy = ns
       mjx = ew
       njx = ew
       niy = ns

       do j = 1,mjx
          do i = 1,miy
               tempout(j,i) = rmsg
          end do
       end do

      do nreqlvs = 1,numlevels
         if(cvcord .eq. 'z') then
!Convert rlevel to meters from km

            rlevel = interp_levels(nreqlvs) * 1000.d0
            vlev = exp(-rlevel/sclht)
         else if(cvcord .eq. 'p') then
            vlev = interp_levels(nreqlvs)
         else if(cvcord .eq. 't') then
            vlev = interp_levels(nreqlvs)
         end if

         do j=1,mjx
         do i=1,miy
cGet the interpolated value that is within the model domain
            ifound = 0
            do k = 1,nz-1
               ripk  = nz-k+1
               vcp1  = vcarray(j,i,ripk-1)
               vcp0  = vcarray(j,i,ripk)
               valp0 = datain(j,i,ripk)
               valp1 = datain(j,i,ripk-1)
               if ((vlev.ge.vcp0.and.vlev.le.vcp1) .or.
     &            (vlev.le.vcp0.and.vlev.ge.vcp1)) then
c                  print *,i,j,valp0,valp1
                  if((valp0 .eq. rmsg).or.(valp1 .eq. rmsg)) then
                     tempout(j,i) = rmsg
                     ifound=1
                  else
                     if(logp .eq. 1) then
                       vcp1  = log(vcp1)
                       vcp0  = log(vcp0)
                       if(vlev .eq. 0.0d0) then
                         print *,"Pressure value = 0"
                         print *,"Unable to take log of 0"
                         stop
                       end if
                       tmpvlev  = log(vlev)
                     else
                       tmpvlev = vlev
                     end if
                     tempout(j,i) = wrf_intrp_value(valp0,valp1,
     &                                  tmpvlev,vcp0,vcp1,icase)
c                     print *,"one ",i,j,tempout(j,i)
                     ifound=1
                  end if
                  goto 115
               end if
             end do !end for the k loop
 115  continue

      if (ifound.eq.1) then !Grid point is in the model domain
          goto 333
      end if

cIf the user has requested no extrapolatin then just assign
call values above or below the model level to rmsg.
      if(extrap .eq. 0) then
         tempout(j,i) = rmsg
         goto 333
      end if

c The grid point is either above or below the model domain
c
c First we will check to see if the grid point is above the
c model domain.
      vclhsl   = vcarray(j,i,1) !lowest model level
      vctophsl = vcarray(j,i,nz)!highest model level
      diff     = vctophsl-vclhsl
      isign    = nint(diff/abs(diff))
C
      if(isign*vlev.ge.isign*vctophsl) then
C Assign the highest model level to the out array
         tempout(j,i)=datain(j,i,nz)
C         print *,"at warn",j,i,tempout(j,i)
         goto 333
      endif

c
c   Only remaining possibility is that the specified level is below
c   lowest model level.  If lowest model level value is missing,
c   set interpolated value to missing.
c
      if (datain(i,j,1) .eq. rmsg) then
          tempout(j,i) = rmsg
          goto 333
      endif

c
c   If the field comming in is not a pressure,temperature or height
C   field we can set the output array to the value at the lowest
c   model level.
c
      tempout(j,i) = datain(j,i,1)
c
c   For the special cases of pressure on height levels or height on
c   pressure levels, or temperature-related variables on pressure or
c   height levels, perform a special extrapolation based on
c   US Standard Atmosphere.  Here we calcualate the surface pressure
c   with the altimeter equation.  This is how RIP calculates the
c   surface pressure.
c
       if (icase.gt.0) then
           plhsl  = pres(j,i,1) * 0.01d0  !pressure at lowest model level
           zlhsl  = ght(j,i,1)            !grid point height a lowest model level
           ezlhsl = exp(-zlhsl/sclht)
           tlhsl  = tk(j,i,1)             !temperature in K at lowest model level
           zsurf  = terrain(j,i)
           qvapor = max((qvp(j,i,1)*.001d0),1.e-15)
c virtual temperature
c          vt     = tlhsl * (eps + qvapor)/(eps*(1.0 + qvapor))
c           psurf  = plhsl * (vt/(vt+ussalr * (zlhsl-zsurf)))**rconst
           psurf    = sfp(j,i)
           psurfsm  = smsfp(j,i)
           ezsurf   = exp(-zsurf/sclht)

cThe if for checking above ground
           if ((cvcord.eq.'z'.and.vlev.lt.ezsurf).or.
     &         (cvcord.eq.'p'.and.vlev.lt.psurf)) then
c
c      We are below the lowest data level but above the ground.
c      Use linear interpolation (linear in prs and exp-height).
c
               if (cvcord.eq.'p') then
                   plev=vlev
                   ezlev=((plev-plhsl)*ezsurf+(psurf-plev)*ezlhsl)/
     &                     (psurf-plhsl)
                   zlev=-sclht*log(ezlev)
                   if (icase .eq. 2) then
                       tempout(j,i)=zlev
                       goto 333
                   endif

                elseif (cvcord.eq.'z') then
                   ezlev=vlev
                   zlev=-sclht*log(ezlev)
                   plev=((ezlev-ezlhsl)*psurf+(ezsurf-ezlev)*plhsl)/
     &                    (ezsurf-ezlhsl)
                   if (icase .eq. 1) then
                      tempout(j,i)=plev
                      goto 333
                  endif
                endif

            else   !else for checking above ground
                ptarget=psurfsm-150.d0
                dpmin=1.e4
                do k=1,nz
                   ripk = nz-k+1
                   dp=abs((pres(j,i,ripk) * 0.01d0)-ptarget)
                   if (dp.gt.dpmin) goto 334
                   dpmin=min(dpmin,dp)
                enddo
 334            kupper=k-1

                ripk       = nz - kupper + 1
                pbot       = max(plhsl,psurf)
                zbot       = min(zlhsl,zsurf)
                pratio     = pbot/(pres(j,i,ripk) * 0.01d0)
                tbotextrap = tk(j,i,ripk)*(pratio)**expon
c virtual temperature
                vt = tbotextrap * (eps + qvapor)/(eps*(1.0d0+qvapor))
                if (cvcord.eq.'p') then
                   plev=vlev
                   zlev=zbot+vt/ussalr*(1.-(vlev/pbot)**expon)
                   if(icase .eq. 2) then
                      tempout(j,i)=zlev
                      goto 333
                   endif
                elseif (cvcord.eq.'z') then
                        zlev=-sclht*log(vlev)
                        plev=pbot*(1.+ussalr/vt*(zbot-zlev))**exponi
                        if (icase .eq. 1) then
                            tempout(j,i)=plev
                            goto 333
                        endif
                 endif
            end if !end if for checking above ground
       end if !for icase gt 0

       if(icase .gt. 2) then !extrapolation for temperature
           tlev=tlhsl+(zlhsl-zlev)*ussalr
           qvapor = max(qvp(j,i,1),1.e-15)
           gammam = gamma*(1.+gammamd*qvapor)
           if(icase .eq. 3) then
              tempout(j,i) = tlev - 273.16d0
           else if(icase .eq. 4) then
              tempout(j,i) = tlev
C Potential temperature - theta
           else if (icase. eq. 5) then
               tempout(j,i)=tlev*(1000.d0/plev)**gammam
C extraolation for equivalent potential temperature
           else if (icase .eq. 6) then
               e    = qvapor*plev/(eps+qvapor)
               tlcl = tlclc1/(log(tlev**tlclc2/e)-tlclc3)+tlclc4
               tempout(j,i)=tlev*(1000.d0/plev)**(gammam)*
     &         exp((thtecon1/tlcl-thtecon2)*qvapor*
     &             (1.+thtecon3*qvapor))
           end if
       end if

 333  continue

         end do
         end do
 !        print *,"----done----",interp_levels(nreqlvs)
          do i = 1,njx
             do j = 1,niy
                dataout(i,j,nreqlvs) = tempout(i,j)
           end do
          end do
       end do !end for the nreqlvs
       return
       end  !wrf_vinterp

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl