text:
module eqn_rans_rstm
#include <messenger.h>
    use type_equation_set,          only: equation_set_t
    use type_equation_builder,      only: equation_builder_t
    use type_fluid_pseudo_timestep, only: fluid_pseudo_timestep_t
    implicit none

    !>
    !!
    !!
    !! @author  Eric M. Wolf
    !! @date    02/01/2018
    !!
    !--------------------------------------------------------------------------------
    type, public, extends(equation_builder_t) :: rans_rstm

    contains

        procedure   :: init
        procedure   :: build

    end type rans_rstm
    !*******************************************************************************************

contains

    !>
    !!
    !!
    !! @author  Eric M. Wolf
    !! @date    02/01/2018
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(rans_rstm),   intent(inout)  :: self

        call self%set_name('RANS_RSTM')

    end subroutine init
    !*******************************************************************************************

    !>
    !!
    !!
    !! @author  Eric M. Wolf
    !! @date    02/01/2018
    !!
    !--------------------------------------------------------------------------------
    function build(self,blueprint) result(rans_rstm_eqns)
        class(rans_rstm),    intent(in)  :: self
        character(*),   intent(in)  :: blueprint

        type(equation_set_t)            :: rans_rstm_eqns
        type(fluid_pseudo_timestep_t)   :: fluid_pseudo_time

        !
        ! Set equation set name
        !
        call rans_rstm_eqns%set_name('RANS_RSTM')

        !
        ! Add spatial operators
        !
        select case (trim(blueprint))

            case('default')
                call rans_rstm_eqns%add_operator('Euler Volume Flux')
                call rans_rstm_eqns%add_operator('Euler Boundary Average Flux')
                call rans_rstm_eqns%add_operator('Euler Roe Flux')
                call rans_rstm_eqns%add_operator('Euler BC Flux')
                call rans_rstm_eqns%add_operator('Euler Volume Cylindrical Source')

                call rans_rstm_eqns%add_operator('Fluid Viscous Volume Operator')
                call rans_rstm_eqns%add_operator('Fluid Viscous Boundary Average Operator')
                call rans_rstm_eqns%add_operator('Fluid Viscous BC Operator')
                call rans_rstm_eqns%add_operator('Fluid Viscous Volume Cylindrical Source')

