! ##################################################################################################################################
! 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 BUILD_G_NM

! For one subcase:

!   1) Calcs UM displacements: UM = GMN*UN where:

!          UN  = Displs calc'd in subr BUILD_N_FS

!   2) Merge UM and UN to get UG

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  NDOFG, NDOFM, NDOFN, NTERM_GMN, BLNK_SUB_NAM
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BUILD_G_NM_BEGEND
      USE CONSTANTS_1, ONLY           :  ONE
      USE PARAMS, ONLY                :  PRTDISP
      USE SPARSE_MATRICES, ONLY       :  I_GMN, J_GMN, GMN, SYM_GMN
      USE COL_VECS, ONLY              :  UG_COL, UM_COL, UN_COL

      USE BUILD_G_NM_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME   = 'BUILD_G_NM'

      INTEGER(LONG)                   :: I                 ! DO loop index
      INTEGER(LONG)                   :: G_SET_COL         ! Col no. in TDOF for G  displ set definition
      INTEGER(LONG)                   :: N_SET_COL         ! Col no. in TDOF for N  displ set definition
      INTEGER(LONG)                   :: M_SET_COL         ! Col no. in TDOF for M  displ set definition
      INTEGER(LONG), PARAMETER        :: NUMCOLS     = 1   ! Variable for number of cols of an array
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BUILD_G_NM_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! Recover UM from GMN x UN

      IF (NDOFM > 0) THEN

         CALL MATMULT_SFF ( 'GMN', NDOFM, NDOFN, NTERM_GMN, SYM_GMN, I_GMN, J_GMN, GMN, 'UN', NDOFN, NUMCOLS, UN_COL, 'Y',         &
                            'UM', ONE, UM_COL )
! Merge UN and UM to get UG

         CALL TDOF_COL_NUM ( 'G ', G_SET_COL )
         CALL TDOF_COL_NUM ( 'N ', N_SET_COL )
         CALL TDOF_COL_NUM ( 'M ', M_SET_COL )

         CALL MERGE_COL_VECS ( N_SET_COL, NDOFN, UN_COL, M_SET_COL, NDOFM, UM_COL,G_SET_COL, NDOFG, UG_COL )

      ELSE

         DO I=1,NDOFG
            UG_COL(I) = UN_COL(I)
         ENDDO

      ENDIF

! Print out displ matrices if PRTDISP says to

      IF (PRTDISP(1) == 1) THEN
         IF (NDOFG  > 0) THEN
            CALL WRITE_VECTOR ( '   G-SET DISPL VECTOR   ', 'DISPL', NDOFG, UG_COL)
         ENDIF
      ENDIF

      IF ((PRTDISP(2) == 1) .OR. (PRTDISP(2) == 3)) THEN
         IF (NDOFN  > 0) THEN
            CALL WRITE_VECTOR ( '   N-SET DISPL VECTOR   ', 'DISPL', NDOFN, UN_COL)
         ENDIF
      ENDIF

      IF ((PRTDISP(2) == 2) .OR. (PRTDISP(2) == 3)) THEN
         IF (NDOFM  > 0) THEN
            CALL WRITE_VECTOR ( '   M-SET DISPL VECTOR   ', 'DISPL', NDOFM, UM_COL)
         ENDIF
      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

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

      END SUBROUTINE BUILD_G_NM


text:
      SUBROUTINE CLOUD0
*
*
*       Cloud initialization.
*       ---------------------
*
      INCLUDE 'common6.h'
      COMMON/CLOUDS/  XCL(3,MCL),XDOTCL(3,MCL),BODYCL(MCL),RCL2(MCL),
     &                CLM(MCL),CLMDOT(MCL),CLDOT,VCL,SIGMA,RB2,PCL2,
     &                TCL,STEPCL,NCL,NEWCL
*
*
*       Initialize cloud variables.
      NCL = 0
      TCL = 0.0D0
      NEWCL = 0
      PCL2 = 0.0
*
*       Read the cloud parameters (pc, km/sec, Msun & pc).
      if (amusein.eq.0) READ (5,*)  NCL, RB2, VCL, SIGMA,
     &            (CLM(J),J=1,NCL), (RCL2(J),J=1,NCL)
      if(rank.eq.0)then
      WRITE (6,1)  NCL, RB2, VCL, SIGMA
    1 FORMAT (/,12X,'CLOUDS:    NCL =',I4,'  RB =',F5.1,
     &                      '   MEAN VELOC =',F5.1,'  DISP =',F5.1)
      WRITE (6,2)  (CLM(J),J=1,NCL)
    2 FORMAT (/,12X,'CLOUD MASSES:   ',1P,10E9.1)
      WRITE (6,3)  (RCL2(J),J=1,NCL)
    3 FORMAT (/,12X,'PLUMMER RADII:  ',1P,10E9.1)
      end if
      RBAR1 = RBAR
      IF (RBAR.EQ.0.0) RBAR1 = 1.0
*       Set cloud parameters in scaled units.
      RB2 = RB2/RBAR1
      A1 = 0.047*SQRT(ZMASS*ZMBAR/RBAR1)
*       Rms velocity of cluster members in km/sec.
      A2 = A1/SQRT(0.5D0*ZMASS)
*       Velocity unit.
      VCL = VCL/A2
*       Cloud velocity in scaled units.
      SIGMA = SIGMA/A2
*       Specify conservative cloud integration step using crossing time.
      STEPCL = 0.002*TCR*RB2/VCL
*
*       Adopt a quantized value.
      CALL STEPK(STEPCL,DTN)
      STEPCL = DTN
*
*       Scale radii & masses to model units.
      DO 10 J = 1,NCL
          RCL2(J) = RCL2(J)/RBAR1
          CLM(J) = CLM(J)/ZMBAR
   10 CONTINUE
*
      if(rank.eq.0)then
      WRITE (6,15)  RB2, VCL, SIGMA, STEPCL
   15 FORMAT (/,12X,'SCALING:     RB =',F6.1,'  VCL =',F5.1,
     &                        '  SIGMA =',F5.1,'  STEP =',1P,E10.2)
      WRITE (6,20)  (CLM(J),J=1,NCL)
   20 FORMAT (/,12X,'SCALED MASSES:  ',10F7.2)
      WRITE (6,25)  (RCL2(J),J=1,NCL)
   25 FORMAT (/,12X,'SCALED RADII:   ',10F7.2)
      end if
      CLDOT = 0.1*RB2/VCL
*       Time scale for 'sun-rise' is 0.05 of the cloud crossing time.
      CLDOT = 1.0/CLDOT
*
*       Define the square of cloud half-mass radii & growth times.
      DO 30 J = 1,NCL
          RCL2(J) = RCL2(J)**2
          CLMDOT(J) = CLM(J)*CLDOT
   30 CONTINUE
*
*       Set square boundary radius & impact parameter.
      RB2 = RB2**2
      PCL2 = RB2
*       Define density centre for routine CLOUD.
      DO 40 K = 1,3
          RDENS(K) = 0.0
   40 CONTINUE
*
*       Initialize new clouds on the boundary.
      DO 50 ICL = 1,NCL
          CALL CLOUD(ICL)
   50 CONTINUE
*
      RETURN
*
      END

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


! @Copyright 2007 Kristjan Haule
SUBROUTINE OUTWIN(A, B, REL, V, Rnet, Nr, dh, ER, FL, Z, VAL,SLO,Nodes)
  !         Integration of skalarrel. Schroedinger Equation
  ! ----------------------------------------------------------------
  !  Input:
  !    ER    Energy in Rydbergs
  !    FL    Angular momentum
  !    Z     Nuclear charge
  !    V     rad.sym. Potential in Hartree
  !    RNOT  first radial grid point
  !    DH    log. increment
  !    Nr    Number of radial mesh points
  !  Output:
  !    VAL,SLO:  Wave function and gradient at the spherical boundary
  !    Nodes:    number of nodes
  !
  !    Rydberg units
  !
  ! ----------------------------------------------------------------
  IMPLICIT NONE
  REAL*8, intent(out) :: A(Nr), B(Nr)    ! solution with large and small component
  LOGICAL, intent(in) :: rel             ! real of complex
  REAL*8, intent(in)  :: V(Nr)           ! radial symmetric potential
  REAL*8, intent(in)  :: Rnet(Nr), dh, ER, Z ! radial grid-points, log. increment, Energy in Rydbergs, nuclear charge
  INTEGER, intent(in) :: Nr, fl          ! Number of radial mesh points, Angular momentum l
  INTEGER, intent(out):: nodes           ! number of nodes
  REAL*8, intent(out) :: val, slo        ! value at the spherical boundary:
  !    for large component u(r)=g(r,1) with psi=(u/r)*ylm
  !           val is u(rmax), i.e. (rmax * radial w.f.)
  !           slo is radial derivative of g at rmax (rmax * radial w.f.)'
  !
  ! Note: if r = b(exp(a*z)-1) then Jacobian dr=a(r+b) dz
  !
  ! locals
  INTEGER :: K, ij
  REAL*8  :: D(2,3), C, FLLP1, R83SQ, R1, R2, R3, H83, G0, AA, X, Y, U, ZZ, DRDI
  REAL*8  :: DG1, DG2, DG3, DF1, DF2, DF3, F0, S, SF, B1, B2, DET, PHI, R
  !
  !     Hartree in Ryd
  !E=EH*2.d0
  !
  A(:)=0.0
  B(:)=0.0
  Nodes = 0
  ZZ = Z + Z
  C = 2.d0*137.0359895d0
  !
  if(.not.rel) C=1e10  ! set to infinity
  ! angular momentum part
  FLLP1 = FL*(FL + 1.d0)
  ! some coefficients for integration
  R83SQ = 64.D0/9.D0
  R1 = 1.D0/9.D0
  R2 = -5.D0*R1
  R3 = 19.D0*R1
  H83 = 8.D0/3.D0
  !
  G0 = 1.d0
  IF (Z .LT. 0.9D0) THEN
     S = FL+1.d0
     SF = FL
     F0 = FL/C
  ELSE
     AA = ZZ/C
     S = DSQRT(FLLP1 + 1.D0 - AA*AA)
     SF = S
     F0 = G0*(S - 1.D0)/AA
  ENDIF
  DO  K = 1,3 !
     R = RNET(K)
     DRDI = DH*R
     A(K) = (R**S)*G0
     B(K) = (R**SF)*F0
     D(1,K) = DRDI*A(K)*S/R
     D(2,K) = DRDI*B(K)*SF/R
  ENDDO
  !
  DG1 = D(1,1)
  DG2 = D(1,2)
  DG3 = D(1,3)
  DF1 = D(2,1)
  DF2 = D(2,2)
  DF3 = D(2,3)
  DO K = 4, Nr
     R = RNET(K)
     DRDI = DH*R
     !       Factor two because V is in Hartree-Rydberg !
     PHI = (ER - 2.d0*V(K)/R)*DRDI/C
     U = DRDI*C + PHI
     X = -DRDI/R
     Y = -FLLP1*X*X/U + PHI
     DET = R83SQ - X*X + U*Y
     B1 = A(K-1)*H83 + R1*DG1 + R2*DG2 + R3*DG3
     B2 = B(K-1)*H83 + R1*DF1 + R2*DF2 + R3*DF3
     A(K) = (B1*(H83-X) + B2*U)/DET
     B(K) = (B2*(H83+X) - B1*Y)/DET
     IF (A(K)*A(K-1) .LT. 0D0) Nodes = Nodes + 1
     DG1 = DG2
     DG2 = DG3
     DG3 = U*B(K) - X*A(K)
     DF1 = DF2
     DF2 = DF3
     DF3 = X*B(K) - Y*A(K)
  ENDDO
  !
  DO ij=1,Nr
     B(ij)=B(ij)*c/2.d0
  ENDDO
  !
  VAL = A(Nr)/RNET(Nr)
  SLO = DG3/(DH*RNET(Nr))
  SLO = (SLO-VAL)/RNET(Nr)
  RETURN
END SUBROUTINE OUTWIN

REAL*8 FUNCTION RINT13(rel,A1,B1,A2,B2,Nr,DX,R0)
  !PERFORMS RADIAL INTEGRALS
  !coded by D.D.KOELLING
  IMPLICIT NONE
  logical, intent(in):: rel
  REAL*8, intent(in) :: A1(Nr),B1(Nr),A2(Nr),B2(Nr)
  REAL*8, intent(in) :: DX, R0
  INTEGER, intent(in):: Nr
  !
  ! temporary
  REAL*8 :: D, S, CIN, R, R1, Z2, Z4, P1, P2
  INTEGER:: J, J1
  !
  CIN=1.d0/137.0359895d0**2
  IF(.NOT.REL) CIN=1E-22
  !
  D=EXP(DX)
  J=3-MOD(Nr,2)
  R=r0*(D**(J-1))
  J1=J-1
  R1=R/D
  Z4=0
  Z2=0
  DO
     Z4=Z4+R*(A1(J)*A2(J)+CIN*B1(J)*B2(J))
     R=R*D
     J=J+1
     IF(J.GE.Nr) EXIT
     Z2=Z2+R*(A1(J)*A2(J)+CIN*B1(J)*B2(J))
     R=R*D
     J=J+1
  ENDDO
  P1=r0*(A1(1)*A2(1)+CIN*B1(1)*B2(1))
  P2=R1*(A1(J1)*A2(J1)+CIN*B1(J1)*B2(J1))
  S=2*Z2+4*Z4+R*(A1(J)*A2(J)+CIN*B1(J)*B2(J))+P2
  S=(DX*S+P1)/3.0D0
  IF(J1.GT.1) S=S+0.5D0*DX*(P1+P2)
  RINT13=S
  RETURN
END FUNCTION RINT13

REAL*8 FUNCTION RINT13g(C1,C2,A1,B1,A2,B2,Nr,DX_,R0_)
  ! C1=1, C2=[1.d0/137.0359895d0**2 or 0]
  IMPLICIT NONE
  !        but here C1 and C2 are arbitrary
  REAL*8, intent(in) :: C1, C2
  REAL*8, intent(in) :: A1(Nr), B1(Nr), A2(Nr), B2(Nr)
  REAL*8, intent(in) :: DX_, R0_
  INTEGER, intent(in):: Nr
  ! locals
  INTEGER::            J, J1
  REAL*8 ::   D, S, P1, P2, R, R1, Z2, Z4
  D = EXP(DX_)
  J = 3-MOD(Nr,2)
  J1=J-1
  R = r0_*(D**(J-1))
  R1 = R/D
  Z4 = 0
  Z2 = 0
  DO
     Z4=Z4+R*(C1*A1(J)*A2(J) + C2*B1(J)*B2(J))     ! Z4 ~ R*A*X (r)
     R=R*D
     J=J+1
     IF (J .GE. Nr) EXIT
     Z2 = Z2 + R*(C1*A1(J)*A2(J) + C2*B1(J)*B2(J))  ! Z2 ~ R*A*X (r+dr)
     R=R*D
     J=J+1
  enddo
  P1 = R0_*(C1*A1(1)*A2(1) + C2*B1(1)*B2(1))
  P2 = R1*(C1*A1(J1)*A2(J1) + C2*B1(J1)*B2(J1))
  S = 2*Z2 + 4*Z4 + R*(C1*A1(J)*A2(J) + C2*B1(J)*B2(J)) + P2
  S = (DX_*S+P1)/3.0D+0
  IF (J1.GT.1) S=S+0.5D0*DX_*(P1+P2)
  RINT13g=S
  RETURN
END FUNCTION RINT13g

! Nr -> jrj(jatom)
SUBROUTINE ReadPotential(VR, filename, maxNr, Nr, nat)
!!! Reading potential
  REAL*8, intent(out)       :: VR(maxNr,nat)
  CHARACTER*200, intent(in) :: filename
  INTEGER, intent(in)       :: nat, maxNr, Nr(nat)
  !! temporary
  INTEGER :: jtape, iscf, idummy, j, jatom
  !
  jtape=18
  open(jtape,file=filename,status='old')
  READ(jtape, '(50X,I2,//)') ISCF

  do jatom=1,nat
     READ(jtape,1980)
     READ(jtape,2000) idummy
     READ(jtape,2031)
     READ(jtape,2022)(VR(J,jatom),J=1,Nr(jatom))
     READ(jtape,2031)
     READ(jtape,2030)
  enddo

  VR(:,:)=VR(:,:)/2.0D0  ! To Hartree's?
  close(jtape)
2022 FORMAT(3X,4E19.12)
2000 FORMAT(16X,I2//)
1980 FORMAT(3X)
2031 FORMAT(/)
2030 FORMAT(///)
END SUBROUTINE ReadPotential

SUBROUTINE ReadLinearizationE(E, filename, nat, lmax2)
!!! Reading linearization energies
  CHARACTER*200, intent(in) :: filename
  INTEGER, intent(in) :: LMAX2
  REAL*8, intent(out) :: E(0:LMAX2,nat)
  ! locals
  INTEGER :: itape, i
  REAL*8  :: elo
  !jatom = jatom+1  ! to convert from python to fortran
  itape=9
  open(itape,file=filename,status='old',form='unformatted')
  DO jatom=1,nat
     READ(itape) E(:,jatom)     ! linearization energies
     READ(itape) elo            ! local
  ENDDO
  close(itape)
END SUBROUTINE ReadLinearizationE

  ! NORM OF V(0,0)=V00(R)*R/SQRT(4.D0*PI)
  !DO i=1,jri
  !   Rx(i)=r0*(exp(dx*(i-1.d0)))
  !ENDDO
  ! dx <- dx(jatom)
  ! El <- E(l)

  ! DX <- DX(jatom)
  ! jrj <- jrj(jatom)
  ! r0 <- r0(jatom)
  ! Nrad, JRJ

SUBROUTINE ATPAR(A,B,Ae,Be,Aee,Bee,Pei,rel,l,VR,Rx,El,r0,dx,Zn,Nr)
  IMPLICIT NONE
  REAL*8, intent(out)  :: A(Nr), B(Nr)   ! large and small component of Sch. equation: u(r)
  REAL*8, intent(out)  :: Ae(Nr), Be(Nr) ! large and small component of derivative of Sch. equation: udot(r)
  REAL*8, intent(out)  :: Aee(Nr), Bee(Nr) ! large and small component of the second derivative of Sch. equation: udot(r)
  REAL*8, intent(out)  :: Pei            ! <udot|udot>
  !INTEGER, intent(out) :: Nodel          ! number of nodes of u(r)
  !REAL*8, intent(out)  :: dUv, dUve, dUvee  ! derivatives at the Rmt: dUv,dUve,dUvee of ul, dot{ul}, dot_dot{ul}
  LOGICAL, intent(in) :: rel             ! relativistic or not
  INTEGER, intent(in) :: l               ! orbital momentum
  REAL*8, intent(in)  :: VR(Nr),Rx(Nr)   ! radial potenial, radial mesh
  REAL*8, intent(in)  :: El              ! linearization energy for this l
  REAL*8, intent(in)  :: r0, dx, Zn      ! first point of the radial grid, dr for the grid, nuclear charge
  INTEGER, intent(in) :: Nr              ! number of radial points
  ! External
  INTERFACE
     FUNCTION RINT13(rel,A1,B1,A2,B2,Nr,DX,R0)
       REAL*8 :: RINT13
       logical, intent(in):: rel
       REAL*8, intent(in) :: A1(Nr),B1(Nr),A2(Nr),B2(Nr)
       REAL*8, intent(in) :: DX, R0
       INTEGER, intent(in):: Nr
     END FUNCTION RINT13
  END INTERFACE
  ! temporary
  REAL*8  :: A1(Nr), B1(Nr), A2(Nr), B2(Nr)
  REAL*8  :: Uv1, dUv1, Uv2, dUv2, Uvee, Uv, dElR, ovlp, cross, Uve, dUv, dUve, dUvee
  INTEGER :: Nodel, i
  !
  !  CALCULATE ENERGY-DERIVATIVE BY FINITE DIFFERENCE
  !  dElR IS THE UPWARD AND DOWNWARD ENERGY SHIFT IN Rydbers
  dElR=4.0D-3

  CALL OUTWIN(A1, B1, rel, VR, Rx, Nr, dx, El-dElR, l, Zn, Uv1, dUv1, Nodel)
  ! Normalization
  ovlp = RINT13(rel, A1, B1, A1, B1, Nr, dx, r0)
  A1(:) = A1(:)/SQRT(ovlp)
  B1(:) = B1(:)/SQRT(ovlp)
  Uv1=Uv1/SQRT(ovlp)         ! Normalized value at Rmt and energy El-dElR
  dUv1=dUv1/SQRT(ovlp)

  CALL OUTWIN(A2, B2, rel, VR, Rx, Nr, dx, El+dElR, l, Zn, Uv2, dUv2, Nodel)
  ! Normalization
  ovlp = RINT13(rel, A2, B2, A2, B2, Nr, dx, r0)
  A2(:) = A2(:)/SQRT(ovlp)
  B2(:) = B2(:)/SQRT(ovlp)
  Uv2=Uv2/SQRT(ovlp)         ! Normalized value at Rmt and energy El+dElR
  dUv2=dUv2/SQRT(ovlp)

  ! Finite difference
  UvE = (Uv2-Uv1)/(2.*dElR)       ! energy derivative of the value at Rmt, i.e., udot(Rmt)
  dUvE= (dUv2-dUv1)/(2.*dElR)     ! d/dr udot(Rmt)
  Ae(:) = (A2(:)-A1(:))/(2.*dElR)
  Be(:) = (B2(:)-B1(:))/(2.*dElR)

  ! Schroedinger equation solver at El
  CALL OUTWIN(A, B, rel, VR, Rx, Nr, dx, El, l, Zn, Uv, dUv, Nodel)
  ovlp = RINT13(rel, A, B, A, B, Nr, dx, r0)
  A(:) = A(:)/SQRT(ovlp)
  B(:) = B(:)/SQRT(ovlp)
  Uv=Uv/SQRT(ovlp)         ! Normalized value at Rmt and energy El+dElR
  dUv=dUv/SQRT(ovlp)

  ! Orthogonalization between u(r) and udot(r)
  CROSS = RINT13(rel, A, B, Ae, Be, Nr, dx, r0)
  Ae(:) = Ae(:)-CROSS*A(:)
  Be(:) = Be(:)-CROSS*B(:)
  UvE = UvE-CROSS*Uv
  dUvE=dUvE-CROSS*dUv

  ! PEI = <udor(r)|udot(r)>
  PEI = RINT13(rel, Ae, Be, Ae, Be, Nr, dx, r0)

  ! Finite difference
  Aee(:) = (A2(:)-2*A(:)+A1(:))/(2.*dElR**2)
  Bee(:) = (B2(:)-2*B(:)+B1(:))/(2.*dElR**2)
  Uvee = (Uv2-2*Uv+Uv1)/(2.*dElR**2)
  dUvee= (dUv2-2*dUv+dUv1)/(2.*dElR**2)
  ! Orthogonalization between udotdot(r) and u(r)
  CROSS = RINT13(rel, Aee, Bee, A, B, Nr, dx, r0)
  Aee(:) = Aee(:)-CROSS*A(:)
  Bee(:) = Bee(:)-CROSS*B(:)
  Uvee = Uvee - CROSS*Uv
  dUvee=dUvee - CROSS*dUv
  ! Orthogonalization between udotdot(r) and udot(r)
  CROSS = RINT13(rel, Aee, Bee, Ae, Be, Nr, dx, r0)
  CROSS = CROSS/PEI
  Aee(:) = Aee(:)-CROSS*Ae(:)
  Bee(:) = Bee(:)-CROSS*Be(:)
  Uvee = Uvee - CROSS*UvE
  dUvee=dUvee - CROSS*dUvE

  !print *, RINT13(rel, A, B, Ae, Be, Nr, dx, r0), RINT13(rel, A, B, Aee, Bee, Nr, dx, r0), RINT13(rel, Ae, Be, Aee, Bee, Nr, dx, r0)

  RETURN
END SUBROUTINE ATPAR


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Copyright (c) 1997 Neil N. Carlson
!!
!! This file is part of MFE1 which is released under the MIT license.  See the
!! file LICENSE or visit http://opensource.org/licenses/MIT for details.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module problem_data

  use mfe_constants, only: wp
  implicit none
  private

  real(kind=wp), save, public :: visc

end module problem_data

module problem_init

  use problem_data
  use common_io
  implicit none
  private

  public :: read_problem_data

  contains

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!
   !!  READ_PROBLEM_DATA
   !!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    subroutine read_problem_data ()

      call read_tagged_data (visc, "PDE viscosity coefficient")

    end subroutine read_problem_data

end module problem_init

module problem_pde

  use problem_data
  use mfe_constants, only: wp
  use mfe_data, only: eqw
  use local_arrays
  use local_laplacian
  implicit none
  private

  public  :: pde_rhs
  private :: flux

  contains

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!
   !!  LOAD_PDE_RHS -- Inner products for gasdynamics.
   !!
   !!        Equation 1 == mass density,
   !!        Equation 2 == momentum density,
   !!        Equation 3 == total energy density (internal plus kinetic).
   !!
   !!     Each equation is in conservation law form
   !!
   !!        du/dt = -df/dx + visc* d2u/dx2.
   !!
   !!  FLUX -- Gas dynamics flux f.
   !!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     subroutine pde_rhs (t)

       real(kind=wp), intent(in) :: t

       ! Local variables
       integer :: i, k
       real(kind=wp) :: rx1, rx2, term
       real(kind=wp), dimension(3) :: umid, f1, f2, favg
       real(kind=wp), parameter :: c1 = 1.0_wp / 6.0_wp, c2 = 4.0_wp / 6.0_wp

       do i = 1, ncell

         f1 = flux (u(1,i) % u)                      ! Flux at left endpoint.
         f2 = flux (u(2,i) % u)                      ! Flux at right endpoint.
         umid = 0.5_wp * (u(1,i) % u + u(2,i) % u)   ! Variables at midpoint.
         favg = c1 * (f1 + f2) + c2 * flux (umid)    ! Average flux (Simpson).

         rx1 = 0.0_wp
         rx2 = 0.0_wp

         do k = 1, 3

           term = (eqw(k) * (favg(k) - f1(k)))
           rx1 = rx1 - term * n(k,i) % x
           r(1,i) % u(k) = - term * n(k,i) % u

           term = (eqw(k) * (f2(k) - favg(k)))
           rx2 = rx2 - term * n(k,i) % x
           r(2,i) % u(k) = - term * n(k,i) % u

         end do

         r(1,i) % x = rx1
         r(2,i) % x = rx2

       end do

       call laplacian (eqno=1, coef=visc)
       call laplacian (eqno=2, coef=visc)
       call laplacian (eqno=3, coef=visc)

     end subroutine pde_rhs

     function flux (u) result (f)

       real(kind=wp), dimension(:), intent(in) :: u
       real(kind=wp), dimension(3) :: f

       real(kind=wp), parameter :: gamma = 1.4_wp,                     &
                                   c1 = 0.5_wp * (3.0_wp - gamma),     &
                                   c2 = gamma - 1.0_wp,                &
                                   c3 = 0.5_wp * (1.0_wp - gamma),     &
                                   c4 = gamma

       ! Mass flux (equation 1).
       f(1) = u(2)

       ! Momentum flux (equation 2).
       f(2) = c1 * (u(2) * (u(2) / u(1))) + c2 * u(3)

       ! Energy flux (equation 3).
       f(3) = (c3 * (u(2) * (u(2) / u(1))) + c4 * u(3)) * (u(2) / u(1))

     end function flux

end module problem_pde


text:
!+ gaseous absorption after Rosenkranz 98 model
subroutine rosen98_gasabs &
 (errorstatus,&     ! out
  freq,  &      ! in
  tempK, &      ! in
  rhoWv, &      ! in
  pres,  &      ! in
  absAir,&      ! out
  absWv  )      ! out

! Description:
!  Based on frequency, temperature, water vapor density, and pressure, this routine
!  calculates the absorption due to air (N2 and O2) and water vapor in the frequency
!  range from 0.1 to 800.0 GHz, pressure range from 10 to 1.2e5 Pa, and absolute
!  temperatures larger than 100 K.
!
! Method:
!  ROSENKRANZ (1998) model  -- reference "Water vapor microwave
!  continuum absorption: a comparison of measurements and results"
!  To appear in Radio Science
!
! Owner: IGMK
!
! History:
!
! Version Date         Comment
! ------- ----         -------
! 0.01    02/08/2001   Fixed division by zero when rhowv = 0 - G. Petty
! 0.1     21/09/2009   Code adaption from G. Petty - M. Mech
! 0.2     27/10/2012   Corrections to WV bands and continuum - E. Orlandi
! 0.3     13/11/2012   Application of European Standards for Writing and
!                      Documenting Exchangeable Fortran 90 Code - M. Mech
! 0.4     27/02/2013   functions in module gasabs_module
!
! Code Description:
!  Language: Fortran 90.
!  Software Standards: "European Standards for Writing and
!   Documenting Exchangeable Fortran 90 Code".
!
! Parent Module: get_gasabs
!
! Declarations:

! Modules used:

  use kinds, only: dbl, &     ! integer parameter specifying double precision
                    long       ! integer parameter specifying long integer

  use gasabs_module            ! functions for calculating absorption by gases

!  use settings, only: verbose

!!  use report_module
! Imported Scalar Variables with intent (in):

  implicit none

!- End of header ---------------------------------------------------------------

! Subroutine arguments
! Scalar arguments with intent(in):

  real(kind=dbl), intent(in) :: freq     ! frequency [GHz]
  real(kind=dbl), intent(in) :: tempK    ! temperature [K]
  real(kind=dbl), intent(in) :: pres     ! pressure [Pa]
  real(kind=dbl), intent(in) :: rhoWv    ! water vapor density [kg/m**3]

! Scalar arguments with intent(out):

  real(kind=dbl), intent(out) :: absAir  ! extinction by dry air [Np/km]
  real(kind=dbl), intent(out) :: absWv   ! extinction by water vapor [Np/km]

! End of Subroutine arguments

! Local scalars:

  real(kind=dbl) :: pmb         ! pressure [mb]
  real(kind=dbl) :: vapden      ! water vpor density [g/m**3]
  real(kind=dbl) :: e           ! water vapor pressure [Pa]
  real(kind=dbl) :: q           ! specific humidity
  real(kind=dbl) :: Tv          ! virtuel temperature [K]
  real(kind=dbl) :: rhoair      ! moist air density [kg/m**3]

! Used Functions

!  real(kind=dbl) :: absn2       ! function to calculate extinction by n_2
!  real(kind=dbl) :: o2abs       ! function to calculate extinction by o2
!  real(kind=dbl) :: abh2o       ! function to calculate extinction by h_2o

! Error handling

  integer(kind=long), intent(out) :: errorstatus
  integer(kind=long) :: err = 0
  character(len=80) :: msg
  character(len=14) :: nameOfRoutine = 'rosen98_gasabs'

!   if (verbose >= 2) call report(info,'Start of ', nameOfRoutine)

!   ! check for "reasonable" input values

! !  print*, freq
!   if ((freq <= 0.0_dbl) .or. (freq > 800.0_dbl)) then
!      errorstatus = fatal
!      msg = 'Frequency not between 0 and 800 GHz in rosen98_gasabs!'
!      call report(errorstatus, msg, nameOfRoutine)
!      return
!   elseif (tempK <= 100.0_dbl) then
!      errorstatus = fatal
!      msg = 'Temperature lower than 100 K in rosen98_gasabs!'
!      call report(errorstatus, msg, nameOfRoutine)
!      return
!   elseif ((pres < 10.0_dbl) .or. (pres > 1.2d5)) then
!      print*, pres
!      errorstatus = fatal
!      msg = 'Pressure not between 10 and 1.2d5 Pa in rosen98_gasabs!'
!      call report(errorstatus, msg, nameOfRoutine)
!      return
!   else
!      err = success
!   end if
  print*,"checking input consistency"
  print*,"freq     ","temp     ","abs_hum      ","press"
  print*,freq, tempK, rhoWv,  pres
  ! convert pressure from Pa to Mb
  pmb = pres / 100.0_dbl

  ! convert vapor density from kg/m**3 to g/m**3

  vapden = rhoWv * 1000.0_dbl

  ! get volume extinction coefficients
  absair = absn2(tempK,pmb,freq) + o2abs(tempK,pmb,vapden,freq)
  abswv = abh2o(tempK,pmb,vapden,freq)

  ! convert vapor density to vapor pressure
  e = rhoWv * (tempK * 461.5_dbl)
  ! calculate specific humidity
  q = 0.622_dbl * e / pres
  ! calculate virtual temperature
  Tv = (1._dbl + 0.61_dbl * q) * tempK
  ! moist air density
  rhoair = pres / (Tv * 287.06_dbl)

!!  errorstatus = err
errorstatus = 1

!!  if (verbose >= 2) call report(info,'End of ', nameOfRoutine)

  return

end subroutine rosen98_gasabs

max_stars_repo_name: DaveOri/pyGasAbs

max_issues_repo_name: DaveOri/pyGasAbs

max_forks_repo_name: DaveOri/pyGasAbs


text:
**==dcint.spg  processed by SPAG 5.00N  at 11:12 on 23 Mar 1999
      SUBROUTINE DCINT(Num)
      IMPLICIT NONE
C
C ======================================================================
C
C  PURPOSE:
C       Decode one integer number.
C
C  METHOD:
C       - CALL MOVEIP to move line pointer IP to first non-blank character
C         in line buffer CHBUFF.
C       - Check sign.
C       - Move IP to next char. after sign.
C       - Decode integer NUM from CHBUFF.
C       - After decoding of NUM, IP will point at the first blank char,
C         after NUM.
C       - If any error, IP will point at char. that caused the error,
C         and error-flag will be set to a negative value.
C
C  ERROR HANDLING:
C       - If call to MOVEIP return error flag on, then error message
C         is printed.
C       - If invalid char. then error message is printed. Also
C         CHBUFF is printed with a ? at position where invalid
C         char. was found.
C       - If invalid char. IP is moved to the next blank character in
C         CHBUFF.
C
C  INPUT:
C       None.
C
C  OUTPUT:
C       NUM  - Decoded integer number, if no error.
C
C  INTERNAL VARIABLES:
C
C  CALL TO:
C       MOVEIP
C
C  CALL FROM:
C       General purpose routine.
C
C  LIBRARY:
C       FFDLIB
C
C  MACHINE DEPENDENT PART:
C       INCLUDE statement
C
C  PROGRAMMED BY:
C       Per Ivar Loekstad / ETOE / Aker Engineering a.s.
C
C  DATE/VERSION:
C       19.02.92/1.0
C
C ======================================================================
C
      INTEGER INDEX
      INTEGER Num , isign , iend , m
      CHARACTER ch*1
C
      INCLUDE 'ffd.inc'
C
C ----------------------------------------------------------------------
C
C Initialize.
C
      Num = 0
      isign = 0
C
C Move IP to next non-blank in line buffer CHBUFF.
C
      CALL MOVEIP(iend)
      IF ( IDCerr.LT.0 ) THEN
         WRITE (LPU,99002)
         IDCerr = -1
         GOTO 300
      ELSE
         IF ( iend.NE.0 ) GOTO 300
C
C Check sign.
C
         ch = CHBuff(IP:IP)
         IF ( ch.NE.'+' ) THEN
            IF ( ch.NE.'-' ) GOTO 200
            isign = 1
         ENDIF
      ENDIF
C
C Move IP to next char. after + or -
C
 100  IP = IP + 1
      IF ( IP.GT.NCOl ) THEN
C
C Shift sign.
C
         IF ( isign.NE.0 ) Num = -Num
         GOTO 300
      ELSE
         ch = CHBuff(IP:IP)
         IF ( ch.EQ.' ' ) THEN
            IF ( isign.NE.0 ) Num = -Num
            GOTO 300
         ENDIF
      ENDIF
C
C Decode integer number.
C
 200  m = INDEX('0123456789',ch) - 1
      IF ( m.LT.0 ) THEN
C
C ----------------------------------------------------------------------
C
C Error section.
C
         m = IP - 1
         CALL PRIERR()
C
C If invalid character, move IP to next blank character in CHBUFF.
C
 250     IF ( IP.GT.NCOl ) THEN
            IDCerr = -1
         ELSEIF ( CHBuff(IP:IP).EQ.' ' ) THEN
            IDCerr = -1
         ELSE
            IP = IP + 1
            GOTO 250
         ENDIF
      ELSE
         Num = Num*10 + m
         GOTO 100
      ENDIF
C
C ----------------------------------------------------------------------
 300  RETURN
C ----------------------------------------------------------------------
C
99001 FORMAT (/,' *** FFDLIB ERROR - Invalid character ',
     &        '- (MODULE: DCINT)',/)
99002 FORMAT (/,' *** FFDLIB ERROR - Call to MOVEIP - (MODULE: DCINT)')
C
      END


text:
Cell coverage with most carriers on campus is very poor, although it is reasonably good throughout most of the town of Davis.  If many people on campus complain about the cell coverage, it is more likely that the carriers will improve coverage.  Here is what to do to pass this information on to the carriers:

Deadcellzones:
Go to: http://www.deadcellzones.com mashup of cell phone coverage complaints
AT&T:
Go to: https://www.wireless.att.com/support/chatLanding.do

Click on technical support chat, then click on go to chat.  Give the tech support person a physical address for a location on campus and explain to them that you get a very poor signal on campus.

AT&T says that the university would actually have to arrange for a cell tower to be put up on campus.  Anyone know who at the university handles such things?
 I was under the impression that AT&T has wanted to put more towers on campus for a long time, but the powers that be on campus have made the process nearly impossible. Verizon has tried for years to get permission to build something and they finally got it. However, the work isnt scheduled to be finished until late 2009. wl

Just get a nextel, UCD has a service contract with them so they have full reception everywhere, even in the MUGA! Or you could just not have a cellphone... Also, they recently replaced one of the light towers next to Toomey with a light / cell tower combo. Users/StevenDaubert
 Nextel doesnt carry the iPhone. No point. :P wl

 Like Will said  no iPhone, no point.  :)  Ive tried campus with a Verizon phone, and got very poor reception.  Of course, one might wonder how someone who doesnt exist gets any reception at all.  In my case, Verizons reception didnt exist!  Users/IDoNotExist

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine PSZHHamps(pa,pb,pc,pd,amp)
C--    Formula taken from Plehn, Spira and Zerwas
C--    Nucl. Phys. B479 (1996) 46
      implicit none
      include 'constants.f'
      include 'masses.f'
      include 'scale.f'
      include 'first.f'
      double precision pa(4),pb(4),pc(4),pd(4),pcsq,pdsq,
     & ss,tt,uu,S,T,U,T1,T2,U1,U2,rhoc,rhod,tauQ,mQsq,lambdaHHH
      double complex Ftriangle,Fbox,Gbox,
     & Dabc,Dbac,Dacb,Cab,Cbc,Cac,Ccd,Cbd,Cad,qlI4,qlI3,
     & Cdelta,Cbox,amp(2)

c--- initialize QCDLoop, if necessary
      if (first) then
        call qlinit
        first=.false.
      endif

      mQsq=mt**2
      ss=(pa(4)+pb(4))**2
     &  -(pa(1)+pb(1))**2-(pa(2)+pb(2))**2-(pa(3)+pb(3))**2
      tt=(pa(4)+pc(4))**2
     &  -(pa(1)+pc(1))**2-(pa(2)+pc(2))**2-(pa(3)+pc(3))**2
      uu=(pb(4)+pc(4))**2
     &  -(pb(1)+pc(1))**2-(pb(2)+pc(2))**2-(pb(3)+pc(3))**2
      S=ss/mQsq
      T=tt/mQsq
      U=uu/mQsq
      lambdaHHH=3d0*(hmass/zmass)**2
      rhoc=hmass**2/mQsq
      rhod=hmass**2/mQsq
      T1=T-rhoc
      U1=U-rhoc
      T2=T-rhod
      U2=U-rhod
      tauQ=4d0/S
C      ss=(pa+pb)**2
C      tt=(pa+pc)**2
C      uu=(pb+pc)**2
      pcsq=pc(4)**2-pc(1)**2-pc(2)**2-pc(3)**2
      pdsq=pd(4)**2-pd(1)**2-pd(2)**2-pd(3)**2
      Cab=qlI3(0d0,0d0,ss,mQsq,mQsq,mQsq,musq,0)
      Cbc=qlI3(0d0,pcsq,uu,mQsq,mQsq,mQsq,musq,0)
      Cac=qlI3(0d0,pcsq,tt,mQsq,mQsq,mQsq,musq,0)
      Ccd=qlI3(pcsq,pdsq,ss,mQsq,mQsq,mQsq,musq,0)
      Cbd=qlI3(0d0,pdsq,tt,mQsq,mQsq,mQsq,musq,0)
      Cad=qlI3(0d0,pdsq,uu,mQsq,mQsq,mQsq,musq,0)
      Dabc=qlI4(0d0,0d0,pcsq,pdsq,ss,uu,mQsq,mQsq,mQsq,mQsq,musq,0)
      Dbac=qlI4(0d0,0d0,pcsq,pdsq,ss,tt,mQsq,mQsq,mQsq,mQsq,musq,0)
      Dacb=qlI4(0d0,pcsq,0d0,pdsq,tt,uu,mQsq,mQsq,mQsq,mQsq,musq,0)

      Ftriangle=2d0/S*(2d0+(4d0-S)*mQsq*Cab)

      Fbox=(4d0*S+8d0*S*mQsq*Cab-2d0*S*(S+rhoc+rhod-8d0)
     & *mQsq**2*(Dabc+Dbac+Dacb)
     & +(rhoc+rhod-8d0)*mQsq*(T1*Cac+U1*Cbc+U2*Cad+T2*Cbd
     & -(T*U-rhoc*rhod)*mQsq*Dacb))/S**2
      Gbox=((T**2+rhoc*rhod-8d0*T)*mQsq
     & *(S*Cab+T1*Cac+T2*Cbd-S*T*mQsq*Dbac)
     & +(U**2+rhoc*rhod-8d0*U)*mQsq
     & *(S*Cab+U1*Cbc+U2*Cad-S*U*mQsq*Dabc)
     & -(T**2+U**2-2d0*rhoc*rhod)*(T+U-8d0)*mQsq*Ccd
     & -2d0*(T+U-8d0)*(T*U-rhoc*rhod)*mQsq**2
     & *(Dabc+Dbac+Dacb))/(S*(T*U-rhoc*rhod))

      Cdelta=lambdaHHH*zmass**2/dcmplx(ss-hmass**2,hmass*hwidth)
      Cbox=cone

      amp(1)=Cdelta*Ftriangle+Cbox*Fbox
      amp(2)=Cbox*Gbox
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


! chk_backslash.f90 --
!     Check if the compiler treats backslashes in literal strings as escape characters
!
!     Background:
!     In the past some compilers treated backslashes in the same way as C compilers would,
!     to introduce special characters. Somewhat annoying.
!
program chk_backslash
    implicit none

    write( *, '(a)' ) 'If a backslash is not treated as an escape character, ', &
                      'this program simply prints a backslash:'
    write( *, '(a)' ) '\'

end program chk_backslash