text:
 SUBROUTINE rubber1d(pr,lb,cmd,model,bet,elas)

 ! calculates stiffness and stresses for rubber and elastomeric
 ! foam materials

 IMPLICIT NONE
 ! DUMMY arguments
 LOGICAL, INTENT(IN) :: cmd ! .FALSE.: compute stress at mechanical strain
                            !  .TRUE.: compute stress and tangent modulus at mechanical strain
 REAL(kind=8), INTENT(IN):: pr(12), & !(12) material constants
                            lb        ! Principal stretch
 INTEGER(kind=4), INTENT(IN) :: model   !Model
                                !1:      Arruda-Boyce
                                !2:      Mooney-Rivlin
                                !3:      Neo-Hooke
                                !4:      Ogden (N=1)
                                !5:      Ogden (N=2)
                                !6:      Ogden (N=3)
                                !7:      Polynomial (N=1)
                                !8:      Polynomial (N=2)
                                !9:      Polynomial (N=3)
                                !10:     Reduced Polynomial (N=1)
                                !11:     Reduced Polynomial (N=2)
                                !12:     Reduced Polynomial (N=3)
                                !13:     Van der Waals (not implemented yet)
                                !14:     Yeoh
                                !15:     Hyperfoam (N=1)
                                !16:     Hyperfoam (N=2)
                                !17:     Hyperfoam (N=3)
 REAL(kind=8), INTENT(OUT) ::      &
      elas,             & !  elasticity tangent modulus
      bet                 !  If CMD Kirchhoff stress tau  Else Second Piola-Kirchhoff

 LOGICAL :: ci1, ci2     !  use first and second invariants

 INTEGER(kind=4) nterm,     & ! number of terms (sum) in the model
                 kode,      & ! model
                 i,j,k        ! index and pointers
 INTEGER(kind=4), PARAMETER :: ii(3) = (/ 1,3,6 /)

 REAL(kind=8) i1,i2,a1,a2,a3,a4,a5,lt,d1,d2,dd1,dd2,lb2  ! auxiliar terms

 !initializes number of terms
 ci1 = .TRUE.       !compute I1 invariant derivatives
 ci2 = .FALSE.      !do NOT compute I2 derivatives
 kode        = model
 SELECT CASE (kode)
 CASE (2)  !Mooney-Rivlin
    ci2 = .TRUE.      !compute I2 derivatives
 CASE (4:6)  !Ogden
    nterm = kode - 3 !number of terms to be considered
 ci1 = .FALSE.      !compute I1 invariant derivatives
 CASE (7:9)  !Polynomial
    nterm = kode - 6 !number of terms to be considered
    ci2 = .TRUE.      !compute I2 derivatives
 CASE (10:12)!reduced polynomial
    nterm = kode - 9 !number of terms to be considered
 CASE (14)   !yeoh
    kode = 12        !reduced polynomial
    nterm = 3        !number of terms to be considered
 CASE (15:17)           !NOTE that isochoric approach is considered
    CALL runen2('Hyperfoam not implemented for 1-D problems')
    !nterm = kode - 14 !number of terms to be considered
 END SELECT

 ! compute stresses (plane stress isochoric)
 lt = 1d0/SQRT(lb)  !isochoric
 lb2 = lb*lb         !lb squared

 IF( ci1 )THEN
   !   first derivatives of I1
   d1 = 2d0*(lb-1d0/lb2)      !I1,L
   i1 = lb2+2d0/lb            !I1
   IF( ci2 )THEN
     ! first derivatives of I2
     d2 = 2d0*(1d0-1d0/lb/lb2)     !I2,L
     i2 = 2d0*lb+1d0/lb2        !I2
   END IF
 END IF

 SELECT CASE (kode)
 CASE (1) !Arruda-Boyce
   ! pr(5) = Mu/5/lm^2,  pr(6) = Mu 11/350/lm^4,  pr(7) =Mu 19/1750/lm^6, ! pr(8) =Mu 519/134750/lm^8
   a1 = pr(4)+i1*(pr(5)+i1*(pr(6)+i1*(pr(7)+i1*pr(8)))) !Phi,I1
   bet = a1*d1

 CASE (2) !Mooney-Rivlin

   bet = pr(1)*d1 + pr(2)*d2

 CASE (3) !Neo-Hookean
   bet = pr(1)*d1

 CASE (4:6) !ogden
   k = 1                                   !pointer to properties array
   bet = 0d0                               !Initializes Kirchhoff stress
   DO j=1,nterm                            !for each term
     a1 = pr(k+1)                          !exponent alpha
     bet = bet + pr(k)*(lb**a1 - lt**a1)
     k = k+2
   END DO
   bet = bet/lb

 CASE ( 7: 9) !polynomial model
   i1 = i1 - 3d0   !I1 - 3
   i2 = i2 - 3d0   !I2 - 3
   a1 = pr(1)                  !initializes Phi,I1
   a2 = pr(2)                  !initializes Phi,I2
   IF( nterm > 1 )THEN
     a1 = a1+2d0*pr(3)*i1 +     pr(4)*i2
     a2 = a2+    pr(4)*i1 + 2d0*pr(5)*i2
   END IF
   IF( nterm > 2 )THEN
     a1 = a1+3d0*pr(6)*i1*i1 +2d0*pr(7)*i1*i2+    pr(8)*i2*i2
     a2 = a2+    pr(7)*i1*i1 +    pr(8)*i1*i2+3d0*pr(9)*i2*i2
   END IF

   bet = a1*d1 + a2*d2

 CASE (10:12) !Reduced polynomial
   i1 = i1 - 3d0        !I1 - 3
   a1 = 0d0             !Initializes Phi,I1
   DO i=1,nterm
     a1 = a1 + pr(ii(i))*i*i1**(i-1)    !Phi,I1
   END DO
   bet = a1*d1

 END SELECT

 IF( cmd ) THEN  !to compute stiffness

   IF( ci1 )THEN
     ! Second Derivatives of First Invariant
     dd1= (2d0+4d0/lb2/lb)         !I1,L,L
     IF( ci2 )THEN
       !  Second Derivatives of Second Invariant
       dd2= 6d0/lb2**2             !I2,L,L
     END IF
   END IF

   SELECT CASE (kode)

   CASE (1) !Arruda-Boyce
     a2 = pr(5)+i1*(2d0*pr(6)+i1*(3d0*pr(7)+i1*4d0*pr(8))) !Phi,I1,i2
     elas = a1*dd1 + a2*d1**2

   CASE (2) !Mooney-Rivlin
     elas = pr(1)*dd1 + pr(2)*dd2

   CASE (3) !Neo-Hookean
     elas = pr(1)*dd1

   CASE (4:6) !ogden
     elas = 0d0                              !initializes
     k = 1
     DO j=1,nterm                            !for each term
       a1 = pr(k+1)                          !exponent
       a2 = a1-1
       a3 = a1/2d0+1d0
       elas = elas + pr(k)*(a2*lb**a1 + a3*lt**a1 ) !DS1/L1
       k = k+2
     END DO                                      !DS1/L1
     elas = elas/lb2

   CASE ( 7: 9) !polynomial model
     a3 = 0           !a1,I1
     a4 = 0           !a1,I2 = a2,I1
     a5 = 0           !a2,I2
     IF( nterm > 1 )THEN
       a3 = 2d0*pr(3)    !a1,I1
       a4 =     pr(4)    !a1,I2
       a5 = 2d0*pr(5)    !a2,I1
     END IF
     IF( nterm > 2 )THEN
       a3 = a3+6d0*pr(6)*i1 +2d0*pr(7)*i2  !Phi,I1,I1
       a4 = a4+2d0*pr(7)*i1 +2d0*pr(8)*i2  !Phi,I1,I2
       a5 = a5+2d0*pr(8)*i1 +6d0*pr(9)*i2  !Phi,I2,I2
     END IF

     elas = a1*dd1 +a2*dd2 +a3*d1*d1 +2d0*a4*d1*d2 +a5*d2*d2

   CASE (10:12) !Reduced polynomial
     a2 = 0d0        !initializes Phi,I1,I1
     IF( nterm > 1 ) a2 = pr(3)*2d0            !Phi,I1,I1
     IF( nterm > 3 ) a2 = a2 + pr(6)*6d0*i1    !Phi,I1,I1
     elas = a1*dd1 + a2*d1**2

   END SELECT

   ! Kirchhoff stress
   bet = bet*lb
   elas  = elas * lb2 + bet             !Dt/Ln2

 ELSE
   !second Piola Kirchhoff
   bet = bet/lb

 END IF

 RETURN
 END SUBROUTINE rubber1d

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
       SUBROUTINE TSDRDATE()
       INTEGER*4 SYEAR,HR,HOUR,I,YEAR
       INTEGER*4 FRAC_SEC,SFDUT1,SFDUT2,SS_1950,SECOND,MSEC
       INTEGER*2 SDR(3),SDR2(3)
       CHARACTER*100 MSG

       CALL XVMESSAGE('Gross check through entire century',' ')
       CALL XVMESSAGE(' ',' ')
       DO 25 I=0,2
        IF (I.EQ.0) HR=0
        IF (I.EQ.1) HR=1
        IF ((I.NE.0).AND.(I.NE.1)) HR=8759
        CALL XVMESSAGE(' ',' ')
        WRITE(MSG,1000) HR
        CALL XVMESSAGE(MSG,' ')
 1000      FORMAT('Hour of year = ', I)
        CALL XVMESSAGE(' ',' ')
        CALL XVMESSAGE('ss_sec,fsec,year,hour,sec,msec',' ')
        CALL XVMESSAGE(' ',' ')

        DO 50 YEAR=0,99
         SDR(1)=HR+24
         SDR(2)=3500
         SDR(3)=890
         CALL SFDUDATE(SFDUT1,SFDUT2,SDR,YEAR)
         CALL SDRDATE(SFDUT1,SFDUT2,SDR2,SYEAR)
         SS_1950 = SFDUT1
         FRAC_SEC = SFDUT2
         HOUR = SDR(1)
         SECOND = SDR(2)
         MSEC = SDR(3)
         WRITE(MSG,1500)SS_1950,FRAC_SEC,YEAR,HOUR,SECOND,MSEC
         CALL XVMESSAGE(MSG,' ')
 1500          FORMAT(6I12)
         IF (YEAR.NE.SYEAR.OR.HOUR.NE.SDR2(1).OR.
     c   SECOND.NE.SDR2(2).OR.MSEC.NE.SDR2(3)) THEN
          HOUR = SDR2(1)
          SECOND = SDR2(2)
          MSEC = SDR2(3)
          WRITE(MSG,1500)SS_1950,FRAC_SEC,SYEAR,HOUR,SECOND,MSEC
          CALL XVMESSAGE(MSG,' ')
         ENDIF
   50  CONTINUE
   25  CONTINUE

       CALL XVMESSAGE('Check through each day of VGR Uranus period',' ')
       CALL XVMESSAGE(' ',' ')
       DO 125 YEAR=85,88
       WRITE(MSG,2000)YEAR
       CALL XVMESSAGE(' ',' ')
       CALL XVMESSAGE(MSG,' ')
       CALL XVMESSAGE(' ',' ')
 2000  FORMAT('Year = ', I)
       DO 150 HR=0,8783,24
       SDR(1)=HR+24
       SDR(2)=3500
       SDR(3)=890
       IF (YEAR.NE.88.AND.HR.EQ.8760) GO TO 125
       CALL SFDUDATE(SFDUT1,SFDUT2,SDR,YEAR)
       CALL SDRDATE(SFDUT1,SFDUT2,SDR2,SYEAR)
       SS_1950 = SFDUT1
       FRAC_SEC = SFDUT2
       HOUR = SDR(1)
       SECOND = SDR(2)
       MSEC = SDR(3)
       WRITE(MSG,1500)SS_1950,FRAC_SEC,YEAR,HOUR,SECOND,MSEC
       CALL XVMESSAGE(MSG,' ')
       IF (YEAR.NE.SYEAR.OR.HOUR.NE.SDR2(1).OR.
     c  SECOND.NE.SDR2(2).OR.MSEC.NE.SDR2(3)) THEN
         HOUR = SDR2(1)
         SECOND = SDR2(2)
         MSEC = SDR2(3)
         WRITE(MSG,1500)SS_1950,FRAC_SEC,SYEAR,HOUR,SECOND,MSEC
         CALL XVMESSAGE(MSG,' ')
       ENDIF
  150  CONTINUE
  125  CONTINUE
       END

ext: f

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


program iprogram
    use utilities
    use percolation
    implicit none

    integer :: num_samples, i, j, spanning_label, num_Ls, num_clusters, fileunit
    logical, dimension(:,:), allocatable :: matrix
    integer, dimension(:,:), allocatable :: label_matrix
    integer, dimension(:), allocatable :: Ls
    real(kind=dp) :: D, const, tmpM
    real(kind=dp), dimension(:), allocatable :: Ms_approx, Ms

    num_samples = 200
    Ls = 2**[(i, i=4, 11)]
    num_Ls = size(Ls)
    allocate(Ms(num_Ls))
    Ms = 0

    do i = 1, num_Ls
        associate(L => Ls(i), M => Ms(i))
            tmpM = 0
            !$omp  parallel do &
            !$omp& private(matrix, label_matrix, num_clusters, spanning_label) &
            !$omp& reduction(+:tmpM)
            do j = 1, num_samples
                matrix = create_binary_matrix(pc, L)
                call label(matrix, label_matrix, num_clusters)
                spanning_label = find_spanning_cluster(label_matrix, num_clusters)
                tmpM = tmpM + count(label_matrix == spanning_label)
            end do
            !$omp end parallel do
            M = tmpM
        end associate
    end do

    Ms = Ms/num_samples

    call linfit(log(1.0d0*Ls), log(Ms), D, const)
    Ms_approx = exp(const) * Ls**D

    open(newunit=fileunit, file="tmp/i.dat", status="replace")
    do i = 1, num_Ls
        write(unit=fileunit, fmt=*) Ls(i), Ms(i), Ms_approx(i)
    end do
    close(unit=fileunit)

    open(newunit=fileunit, file="tmp/i_D.dat", status="replace")
    write(unit=fileunit, fmt="(f0.3)") D
    close(unit=fileunit)

end program iprogram


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:34:48 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE MKSTRESS__genmod
          INTERFACE
            SUBROUTINE MKSTRESS
            END SUBROUTINE MKSTRESS
          END INTERFACE
        END MODULE MKSTRESS__genmod


SUBROUTINE LA_TEST_DPTSV( N, NRHS, D, E, B, LDB, INFO )
!
!  -- LAPACK95 interface driver routine (version 1.1) --
!     UNI-C, Denmark,
!     April 7, 1999
!
!  .. Use Statements ..
   USE LA_PRECISION, ONLY: WP => DP
   USE F95_LAPACK, ONLY: LA_PTSV
!  .. Implicit Statement ..
   IMPLICIT NONE
!  .. Scalar Arguments ..
   INTEGER, INTENT(IN) :: N, NRHS, LDB
   INTEGER, INTENT(INOUT) :: INFO
!  .. Array Arguments ..
      REAL(WP), INTENT(INOUT) :: D(1:N), E(1:N-1), B(1:LDB,1:N)
!  .. Parameters ..
   CHARACTER(LEN=8),  PARAMETER :: SRNAME = 'LA_PTSV '
   CHARACTER(LEN=14), PARAMETER :: SRNAMT = 'LA_TEST_DPTSV '
!  .. Common blocks ..
   INTEGER :: INFOTC
   COMMON /LINFO95/ INFOTC
!  .. Local Scalars ..
   INTEGER :: I, J, ID, IE, IB1, IB2
!  .. Local Arrays ..
   LOGICAL, SAVE :: CTEST = .TRUE., ETEST = .TRUE.
!  .. Executable Statements ..
   ID = N; IE = N-1; IB1 = N; IB2 = NRHS
   I = INFO / 100; J = INFO - I*100
   SELECT CASE(I)
   CASE(0)
      IF (NRHS == 1) THEN
        CALL LA_PTSV( D(1:ID), E(1:IE),  B(1:IB1,1), INFO )
        INFO = INFOTC
      ELSE
        CALL LA_PTSV( D(1:ID),  E(1:IE),  B(1:IB1,1:IB2), INFO )
        INFO = INFOTC
      END IF
   CASE (2)
      IE = ID - 2
   CASE(3)
      IB1 = ID - 1
   CASE(:-1, 1, 4:)
      CALL UESTOP(SRNAMT)
   END SELECT
      IF( I /= 0 ) THEN
      SELECT CASE (NRHS)
      CASE (2:)
         CALL LA_PTSV( D(1:ID), E(1 :IE), B(1:IB1,1:IB2),  INFO )
      CASE(1)
         CALL LA_PTSV( D(1:ID), E(1:IE), B(1:IB1,1), INFO )
      CASE(:-1)
         CALL UESTOP(SRNAMT)
      END SELECT
   END IF
   CALL LA_AUX_AA01( I, CTEST, ETEST, SRNAMT )
END SUBROUTINE LA_TEST_DPTSV


text:
! a simple code to check the analytic Jacobian via numerical
! differencing

subroutine test_jacobian() bind(C)

  use network
  use eos_module
  use burner_module
  use actual_burner_module
  use actual_rhs_module
  use burn_type_module
  use numerical_jac_module

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  type (burn_t) :: state_ana, state_num

  type (eos_t) :: eos_state

  character (len=32) :: probin_file
  integer :: probin_pass(32)
  integer :: i, j

  probin_file = "probin"
  do i = 1, len(trim(probin_file))
     probin_pass(i) = ichar(probin_file(i:i))
  enddo

  call runtime_init(probin_pass(1:len(trim(probin_file))), len(trim(probin_file)))

  call network_init()
  call actual_rhs_init()
  call burner_init()
  call eos_init()

  state_ana % rho   = 2.0e7_rt
  state_ana % T     = 8.0e9_rt

  state_ana % xn(:) = ONE / nspec

  call test_numerical_jac(state_ana)

end subroutine test_jacobian

max_stars_repo_name: doreenfan/Microphysics

max_issues_repo_name: doreenfan/Microphysics

max_forks_repo_name: doreenfan/Microphysics


text:
module twinOne
    implicit none
    use twinTwo
    write (*,*) <caret>x
endmodule

max_stars_repo_name: LChernigovskaya/fortran-plugin

max_issues_repo_name: LChernigovskaya/fortran-plugin

max_forks_repo_name: LChernigovskaya/fortran-plugin


text:
      SUBROUTINE DTSYS (ITASK,DATEA,FSEC,DPDTTM,ERR,MESSAG)
      IMPLICIT NONE

*     FORMAL_PARAMETERS:
      INTEGER ITASK, DATEA(6)
      LOGICAL ERR
      REAL FSEC
      DOUBLE PRECISION DPDTTM
      CHARACTER*(*) MESSAG

**    Local variables
      INTEGER I1, NZERO, DAYMAX, NDAYS
      INTEGER IMOTB1(12),IMOTB2(12)
      DOUBLE PRECISION DTMP, DTMAX
      LOGICAL DTLEAP
      INTEGER INDAY, ISECS
      DOUBLE PRECISION FRDAY, DSECS, FSECS

*     The absolute day number belonging to a certain day is
*     the number of days passed since 1-jan-0001. Hence, the
*     absolute day number of 1-jan-0001 is 0. The absolute
*     day number belonging to 1-jan-0002 is 365.
*     The number of days N calculated for a certain year Y is
*     given by N = Y*365 + Y/4 - Y/100 + Y/400 - 1.
*     The 1 subtracted simplifies the use of DOY numbers.
*     For 1-jan-0002 for instance DOY=1. If this number is added
*     to the N = 364 calculated for the passed year, year 1, we
*     have again the absolute day number belonging to 1-jan-0002.
*     The offset value used in the TTUTIL day counting should
*     give the number 1 to 1-jan-1900, which implies:
*     ttutiloffset = 1899 * 365 + 1899/4 - 1899/100 + 1899/400 - 1.
*     Dates below 1-jan-0001 and above 31-dec-9999 are illegal.

      INTEGER OFFSET
      PARAMETER (OFFSET = 693594)
      INTEGER N400, N100, N4, N1, LEFT, LASTY, DLAST
      SAVE

      DATA IMOTB1 /0,31,59,90,120,151,181,212,243,273,304,334/
      DATA IMOTB2 /31,28,31,30,31,30,31,31,30,31,30,31/

      ERR    = .FALSE.
      MESSAG = ' '

*     value corresponding to 31-dec-9999 24:00
      DTMAX = DBLE (9999*365 + 9999/4 - 9999/100 + 9999/400 - OFFSET)

      IF (ITASK.EQ.1) THEN

*        convert from array to double precision real
         NZERO  = 0
         DO 10 I1=1,3
            IF (DATEA(I1).EQ.0) NZERO = NZERO+1
10       CONTINUE

         IF (NZERO.GT.0.AND.NZERO.LT.3) THEN
            ERR    = .TRUE.
            MESSAG = 'year 0 or partial date not allowed'
            RETURN
         END IF

*        days for the year
         IF (DATEA(1).GE.1) THEN
*           absolute day number belonging to 31 december previous year
            LASTY = DATEA(1) - 1
            INDAY = LASTY*365 + LASTY/4 - LASTY/100 + LASTY/400 - 1
         ELSE IF (DATEA(1).EQ.0) THEN
*           set date to TTUTIL offset value
            INDAY = OFFSET
         ELSE
            ERR    = .TRUE.
            MESSAG = 'illegal year'
            RETURN
         END IF

*        count days in month
         IF (DATEA(2).NE.0) THEN
            IF (DATEA(2).GE.1.AND.DATEA(2).LE.12) THEN
               IF (DATEA(2).GE.2) THEN
                  INDAY = INDAY+IMOTB1(DATEA(2))
                  IF (DTLEAP(DATEA(1)).AND.DATEA(2).GT.2)
     &                INDAY = INDAY + 1
               END IF
            ELSE
               ERR    = .TRUE.
               MESSAG = 'month invalid'
               RETURN
            END IF
         END IF

*        count days in day
         IF (DATEA(3).NE.0) THEN
            DAYMAX = IMOTB2(DATEA(2))
            IF (DTLEAP (DATEA(1)).AND.DATEA(2).EQ.2) DAYMAX = DAYMAX+1

            IF (DATEA(3).GE.1.AND.DATEA(3).LE.DAYMAX) THEN
               INDAY = INDAY + DATEA(3)
            ELSE
               ERR    = .TRUE.
               MESSAG = 'day invalid'
               RETURN
            END IF
         END IF

*        subtract TTUTIL offset value and convert
         DTMP = DBLE(INDAY-OFFSET)

         IF (DATEA(4).GE.0) THEN
*           important: hours can be more than 24 !!
            DTMP = DTMP+DBLE (DATEA(4))/24.D0
         ELSE
            ERR    = .TRUE.
            MESSAG = 'hour invalid'
            RETURN
         END IF

         IF (DATEA(5).GE.0.AND.DATEA(5).LE.59) THEN
            DTMP = DTMP+DBLE (DATEA(5))/1440.D0
         ELSE
            ERR    = .TRUE.
            MESSAG = 'minute invalid'
            RETURN
         END IF

         IF (DATEA(6).GE.0.AND.DATEA(6).LE.59) THEN
            DTMP = DTMP+DBLE (DATEA(6))/86400.D0
         ELSE
            ERR    = .TRUE.
            MESSAG = 'seconds invalid'
            RETURN
         END IF

         IF (FSEC.GE.0.AND.FSEC.LT.1.) THEN
            DTMP = DTMP+DBLE (FSEC)/86400.D0
         ELSE
            ERR    = .TRUE.
            MESSAG = 'fractional seconds invalid'
            RETURN
         END IF

         IF (DTMP .GE. DTMAX) THEN
            ERR    = .TRUE.
            MESSAG = 'year above 9999'
            RETURN
         END IF

         DPDTTM = DTMP

      ELSE IF (ITASK.EQ.2) THEN

*        convert from double precision real to date array and
*        fractional seconds

         IF (DPDTTM .GE. DTMAX) THEN
            ERR    = .TRUE.
            MESSAG = 'year above 9999'
            RETURN
         END IF

*        get integer and fractional day
         IF (DPDTTM .GE. 0.0D0) THEN
            INDAY = INT(DPDTTM)
            FRDAY = DPDTTM - DBLE(INDAY)
         ELSE
            INDAY = -INT(ABS(DPDTTM))
            FRDAY = ABS(DPDTTM) - DBLE(ABS(INDAY))
            IF (FRDAY .GT. 0.0D0) THEN
*              there is a fractional day ; get it right
               FRDAY = 1.0D0 - FRDAY
               INDAY = INDAY - 1
            END IF
         END IF

*        get fractional day in seconds
         DSECS = 86400.D0 * FRDAY
         ISECS = INT(DSECS)
         FSECS = DSECS - DBLE(ISECS)

*        if the fractional seconds are close to 1.0 ?
         IF (FSECS .GT. 0.999999D0) THEN
*           regard as 1.0
            FSECS = 0.0D0
            ISECS = ISECS + 1
            IF (ISECS.EQ.86400) THEN
*              regard as a day
               ISECS = 0
               INDAY = INDAY + 1
            END IF
         END IF

*        now we have integer days and integer seconds for the day
*        corrected for a nearly finished second
*        now the date calculation can be done in integer arithmetic

*        alogorithm from FSE 4 calendar routine
         INDAY = INDAY + OFFSET

         IF (INDAY.LT.0) THEN
            ERR    = .TRUE.
            MESSAG = 'date before 1-jan-0001'
            RETURN
         END IF

*        number of 400 year periods and number of days left
         N400 = INDAY / 146097
         LEFT = INDAY - N400 * 146097

*        number of 100 year periods (max 3) and the number of days left
         N100 = MIN (3, LEFT / 36524)
         LEFT = LEFT - N100 * 36524

*        number of 4 year periods and the number of days left
         N4   = LEFT / 1461
         LEFT = LEFT - N4 * 1461

*        number of 1 year periods (max 3) and the number of days left
         N1   = MIN (3, LEFT / 365)
         LEFT = LEFT - N1 * 365

*        current year number is all years passed plus 1
         DATEA(1) = n400*400 + n100*100 + n4*4 + n1 + 1

*        find day number in year (DOY) by subtracting last years
         LASTY = DATEA(1) - 1
         DLAST = LASTY*365 + LASTY/4 - LASTY/100 + LASTY/400 - 1
         INDAY = INDAY - DLAST

*        find out month number
         DATEA(2) = 1
         NDAYS = IMOTB2(DATEA(2))
         IF (DTLEAP (DATEA(1)).AND.DATEA(2).EQ.2) NDAYS = NDAYS+1

20       IF (INDAY.GT.NDAYS) THEN
            INDAY    = INDAY-NDAYS
            DATEA(2) = DATEA(2)+1
            NDAYS    = IMOTB2(DATEA(2))
            IF (DTLEAP (DATEA(1)).AND.DATEA(2).EQ.2) NDAYS = NDAYS+1
         GOTO 20
         END IF

*        day number in month
         DATEA(3) = INDAY

*        hh:mm:ss
         DATEA(4) = MOD(ISECS,86400) / 3600
         DATEA(5) = MOD(ISECS,3600) / 60
         DATEA(6) = MOD(ISECS,60)

*        fractional seconds
         FSEC = REAL (FSECS)

      ELSE
         CALL FATALERR ('DTSYS','unknown task')
      END IF

      RETURN
      END

ext: for


text:
*DECK DSRCMA
      SUBROUTINE DSRCMA (RSAV, ISAV, JOB)
C-----------------------------------------------------------------------
C This routine saves or restores (depending on JOB) the contents of
C the Common blocks DLS001, DLSA01, which are used
C internally by one or more ODEPACK solvers.
C
C RSAV = real array of length 240 or more.
C ISAV = integer array of length 46 or more.
C JOB  = flag indicating to save or restore the Common blocks:
C        JOB  = 1 if Common is to be saved (written to RSAV/ISAV)
C        JOB  = 2 if Common is to be restored (read from RSAV/ISAV)
C        A call with JOB = 2 presumes a prior call with JOB = 1.
C-----------------------------------------------------------------------
      INTEGER ISAV, JOB
      INTEGER ILS, ILSA
      INTEGER I, LENRLS, LENILS, LENRLA, LENILA
      DOUBLE PRECISION RSAV
      DOUBLE PRECISION RLS, RLSA
      DIMENSION RSAV(*), ISAV(*)
      SAVE LENRLS, LENILS, LENRLA, LENILA
      COMMON /DLS001/ RLS(218), ILS(37)
      COMMON /DLSA01/ RLSA(22), ILSA(9)
      DATA LENRLS/218/, LENILS/37/, LENRLA/22/, LENILA/9/
C
      IF (JOB .EQ. 2) GO TO 100
      DO 10 I = 1,LENRLS
 10     RSAV(I) = RLS(I)
      DO 15 I = 1,LENRLA
 15     RSAV(LENRLS+I) = RLSA(I)
C
      DO 20 I = 1,LENILS
 20     ISAV(I) = ILS(I)
      DO 25 I = 1,LENILA
 25     ISAV(LENILS+I) = ILSA(I)
C
      RETURN
C
 100  CONTINUE
      DO 110 I = 1,LENRLS
 110     RLS(I) = RSAV(I)
      DO 115 I = 1,LENRLA
 115     RLSA(I) = RSAV(LENRLS+I)
C
      DO 120 I = 1,LENILS
 120     ILS(I) = ISAV(I)
      DO 125 I = 1,LENILA
 125     ILSA(I) = ISAV(LENILS+I)
C
      RETURN
C----------------------- End of Subroutine DSRCMA ----------------------
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
      MODULE WW1_I
      INTERFACE
!
      SUBROUTINE WW1(IK,BK,ID,BD,K2,QM1,QM2,QM3,QM4,WW)
      USE vast_kind_param, ONLY:  DOUBLE
      INTEGER,      INTENT(IN)               :: K2
      INTEGER,      INTENT(IN), DIMENSION(7) :: IK, ID
      REAL(DOUBLE), INTENT(IN)               :: QM1, QM2, QM3, QM4
      REAL(DOUBLE), INTENT(IN), DIMENSION(3) :: BK, BD
      REAL(DOUBLE), INTENT(OUT)              :: WW
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! 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 and documentation.

! 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.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE PRINT_ORDER

! Writes abcissa's and weights from subroutines ORDER_GAUSS and ORDER_TRIA used in isoparametric element matrix generation subr's

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, MAX_ORDER_GAUSS, MAX_ORDER_TRIA, NUM_TRIA_ORDERS, TRIA_ORDER_NUMS

      USE PRINT_ORDER_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'PRINT_ORDER'
      CHARACTER(18*BYTE)              :: SSC, HHC              ! Char var used in writing to internal file for formatting output
      CHARACTER(18*BYTE)              :: SSC_I, SSC_J          ! Char var used in writing to internal file for formatting output

      INTEGER(LONG)                   :: J, KORDER             ! DO loop indices
      INTEGER(LONG)                   :: TORDER                ! Trianular integration order number

      REAL(DOUBLE)                    :: SSS(MAX_ORDER_GAUSS)  ! Gauss abscissa's
      REAL(DOUBLE)                    :: SS_I(MAX_ORDER_TRIA)  ! Tria X abscissa's
      REAL(DOUBLE)                    :: SS_J(MAX_ORDER_TRIA)  ! Tria Y abscissa's
      REAL(DOUBLE)                    :: HHH(MAX_ORDER_GAUSS)  ! Weights

! **********************************************************************************************************************************
! Write header

      WRITE(F06,1101)

! Write Gaussian orders

      WRITE(F06,1102)
      WRITE(F06,*)

      DO KORDER=1,MAX_ORDER_GAUSS
         CALL ORDER_GAUSS ( KORDER, SSS, HHH )
         WRITE(F06,1103) KORDER
         DO J=1,KORDER
            WRITE(SSC,1104) SSS(J)
            WRITE(HHC,1104) HHH(J)
            WRITE(F06,1105) J, SSC(1:8), SSC(9:13), SSC(14:18), HHC(1:8), HHC(9:13), HHC(14:18)
         ENDDO
         WRITE(F06,*)
         WRITE(F06,*)
      ENDDO

! Write Triangular orders

      WRITE(F06,1112)
      WRITE(F06,*)

      DO KORDER=1,NUM_TRIA_ORDERS
         TORDER = TRIA_ORDER_NUMS(KORDER)
         IF (TORDER <= MAX_ORDER_TRIA) THEN
            CALL ORDER_TRIA ( TORDER, SS_I, SS_J, HHH )
            WRITE(F06,1113) TORDER
            DO J=1,TORDER
               WRITE(SSC_I,1114) SS_I(J)
               WRITE(SSC_J,1114) SS_J(J)
               WRITE(HHC,1114) HHH(J)
               WRITE(F06,1115) J, SSC_I(1:8), SSC_I(9:13), SSC_I(14:18), SSC_J(1:8), SSC_J(9:13), SSC_J(14:18),                    &
                               HHC(1:8), HHC(9:13), HHC(14:18)
            ENDDO
            WRITE(F06,*)
            WRITE(F06,*)
         ELSE
            WRITE(F06,2001) TORDER, MAX_ORDER_TRIA
            EXIT
         ENDIF
      ENDDO

! Write trailer

      WRITE(F06,1106)

! **********************************************************************************************************************************
 1101 FORMAT(' __________________________________________________________________________________________________________________',&
             '_________________'                                                                                ,//,&
             ' ::::::::::::::::::::::::::::::::::::::::::START DEBUG( 5) OUTPUT FROM SUBROUTINE ORDER::::::::::::::::::::::::::::',&
              ':::::::::::::::::',/)

 1102 FORMAT(32X,'Abcissa and Weight Coefficients of the Gaussian Quadrature Formula')

 1103 FORMAT(65X,'n =',I3,/,65x,'------',//,41X,'J            SSS(J)                   HHH(J)')

 1104 FORMAT(F18.15)

 1105 FORMAT(39X, I3, 3X, A8,' ',A5,' ',A5,5X,A8,' ',A5,' ',A5)

 1112 FORMAT(32X,'Abcissa and Weight Coefficients of the Triangulur Quadrature Formula')

 1113 FORMAT(65X,'n =',I3,/,65X,'------',//,29X,'J          SS_I(J)                  SS_J(J)                   HHH(J)')

 1114 FORMAT(F18.15)

 1115 FORMAT(27X, I3, 3X, A8,' ',A5,' ',A5,5X,A8,' ',A5,' ',A5,5X,A8,' ',A5,' ',A5)

 1106 FORMAT(' :::::::::::::::::::::::::::::::::::::::::::END DEBUG( 5) OUTPUT FROM SUBROUTINE ORDER:::::::::::::::::::::::::::::',&
              ':::::::::::::::::'                                                                                ,/,&
             ' __________________________________________________________________________________________________________________',&
             '_________________',/)

 2001 FORMAT(' *INFORMATION: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' TRIANGULAR INTEGRATION ORDER ',I3,' IS GREATER THAM MAX_ORDER_TRIA = ',I3                             &
                    ,/,14X,' CANNOT COMPLETE PRINTOUT OF TRIANGULAR INTEGRATION ORDERS. THIS IS NOT A FATAL ERROR')

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

      END SUBROUTINE PRINT_ORDER


text:
      subroutine ptrim(word,latt,twiss,imon,emon,nmon,lfno)
      use tfstk
      use ffs
c      trim data greater than cut*sigma
      use tffitcode
      implicit real*8(a-h,o-z)
      integer*8 latt(nlat)
      dimension twiss(nlat,-ndim:ndim,ntwissfun),
     1          imon(nmona,4),emon(nmona,4),rms(2)
      external pack
      logical exist
      character*(*) word
      include 'inc/common.inc'
c
      call getwdl(word)
      if(word.eq.'S') then
      elseif(word.eq.'O') then
      elseif(word.eq.'M') then
        cut=getva(exist)
        if(.not.exist) then
          call getwdl(word)
          return
        endif
        ia=italoc(2*nmon)
        ia1=ia-1
        ia2=ia1+nmon
        if(simulate) then
          do 10 i=1,nmon
            j=imon(i,2)
            nq=imon(j,4)
            rlist(ia1+i)=twiss(imon(j,1),0,mfitdx)
     1           -twiss(imon(j,1),ndim,mfitdx)-rlist(latt(nq)+5)
     1           -emon(j,1)
            rlist(ia2+i)=twiss(imon(j,1),0,mfitdy)
     1           -twiss(imon(j,1),ndim,mfitdy)-rlist(latt(nq)+6)
     1           -emon(j,2)
 10       continue
        else
          do 11 i=1,nmon
            j=imon(i,2)
            rlist(ia1+i)=twiss(imon(j,1),0,mfitdx)
            rlist(ia2+i)=twiss(imon(j,1),0,mfitdy)
 11     continue
        endif
        rms(1)=0d0
        rms(2)=0d0
        nm=0
        do 20 i=1,nmon
          nm=nm+1
          rms(1)=rms(1)+rlist(ia1+i)**2
          rms(2)=rms(2)+rlist(ia2+i)**2
   20   continue
        rms(1)=sqrt(rms(1)/max(nm,1))
        rms(2)=sqrt(rms(2)/max(nm,1))
        fcut1=cut*rms(1)
        fcut2=cut*rms(2)
        do 30 i=1,nmon
          j=imon(i,2)
          if(abs(rlist(ia1+i)).gt.fcut1.or.
     1       abs(rlist(ia2+i)).gt.fcut2) then
            imon(j,3)=1
          endif
   30   continue
        call pack(imon(1,2),imon(1,3),nmon,nmon)
        write(lfno,'(2(A,I4))') '  BPM available :',nmon,' in ',nmonact
        nmonact=nmon
        call tfree(int8(ia))
        call getwdl(word)
      else
        return
      endif
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


subroutine compute_fact10(u,v,t,q,ps,prsi1,prsi2,skint,z0rl,islimsk,f10m)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    compute_fact10                compute 10m wind factor
!   prgmmr: treadon          org: np23                date: 2006-09-28
!
! abstract: Use GFS surface physics routines to compute 10m wind factor
!
! program history log:
!   2006-09-28 treadon - initial routine
!   2008-06-05 safford - rm unused vars and uses, comment out unused params
!
!   input argument list:
!      u       - u wind component (2d field, 1st model layer)
!      v       - v wind component
!      t       - sensible temperature
!      q       - specific humidity
!      ps      - surface pressure
!      prsi1   - surface pressure pressure
!      prsi2   - interface pressure at top of first layer
!      skint   - skin temperature
!      z0rl    - surface roughness
!      islimsk - land/sea/ice mask
!
!   output argument list:
!     f10m     - 10m wind factor
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use constants, only: grav,zero,quarter,half,one,two,four,five,&
       fv,rd,rd_over_cp,r1000
  implicit none

! Passed Variables
  real(r_kind)   ,intent(in   ) :: u,v,t,q,ps,skint,z0rl
  integer(i_kind),intent(in   ) :: islimsk
  real(r_kind)   ,intent(  out) :: f10m
  real(r_kind)   ,intent(in   ) :: prsi1,prsi2

! Local Variables
  real(r_kind):: prsl,prkl
  real(r_kind):: prki1,prki2
  real(r_kind):: q0,tem,del,rkap,rkapi,rkapp1
  real(r_kind)::ustar,wind,z0,z0max,ztmax,tv1,z1, &
       rat,thv1,theta1,tvs,dtv,rb,fm,fh,hlinf, &
       hl1,pm,ph,pm10,hl12,ph2,fm10
  real(r_kind):: psurf,ps1
  real(r_kind) restar,aa0,bb0,fhs,fms,hl0,hlt,adtv,bb,aa,hl0inf,hltinf,&
       hl110,olinf

! Local Parameters
  real(r_kind),parameter::  charnok=0.014_r_kind
! real(r_kind),parameter::  ca=0.4_r_kind
  real(r_kind),parameter::  alpha=five
  real(r_kind),parameter::  a0=-3.975_r_kind
  real(r_kind),parameter::  a1=12.32_r_kind
  real(r_kind),parameter::  b1=-7.755_r_kind
  real(r_kind),parameter::  b2=6.041_r_kind
  real(r_kind),parameter::  a0p=-7.941_r_kind
  real(r_kind),parameter::  a1p=24.75_r_kind
  real(r_kind),parameter::  b1p=-8.705_r_kind
  real(r_kind),parameter::  b2p=7.899_r_kind
  real(r_kind),parameter::  vis=1.4e-5_r_kind
! real(r_kind),parameter::  aa1=-1.076_r_kind
! real(r_kind),parameter::  bb1=0.7045_r_kind
! real(r_kind),parameter::  cc1=-0.05808_r_kind
! real(r_kind),parameter::  bb2=-0.1954_r_kind
! real(r_kind),parameter::  cc2=0.009999_r_kind
! real(r_kind),parameter::  rnu=1.51e-5_r_kind
! real(r_kind),parameter::  arnu=0.135_r_kind*rnu
  real(r_kind),parameter::  ten=10.0_r_kind

!  INITIALIZE VARIABLES. ALL UNITS ARE SUPPOSEDLY M.K.S. UNLESS SPECIFIED
!  PSURF IS IN PASCALS
!  WIND IS WIND SPEED, THETA1 IS ADIABATIC SURFACE TEMP FROM LEVEL 1
!  SURFACE ROUGHNESS LENGTH IS CONVERTED TO M FROM CM

   rkap = rd_over_cp
   RKAPI  = one / RKAP
   RKAPP1 = one + RKAP

   rat=zero ; restar=zero ; ustar=zero
   del = prsi1-prsi2
   tem   = rkapp1*del
   prki1 = (prsi1*0.01_r_kind)**rkap
   prki2 = (prsi2*0.01_r_kind)**rkap
   prkl  = (prki1*prsi1-prki2*prsi2)/tem
   prsl  = 100.0_r_kind*prkl**rkapi

   psurf=ps*r1000
   ps1=prsl*r1000
   wind= sqrt( u*u + v*v )
   wind=max(wind,one)
   q0=max(q,1.e-8_r_kind)
   theta1=t*(prki1/prkl)
   tv1 =t*(one+fv*q0)
   thv1=theta1 *(one+fv*q0)
   tvs =skint*(one+fv*q0)
   z0=0.01_r_kind*z0rl
   z1=-rd*tv1*log(ps1/psurf)/grav

!  COMPUTE STABILITY DEPENDENT EXCHANGE COEFFICIENTS
   if (islimsk == 0) then
      ustar = sqrt(grav * z0 / charnok)
   end if
!  COMPUTE STABILITY INDICES (RB AND HLINF)
   z0max = min(z0,one*z1)
   ztmax = z0max
   if (islimsk == 0) then
      restar=ustar*z0max/vis
      restar=max(restar,1.e-6_r_kind)

!     Rat taken from Zeng, Zhao and Dickinson 1997
      rat = 2.67_r_kind * restar**quarter - 2.57_r_kind
      rat = min(rat,7.0_r_kind)
      ztmax = z0max * exp(-rat)
   end if

   dtv = thv1 - tvs
   adtv = abs(dtv)
   adtv = max(adtv,1.e-3_r_kind)
   dtv = sign(one,dtv)*adtv
   rb = grav*dtv*z1 / (half*(thv1+tvs)*wind*wind)
   rb=max(rb,-5.e3_r_kind)
   fm=log((z0max+z1) / z0max)
   fh=log((ztmax+z1) / ztmax)
   hlinf=rb*fm*fm / fh
   fm10=log((z0max+ten) / z0max)

!  STABLE CASE
   if (dtv >= zero) then
      hl1=hlinf
   end if
   if ((dtv >= zero) .AND. (hlinf > quarter)) then
      hl0inf=z0max*hlinf/z1
      hltinf=ztmax*hlinf/z1
      aa=sqrt(one + four*alpha*hlinf)
      aa0=sqrt(one + four*alpha*hl0inf)
      bb=aa
      bb0=sqrt(one + four*alpha*hltinf)
      pm=aa0 - aa + log((aa+one)/(aa0+one))
      ph=bb0 - bb + log((bb+one)/(bb0+one))
      fms=fm-pm
      fhs=fh-ph
      hl1=fms*fms*rb/fhs
   end if

!  SECOND ITERATION
   if (dtv >= zero) then
      hl0=z0max*hl1/z1
      hlt=ztmax*hl1/z1
      aa=sqrt(one + four*alpha*hl1)
      aa0=sqrt(one + four*alpha*hl0)
      bb=aa
      bb0=sqrt(one + four*alpha*hlt)
      pm=aa0 - aa + log((aa+one)/(aa0+one))
      ph=bb0 - bb + log((bb+one)/(bb0+one))
      hl110=hl1*ten/z1
      aa=sqrt(one + four*alpha*hl110)
      pm10=aa0 - aa + log((aa+one)/(aa0+one))
      hl12=hl1*two/z1
      bb=sqrt(one + four*alpha*hl12)
      ph2=bb0 - bb + log((bb+one)/(bb0+one))
   end if

!  UNSTABLE CASE
!  CHECK FOR UNPHYSICAL OBUKHOV LENGTH
   if (dtv < zero) then
      olinf = z1/hlinf
      if ( abs(olinf) <= z0max*50.0_r_kind ) then
         hlinf = -z1/(50.0_r_kind*z0max)
      end if
   end if

!  GET PM AND PH
   if (dtv < zero .AND. hlinf >= (-half)) then
      hl1=hlinf
      pm=(a0 + a1*hl1)*hl1/(one + b1*hl1 + b2*hl1*hl1)
      ph=(a0p + a1p*hl1)*hl1/(one + b1p*hl1 + b2*hl1*hl1)
      hl110=hl1*ten/z1
      pm10=(a0 + a1*hl110)*hl110/(one + b1*hl110 + b2*hl110*hl110)
      hl12=hl1*two/z1
      ph2=(a0p + a1p*hl12)*hl12/(one + b1p*hl12 + b2p*hl12*hl12)
   end if
   if (dtv < zero .AND. hlinf < (-half)) then
      hl1=-hlinf
      pm=log(hl1) + two*hl1**(-quarter) - 0.8776_r_kind
      ph=log(hl1) + half*hl1**(-half) + 1.386_r_kind
      hl110=hl1*ten/z1
      pm10=log(hl110) + two*hl110**(-quarter) - 0.8776_r_kind
      hl12=hl1*two/z1
      ph2=log(hl12) + half*hl12**(-half) + 1.386_r_kind
   end if

!  FINISH THE EXCHANGE COEFFICIENT COMPUTATION TO PROVIDE FM AND FH
   fm=fm-pm
   fh=fh-ph
   fm10=fm10-pm10
   f10m=fm10/fm

   return
  end subroutine compute_fact10