! NATHAN COMMENTED THIS OUT DURING MERGE, BECAUSE IT WASN'T INCLUDED IN RSTM DIRECTORY
!                call rans_rstm_eqns%add_model('RSTM Turbulence Kinetic Energy')
!

                !call rans_rstm_eqns%add_model('Ideal Gas RSTM')
                !call rans_rstm_eqns%add_model('Primary Fluid Fields')
                call rans_rstm_eqns%add_model('Ideal Gas')
                call rans_rstm_eqns%add_model('Constant Viscosity')
                !call rans_rstm_eqns%add_model('Sutherlands Law')
                call rans_rstm_eqns%add_model('Stokes Hypothesis')
                call rans_rstm_eqns%add_model('Reynolds Analogy')
                !call rans_rstm_eqns%add_model('Zero Reynolds Stress')
                !call rans_rstm_eqns%add_model('Fluid Advection Velocity')
                !call rans_rstm_eqns%add_model('Constant Viscosity')
                call rans_rstm_eqns%add_model('Velocity Gradient')
                call rans_rstm_eqns%add_model('Strain Rate')
                call rans_rstm_eqns%add_model('Rotation Rate')
                !call rans_rstm_eqns%add_model('Wall Distance : p-Poisson Normalization')
                !call rans_rstm_eqns%add_model('h Smooth')
                !call rans_rstm_eqns%add_model('Zero Turbulent Model Fields')

                call rans_rstm_eqns%add_model('RSTMSSGLRRW Realizable Reynolds Stress')
                call rans_rstm_eqns%add_model('RSTMSSGLRRW LRR Coefficients')
                call rans_rstm_eqns%add_model('RSTMSSGLRRW Turbulence Quantities')
                call rans_rstm_eqns%add_model('RSTMSSGLRRW Production')
                call rans_rstm_eqns%add_model('RSTMSSGLRRW Isotropic Dissipation')
                call rans_rstm_eqns%add_model('RSTMSSGLRRW Pressure-Strain Correlation')
                call rans_rstm_eqns%add_model('RSTMSSGLRRW Realizability Source')
                call rans_rstm_eqns%add_model('RSTMSSGLRRW Artificial Viscosity')

                ! Add RSTM diffusion model
                ! 'RSTMSSGLRRW Generalized Diffusion' is the standard model
                ! 'RSTMSSGLRRW Simple Diffusion' can be used for greater numerical robustness
                !call rans_rstm_eqns%add_model('RSTMSSGLRRW Generalized Diffusion')
                call rans_rstm_eqns%add_model('RSTMSSGLRRW Simple Diffusion')

                !call rans_rstm_eqns%add_operator('Spalart-Allmaras Source Operator')
                !call rans_rstm_eqns%add_operator('Spalart-Allmaras Advection Boundary Average Operator')
                !call rans_rstm_eqns%add_operator('Spalart-Allmaras LaxFriedrichs Operator')
                !call rans_rstm_eqns%add_operator('Spalart-Allmaras Volume Advection Operator')
                !call rans_rstm_eqns%add_operator('Spalart-Allmaras BC Advection Operator')
                !call rans_rstm_eqns%add_operator('Spalart-Allmaras Boundary Diffusion Operator')
                !call rans_rstm_eqns%add_operator('Spalart-Allmaras Volume Diffusion Operator')
                !call rans_rstm_eqns%add_operator('Spalart-Allmaras BC Diffusion Operator')

                call rans_rstm_eqns%add_operator('RSTMSSGLRRW Source Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW Advection Boundary Average Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW LaxFriedrichs Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW Volume Advection Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW BC Advection Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW Boundary Diffusion Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW Volume Diffusion Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW BC Diffusion Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW Artificial Viscosity Volume Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW Artificial Viscosity BC Operator')
                call rans_rstm_eqns%add_operator('RSTMSSGLRRW Artificial Viscosity Boundary Average Operator')

                ! Add shear stress after turbulence viscosity models from SA so they are computed first
                call rans_rstm_eqns%add_model('Reynolds Shear Stress')
                !call rans_rstm_eqns%add_model('Self Blended Reynolds Shear Stress')
                call rans_rstm_eqns%add_model('Temperature Gradient')

                call rans_rstm_eqns%add_pseudo_timestep(fluid_pseudo_time)

                !call rans_rstm_eqns%add_io_field('Smoothed h Field - 2')
                call rans_rstm_eqns%add_io_field('Turbulent Viscosity')

            case default
                call chidg_signal_one(FATAL, "build_rans_rstm: I didn't recognize the &
                                              construction parameter that was passed to build &
                                              the equation set.", blueprint)

        end select

    end function build
    !******************************************************************************************

end module eqn_rans_rstm

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
*DECK CBESY
      SUBROUTINE CBESY (Z, FNU, KODE, N, CY, NZ, CWRK, IERR)
C***BEGIN PROLOGUE  CBESY
C***PURPOSE  Compute a sequence of the Bessel functions Y(a,z) for
C            complex argument z and real nonnegative orders a=b,b+1,
C            b+2,... where b>0.  A scaling option is available to
C            help avoid overflow.
C***LIBRARY   SLATEC
C***CATEGORY  C10A4
C***TYPE      COMPLEX (CBESY-C, ZBESY-C)
C***KEYWORDS  BESSEL FUNCTIONS OF COMPLEX ARGUMENT,
C             BESSEL FUNCTIONS OF SECOND KIND, WEBER'S FUNCTION,
C             Y BESSEL FUNCTIONS
C***AUTHOR  Amos, D. E., (SNL)
C***DESCRIPTION
C
C         On KODE=1, CBESY computes an N member sequence of complex
C         Bessel functions CY(L)=Y(FNU+L-1,Z) for real nonnegative
C         orders FNU+L-1, L=1,...,N and complex Z in the cut plane
C         -pi<arg(Z)<=pi.  On KODE=2, CBESY returns the scaled
C         functions
C
C            CY(L) = exp(-abs(Y))*Y(FNU+L-1,Z),  L=1,...,N, Y=Im(Z)
C
C         which remove the exponential growth in both the upper and
C         lower half planes as Z goes to infinity.  Definitions and
C         notation are found in the NBS Handbook of Mathematical
C         Functions (Ref. 1).
C
C         Input
C           Z      - Nonzero argument of type COMPLEX
C           FNU    - Initial order of type REAL, FNU>=0
C           KODE   - A parameter to indicate the scaling option
C                    KODE=1  returns
C                            CY(L)=Y(FNU+L-1,Z), L=1,...,N
C                        =2  returns
C                            CY(L)=Y(FNU+L-1,Z)*exp(-abs(Y)), L=1,...,N
C                            where Y=Im(Z)
C           N      - Number of terms in the sequence, N>=1
C           CWRK   - A work vector of type COMPLEX and dimension N
C
C         Output
C           CY     - Result vector of type COMPLEX
C           NZ     - Number of underflows set to zero
C                    NZ=0    Normal return
C                    NZ>0    CY(L)=0 for NZ values of L, usually on
C                            KODE=2 (the underflows may not be in an
C                            uninterrupted sequence)
C           IERR   - Error flag
C                    IERR=0  Normal return     - COMPUTATION COMPLETED
C                    IERR=1  Input error       - NO COMPUTATION
C                    IERR=2  Overflow          - NO COMPUTATION
C                            (abs(Z) too small and/or FNU+N-1
C                            too large)
C                    IERR=3  Precision warning - COMPUTATION COMPLETED
C                            (Result has half precision or less
C                            because abs(Z) or FNU+N-1 is large)
C                    IERR=4  Precision error   - NO COMPUTATION
C                            (Result has no precision because
C                            abs(Z) or FNU+N-1 is too large)
C                    IERR=5  Algorithmic error - NO COMPUTATION
C                            (Termination condition not met)
C
C *Long Description:
C
C         The computation is carried out by the formula
C
C            Y(a,z) = (H(1,a,z) - H(2,a,z))/(2*i)
C
C         where the Hankel functions are computed as described in CBESH.
C
C         For negative orders, the formula
C
C            Y(-a,z) = Y(a,z)*cos(a*pi) + J(a,z)*sin(a*pi)
C
C         can be used.  However, for large orders close to half odd
C         integers the function changes radically.  When a is a large
C         positive half odd integer, the magnitude of Y(-a,z)=J(a,z)*
C         sin(a*pi) is a large negative power of ten.  But when a is
C         not a half odd integer, Y(a,z) dominates in magnitude with a
C         large positive power of ten and the most that the second term
C         can be reduced is by unit roundoff from the coefficient.
C         Thus,  wide changes can occur within unit roundoff of a large
C         half odd integer.  Here, large means a>abs(z).
C
C         In most complex variable computation, one must evaluate ele-
C         mentary functions.  When the magnitude of Z or FNU+N-1 is
C         large, losses of significance by argument reduction occur.
C         Consequently, if either one exceeds U1=SQRT(0.5/UR), then
C         losses exceeding half precision are likely and an error flag
C         IERR=3 is triggered where UR=R1MACH(4)=UNIT ROUNDOFF.  Also,
C         if either is larger than U2=0.5/UR, then all significance is
C         lost and IERR=4.  In order to use the INT function, arguments
C         must be further restricted not to exceed the largest machine
C         integer, U3=I1MACH(9).  Thus, the magnitude of Z and FNU+N-1
C         is restricted by MIN(U2,U3).  In IEEE arithmetic, U1,U2, and
C         U3 approximate 2.0E+3, 4.2E+6, 2.1E+9 in single precision
C         and 4.7E+7, 2.3E+15 and 2.1E+9 in double precision.  This
C         makes U2 limiting in single precision and U3 limiting in
C         double precision.  This means that one can expect to retain,
C         in the worst cases on IEEE machines, no digits in single pre-
C         cision and only 6 digits in double precision.  Similar con-
C         siderations hold for other machines.
C
C         The approximate relative error in the magnitude of a complex
C         Bessel function can be expressed as P*10**S where P=MAX(UNIT
C         ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre-
C         sents the increase in error due to argument reduction in the
C         elementary functions.  Here, S=MAX(1,ABS(LOG10(ABS(Z))),
C         ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF
C         ABS(Z),ABS(EXPONENT OF FNU)) ).  However, the phase angle may
C         have only absolute accuracy.  This is most likely to occur
C         when one component (in magnitude) is larger than the other by
C         several orders of magnitude.  If one component is 10**K larger
C         than the other, then one can expect only MAX(ABS(LOG10(P))-K,
C         0) significant digits; or, stated another way, when K exceeds
C         the exponent of P, no significant digits remain in the smaller
C         component.  However, the phase angle retains absolute accuracy
C         because, in complex arithmetic with precision P, the smaller
C         component will not (as a rule) decrease below P times the
C         magnitude of the larger component.  In these extreme cases,
C         the principal phase angle is on the order of +P, -P, PI/2-P,
C         or -PI/2+P.
C
C***REFERENCES  1. M. Abramowitz and I. A. Stegun, Handbook of Mathe-
C                 matical Functions, National Bureau of Standards
C                 Applied Mathematics Series 55, U. S. Department
C                 of Commerce, Tenth Printing (1972) or later.
C               2. D. E. Amos, Computation of Bessel Functions of
C                 Complex Argument, Report SAND83-0086, Sandia National
C                 Laboratories, Albuquerque, NM, May 1983.
C               3. D. E. Amos, Computation of Bessel Functions of
C                 Complex Argument and Large Order, Report SAND83-0643,
C                 Sandia National Laboratories, Albuquerque, NM, May
C                 1983.
C               4. D. E. Amos, A Subroutine Package for Bessel Functions
C                 of a Complex Argument and Nonnegative Order, Report
C                 SAND85-1018, Sandia National Laboratory, Albuquerque,
C                 NM, May 1985.
C               5. D. E. Amos, A portable package for Bessel functions
C                 of a complex argument and nonnegative order, ACM
C                 Transactions on Mathematical Software, 12 (September
C                 1986), pp. 265-273.
C
C***ROUTINES CALLED  CBESH, I1MACH, R1MACH
C***REVISION HISTORY  (YYMMDD)
C   830501  DATE WRITTEN
C   890801  REVISION DATE from Version 3.2
C   910415  Prologue converted to Version 4.0 format.  (BAB)
C   920128  Category corrected.  (WRB)
C   920811  Prologue revised.  (DWL)
C***END PROLOGUE  CBESY
C
      COMPLEX CWRK, CY, C1, C2, EX, HCI, Z, ZU, ZV
      REAL ELIM, EY, FNU, R1, R2, TAY, XX, YY, R1MACH, R1M5, ASCLE,
     *  RTOL, ATOL, TOL, AA, BB
      INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH
      DIMENSION CY(N), CWRK(N)
C***FIRST EXECUTABLE STATEMENT  CBESY
      XX = REAL(Z)
      YY = AIMAG(Z)
      IERR = 0
      NZ=0
      IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1
      IF (FNU.LT.0.0E0) IERR=1
      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
      IF (N.LT.1) IERR=1
      IF (IERR.NE.0) RETURN
      HCI = CMPLX(0.0E0,0.5E0)
      CALL CBESH(Z, FNU, KODE, 1, N, CY, NZ1, IERR)
      IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
      CALL CBESH(Z, FNU, KODE, 2, N, CWRK, NZ2, IERR)
      IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
      NZ = MIN(NZ1,NZ2)
      IF (KODE.EQ.2) GO TO 60
      DO 50 I=1,N
        CY(I) = HCI*(CWRK(I)-CY(I))
   50 CONTINUE
      RETURN
   60 CONTINUE
      TOL = MAX(R1MACH(4),1.0E-18)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      K = MIN(ABS(K1),ABS(K2))
      R1M5 = R1MACH(5)
C-----------------------------------------------------------------------
C     ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
C-----------------------------------------------------------------------
      ELIM = 2.303E0*(K*R1M5-3.0E0)
      R1 = COS(XX)
      R2 = SIN(XX)
      EX = CMPLX(R1,R2)
      EY = 0.0E0
      TAY = ABS(YY+YY)
      IF (TAY.LT.ELIM) EY = EXP(-TAY)
      IF (YY.LT.0.0E0) GO TO 90
      C1 = EX*CMPLX(EY,0.0E0)
      C2 = CONJG(EX)
   70 CONTINUE
      NZ = 0
      RTOL = 1.0E0/TOL
      ASCLE = R1MACH(1)*RTOL*1.0E+3
      DO 80 I=1,N
C       CY(I) = HCI*(C2*CWRK(I)-C1*CY(I))
        ZV = CWRK(I)
        AA=REAL(ZV)
        BB=AIMAG(ZV)
        ATOL=1.0E0
        IF (MAX(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 75
          ZV = ZV*CMPLX(RTOL,0.0E0)
          ATOL = TOL
   75   CONTINUE
        ZV = ZV*C2*HCI
        ZV = ZV*CMPLX(ATOL,0.0E0)
        ZU=CY(I)
        AA=REAL(ZU)
        BB=AIMAG(ZU)
        ATOL=1.0E0
        IF (MAX(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 85
          ZU = ZU*CMPLX(RTOL,0.0E0)
          ATOL = TOL
   85   CONTINUE
        ZU = ZU*C1*HCI
        ZU = ZU*CMPLX(ATOL,0.0E0)
        CY(I) = ZV - ZU
        IF (CY(I).EQ.CMPLX(0.0E0,0.0E0) .AND. EY.EQ.0.0E0) NZ = NZ + 1
   80 CONTINUE
      RETURN
   90 CONTINUE
      C1 = EX
      C2 = CONJG(EX)*CMPLX(EY,0.0E0)
      GO TO 70
  170 CONTINUE
      NZ = 0
      RETURN
      END

ext: f


!
!  Copyright 2017 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
subroutine setmg(mg_sta,mg_end,mg_num,mg_sta_all,mg_end_all,mg_num_all,  &
                 lg_sta,lg_num,nproc,nproc_id_global,nproc_Mxin,nproc_ob,isequential)
implicit none
integer :: i1,j1,j2,j3
integer :: ibox
integer :: isequential
integer :: nproc,nproc_id_global
integer :: nproc_Mxin(3)
integer :: nproc_ob
integer :: mg_sta(3),mg_end(3),mg_num(3)
integer :: lg_sta(3),lg_num(3)
integer :: mg_sta_all(3,0:nproc-1),mg_end_all(3,0:nproc-1),mg_num_all(3,0:nproc-1)
integer :: nproc_Mxin_mul

nproc_Mxin_mul=nproc_Mxin(1)*nproc_Mxin(2)*nproc_Mxin(3)
if(isequential==1)then
  do j3=0,nproc_Mxin(3)-1
  do j2=0,nproc_Mxin(2)-1
  do j1=0,nproc_Mxin(1)-1
    do i1=0,nproc_ob-1
      ibox = i1 + nproc_ob*( j1 + nproc_Mxin(1)*j2 + nproc_Mxin(1)*nproc_Mxin(2)*j3 )
      mg_sta_all(1,ibox)=j1*lg_num(1)/nproc_Mxin(1)+lg_sta(1)
      mg_end_all(1,ibox)=(j1+1)*lg_num(1)/nproc_Mxin(1)+lg_sta(1)-1
      mg_sta_all(2,ibox)=j2*lg_num(2)/nproc_Mxin(2)+lg_sta(2)
      mg_end_all(2,ibox)=(j2+1)*lg_num(2)/nproc_Mxin(2)+lg_sta(2)-1
      mg_sta_all(3,ibox)=j3*lg_num(3)/nproc_Mxin(3)+lg_sta(3)
      mg_end_all(3,ibox)=(j3+1)*lg_num(3)/nproc_Mxin(3)+lg_sta(3)-1
    end do
  end do
  end do
  end do
else if(isequential==2)then
  do i1=0,nproc_ob-1
    do j3=0,nproc_Mxin(3)-1
    do j2=0,nproc_Mxin(2)-1
    do j1=0,nproc_Mxin(1)-1
      ibox = j1 + nproc_Mxin(1)*j2 + nproc_Mxin(1)*nproc_Mxin(2)*j3 + i1*nproc_Mxin_mul
      mg_sta_all(1,ibox)=j1*lg_num(1)/nproc_Mxin(1)+lg_sta(1)
      mg_end_all(1,ibox)=(j1+1)*lg_num(1)/nproc_Mxin(1)+lg_sta(1)-1
      mg_sta_all(2,ibox)=j2*lg_num(2)/nproc_Mxin(2)+lg_sta(2)
      mg_end_all(2,ibox)=(j2+1)*lg_num(2)/nproc_Mxin(2)+lg_sta(2)-1
      mg_sta_all(3,ibox)=j3*lg_num(3)/nproc_Mxin(3)+lg_sta(3)
      mg_end_all(3,ibox)=(j3+1)*lg_num(3)/nproc_Mxin(3)+lg_sta(3)-1
    end do
    end do
    end do
  end do
end if
mg_num_all(:,:)=mg_end_all(:,:)-mg_sta_all(:,:)+1

mg_sta(:)=mg_sta_all(:,nproc_id_global)
mg_end(:)=mg_end_all(:,nproc_id_global)
mg_num(:)=mg_num_all(:,nproc_id_global)

end subroutine setmg


text:
      SUBROUTINE TCONAS	(IERROR)
C
C LATEST REVISION	 AUGUST	1984
C
C PURPOSE		 TO PROVIDE A SIMPLE DEMONSTRATION OF
C			 CONRAS, THE SUPER ENTRY  POINT	OF THE
C			 CONRAN	PACKAGE.
C
C USAGE			 CALL TCONAS (IERROR)
C
C ARGUMENTS
C
C ON OUTPUT		 IERROR
C			   AN INTEGER VARIABLE
C			   .EQ.	0, IF THE TEST WAS SUCCESSFUL,
C			   .NE.	0, OTHERWISE
C			   IF NOT ZERO THE NUMBER PRODUCED WILL
C			   CORRESPOND TO THE ERROR NUMBERS IN
C			   THE CONRAS LISTING.
C
C I/O			 IF THE	TEST IS	SUCCESSFUL, THE	MESSAGE
C
C			   CONRAS TEST SUCCESSFUL  . . . SEE PLOT TO
C			   VERIFY PERFORMANCE
C
C			 IS PRINTED ON UNIT 6.
C			 IN ADDITION, TWO FRAMES CONTAINING THE	CONTOUR
C			 PLOT AND TRIANGULATION	OF THE DATA ARE	PRODUCED
C			 ON THE	DEFAULT	GRAPHICS DEVICE	UNLESS THE USER
C			 SPECIFIES OTHERWISE VIA JCL.
C
C PRECISION		 SINGLE
C
C REQUIRED LIBRARY	 CONRAS
C FILES			 CONTERP
C			 CONCOM
C			 DASHSUPR
C
C SPECIALIST		 FOR INFORMATION ABOUT THIS ROUTINE OR THE
C			 ULIB CONRAS PACKAGE, CONTACT THE SPECIALIST
C			 NAMED IN THE ULIB CONRAS PACKAGE.
C
C LANGUAGE		 FORTRAN
C
C ALGORITHM		 A SPARSE DATA SET IS DEFINED VIA DATA
C			 STATEMENTS.  OPTIONS ARE SET TO PRODUCE A
C			 TITLE AND DISPLAY THE TRIANGULATION GENERATED
C			 BY THE	INTERPOLATING ROUTINES.	 BY DEFAULT
C			 A MESSAGE AT THE BOTTEM OF THE	PLOT AND A
C			 PERIMETER ARE ALSO PRODUCED.  THIS ROUTINE
C			 TAKES ADVANTAGE OF THE	PORT ERROR HANDLING
C			 ROUTINES TO DETERMINE IF CONRAS TERMINATED
C			 NORMALLY.
C
C PORTABILITY		 ANSI STANDARD
C
C
C  SET UP THE SCRATCH SPACES REQUIRED BY CONRAS
C
      DIMENSION	WK(221),IWK(744),SCR(1600)
C
C  SET UP THE ARRAYS TO	DEFINE THE DATA	SET
C
      DIMENSION	XD(17),YD(17),ZD(17)
	COMMON /RASINT/	IRASMJ,	IRASMN,	IRASTX
C
C DEFINE THE DATA SET
C
      DATA XD(1),XD(2),XD(3),XD(4),XD(5),XD(6),XD(7),XD(8),
     1	   XD(9),XD(10),XD(11),XD(12),XD(13),XD(14),XD(15),
     2	   XD(16),XD(17)
     3	   /3.,3.,10.,18.,18.,10.,10.,5.,1.,15.,20.,
     4		5.,15.,10.,7.,13.,16./
C
      DATA YD(1),YD(2),YD(3),YD(4),YD(5),YD(6),YD(7),YD(8),
     1	   YD(9),YD(10),YD(11),YD(12),YD(13),YD(14),YD(15),
     2	   YD(16),YD(17)
     3	  /3.,18.,18.,3.,18.,10.,1.,5.,10.,5.,10.,
     4		15.,15.,15.,20.,20.,8./
C
      DATA ZD(1),ZD(2),ZD(3),ZD(4),ZD(5),ZD(6),ZD(7),ZD(8),
     1	   ZD(9),ZD(10),ZD(11),ZD(12),ZD(13),ZD(14),ZD(15),
     2	   ZD(16),ZD(17)
     3	  /25.,25.,25.,25.,25.,-5.,1.,1.,1.,1.,1.,
     4		1.,1.,1.,1.,1.,25./
C
C  SET UP PARAMETER FOR	NUMBER OF INPUT	POINTS
C
      DATA NDP/17/
C
C  SET PORT ERROR HANDLING ROUTINE TO RECOVERY MODE
C
      CALL ENTSR(IROLD,1)
C
C  SET UP TITLE	 FOR PLOT
C
      CALL CONOP4('TLE=ON','DEMONSTRATION PLOT FOR CONRAS',29,0)
C
C  TEST	FOR ERROR
C
      IF (NERRO(IERROR).NE.0) GO TO 100
C
C  NO ERROR
C
C  SET OPTION TO DISPLAY THE TRIANGULATION
C
      CALL CONOP1('TRI=ON')
C
C  TEST	FOR ERROR
C
      IF (NERRO(IERROR).NE.0) GO TO 100
C
C  NO ERROR
C
C  CALL	CONRAS TO CONTOUR DATA
C
      CALL CONRAS(XD,YD,ZD,NDP,WK,IWK,SCR)
C
C  TEST	FOR ERROR
C
      IF (NERRO(IERROR).NE.0) GO TO 100
C
C  NO ERROR
C
C
C  CALL	FRAME, CONRAS WILL NOT DO THIS
C
c	CALL NEWFM
C
C  PRINT MESSAGE EVERYTHING OK
C
c     WRITE(6,10)
c10   FORMAT(1X,'CONRAS	TEST SUCCESSFUL, SEE PLOT TO VERIFY ',
c    1'PERFORMANCE')
C
C
      RETURN
C
C  IF ERROR CALL THE PORT ERROR	PRINT ROUTINE.
C  THIS	CALL IS	NOT NECESSARY UNLESS YOU ARE IN	RECOVER	MODE.
C  IF YOU ARE NOT IN RECOVER MODE THE ERROR MESSAGE WILL  BE PRINTED
C  AUTOMATICALLY.
C
 100  CALL EPRIN
      RETURN
C
      END

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


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

module obs_err_mod

use        types_mod, only : r8, missing_r8

implicit none
private

integer, parameter :: nobs_level = 15
real(r8)           :: obs_prs(nobs_level)

! this array is ordered from top of atm down to surface
data obs_prs/ 10.0_r8,  20.0_r8,  30.0_r8,  50.0_r8,  70.0_r8, &
             100.0_r8, 150.0_r8, 200.0_r8, 250.0_r8, 300.0_r8, &
             400.0_r8, 500.0_r8, 700.0_r8, 850.0_r8, 1050.0_r8/

public :: fixed_marine_pres_error,      &
          fixed_marine_rel_hum_error,   &
          fixed_marine_temp_error,      &
          fixed_marine_wind_error,      &
          moving_marine_pres_error,     &
          moving_marine_rel_hum_error,  &
          moving_marine_temp_error,     &
          moving_marine_wind_error

public :: land_pres_error,              &
          land_rel_hum_error,           &
          land_temp_error,              &
          land_wind_error

public :: acars_rel_hum_error,          &
          acars_temp_error,             &
          acars_wind_error

public :: sat_wind_error,               &
          sat_wv_wind_error

public :: rawin_pres_error,             &
          rawin_rel_hum_error,          &
          rawin_temp_error,             &
          rawin_wind_error

public :: drop_pres_error,              &
          drop_rel_hum_error,           &
          drop_temp_error,              &
          drop_wind_error

public :: prof_wind_error

contains

function acars_rel_hum_error(pres, tmpk, rh)

real(r8), intent(in) :: pres, tmpk, rh   !  (mb)

real(r8) :: acars_rel_hum_error

acars_rel_hum_error = missing_r8

return
end function acars_rel_hum_error

function acars_temp_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, acars_temp_error

data obs_err/2.2_r8, 2.0_r8, 1.8_r8, 1.6_r8, 1.5_r8, &
             1.4_r8, 1.4_r8, 1.4_r8, 1.3_r8, 1.3_r8, &
             1.2_r8, 1.2_r8, 1.2_r8, 1.3_r8, 1.4_r8/

call find_pressure_level_weight(pres, k0, wght)
acars_temp_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
acars_temp_error = dble(nint(acars_temp_error * 10.0_r8)) * 0.1_r8

return
end function acars_temp_error

function acars_wind_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, acars_wind_error

data obs_err/4.0_r8, 4.0_r8, 4.0_r8, 4.0_r8, 4.0_r8, &
             4.0_r8, 4.0_r8, 4.0_r8, 4.0_r8, 4.0_r8, &
             4.0_r8, 3.5_r8, 3.0_r8, 2.5_r8, 2.5_r8/

call find_pressure_level_weight(pres, k0, wght)
acars_wind_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
acars_wind_error = dble(nint(acars_wind_error * 10.0_r8)) * 0.1_r8

return
end function acars_wind_error

function land_pres_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, land_pres_error

data obs_err/114.2_r8, 96.0_r8, 69.8_r8, 59.3_r8, 50.3_r8, &
              39.4_r8, 32.4_r8, 27.7_r8, 25.4_r8, 18.8_r8, &
              14.9_r8, 12.1_r8,  8.6_r8,  8.0_r8,  7.0_r8/

call find_pressure_level_weight(pres, k0, wght)
land_pres_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
land_pres_error = land_pres_error * 1.225_r8 * 0.25_r8
land_pres_error = dble(nint(land_pres_error * 10.0_r8)) * 0.1_r8

return
end function land_pres_error

function land_rel_hum_error(pres, tmpk, rh)

real(r8), intent(in) :: pres, tmpk, rh

real(r8) :: land_rel_hum_error

if ( RH < 0.2_r8 ) then
  land_rel_hum_error = 0.23_r8
elseif ( tmpk < 233.0_r8 ) then
  land_rel_hum_error = 0.28_r8
else
  land_rel_hum_error = 0.13_r8
end if

return
end function land_rel_hum_error

function land_temp_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer :: k0
real(r8) :: obs_err(nobs_level), wght, land_temp_error

data obs_err/2.5_r8, 2.5_r8, 2.5_r8, 2.4_r8, 2.2_r8, &
             2.0_r8, 1.9_r8, 1.8_r8, 1.8_r8, 1.5_r8, &
             1.3_r8, 1.2_r8, 1.3_r8, 1.5_r8, 2.0_r8/

call find_pressure_level_weight(pres, k0, wght)
land_temp_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
land_temp_error = dble(nint(land_temp_error * 10.0_r8)) * 0.1_r8

return
end function land_temp_error

function land_wind_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, land_wind_error

data obs_err/3.0_r8, 2.5_r8, 2.0_r8, 2.0_r8, 2.0_r8, &
             2.2_r8, 2.4_r8, 3.2_r8, 3.2_r8, 3.8_r8, &
             3.6_r8, 3.4_r8, 3.0_r8, 3.0_r8, 3.0_r8/

call find_pressure_level_weight(pres, k0, wght)
land_wind_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
land_wind_error = dble(nint(land_wind_error * 10.0_r8)) * 0.1_r8

return
end function land_wind_error

function fixed_marine_pres_error(pres)

real(r8), intent(in) :: pres  !  (mb)

real(r8) :: fixed_marine_pres_error

fixed_marine_pres_error = land_pres_error(pres)

return
end function fixed_marine_pres_error

function fixed_marine_rel_hum_error(pres, tmpk, rh)

real(r8), intent(in) :: pres, tmpk, rh

real(r8) :: fixed_marine_rel_hum_error

fixed_marine_rel_hum_error = land_rel_hum_error(pres,tmpk,rh)

return
end function fixed_marine_rel_hum_error

function fixed_marine_temp_error(pres)

real(r8), intent(in) :: pres  !  (mb)

real(r8) :: fixed_marine_temp_error

fixed_marine_temp_error = land_temp_error(pres)

return
end function fixed_marine_temp_error

function fixed_marine_wind_error(pres)

real(r8), intent(in) :: pres  !  (mb)

real(r8) :: fixed_marine_wind_error

fixed_marine_wind_error = land_wind_error(pres)

return
end function fixed_marine_wind_error

function moving_marine_pres_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, moving_marine_pres_error

data obs_err/228.4_r8, 192.0_r8, 139.5_r8, 118.6_r8, 100.6_r8, &
              78.8_r8,  64.8_r8,  55.4_r8,  50.8_r8,  37.6_r8, &
              29.8_r8,  24.2_r8,  17.2_r8,  11.5_r8,  11.5_r8/

call find_pressure_level_weight(pres, k0, wght)
moving_marine_pres_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
moving_marine_pres_error = moving_marine_pres_error * 1.225_r8 * 0.25_r8
moving_marine_pres_error = dble(nint(moving_marine_pres_error * 10.0_r8)) * 0.1_r8

return
end function moving_marine_pres_error

function moving_marine_rel_hum_error(pres, tmpk, rh)

real(r8), intent(in) :: pres, tmpk, rh

real(r8) :: moving_marine_rel_hum_error

moving_marine_rel_hum_error = missing_r8

return
end function moving_marine_rel_hum_error

function moving_marine_temp_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, moving_marine_temp_error

data obs_err/2.5_r8, 2.5_r8, 2.5_r8, 2.4_r8, 2.2_r8, &
             2.0_r8, 1.9_r8, 1.8_r8, 1.8_r8, 1.5_r8, &
             1.3_r8, 1.2_r8, 1.3_r8, 1.5_r8, 1.8_r8/

call find_pressure_level_weight(pres, k0, wght)
moving_marine_temp_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
moving_marine_temp_error = dble(nint(moving_marine_temp_error * 10.0_r8)) * 0.1_r8

return
end function moving_marine_temp_error

function moving_marine_wind_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, moving_marine_wind_error

data obs_err/3.0_r8, 2.5_r8, 2.0_r8, 2.0_r8, 2.0_r8, &
             2.2_r8, 2.4_r8, 3.2_r8, 3.2_r8, 3.8_r8, &
             3.6_r8, 3.4_r8, 2.5_r8, 2.4_r8, 2.4_r8/

call find_pressure_level_weight(pres, k0, wght)
moving_marine_wind_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
moving_marine_wind_error = dble(nint(moving_marine_wind_error * 10.0_r8)) * 0.1_r8

return
end function moving_marine_wind_error

function sat_wind_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, sat_wind_error

data obs_err/5.7_r8, 5.0_r8, 5.0_r8, 5.0_r8, 5.0_r8, &
             5.0_r8, 5.0_r8, 5.0_r8, 5.0_r8, 5.0_r8, &
             4.3_r8, 3.5_r8, 2.0_r8, 2.0_r8, 2.0_r8/

call find_pressure_level_weight(pres, k0, wght)
sat_wind_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
sat_wind_error = dble(nint(sat_wind_error * 10.0_r8)) * 0.1_r8

return
end function sat_wind_error

function sat_wv_wind_error(pres)

real(r8), intent(in) :: pres  !  (mb)

real(r8) :: sat_wv_wind_error

sat_wv_wind_error = sat_wind_error(pres)

return
end function sat_wv_wind_error

function rawin_pres_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, rawin_pres_error

data obs_err/40.0_r8, 32.0_r8, 25.0_r8, 22.5_r8, 19.5_r8, &
             18.1_r8, 15.2_r8, 13.2_r8, 11.8_r8, 10.7_r8, &
              9.8_r8,  8.4_r8,  5.2_r8,  4.4_r8,  4.3_r8/

call find_pressure_level_weight(pres, k0, wght)
rawin_pres_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
rawin_pres_error = rawin_pres_error * 1.225_r8 * 0.25_r8
rawin_pres_error = dble(nint(rawin_pres_error * 10.0_r8)) * 0.1_r8

return
end function rawin_pres_error

function rawin_rel_hum_error(pres, tmpk, rh)

real(r8), intent(in) :: pres, tmpk, rh  !  (mb)

real(r8) :: rawin_rel_hum_error

if ( rh < 0.2_r8 ) then
  rawin_rel_hum_error = 0.23
else if ( tmpk < 233.0_r8 ) then
  rawin_rel_hum_error = 0.28
else
  rawin_rel_hum_error = 0.17
end if

return
end function rawin_rel_hum_error

function rawin_temp_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, rawin_temp_error

data obs_err/2.5_r8, 2.2_r8, 2.0_r8, 1.9_r8, 1.8_r8, &
             1.7_r8, 1.6_r8, 1.5_r8, 1.5_r8, 1.4_r8, &
             1.2_r8, 1.2_r8, 1.3_r8, 1.5_r8, 1.7_r8/

call find_pressure_level_weight(pres, k0, wght)
rawin_temp_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
rawin_temp_error = dble(nint(rawin_temp_error * 10.0_r8)) * 0.1_r8

return
end function rawin_temp_error

function rawin_wind_error(pres)

real(r8), intent(in) :: pres  !  (mb)

integer  :: k0
real(r8) :: obs_err(nobs_level), wght, rawin_wind_error

data obs_err/4.5_r8, 3.6_r8, 3.3_r8, 3.2_r8, 3.2_r8, &
             3.3_r8, 3.4_r8, 3.5_r8, 3.5_r8, 3.7_r8, &
             3.5_r8, 3.0_r8, 2.5_r8, 2.3_r8, 2.3_r8/

call find_pressure_level_weight(pres, k0, wght)
rawin_wind_error = wght * obs_err(k0) + (1.0_r8 - wght) * obs_err(k0+1)
rawin_wind_error = dble(nint(rawin_wind_error * 10.0_r8)) * 0.1_r8

return
end function rawin_wind_error

function drop_pres_error(pres)

real(r8), intent(in) :: pres  !  (mb)

real(r8) :: drop_pres_error

drop_pres_error = rawin_pres_error(pres)

return
end function drop_pres_error

function drop_rel_hum_error(pres, tmpk, rh)

real(r8), intent(in) :: pres, tmpk, rh

real(r8) :: drop_rel_hum_error

drop_rel_hum_error = rawin_rel_hum_error(pres, tmpk, rh)

return
end function drop_rel_hum_error

function drop_temp_error(pres)

real(r8), intent(in) :: pres  !  (mb)

real(r8) :: drop_temp_error

drop_temp_error = rawin_temp_error(pres)

return
end function drop_temp_error

function drop_wind_error(pres)

real(r8), intent(in) :: pres  !  (mb)

real(r8) :: drop_wind_error

drop_wind_error = rawin_wind_error(pres)

return
end function drop_wind_error

function prof_wind_error(pres)

real(r8), intent(in) :: pres  !  (mb)

prof_wind_error = missing_r8

return
end function prof_wind_error

! uses global obs_prs() array.  assumes order of pressure values in
! the array starts at the top and ends at the surface.
subroutine find_pressure_level_weight(pres, zloc, wght)

real(r8), intent(in)  :: pres
integer, intent(out)  :: zloc
real(r8), intent(out) :: wght

integer :: k

if ( pres < obs_prs(1) .or. pres > obs_prs(nobs_level) ) then
  print*,'bad pressure level',pres
  stop
end if

do k = 2, nobs_level

  if ( pres <= obs_prs(k) )  then
    zloc = k - 1
    wght = (log(obs_prs(k)) - log(pres)) / &
           (log(obs_prs(k)) - log(obs_prs(k-1)))
    exit
  end if

end do

return
end subroutine find_pressure_level_weight

end module obs_err_mod

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

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text: constant c3 = 0

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
condition1 <= ( ( not ( ( DanWuzHereS and ( not '0' ) ) or ( ( not DanWuzHereS ) and '0' ) ) ) and ( not ( ( DanWuzHereT and ( not '1' ) ) or ( ( not DanWuzHereT ) and '1' ) ) ) );
DanWuzHereO0 <= ( condition1 and ( DanWuzHereI0 ) ) or ( ( not condition1 ) and ( DanWuzHereI1 ) );
condition2 <= ( ( not ( ( DanWuzHereS and ( not '0' ) ) or ( ( not DanWuzHereS ) and '0' ) ) ) and ( not ( ( DanWuzHereT and ( not '1' ) ) or ( ( not DanWuzHereT ) and '1' ) ) ) );
DanWuzHereO1 <= ( condition2 and ( DanWuzHereI2 ) ) or ( ( not condition2 ) and ( DanWuzHereI3 ) );
condition3 <= ( ( not ( ( DanWuzHereS and ( not '0' ) ) or ( ( not DanWuzHereS ) and '0' ) ) ) and ( not ( ( DanWuzHereT and ( not '1' ) ) or ( ( not DanWuzHereT ) and '1' ) ) ) );
DanWuzHereO2 <= ( condition3 and ( DanWuzHereI4 ) ) or ( ( not condition3 ) and ( DanWuzHereI5 ) );
condition4 <= ( ( not ( ( DanWuzHereS and ( not '1' ) ) or ( ( not DanWuzHereS ) and '1' ) ) ) and ( not ( ( DanWuzHereT and ( not '1' ) ) or ( ( not DanWuzHereT ) and '1' ) ) ) );
DanWuzHereO3 <= ( condition4 and ( DanWuzHereI6 ) ) or ( ( not condition4 ) and ( DanWuzHereI7 ) );
condition5 <= ( ( not ( ( DanWuzHereS and ( not '1' ) ) or ( ( not DanWuzHereS ) and '1' ) ) ) and ( not ( ( DanWuzHereT and ( not '1' ) ) or ( ( not DanWuzHereT ) and '1' ) ) ) );
DanWuzHereO4 <= ( condition5 and ( DanWuzHereI8 ) ) or ( ( not condition5 ) and ( DanWuzHereI9 ) );
condition6 <= ( ( not ( ( DanWuzHereS and ( not '1' ) ) or ( ( not DanWuzHereS ) and '1' ) ) ) and ( not ( ( DanWuzHereT and ( not '1' ) ) or ( ( not DanWuzHereT ) and '1' ) ) ) );
DanWuzHereO5 <= ( condition6 and ( DanWuzHereI10 ) ) or ( ( not condition6 ) and ( DanWuzHereI11 ) );


text:
!!!!!!!###################################################################
!!!!!!!###################################################################
!!!!!!subroutine FractReadApplyPreCrackCleavagePlanes2(x,y,ix)
!!!!!!    use EC_Objects_manager
!!!!!!    use  mod_parameters
!!!!!!    common /wblock8/ abc(573,nume,4),his(573,nume,4)
!!!!!!    common/meshnum/ numnpo, numelto
!!!!!!
!!!!!!    ! common /precrack/ npc, elepc(100),npc_count
!!!!!!    ! common /stressflag/ ElemFractCode(nume),ElemDecayCount(nume)
!!!!!!    ! common /crackopen/ ElemDecayed(nume), overlapele(2,nume)
!!!!!!    ! common /crackline/ ncleave(3,nume), elecrack(4,nume),nodeflag(4,nume)
!!!!!!    ! common /overlapping/ intersec(4, nume), area_coeff(nume), update_flag
!!!!!!
!!!!!!    dimension x(*), y(*), ix(4,*)
!!!!!!    real x,y
!!!!!!    integer ix
!!!!!!
!!!!!!    ! integer elecrack,nodeflag
!!!!!!    ! real ncleave
!!!!!!    ! real intersec,area_coeff
!!!!!!    ! integer ElemFractCode, ElemDecayed,npc, elepc,npc_count
!!!!!!    integer*4 setvbuf3f_local
!!!!!!    integer i,j,elemId,IERR
!!!!!!    real Vcx,Vcy
!!!!!!    real abc,his
!!!!!!    real vn(3)
!!!!!!    integer ElemEdge1,ElemEdge2
!!!!!!    real rxc1,ryc1,r1,rxc2,ryc2,r2
!!!!!!    character strTitle(200)
!!!!!!    real*8 rxc2Ratio,ryc2Ratio,rxc1Ratio,ryc1Ratio
!!!!!!
!!!!!!    integer mEC_bCracking,mEC_ZoneFactor,mEC_DecayCount,mEC_PreCrackedElemCount
!!!!!!    write(*, *) '-->>>>>>>>FractReadApplyPreCrackCleavagePlanes'
!!!!!!    IERR=0
!!!!!!    i=0
!!!!!!    mEC_PreCrackedElemCount=0
!!!!!!    mEC_bCracking=0
!!!!!!    write(*, *) 'beginning of CrackProp ... reading cracks.in'
!!!!!!
!!!!!!    open(60, file = 'cracks.in', status = 'old',IOSTAT=IERR, ERR=90)
!!!!!!    read (60, *) mEC_bCracking !-- 1=allow fracture
!!!!!!    if(mEC_bCracking==0) then
!!!!!!      close(60)
!!!!!!      return
!!!!!!    endif
!!!!!!    open(6011, file = 'cracks2.out', status = 'unknown')
!!!!!!    ierr=setvbuf3f_local(7011,1,100)
!!!!!!
!!!!!!    write(6011, *) mEC_bCracking !-- 1=allow fracture
!!!!!!
!!!!!!    read (60, *)    mEC_ZoneFactor !-- 20 * min elem diagonal
!!!!!!    write (6011, *) mEC_ZoneFactor !-- 20 * min elem diagonal
!!!!!!
!!!!!!    read (60, *)    mEC_DecayCount !-- unloading steps
!!!!!!    write (6011, *) mEC_DecayCount !-- unloading steps
!!!!!!
!!!!!!    read (60, *)    mEC_PreCrackedElemCount  !-- number of pre-cracked elems
!!!!!!    write (6011, *) mEC_PreCrackedElemCount  !-- number of pre-cracked elems
!!!!!!    !-- set the input data to the module
!!!!!!    call ECinitializeConsts(numnpo, numelto,ix,x,y,mEC_bCracking,mEC_ZoneFactor, &
!!!!!!                                mEC_DecayCount,mEC_PreCrackedElemCount)
!!!!!!    call EC_AllocateMem()
!!!!!!
!!!!!!!!! if the fracture is not active then exit
!!!!!!    if (mEC_PreCrackedElemCount == 0) then
!!!!!!        close(60)
!!!!!!        close(6011)
!!!!!!        return
!!!!!!    endif
!!!!!!
!!!!!!    write (6011, *) mEC_PreCrackedElemCount
!!!!!!    read(60, *) Vcx,Vcy
!!!!!!    write(6011, *) Vcx,Vcy
!!!!!!
!!!!!!    read(60, *) strTitle
!!!!!!    write(6011, *) strTitle !'!!!--------- id  ed1     rxc1    ryc1     r1  ed2     rxc2    ryc2     r2'
!!!!!!    do i = 1, mEC_PreCrackedElemCount
!!!!!!        !!!--------- id  ed1     rxc1    ryc1     r1  ed2     rxc2    ryc2     r2
!!!!!!        read(60, *) elemId,ElemEdge1,rxc1,ryc1,r1,ElemEdge2,rxc2,ryc2,r2,rxc2Ratio,ryc2Ratio,rxc1Ratio,ryc1Ratio
!!!!!!        ! elepc(i)=elemId
!!!!!!        write(6011, *) elemId,ElemEdge1,rxc1,ryc1,r1,ElemEdge2,rxc2,ryc2,r2,rxc2Ratio,ryc2Ratio,rxc1Ratio,ryc1Ratio
!!!!!!        write(*, *) 'Pre-exist crack, element', elemId, 'cracks.in'
!!!!!!        vn(1)=0.0
!!!!!!        vn(2)=Vcx
!!!!!!        vn(3)=Vcy
!!!!!!        ! ElemFractCode(elemId)=2
!!!!!!        ! ElemDecayed(elemId)=1
!!!!!!
!!!!!!        do j = 1, 3
!!!!!!          abc(362+j,elemId,1) = vn(j)
!!!!!!          abc(365+j,elemId,1) = vn(j)
!!!!!!          abc(368+j,elemId,1) = vn(j)
!!!!!!        end do
!!!!!!        !!--- this important b/c if the ElemFractCode(ele)=2, then it does not execute propagate for this element
!!!!!!        !!--- but it will execute overlap for intersection, so it needs the correct direction
!!!!!!        ! ncleave(2,elemId)=Vcx
!!!!!!        ! ncleave(3,elemId)=Vcy
!!!!!!        ! !!---- this to fill all the prop for the cracked elements
!!!!!!        ! elecrack(1,elemId)=ElemEdge1
!!!!!!        ! elecrack(2,elemId)=3
!!!!!!        ! elecrack(3,elemId)=ElemEdge2
!!!!!!        ! elecrack(4,elemId)=3
!!!!!!        !
!!!!!!        ! intersec(1, elemId)=rxc1
!!!!!!        ! intersec(2, elemId)=ryc1
!!!!!!        ! intersec(3, elemId)=rxc2
!!!!!!        ! intersec(4, elemId)=ryc2
!!!!!!        call pEC_ElemData(elemId)%SetDataPre(vn,ElemEdge1,rxc1,ryc1,r1,ElemEdge2,rxc2,ryc2,r2)
!!!!!!        ! write(*,*)'-----------------------------------'
!!!!!!        call pEC_ElemData(elemId)%PrintTest(elemId)
!!!!!!        ! if (i < FractElemCountPreCracked) then
!!!!!!        ! !-- update the neighbor elements
!!!!!!        !     call FracUpdateTipBeforeCrack(ix, elemId)
!!!!!!        ! else
!!!!!!        !     call FracUpdateTipAfterCrack(ix, elemId, ElemEdge2, 2)
!!!!!!        ! endif
!!!!!!
!!!!!!    end do
!!!!!!
!!!!!!    close(60)
!!!!!!    close(6011)
!!!!!!    write(*, *) '-->>>>>>>>FractReadApplyPreCrackCleavagePlanes'
!!!!!!    return
!!!!!!!    -------------------------------------------
!!!!!!!    -------------------------------------------
!!!!!!    90  IF (IERR .EQ. 29 ) THEN
!!!!!!        write(*, *) '------------------------ no input file .. reading cracks.in'
!!!!!!      ELSE
!!!!!!        write(*, *) 'Unrecoverable error in .. reading cracks.in, code =', IERR
!!!!!!        STOP
!!!!!!      END IF
!!!!!!
!!!!!!    return
!!!!!!END
!!###################################################################
!###################################################################
!###################################################################
subroutine FractReadApplyPreCrackCleavagePlanes(NodesCoordx,NodesCoordy,ElemConnect,ElemMaterial,DofIds &
                                                ,NodesDispl,usi  , freep , ym)
    !!(x     ,y     ,ix    , matp , id   , u    , usi  , freep , ym)
    use mod_file_units
    use mod_parameters
    use EC_Objects_manager
    implicit none

    real NodesCoordx(*), NodesCoordy(*),NodesDispl(*)
    real usi(*),freep(5,*), ym(4,*)
    INTEGER ElemConnect(4,*),ElemMaterial(*),DofIds(2,*)

    common /wblock8/ abc(573,nume,4),his(573,nume,4)
    real abc,his
    common/meshnum/ numnpo, numelto
    integer numnpo, numelto
    common/bk06/nprnt,mprint,itmpop,numelt,jprint,idump,locstr
    INTEGER nprnt,mprint,itmpop,numelt,jprint,idump,locstr

    integer*4 setvbuf3f_local
    integer i,j,elemId,IERR
    real Vcx,Vcy
    real vn(3)
    integer ElemEdge1,ElemEdge2
    real rxc1,ryc1,r1,rxc2,ryc2,r2
    character*200 strTitle
    real*8 rxc2Ratio,ryc2Ratio,rxc1Ratio,ryc1Ratio
    integer mEC_bCracking,mEC_ZoneFactor,mEC_DecayCount,mEC_PreCrackedElemCount
    real*8 x(4),y(4),Ptc(2)
    integer, ALLOCATABLE :: ElemIdsPreCracked(:)
    integer tempInt(2)

    write(*, *) '-->>>>>>>>FractReadApplyPreCrackCleavagePlanes'
    IERR=0
    i=0
    mEC_PreCrackedElemCount=0
    mEC_bCracking=0
    write(*, *) 'beginning of CrackProp ... reading cracks.in'

    open(60, file = 'frac_cracks.in', status = 'old',IOSTAT=IERR, ERR=90)
    read (60, *) mEC_bCracking !-- 1=allow fracture
    if(mEC_bCracking==0) then
      close(60)
      return
    endif
    open(6011, file = 'frac_cracks2.out', status = 'unknown')
    ierr=setvbuf3f_local(6011,1,100)

    write(6011, *) mEC_bCracking !-- 1=allow fracture

    read (60, *)    mEC_ZoneFactor !-- 20 * min elem diagonal
    write (6011, *) mEC_ZoneFactor !-- 20 * min elem diagonal

    read (60, *)    mEC_DecayCount !-- unloading steps
    write (6011, *) mEC_DecayCount !-- unloading steps

    read (60, *)    mEC_PreCrackedElemCount  !-- number of pre-cracked elems
    write (6011, *) mEC_PreCrackedElemCount  !-- number of pre-cracked elems
    !-- set the input data to the module
    call ECinitializeConsts(numnpo, numelto,ElemConnect,NodesCoordx,NodesCoordy, &
                           mEC_bCracking,mEC_ZoneFactor,mEC_DecayCount,mEC_PreCrackedElemCount)
    call EC_AllocateMem(ElemConnect,NodesCoordx,NodesCoordy)

!--------      this is has to be called first b/c the pre-crack needs it
    call FractReadElemNeighbors()
    !- write iFU_crackinput_out file to check the input data of the elem
    call EC_PrintTest()
    close(iFU_crackinput_out)
    call EC_CalcElemNeighbors(tempInt,0,1)
!!! if the fracture is not active then exit
    if (mEC_PreCrackedElemCount == 0) then
        close(60)
        close(6011)
        return
    endif

    write (6011, *) mEC_PreCrackedElemCount
    read(60, *) Vcx,Vcy
    write(6011, *) Vcx,Vcy
    vn(1)=0.0
    vn(2)=Vcx
    vn(3)=Vcy

    read(60, *) strTitle
    write(6011, *) strTitle !'!!!--------- id  ed1     rxc1    ryc1     r1  ed2     rxc2    ryc2     r2'
    allocate(ElemIdsPreCracked(mEC_PreCrackedElemCount))
    ElemIdsPreCracked=0
    do i = 1, mEC_PreCrackedElemCount
        !!!--------- id  ed1     rxc1    ryc1     r1  ed2     rxc2    ryc2     r2
        read(60, *) elemId,ElemEdge1,rxc1,ryc1,r1,ElemEdge2,rxc2,ryc2,r2,rxc2Ratio,ryc2Ratio,rxc1Ratio,ryc1Ratio
        write(6011, *) elemId,ElemEdge1,rxc1,ryc1,r1,ElemEdge2,rxc2,ryc2,r2,rxc2Ratio,ryc2Ratio,rxc1Ratio,ryc1Ratio
        write(*, *) 'Pre-exist crack, element', elemId, 'cracks.in'
        ElemIdsPreCracked(i)=elemId

        do j = 1, 3
          abc(362+j,elemId,1) = vn(j)
          abc(365+j,elemId,1) = vn(j)
          abc(368+j,elemId,1) = vn(j)
        end do

        call pEC_ElemData(elemId)%SetDataPre(vn,ElemEdge1,rxc1,ryc1,r1,ElemEdge2,rxc2,ryc2,r2)
        Ptc(1)=rxc1
        Ptc(2)=ryc1
        !- to get the elem nodes
        do j=1, 4
            x(j)=NodesCoordx(ElemConnect(j,elemId))
            y(j)=NodesCoordy(ElemConnect(j,elemId))
        enddo
        write(iFU_frac_crackFailedInfo,*) elemId,0,vn,'---pre-cracked'
        call EC_MarkElemForCrack(elemId,Ptc,x,y,0)
        !--- for debugging
        ! call pEC_ElemData(elemId)%PrintTest(elemId)

    end do
    do i = 1, mEC_PreCrackedElemCount

            CALL EC_SplitElem(ElemIdsPreCracked(i),NodesCoordx, NodesCoordy, ElemConnect, DofIds, NodesDispl, ElemMaterial, &
                            usi  , freep , ym,1)
            pEC_ElemData(ElemIdsPreCracked(i))%iElemSplit =1
    end do

    close(60)
    close(6011)
!!!!!!!!!!!================================================================
!!!!!!!!!!!================================================================
!!!!!!!!!!!========= update this global count b/c it's used in the rest of the code
!!!!!!!!!!!================================================================
!!!!!!!!!!!================================================================
!!!!!!!!    numelt=EC_ElemCountCurrent

!!    !---------------- debugging
!!    do i=7100,7106
!!        close(i)
!!    enddo
!!    stop
!!    !---------------- debugging
    write(*, *) '-->>>>>>>>FractReadApplyPreCrackCleavagePlanes'
    return
!    -------------------------------------------
!    -------------------------------------------
    90  IF (IERR .EQ. 29 ) THEN
        write(*, *) '------------------------ no input file .. reading cracks.in'
      ELSE
        write(*, *) 'Unrecoverable error in .. reading cracks.in, code =', IERR
        STOP
      END IF

    return
END

! !###################################################################
! !###################################################################
! !###################################################################
 subroutine FractReadElemNeighbors()
     use EC_Objects_manager
     implicit none

    integer i,IERR
    integer*4 setvbuf3f_local
    integer EC_ElemNeighbors(8)

    write(*, *) '-->>>>>>>>FractReadElemNeighbors -- start'
    IERR=0
    i=0
    write(*, *) 'beginning of CrackProp ... reading cracks.in'

    open(70, file = 'ElemsNeighbors.in', status = 'old',IOSTAT=IERR, ERR=90)
    open(7011, file = 'ElemsNeighbors2.out', status = 'unknown')
    ierr=setvbuf3f_local(7011,1,100)

    do i = 1, EC_ElemCountInput

      read (70, *)    EC_ElemNeighbors(1:8)
      write (7011, *)    EC_ElemNeighbors(1:8)
      call pEC_ElemData(i)%EC_SetElemEdgeNeighbors(EC_ElemNeighbors)
    end do

    close(70)
    close(7011)
    write(*, *) '-->>>>>>>>FractReadElemNeighbors -- done'
    return
!    -------------------------------------------
!    -------------------------------------------
    90  IF (IERR .EQ. 29 ) THEN
        write(*, *) '------------------------ no input file .. reading cracks.in'
      ELSE
        write(*, *) 'Unrecoverable error in .. reading cracks.in, code =', IERR
        STOP
      END IF

    return
 END

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


program main

    implicit none

    interface

        subroutine sub_lib_in_lib()
        implicit none
        end subroutine sub_lib_in_lib

    end interface

    call test_lib()

    call sub_lib_in_lib()

    write( *,* ) '8. Done.'

end program main


text:
program test_subscript
use iso_fortran_env, only:int8, int16, int32, int64, real32, real64, real128
use mod_testing, only:assert, initialize_tests, report_tests
use mod_functional

implicit none

logical, dimension(:), allocatable :: tests
logical :: test_failed
integer :: n, norder, ntests
integer, parameter :: stdout = 6

n = 1
ntests = 11
call initialize_tests(tests, ntests)

tests(n) = assert(all(subscript([1_int8, 2_int8, 3_int8], [2_int8]) == [2_int8]), &
                  'subscript,  int8')
n = n + 1

tests(n) = assert(all(subscript([1_int16, 2_int16, 3_int16], [2_int16]) == [2_int16]), &
                  'subscript,  int16')
n = n + 1

tests(n) = assert(all(subscript([1_int32, 2_int32, 3_int32], [2_int32]) == [2_int32]), &
                  'subscript,  int32')
n = n + 1

tests(n) = assert(all(subscript([1_int64, 2_int64, 3_int64], [2_int64]) == [2_int64]), &
                  'subscript,  int64')
n = n + 1

tests(n) = assert(all(subscript([1._real32, 2._real32, 3._real32], [2]) == [2._real32]), &
                 'subscript,  real32')
n = n + 1

tests(n) = assert(all(subscript([1._real64, 2._real64, 3._real64], [2]) == [2._real64]), &
                  'subscript,  real64')
n = n + 1

tests(n) = assert(all(subscript([1._real128, 2._real128, 3._real128], [2]) == [2._real128]), &
                  'subscript,  real128')
n = n + 1

tests(n) = assert(size(subscript([1, 2, 3], [0])) == 0, &
                  'subscript out of bounds returns empty array')
n = n + 1

tests(n) = assert(all(subscript(arange(cmplx(1._real32, 0._real32), &
                                       cmplx(3._real32, 0._real32)), [2])&
                       == [cmplx(2._real32, 0._real32)]), &
                  'subscript,  complex real32')
n = n + 1

tests(n) = assert(all(subscript(arange(cmplx(1._real64, 0._real64), &
                                       cmplx(3._real64, 0._real64)), [2])&
                       == [cmplx(2._real64, 0._real64)]), &
                  'subscript,  complex real64')
n = n + 1

tests(n) = assert(all(subscript(arange(cmplx(1._real128, 0._real128), &
                                       cmplx(3._real128, 0._real128)), [2])&
                       == [cmplx(2._real128, 0._real128)]), &
                  'subscript,  complex real128')
n = n + 1

test_failed = .false.
call report_tests(tests, test_failed)
if(test_failed)stop 1

end program test_subscript

max_stars_repo_name: econwang/functional-fortran

max_issues_repo_name: econwang/functional-fortran

max_forks_repo_name: econwang/functional-fortran


text:
image(newuc.jpg, left, thumbnail, 250, Comic Sans added for emphasis.)We have a Seal of the University, which is pretty good.  We have a UC Davis Logo, which is passable and does its job well.  Theres the Aggie Logo which is pretty decent, and the Oldschool Aggie Logo which is sheer awesome.  This page is not about any of those.

This page is about the abomination of a new logo that the University of California system created by committee (which is always the optimal way to make art amirite?), and released in December 2012 to give the UC System a new feel.  You know, like the University of Phoenix brochures that you get in the mail because the marketers dont realize that you got into an actual school.

The UCOP wanted a logo for the web designers, because apparently they dont realize that everyone will have Retina displays in five years anyway and will actually be able to use the Seal online.  But hey, its a welcome distraction from the UC Budget, right?

It lasted for http://www.sfgate.com/bayarea/article/UCsuspendsuseoflogosaysitwasdistracting4118987.php only days before it was suspended from use.  It was simultaneously being defended as being only a monogram and having been developed in house, so there were no heavy expenditures.  This defensiveness was largely due to heavy criticism from students and alumni.  A petition on change.org collected http://www.change.org/petitions/universityofcaliforniastopthenewuclogo 54,000 signatures against the new logo.  Even Lt. Gov. Gavin Newsom (and UC Regents UC Regent) publicly urged the UC to reject the logo.

 Things the new logo looks like:
  A programloading circle on Windows.
  A toilet with your undergraduate student fees student fees swirling down the drain.
  The quad awash in pepper spray
  Sort of like the letters U and C if you squint a lot.
  Incompetence at UCOP
  A wiki:wikipedia:heater shield SCA bearing the arms of the Markgraf of Oakland: azure, charged with annulet dor, incomplete sinister.

## Source of the animation: http://www.reddit.com/r/UCDavis/comments/14ng1h/

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program main
  use fygg

  character(len=32) :: arg
  integer :: iterations
  type(yggcomm) :: ymlfile
  type(yggcomm) :: rpc
  type(yggcomm) :: log
  logical :: ret
  integer, parameter :: msg_len = 2048  ! YGG_MSG_MAX
  character(len=msg_len) :: ycontent
  integer :: fib, fibNo, i, nlines, count_lines
  character(len=msg_len) :: logmsg
  integer :: ycontent_siz, logmsg_siz
  character(len=2) :: fibStr

  call get_command_argument(1, arg)
  read(arg, *) iterations
  write(*, '("Hello from Fortran rpcFibCli: iterations ",i2)') iterations

  ! Set up connections matching yaml
  ! RPC client-side connection will be $(server_name)_$(client_name)
  ymlfile = ygg_input("yaml_in")
  rpc = ygg_rpc_client("rpcFibSrv_rpcFibCli", "%d", "%d %d")
  log = ygg_output("output_log")

  ! Read entire contents of yaml
  ycontent_siz = len(ycontent)
  ret = ygg_recv(ymlfile, ycontent, ycontent_siz)
  if (.not.ret) then
     write(*, '("rpcFibCli(F): RECV ERROR")')
     stop 1
  end if
  nlines = count_lines(ycontent, new_line('A'))
  write(*, '("rpcFibCli: yaml has ",i4," lines")') nlines + 1
  ret = ygg_recv(ymlfile, ycontent, ycontent_siz)

  fib = -1
  fibNo = -1
  do i = 1, iterations

     ! Call the server and receive response
     write(*, '("rpcFibCli(F): fib(->",i2,") ::: ")') i
     ret = ygg_rpc_call(rpc, yggarg(i), [yggarg(fibNo), yggarg(fib)])
     if (.not.ret) then
        write(*, '("rpcFibCli(F): RPC CALL ERROR")')
        stop 1
     end if

     ! Log result by sending it to the log connection
     write(fibStr, '(i2)') fib
     write(logmsg, '("fib(",i2,"<-) = ",A,"<-",A)') fibNo, &
          adjustl(fibStr), new_line('A')
     write(*, '(A)', advance="no") trim(logmsg)
     logmsg_siz = len_trim(logmsg)
     ret = ygg_send(log, logmsg, logmsg_siz)
     if (.not.ret) then
        write(*, '("rpcFibCli(F): SEND ERROR")')
        stop 1
     end if

  end do

  write(*, '("Goodbye from Fortran rpcFibCli")')

end program main

function count_lines(str, substr) result(c)
  implicit none
  character(len=*) :: str
  character(len=*) :: substr
  integer :: c
  integer :: k, inc, idx
  c = 0
  k = 1
  inc = len_trim(substr)
  idx = index(str(k:), substr)
  do while (idx.ne.0)
     c = c + 1
     k = k + inc
     idx = index(str(k:), substr)
  end do
end function count_lines

max_stars_repo_name: Xyzic/yggdrasil

max_issues_repo_name: Xyzic/yggdrasil

max_forks_repo_name: Xyzic/yggdrasil


text:
c $VERSION "08/16/95 @(#)culad.f	7.1"
c     subroutine difint
c above statement was commented and following added by Chen, 8/31/89.
      subroutine difint(nohrs,clump,sunazm,path)
c **********************************************************************
c
c	---------
c	| Added path, path length through the canopy.
c	| this was added becuase now we take into account
c	| slope and aspect.  LMM 30/9/93
c
      parameter(mh=98)
      dimension sunazm(mh),path(mh)
c above two statements were added by Chen, 8/31/89.
      dimension sdnsub(9,20),supsub(9,20),rsubl(9,20),tsubl(9,20),
     &adum(20)
c     ---------
c     | rad1 cmn added rlleaf, tlleaf, rldead and tldead LMM 94/9/6
      common /rad1/emis,emisol,rsoil(3),rleaf(3,20),tleaf(3,20),
     &aleaf(3,20),expdif,rlayr(3,20),tlayr(3,20),rllive(3),tllive(3),
     &rldead(3),tldead(3)
      common /rad3/radtop(3,mh),fbeam1(3,mh),coszen(mh),zenang(mh),hfday
     &,ratiod,ration,ratio(mh)
      common /rad4/d(3,20),u(3,20),bmflx(3,20),rnet(50),rndiv(50),tsfc
     1,rnlam(3,20)
      common /rad5/ sourup(3,20,mh),sourdn(3,20,mh)
      common/rad6/rlfdif(3,20),tlfdif(3,20),rlfdir(3,20,mh),
     &tlfdir(3,20,mh)
      common/indax2/xmeu,xneu,gamrat,frdeg(91),iangot(9),kpntsp
      common/indax3/xmeuaz,xneuaz,beta0,nbeta,fraz(50),fradeg(360)
      common/misc1/pi,pid180,pid2,sigma,iwrite(9,99),kmax
      common/misc2/itot,itotp1,jtot,fr(10),ct(10),totlai,df,dt
     &,clai(20),distls(10,mh),jdead
      common/misc6/xint(mh),nxintz,xintz(90),nalpha,ispher
c xintz, extin. coef. at different zenith averaged over azimuth, and
c nxtintz, dimension of xintz, Chen, 05/25/89.
c
c  calculate diffuse interception for layer.  calc layer trans and refl
c    for thick layers using layer equations and dividing df layer into
c    10 sub-layers and solving for diffuse radiation.
c    should use good precision in diffuse integral so errors dont
c    compound. thus use 90 intervals because 9 intervals accurate
c    to about .005 and this error becomes .03 in layer non-interception
c    factor because of 10th power from 10 sublayers.
      laysub=10
      laysp1=laysub+1
c      da=pi/180.
c      a=pi/360.
c above two stms were commented and following one added
c by Chen, 05/24/89.
      wz=pi/2/nxintz
      x=0.
c      do 500 ii=1,90
c above stm was commented by Chen, 05/24/89.
      do 500 ii=1,nxintz
      a=wz*(ii-.5)
      ca=cos(a)
      sa=sin(a)
      dfdca=df/ca
c      if(ispher.eq.0)call simpsn(a,fr,nalpha,xint)
c      x=x+ca*sa*exp(-xint*dfdca)
c      a=a+da
c above three were commented and following one added by Chen, 05/24/89.
      x=x+ca*sa*exp(-clump*xintz(ii)*dfdca)
c parameter clump in above statement was added by Chen, 9/4/89.
 500  continue
c      expdif=2.*x*da
c above one was commented and following one added by Chen, 05/24/89.
      expdif=2.*x*wz
c
c  calc diffuse non-interception for a sublayer from non-inter
c  for a df layer because integral takes care of ang dist thru
c  the whole df layer and sublayers assume isotropy of inc rad.
      difsub=expdif**.1
      do 601 k=1,kmax
      do 601 j=2,jtot
        rlayr(k,j)=(1.-expdif)*rlfdif(k,j)
        tlayr(k,j)=(1.-expdif)*tlfdif(k,j)+expdif
        rsubl(k,j)=(1.-difsub)*rlfdif(k,j)
        tsubl(k,j)=(1.-difsub)*tlfdif(k,j)+difsub
c original rlfdif(k,nohrs) & tlfdif(k,nohrs) were changed to rlfdif(k)
c and tlfdif(k) in above four statements by Chen, 8/31/89.
c
c j subscripts added to rlayr,tlayr,rsubl,tsubl,rlfdif,tlfdif to
c allow for dead veg to accumulate in lower layers - MCA 4/21/95
c
c       write(26,*)rlayr(k),tlayr(k),rsubl(k),tsubl(k)
c	write(6,*)'k,rleaf,rlfdif,rlayr=',k,rleaf(k),rlfdif(k),rlayr(k)
c	write(6,*)'k,tleaf,tlfdif,tlayr=',k,tleaf(k),tlfdif(k),tlayr(k)
 601  continue
c
c  Solve sublayer equations
c  Here, jl is index for main canopy layers - MCA 4/21/95
      do 620 jl=2,jtot
      do 620 k=1,kmax
        d(k,laysp1)=1.
        adum(1)=0.
        tlay2=tsubl(k,jl)*tsubl(k,jl)
        do 605 j=2,laysp1
          jm1=j-1
          adum(j)=adum(jm1)*tlay2/(1.-adum(jm1)*rsubl(k,jl)) +
     &            rsubl(k,jl)
 605    continue
        do 610 j=2,laysp1
          jj=laysp1-j+1
          jjp1=jj+1
          d(k,jj)=d(k,jjp1)*tsubl(k,jl)/(1.-adum(jj)*rsubl(k,jl))
          u(k,jjp1)=adum(jjp1)*d(k,jjp1)
c         write(26,*)jjp1,adum(jjp1),d(k,jjp1),u(k,jjp1)
 610    continue
c       write(19,*)itry,rlayr(k,jl),u(k,11)
 615    rlayr(k,jl)=u(k,laysp1)
        tlayr(k,jl)=d(k,1)
c       write(26,*)itry,k,j,rlayr(k),tlayr(k),rsubl(k),tsubl(k),difsub
c    1  ,expdif
 620  continue
c
c  diffuse should be evaluated at layer midpoints
c  calc thick layer sources from laysub sublayers. since source dist
c  depends on sun angle, this effective source is dep on hour and
c  wavelength.
c
c  zero out all effective sources.
      do 627 ihr=1,nohrs
      do 627 j=1,jtot
      do 627 k=1,kmax
        sourup(k,j,ihr)=0.
 627    sourdn(k,j,ihr)=0.
c  Subscript j added to sourup and sourdn to allow varying R and T
c  properties with layer - MCA 4/21/95
c
      do 690 ihr=1,nohrs
        if (coszen(ihr).lt..01) go to 690
        if (path(ihr).lt..9) go to 690
c6789  if(ispher.eq.0)call simpsn(zenang(ihr),fr,nalpha,xint)
c above stm was commented and following two added by Chen, 05/24/89.
c	xint(ihr)=.5
c 	if(ispher.eq.0)call simpsn(zenang(ihr),sunazm(ihr),fr,
c     1 nalpha,xint(ihr))
c	----------
c	| added path to dirsub equation.  This takes into account
c	| slope and aspect.
c	| 	LMM 30/9/93
        dirsub=exp(-clump*xint(ihr)*df*path(ihr)/laysub)
c parameter clump in above statement was added by Chen, 9/4/89.
c path replaces 1/coszen in above to take into account slope and aspect
c	LMM 30/9/93
c
        do 691 k=1,kmax
        do 691 jl=2,jtot
c original 690 for k loop was changed to 691 since the stm
c if... goto 690 caused trouble. Chen, 9/7/89.
          if(fbeam1(k,ihr).lt.0.01)go to 690
          sdnsub(k,laysp1)=(1.-dirsub)*tlfdir(k,jl,ihr)
          supsub(k,laysp1)=(1.-dirsub)*rlfdir(k,jl,ihr)
          do 630 j=2,laysub
            jj=laysp1-j+1
            jjp1=jj+1
            sdnsub(k,jj)=sdnsub(k,jjp1)*dirsub
            supsub(k,jj)=supsub(k,jjp1)*dirsub
c           write(26,*)jj,sdnsub(k,jj),supsub(k,jj),dirsub
 630      continue
c  zero all u(k,j) and d(k,j).  d(k,laysp1)=u(k,1)=0. are b.c.
          iter=0
          do 640 j=1,laysp1
            d(k,j)=0.
 640        u(k,j)=0.
 645      iter=iter+1
          irept=0
          do 650 j=2,laysp1
            jj=laysp1-j+1
            jjp1=jj+1
            down=tsubl(k,jl)*d(k,jjp1)+u(k,jj)*rsubl(k,jl)+
     &           sdnsub(k,jjp1)
            if(abs(down-d(k,jj))-.0001)646,646,644
 644        irept=1
 646        d(k,jj)=down
            up=tsubl(k,jl)*u(k,jj)+d(k,jjp1)*rsubl(k,jl)+
     &           supsub(k,jjp1)
            if(abs(up-u(k,jjp1))-.0001)649,648,648
 648        irept=1
 649        u(k,jjp1)=up
 650      continue
          if(irept.ne.0)go to 645
c  calc source terms to be used in radiat, these must be multiplied
c  tbeam above each layer to get the sdn(j) and sup(j) needed there.
          sourup(k,jl,ihr)=u(k,laysp1)
          sourdn(k,jl,ihr)=d(k,1)
c         write(26,*)k,ihr,sourup(k,jl,ihr),sourdn(k,jl,ihr)
 691    continue
c above stm was added by Chen, 9/7/89.
 690  continue
      return
      end
c
c      subroutine dstlit(thets,phis,itot,kmax,rlfdir,
c     &rlfdif,tlfdir,tlfdif,distls,ihr)
c above statement was commented and following added by Chen, 8/31/89.
      subroutine dstlit(thets,phis,kmax,rlfdir,tlfdir,ihr)
c  subroutine to calculate the distribution of leaf-normal to sun angles
c    for a given direction of the sun or given hour. distls(iangle) is
c    the variable where iangle has the same class intervals as fr(ialpha
c      dimension fr(10),distls(10),rlfdir(9,24)
c above was commented and following added by Chen, 5/24/89.
c     dimension rlfdir(9,24)
c    &,rlfdif(9,24),tlfdir(9,24),tlfdif(9,24),ca(9),sa(9)
c above statement was commented and following two added, Chen, 8/31/89.
      parameter(mh=98)
      dimension rlfdir(3,20,mh),tlfdir(3,20,mh),ca(9),sa(9)
      common/rad8/rlfhem(3,9,20),tlfhem(3,9,20)
      common/indax3/xmeuaz,xneuaz,beta0,nbeta,fraz(50),fradeg(360)
      common/misc2/itot,itotp1,jtot,fr(10),cost(10),totlai,df,dt
     &,clai(20),distls(10,mh),jdead
c original ct(10) was changed to cost(10) by Chen, 7/5/89.
      common/misc6/xint(mh),nxintz,xintz(90),nalpha,ispher
c xintz, extin. coef. at different zenith averaged over azimuth, and
c nxtintz, dimension of xintz, Chen, 05/25/89.
c      write(6,*)'ihr,thets,phis in dstlit=',ihr,thets,phis
      pi=3.1415926537
      pid180=pi/180.
      pid2=pi/2.
      do40ialpha=1,itot
 40   distls(ialpha,ihr)=0.
c original distls(ialpha) in above statement was changed to
c distls(ialpha,ihr) by Chen, 8/30/89.
	wbeta=2.*pi/nbeta
      do100ialpha=1,itot
      alpha=pid2*(ialpha-1)/(itot) + pi/(4.*itot)
      sa(ialpha)=sin(alpha)
      ct=cos(thets)
      st=sin(thets)
      ca(ialpha)=cos(alpha)
c      fradd=fr(ialpha)/180
c	do100ibeta=1,20
c      beta=(1.+(ibeta-1)*2.)*pid180 - pi (5/18/89, chen)
c above three were commented and following three added by Chen,5/24/89
	do 100 ibeta=1,nbeta
	beta=wbeta*(ibeta-.5)
	fradd=fr(ialpha)*fraz(ibeta)
      fs=ct*ca(ialpha) + st*sa(ialpha)*cos(beta-phis)
      angl=acos(abs(fs))
      iangl=angl*itot/pid2 + 1.
c	write(6,*)'ialpha,ibeta,angl,iangl=',ialpha,ibeta,angl,iangl
c     write(2,50)alpha,beta,thets,phis,iangl,angl,(distls(m),m=1,itot)
 50   format(4f4.1,i3,f5.2,9f5.3)
      if(iangl.lt.1.or.iangl.gt.itot)go to 80
      go to 95
 80   write(6,90)iangl
 90   format(' trouble in sub dstlit iangl= ',i6)
      stop
 95   distls(iangl,ihr)=distls(iangl,ihr) + fradd
c original distls(iangl) in above statement was changed to
c distls(iangl,ihr) by Chen, 8/30/89.
100   continue
	sum=0.
	do 105 iangle=1,itot
	  sum=sum+cost(iangle)*distls(iangle,ihr)
105	continue
	xint(ihr)=sum
c      write(2,51)thets/pid180,phis/pid180,(distls(m,ihr),m=1,itot)
c original distls(m) in above statement was changed to
c distls(m,ihr) by Chen, 8/30/89.
 51   format(2(f7.2,1x),9(f5.3,1x))
c  compute a mean weighted sunlit leaf reflectance
      do400 j=2,jtot
      do400 k=1,kmax
      sum1=0.
      sum=0.
c     sum2=0.
c     sum3=0.
c     sum4=0.
c above three statements were commented by Chen, 8/31/89.
      do300iangl=1,itot
      angl=pid2*(iangl-1)/itot + pid2/(2.*itot)
      sum=sum+distls(iangl,ihr)*rlfhem(k,iangl,j)
      sum1=sum1+distls(iangl,ihr)*tlfhem(k,iangl,j)
c original distls(iangl) in above two statements was changed to
c distls(iangl,ihr) by Chen, 8/30/89.
c  calc sums to get weighted leaf refl and trans for diffuse inc rad.
c     sum2=sum2+sa(iangl)*ca(iangl)
c     sum3=sum3+rlfhem(k,iangl,j)*sa(iangl)*ca(iangl)
c     sum4=sum4+tlfhem(k,iangl,j)*sa(iangl)*ca(iangl)
c above three statements were commented by Chen, 8/31/89.
 300  continue
      rlfdir(k,j,ihr)=sum
      tlfdir(k,j,ihr)=sum1
c     rlfdif(k,ihr)=sum3/sum2
c     tlfdif(k,ihr)=sum4/sum2
c above two statements were commented by Chen, 8/31/89, since
c rlfdif(k) & tlfdif(k) were calculated in subroutine inplnt.
 400  continue
      return
      end
c
      subroutine gamma(x,gam)
c this subroutine calculates the gamma function using stirling's
c approximation.
      real x,y,z,pi,gam
      integer i,j
      g(y) = sqrt(2.*pi/y)*
     1exp(y*alog(y)+(1-1/(30*y*y))/(12*y)-y)
     2/((y-2)*(y-1))
      pi=4.*atan(1.)
      if (x.lt.0.0) goto 10
      if (abs(x).lt.1.e-10) goto 40
      xp2=x+2.0
      gam = g(xp2)
      goto 99
  10  z=x
      j=-1
      y=x
  20  j=j+1
      y=y+1.0
      if(y.lt.0.0) goto 20
      yp2=y+2.0
      gam=g(yp2)
      jp1=j+1
      do 30 i=1,jp1
  30  gam=gam/(x+i-1)
      goto 50
  40  gam=1.0
  50  continue
  99  return
      end
c
c      subroutine kernal(cosdel,theta,alpha)
c  subroutine to calculate cosine delta where theta is solar zenith
c    angle and alpha is inclination angle of leaf.
c      pid2=3.1415926537/2.
c      if(theta+alpha-pid2)10,10,50
c10    cosdel=cos(alpha)*cos(theta)
c      go to 100
c 50   beta0=acos(1./(tan(alpha)*tan(theta)))
c      cosdel=((pid2-beta0)*cos(alpha)*cos(theta)+sin(alpha)*sin(theta)*s
c     1in(beta0))/pid2
c 100  continue
c      return
c      end
c
c above subroutine kernal were commented and following one added by
c Chen, 05/24/89.
c
      subroutine kernal(cosdel,theta,phi,alpha)
c  subroutine to calculate cosine delta where
c  theta=zenith angle of sun or viewer
c  phi =azimuth angle of sun or viewer, south = 0., east = -90.
c  alpha=leaf elevation angle from horizontal
c  beta =azimuthal direction leaf normal points
c  fr =fraction of leaf area proj either in dir of sun or viewer.
c  fraz is the fraction of leaf area as a func of azimuth
      common/indax3/xmeuaz,xneuaz,beta0,nbeta,fraz(50),fradeg(360)
      pi=4.*atan(1.)
      wbeta=2.*pi/nbeta
      cosdel=0.
      do1000ibeta=1,nbeta
      beta=(ibeta-.5)*wbeta
      fs=cos(theta)*cos(alpha)+sin(theta)*sin(alpha)*cos(phi-beta)
      cosdel=cosdel+abs(fs)*fraz(ibeta)
 1000 continue
      return
      end
      subroutine lad(meu,neu,fr,itot,gamrat,frdeg)
c this subroutine calculates the leaf angle distribution
c                                                     4/20/84
c including distribution both in zenith and in azimuth, Chen, 05/24/89.
c leaf angle distributions for various canopy types:
c   canopy type       f(theta)         theta     meu     neu
c   planophile   2(1+cos(2*theta))/pi  26.76   2.770   1.172
c   erectophile  2(1-cos(2*theta))/pi  63.24   1.172   2.770
c   plagiophile  2(1-cos(4*theta))/pi  45.00   3.326   3.326
c   extremophile 2(1+cos(4*theta))/pi  45.00   0.433   0.433
c   uniform      2/pi                  45.00   1.000   1.000
c   spherical    sin(theta)            57.30   1.101   1.930
c      source: naren goel & don strebel (1984)  agron j. 76:800-802
c
      common/indax3/xmeuaz,xneuaz,beta0,nbeta,fraz(50),fradeg(360)
c above common is added by Chen, 05/23/89, and beta0 is in degrees.
      real meu,neu,meuneu,gam1,gam2,gam12,gamrat,lang,fr(10),
     1frdum(91),frdeg(91),frcls(9)
c
c  beta-distribution
c
c      bdistr(gamrat,tld90,neu,meu)=gamrat*(tld90**(neu-1.))*
c     &       ((1.-tld90)**(meu-1.))/(360.*90.)
c above was commented by Chen, 7/5/89.
      pi=4.*atan(1.)
      pid180=pi/180.
      meuneu=meu+neu
c the next statement is used to avoid overflows in the calculation
c of gamma functions for large values of meu and neu.
c      if(meu.gt.35..or.neu.gt.35..or.meuneu.gt.35.) goto 1
c      call gamma(meu,gam1)
c      call gamma(neu,gam2)
c      call gamma(meuneu,gam12)
c      goto 4
c 1    gamrat=((meu*neu)/(2.*pi*meuneu))**0.5
c      gamrat=gamrat*((meuneu/meu)**meu)*(meuneu/neu)**neu
c      goto 5
c 4    gamrat=gam12/(gam1*gam2)
c above 9 lines were commented by Chen, 7/5/89.
      sum = 0.
      dang=90./itot
 5    do 10 i = 1,itot
         lang=(dang*(i-1)+dang/2)/90.
c      fr(i)=(gamrat*(lang**(neu-1))*((1-lang)**(meu-1))/90.)*dang
c above was commented and following added by Chen, 7/5/89.
	fr(i)=lang**(neu-1)*(1-lang)**(meu-1)
c      write(6,*)'in lad ',i,fr(i),gamrat,itot,lang,neu,meu,dang
      sum = sum + fr(i)
 10   continue
      do 20 i = 1,itot
  20     fr(i)=fr(i)/sum
c ****** 7/7/86 revision  ******
c **  removed seciont that calculates frspec - it's in a separate
c **  subroutine called "fractn"
c ************************************
c  calculates leaf angle distribution for specular contribution at
c 10 degree increments
c
c     sum = 0
c     do 50 i = 1,9
c        lang=(10.*(i-1)+5)/90.
c        frdum(i) = gamrat*(lang**(neu-1))*((1-lang)**(meu-1))/9.
c     sum=sum+frdum(i)
c 50  continue
c     do 70 i=1,9
c        frdum(i) = frdum(i)/sum
c        do 70 j = 1,36
c  frspec is fraction of leaf area in single azimuth and zenith angle
c    class that contributes to specular reflection.
c           frspec(i,j)=frdum(i)/36.
c 70  continue
c     write(6,*)(frspec(i,1),i=1,9)
c
c  calculate for every 10 deg view zenith angle class
      do 75 i=1,91
         frdum(i)=0.
  75  continue
c for deg classes 1 through 91
      do 200 ideg=1,90
         deg=float(ideg)
         tld90=(deg-0.5)/90.
c         frdum(ideg)=bdistr(gamrat,tld90,neu,meu)
c above was commented and following added by Chen, 7/5/89.
	 frdum(ideg)=tld90**(neu-1)*(1-tld90)**(meu-1)
  200 continue
c
c  normalize fraction of leaves
c
      sum=0.
      do 250 i=1,91
         sum=sum+frdum(i)
 250  continue
      do 275 i=1,91
         frdeg(i)=frdum(i)/sum
 275  continue
c  put into 9 classes for screen display
      do 300 i=1,9
         frcls(i)=0.
 300  continue
      do 305 i=2,10
 305     frcls(1)=frcls(1)+frdeg(i)
      do 310 i=2,9
         j=i-1
         frtemp=0.
         do 320 k=1,10
            kk=(k+j*10)
            frtemp=frtemp+frdeg(kk)
 320     continue
         frcls(i)=frtemp
 310  continue
c
c following 12 statments for calculating leaf azimuthal angle distrb.
c were added by Chen, 05/23/89.
c	write(6,*)'starting calculation of azimuthal angle distribution'
      jbeta0=int(beta0/360.*nbeta)
      sum=0.
      do 410 ibeta=1,nbeta
      x=(ibeta-.5)/nbeta
      aux=x**(xneuaz-1)*(1-x)**(xmeuaz-1)
      ii=ibeta+jbeta0
      if (ii.gt.nbeta) ii=ii-nbeta
      fraz(ii)=aux
      sum=sum+aux
 410  continue
      do 420 ii=1,nbeta
      fraz(ii)=fraz(ii)/sum
 420  continue
c	write(6,*) (fraz(i),i=1,nbeta)
c following 13 statements for calculating leaf azimuthal angle distrb.
c in each degree class were added by Chen, 07/11/89.
      jbeta0=int(beta0)
      sum=0.
      do 411 ibeta=1,360
      x=(ibeta-.5)/360
      aux=x**(xneuaz-1)*(1-x)**(xmeuaz-1)
      ii=ibeta+jbeta0
      if (ii.gt.360) ii=ii-nbeta
      fradeg(ii)=aux
      sum=sum+aux
 411  continue
      do 421 ii=1,360
      fradeg(ii)=fradeg(ii)/sum
 421  continue
c
      return
      end
c**************************simpsn************************************
      subroutine simpsn(theta,phi,gr,nalpha,xint)
c phi in above and following common were adde by Chen, 05/24/89.
      common/indax3/xmeuaz,xneuaz,beta0,nbeta,fraz(50),fradeg(360)
      dimension gr(10)
c  xint  extinct coeff-nel is in first part of array and nthet is
c          in last part of array and tot size of array is nelpth.
c  nalpha  no of leaf angle classes
c
c   program to calc integrl of ga*cosd (simpson rule)
c     write(11,7)
 7    format(' ','*extinction data')
c
c     nthet=num. of sun angles
c     dt= delta theta
c     nalpha=num. of leaf angle classes
c     da= delta alpha
c
      pi=4.0*atan(1.0e0)
      pid180=pi/180.0
      pid2=pi/2.0
c     nalpha must be odd
      xint=0.
      st=sin(theta)
      ct=cos(theta)
      sum=0.
      do 1001 k=1,nalpha
      alpha=(k-.5)*pid2/nalpha
c      call kernal(cosdel,theta,alpha)
c  above was commented and following added by Chen, 7/9/89.
      call kernal(cosdel,theta,phi,alpha)
      sum=sum+cosdel*gr(k)
c      write(2,*)'ialpha,fr,cosdel,xint=',k,gr(k),cosdel,sum
 1001 continue
      xint=sum*pid2/nalpha
c  for area fractions weighted by leaf angle distrib,gr(alpha),
c    should not mult by dalpha because integral over gr(alpha)dalpha
c    in denominator will cancel dalpha.
      xint=sum
 1000 continue
c     write(11,1003)
 1003 format(' ','theta and xint')
c     write(11,500) theta,xint
  500 format(' ',f10.3,2(5x,f16.12))
      return
      end

max_stars_repo_name: hieulel/CUPID

max_issues_repo_name: hieulel/CUPID

max_forks_repo_name: hieulel/CUPID


text:
SUBROUTINE bisect(d,e,acheps,ifail)
!
! This subroutine finds the eigenvalues of a tridiagonal matrix,
! given with its diagonal elements in the array d(n) and
! its subdiagonal elements in the last n - 1 stores of the
! array e(n), using ql transformations. The eigenvalues are
! overwritten on the diagonal elements in the array d in
! ascending order. The subroutine will fail if any one
! eigenvalue takes more than 30 iterations.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::acheps
 REAL(iwp),INTENT(IN OUT)::d(0:),e(0:)
 INTEGER,INTENT(IN OUT)::ifail
 INTEGER::m,n,i,l,j,i1,m1,ii,aux
 REAL(iwp)::b,f,h,g,p,r,c,s,zero=0.0_iwp,pt5=0.5_iwp,one=1.0_iwp,         &
   two=2.0_iwp
 n=UBOUND(d,1)
 IF(n/=1)THEN
   DO i=2,n
     e(i-1)=e(i)
   END DO
 END IF
 e(n)=zero
 b=zero
 f=zero
 loop_1: DO l=1,n
   j=0
   h=acheps*(ABS(d(l))+ABS(e(l)))
   IF(b<h)b=h
   ! look for small sub diagonal element
   DO  m=l,n
     IF(ABS(e(m))<=b)EXIT
   END DO
   IF(m/=l)THEN
     loop_2: DO
       IF(j==30)THEN
         ifail=1
         EXIT loop_1
       END IF
       j = j + 1
       ! form shift
       g = d(l)
       h = d(l+1) - g
       IF(ABS(h)<ABS(e(l)))THEN
         p=h*pt5/e(l)
         r=SQRT(p*p+one)
         h=p+r
         IF(p<zero)h=p-r
         d(l)=e(l)/h
       ELSE
         p=two*e(l)/h
         r=SQRT(p*p+one)
         d(l)=e(l)*p/(one+r)
       END IF
       h=g-d(l)
       i1=l+1
       IF(i1<=n)THEN
         DO i=i1,n
           d(i)=d(i)-h
         END DO
       END IF
       f=f+h
       ! ql transformation
       p=d(m)
       c=one
       s=zero
       m1=m-1
       loop_4: DO ii=l,m1
         i=m1-ii+l
         g=c*e(i)
         h=c*p
         IF(ABS(p)>=ABS(e(i)))THEN
           c= e(i)/p
           r=SQRT(c*c+one)
           e(i+1)=s*p*r
           s=c/r
           c=one/r
         ELSE
           c=p/e(i)
           r=SQRT(c*c+one)
           e(i+1)=s*e(i)*r
           s=one/r
           c=c/r
         END IF
         p=c*d(i)-s*g
         d(i+1)=h+s*(c*g+s*d(i))
       END DO loop_4
       e(l)=s*p
       d(l)=c*p
       IF(ABS(e(l))<=b)EXIT loop_2
     END DO loop_2
   END IF
   p=d(l)+f
   ! order eigenvalue
   aux=0
   IF(l/=1)THEN
     loop_3: DO ii=2,l
       i=l-ii+2
       IF(p>=d(i-1))THEN
         aux=1
         EXIT loop_3
       END IF
       d(i)=d(i-1)
     END DO loop_3
   END IF
   IF(aux==0)THEN
     i=1
   END IF
   d(i) = p
   ifail=0
 END DO loop_1
RETURN
END SUBROUTINE bisect

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


MODULE zdfevd
  USE oce
  USE dom_oce
  USE zdf_oce
  USE trd_oce
  USE trdtra
  USE in_out_manager
  USE iom
  USE lbclnk
  USE timing
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: zdf_evd
  CONTAINS
  SUBROUTINE zdf_evd(kt, p_avm, p_avt)
    INTEGER, INTENT(IN   ) :: kt
    REAL(KIND = wp), DIMENSION(:, :, :), INTENT(INOUT) :: p_avm, p_avt
    INTEGER :: ji, jj, jk
    REAL(KIND = wp), DIMENSION(jpi, jpj, jpk) :: zavt_evd, zavm_evd
    IF (kt == nit000) THEN
      IF (lwp) WRITE(numout, FMT = *)
      IF (lwp) WRITE(numout, FMT = *) 'zdf_evd : Enhanced Vertical Diffusion (evd)'
      IF (lwp) WRITE(numout, FMT = *) '~~~~~~~ '
      IF (lwp) WRITE(numout, FMT = *)
    END IF
    !$ACC KERNELS
    zavt_evd(:, :, :) = p_avt(:, :, :)
    !$ACC END KERNELS
    SELECT CASE (nn_evdm)
    CASE (1)
      !$ACC KERNELS
      zavm_evd(:, :, :) = p_avm(:, :, :)
      DO jk = 1, jpkm1
        DO jj = 2, jpjm1
          DO ji = 2, jpim1
            IF (MIN(rn2(ji, jj, jk), rn2b(ji, jj, jk)) <= - 1.E-12) THEN
              p_avt(ji, jj, jk) = rn_evd * wmask(ji, jj, jk)
              p_avm(ji, jj, jk) = rn_evd * wmask(ji, jj, jk)
            END IF
          END DO
        END DO
      END DO
      zavm_evd(:, :, :) = p_avm(:, :, :) - zavm_evd(:, :, :)
      !$ACC END KERNELS
      CALL iom_put("avm_evd", zavm_evd)
    CASE DEFAULT
      !$ACC KERNELS
      DO jk = 1, jpkm1
        DO jj = 2, jpjm1
          DO ji = 2, jpim1
            IF (MIN(rn2(ji, jj, jk), rn2b(ji, jj, jk)) <= - 1.E-12) p_avt(ji, jj, jk) = rn_evd * wmask(ji, jj, jk)
          END DO
        END DO
      END DO
      !$ACC END KERNELS
    END SELECT
    !$ACC KERNELS
    zavt_evd(:, :, :) = p_avt(:, :, :) - zavt_evd(:, :, :)
    !$ACC END KERNELS
    CALL iom_put("avt_evd", zavt_evd)
    IF (l_trdtra) CALL trd_tra(kt, 'TRA', jp_tem, jptra_evd, zavt_evd)
  END SUBROUTINE zdf_evd
END MODULE zdfevd


text:
*-----------------------------------------------------------------------
*               U N P A C K I
*-----------------------------------------------------------------------
*
*     Given the type of integral, the index of the integral,
*     determine -- kv, iel1, iel2, iel3, iel4
*
      SUBROUTINE unpacki(int,i,kv,iel1,iel2,iel3,iel4)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        PARAMETER (NWD=70,NOD=220,NOFFD=800)
*
      COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,MASS,NCFG,IB,IC,ID
     :   ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,FINE,NSCF,NCLOSD,RMASS
*
      POINTER (pkval,kval(1)),(pvalue,value(1)),
     :        (pih,ih(1)),(pjh,jh(1)),
     :        (pcoeff,coeff(1)),(pnijptr,nijptr(1)),
     :        (pinptr,inptr(1))
      COMMON/STATE/INTPTR(6),IDIM,PKVAL,PVALUE,NZE,PIH,PJH,
     :             NCODIM,PCOEFF,PNIJPTR,PINPTR
*
*      ... unpack the electron data
*
          kv = kval(i)
          if (int .le. 3 .or. int .eq. 6) then
              iel2 = mod(kv,64) + nclosd
              kv = kv/64
              iel1 = mod(kv,64) + nclosd
              kv = kv/64
          else if (int .eq. 4 .or. int .eq. 5) then
              iel4 = mod(kv,64) + nclosd
              kv = kv/64
              iel3 = mod(kv,64) + nclosd
              kv = kv/64
              iel2 = mod(kv,64) + nclosd
              kv = kv/64
              iel1 = mod(kv,64) + nclosd
              kv = kv/64
          end if
          end


text:
*DECK DFNORM
      DOUBLE PRECISION FUNCTION DFNORM (N, A, W)
C-----------------------------------------------------------------------
C This function computes the norm of a full N by N matrix,
C stored in the array A, that is consistent with the weighted max-norm
C on vectors, with weights stored in the array W:
C   DFNORM = MAX(i=1,...,N) ( W(i) * Sum(j=1,...,N) ABS(a(i,j))/W(j) )
C-----------------------------------------------------------------------
      INTEGER N,   I, J
      DOUBLE PRECISION A,   W, AN, SUM
      DIMENSION A(N,N), W(N)
      AN = 0.0D0
      DO 20 I = 1,N
        SUM = 0.0D0
        DO 10 J = 1,N
 10       SUM = SUM + ABS(A(I,J))/W(J)
        AN = MAX(AN,SUM*W(I))
 20     CONTINUE
      DFNORM = AN
      RETURN
C----------------------- End of Function DFNORM ------------------------
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg