text:
!=======================================================================
!
!  H2 photodissociation rate taking into
!  account shielding and grain extinction
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  K0  = Unattenuated photodissociation rate (in cm^3/s)
!  G0  = Incident FUV field (in Draine units)
!  AV  = visual extinction (in magnitudes)
!  NH2 = H2 column density (in cm^-2)
!
!  Program variables:
!  RATE   = H2 photodissociation rate taking into
!           account self-shielding and grain extinction
!  LAMBDA = wavelength (in Å) of a typical transition
!  DOPW   = Doppler linewidth (in Hz) of a typical transition
!           (assuming turbulent broadening with b=3 km/s)
!  RADW   = radiative linewidth (in Hz) of a typical transition
!
!  Functions called:
!  SCATTER  = attenuation due to scattering by dust
!  H2SHIELD = H2 shielding function (FGK or LEE)
!
!-----------------------------------------------------------------------
FUNCTION H2_PHOTODISSOCIATION_RATE(K0,G0,AV,NH2) RESULT(RATE)

   USE HEALPIX_TYPES
   USE GLOBAL_MODULE, ONLY : V_TURB
   IMPLICIT NONE

   REAL(KIND=DP) :: RATE
   REAL(KIND=DP), INTENT(IN) :: K0,G0,AV,NH2

   REAL(KIND=DP) :: LAMBDA,DOPW,RADW,SCATTER
   REAL(KIND=DP) :: H2SHIELD_FGK,H2SHIELD_LEE

   LAMBDA=1000.0D0
   DOPW=V_TURB/(LAMBDA*1.0D-8)
   RADW=8.0D7

!  Calculate the H2 photodissociation rate
   RATE=K0*G0*SCATTER(AV,LAMBDA)*H2SHIELD_FGK(NH2,DOPW,RADW)
!   RATE=K0*G0*SCATTER(AV,LAMBDA)*H2SHIELD_LEE(NH2)
!   RATE=K0*G0*H2SHIELD_LEE(NH2)

   RETURN
END FUNCTION H2_PHOTODISSOCIATION_RATE
!=======================================================================

!=======================================================================
!
!  CO photodissociation rate taking into
!  account shielding and grain extinction
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  K0  = Unattenuated photodissociation rate (in cm^3/s)
!  G0  = Incident FUV field (in Draine units)
!  AV  = visual extinction (in magnitudes)
!  NCO = CO column density (in cm^-2)
!  NH2 = H2 column density (in cm^-2)
!
!  Program variables:
!  RATE   = CO photodissociation rate taking into
!           account self-shielding and grain extinction
!  LAMBDA = wavelength (in Å) of a typical transition
!
!  Functions called:
!  LBAR     = function to determine the wavelength
!  SCATTER  = attenuation due to scattering by dust
!  COSHIELD = CO shielding function
!
!-----------------------------------------------------------------------
FUNCTION CO_PHOTODISSOCIATION_RATE(K0,G0,AV,NCO,NH2) RESULT(RATE)

   USE HEALPIX_TYPES
   IMPLICIT NONE

   REAL(KIND=DP) :: RATE
   REAL(KIND=DP), INTENT(IN) :: K0,G0,AV,NCO,NH2

   REAL(KIND=DP) :: LAMBDA,LBAR,SCATTER,COSHIELD

!  Calculate the mean wavelength of the CO dissociating bands
   LAMBDA=LBAR(NCO,NH2)

!  Calculate the CO photodissociation rate
   RATE=K0*G0*SCATTER(AV,LAMBDA)*COSHIELD(NCO,NH2)

   RETURN
END FUNCTION CO_PHOTODISSOCIATION_RATE
!=======================================================================

!=======================================================================
!
!  CI photoionization rate taking into account grain extinction
!  and shielding by CI and H2 lines, adopting the treatment of
!  Kamp & Bertoldi (2000, A&A, 353, 276, Equation 8)
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  K0   = Unattenuated photoionization rate (in cm^3/s)
!  G0   = Incident FUV field (in Draine units)
!  AV   = visual extinction (in magnitudes)
!  KAV  = tau(λ)/tau(V) correction factor
!  NCI  = CI column density (in cm^-2)
!  NH2  = H2 column density (in cm^-2)
!  TGAS = gas temperature (in K)
!
!  Program variables:
!  RATE = CI photoionization rate taking into
!         account shielding and grain extinction
!  TAUC = optical depth in the CI absorption band
!
!-----------------------------------------------------------------------
FUNCTION CI_PHOTOIONIZATION_RATE(K0,G0,AV,KAV,NCI,NH2,TGAS) RESULT(RATE)

   USE HEALPIX_TYPES
   IMPLICIT NONE

   REAL(KIND=DP) :: RATE
   REAL(KIND=DP), INTENT(IN) :: K0,G0,AV,KAV,NCI,NH2,TGAS

   REAL(KIND=DP) :: TAUC

!  Calculate the optical depth in the CI absorption band, accounting
!  for grain extinction and shielding by CI and overlapping H2 lines
   TAUC=KAV*AV+1.1D-17*NCI+(0.9D0*TGAS**0.27D0*(NH2/1.59D21)**0.45D0)
!  Calculate the CI photoionization rate
   RATE=K0*G0*EXP(-TAUC)
   write(79,8000) AV,NH2,NCI
   8000 format(1pe15.5,','1pe15.5,','1pe15.5)
   RETURN
END FUNCTION CI_PHOTOIONIZATION_RATE
!=======================================================================

!=======================================================================
!
!  SI photoionization rate -- needs to be implemented!
!  For now, use the standard expression for photorates
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  K0   = Unattenuated photoionization rate (in cm^3/s)
!  G0   = Incident FUV field (in Draine units)
!  AV   = visual extinction (in magnitudes)
!  KAV  = tau(λ)/tau(V) correction factor
!  NSI  = SI column density (in cm^-2)
!
!  Program variables:
!  RATE = SI photoionization rate taking into
!         account shielding and grain extinction
!  TAUS = optical depth in the SI absorption band
!
!-----------------------------------------------------------------------
FUNCTION SI_PHOTOIONIZATION_RATE(K0,G0,AV,KAV,NSI) RESULT(RATE)

   USE HEALPIX_TYPES
   IMPLICIT NONE

   REAL(KIND=DP) :: RATE
   REAL(KIND=DP), INTENT(IN) :: K0,G0,AV,KAV,NSI

   REAL(KIND=DP) :: TAUS

!  Calculate the optical depth in the SI absorption band, accounting
!  for grain extinction and shielding by ???
   TAUS=KAV*AV

!  Calculate the SI photoionization rate
   RATE=K0*G0*EXP(-TAUS)

   RETURN
END FUNCTION SI_PHOTOIONIZATION_RATE
!=======================================================================

!=======================================================================
!
!  H2 line self-shielding, adopting the treatment of
!  Federman, Glassgold & Kwan (1979, ApJ, 227, 466)
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  NH2  = H2 column density (in cm^-2)
!  DOPW = Doppler linewidth (in Hz)
!  RADW = radiative linewidth (in Hz)
!
!  Program variables:
!  SHIELDING_FACTOR = H2 self-shielding factor containing
!              both Doppler and radiative contributions
!  FPARA     = fraction of H2 in para state: 1/(1+o/p ratio)
!  FOSC      = oscillator strength of a typical transition
!  TAUD      = parameter tauD (eq. A7) in Federman's paper
!              (optical depth at line centre)
!  R         = parameter r  (eq. A2) in Federman's paper
!  T         = parameter t1 (eq. A6) in Federman's paper
!  U         = parameter u1 (eq. A6) in Federman's paper
!  JD        = parameter JD (eq. A8) in Federman's paper
!              (Doppler contribution to self-shielding)
!  JR        = parameter JR (eq. A9) in Federman's paper
!              (radiative contribution to self-shielding)
!
!-----------------------------------------------------------------------
FUNCTION H2SHIELD_FGK(NH2,DOPW,RADW) RESULT(SHIELDING_FACTOR)

   USE HEALPIX_TYPES
   IMPLICIT NONE

   REAL(KIND=DP) :: SHIELDING_FACTOR
   REAL(KIND=DP), INTENT(IN) :: NH2,DOPW,RADW

   REAL(KIND=DP) :: FPARA,FOSC,TAUD
   REAL(KIND=DP) :: R,T,U,JD,JR

!  Calculate the optical depth at line centre = N(H2)*f_para*(πe^2/mc)*f/(√πß) ≈ N(H2)*f_para*(1.5E-2)*f/ß
   FPARA=0.5D0 ! (assume o/p ratio=1)
   FOSC=1.0D-2 ! Typical oscillator strength
   TAUD=NH2*FPARA*(1.497358985D-2)*FOSC/DOPW

!  Calculate the Doppler core contribution to the self-shielding (JD)
   IF(TAUD.EQ.0.0D0) THEN
      JD=1.0D0
   ELSE IF(TAUD.LT.2.0D0) THEN
      JD=EXP(-(0.666666667D0*TAUD))
   ELSE IF(TAUD.LT.10.0D0) THEN
      JD=0.638D0*TAUD**(-1.25D0)
   ELSE IF(TAUD.LT.100.0D0) THEN
      JD=0.505D0*TAUD**(-1.15D0)
   ELSE
      JD=0.344D0*TAUD**(-1.0667D0)
   ENDIF

!  Calculate the radiative wing contribution to self-shielding (JR)
   IF(RADW.EQ.0.0D0) THEN
      JR=0.0D0
   ELSE
      R=RADW/(1.772453851D0*DOPW)
      T=3.02D0*((R*1.0D3)**(-0.064D0))
      U=SQRT(TAUD*R)/T
      JR=R/(T*SQRT(0.785398163D0+U**2))
   ENDIF

!  Calculate the total self-shielding factor
   SHIELDING_FACTOR=JD+JR

   RETURN
END FUNCTION H2SHIELD_FGK
!=======================================================================

!=======================================================================
!
!  H2 line shielding, using the computed values listed in
!  Lee et al. (1996, A&A, 311, 690, Table 10)
!
!  Appropriate shielding factors are determined by performing a
!  1-dimensional spline interpolation over the values listed in
!  Table 10 of Lee et al. which include contributions from line
!  overlap, self-shielding and continuum absorption
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  NH2 = H2 column density (in cm^-2)
!
!  Program variables:
!  SHIELDING_FACTOR = total H2 shielding factor with
!              contributions from both H2 and H lines
!              and dust continuum absorption by using
!              1D spline interpolation over the grid
!  SH2_GRID  = H2 shielding factors from Lee et al. (1996)
!              as a function of H2 column density
!  SH2_DERIV = 2nd derivative of SH2_GRID values from SPLINE
!  COL_GRID  = H2 column densities (in cm^-2)
!  NUM_COL   = number of H2 column densities
!  START     = .TRUE. when H2SHIELD_LEE is first called
!
!  Functions called:
!  SPLINE = second derivative of the supplied 1D function (in spline.f90)
!  SPLINT = 1-dimensional cubic spline interpolated value (in spline.f90)
!
!-----------------------------------------------------------------------
FUNCTION H2SHIELD_LEE(NH2) RESULT(SHIELDING_FACTOR)

   USE HEALPIX_TYPES
   USE SHIELDING_MODULE, ONLY : START,NUM_COL,COL_GRID,SH2_GRID,SH2_DERIV
   IMPLICIT NONE

   REAL(KIND=DP) :: SHIELDING_FACTOR
   REAL(KIND=DP), INTENT(IN) :: NH2

   REAL(KIND=DP) :: NH2_VALUE

   IF(START) THEN
      CALL SPLINE(COL_GRID,SH2_GRID,NUM_COL,1.0D30,1.0D30,SH2_DERIV)
   ENDIF

   NH2_VALUE=NH2
   IF(NH2.LT.COL_GRID(1)) NH2_VALUE=COL_GRID(1)
   IF(NH2.GT.COL_GRID(NUM_COL)) NH2_VALUE=COL_GRID(NUM_COL)

   CALL SPLINT(COL_GRID,SH2_GRID,SH2_DERIV,NUM_COL,NH2_VALUE,SHIELDING_FACTOR)
   IF(SHIELDING_FACTOR.LT.0.0D0) SHIELDING_FACTOR=0.0D0

   RETURN
END FUNCTION H2SHIELD_LEE
!=======================================================================

!=======================================================================
!
!  CO line shielding, using the computed values listed in
!  van Dishoeck & Black (1988, ApJ, 334, 771, Table 5)
!
!  Appropriate shielding factors are determined by performing a
!  2-dimensional spline interpolation over the values listed in
!  Table 5 of van Dishoeck & Black, which include contributions
!  from self-shielding and H2 screening
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  NCO = CO column density (in cm^-2)
!  NH2 = H2 column density (in cm^-2)
!
!  Program variables:
!  SHIELDING_FACTOR = total 12CO shielding factor with
!              contributions from both H2 and CO lines
!              by 2D spline interpolation over the grid
!  SCO_GRID  = log10 values of the 12CO shielding factors
!              from van Dishoeck & Black (1988) as a function
!              of CO column density (1st index) and H2 column
!              density (2nd index)
!  SCO_DERIV = 2nd derivative of SCO_GRID values from SPLIE2
!  NCO_GRID  = log10 values of CO column densities (in cm^-2)
!  NH2_GRID  = log10 values of H2 column densities (in cm^-2)
!  NUM_NCO   = number of CO column densities
!  NUM_NH2   = number of H2 column densities
!  START     = .TRUE. when COSHIELD is first called
!
!  Functions called:
!  SPLIE2 = second derivative of the supplied 2D function (in spline.f90)
!  SPLIN2 = 2-dimensional cubic spline interpolated value (in spline.f90)
!
!-----------------------------------------------------------------------
FUNCTION COSHIELD(NCO,NH2) RESULT(SHIELDING_FACTOR)

   USE HEALPIX_TYPES
   USE SHIELDING_MODULE, ONLY : START,NUM_NCO,NUM_NH2,NCO_GRID,NH2_GRID,SCO_GRID,SCO_DERIV
   IMPLICIT NONE

   REAL(KIND=DP) :: SHIELDING_FACTOR
   REAL(KIND=DP), INTENT(IN) :: NCO,NH2

   REAL(KIND=DP) :: LOG_NCO,LOG_NH2

   IF(START) THEN
      CALL SPLIE2(NCO_GRID,NH2_GRID,SCO_GRID,NUM_NCO,NUM_NH2,SCO_DERIV)
      START=.FALSE.
   ENDIF

   LOG_NCO=DLOG10(NCO+1.0D0)
   LOG_NH2=DLOG10(NH2+1.0D0)

   IF(LOG_NCO.LT.NCO_GRID(1)) LOG_NCO=NCO_GRID(1)
   IF(LOG_NH2.LT.NH2_GRID(1)) LOG_NH2=NH2_GRID(1)
   IF(LOG_NCO.GT.NCO_GRID(NUM_NCO)) LOG_NCO=NCO_GRID(NUM_NCO)
   IF(LOG_NH2.GT.NH2_GRID(NUM_NH2)) LOG_NH2=NH2_GRID(NUM_NH2)

   CALL SPLIN2(NCO_GRID,NH2_GRID,SCO_GRID,SCO_DERIV,NUM_NCO,NUM_NH2,LOG_NCO,LOG_NH2,SHIELDING_FACTOR)
   SHIELDING_FACTOR=10.0D0**SHIELDING_FACTOR

   RETURN
END FUNCTION COSHIELD
!=======================================================================

!=======================================================================
!
!  Scattering by dust grains, adopting the treatment of
!  Wagenblast & Hartquist (1989, MNRAS, 237, 1019) and
!  Flannery, Roberge & Rybicki (1980, ApJ, 236, 598)
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  AV     = visual extinction (in magnitudes)
!  LAMBDA = wavelength (in Å) of incident radiation
!
!  Program variables:
!  SCATTER = attenuation factor describing the influence of
!            grain scattering on the FUV flux, dependening
!            on the total column density and wavelength of
!            light (assuming albedo=0.3 gscat=0.8)
!  TAUV    = optical depth at visual wavelength (λ=5500Å)
!  TAUL    = optical depth at wavelength LAMBDA
!  A(0)    = a(0)*exp(-k(0)*tau)
!          = relative intensity decrease for 0 < tau < 1
!  A(I)    = ∑ a(i)*exp(-k(i)*tau) for i=1,5
!            relative intensity decrease for tau ≥ 1
!  K(0)    = see A0
!  K(I)    = see A(I)
!
!  Functions called:
!  XLAMBDA = function to determine tau(λ)/tau(V)
!
!-----------------------------------------------------------------------
FUNCTION SCATTER(AV,LAMBDA)

   USE HEALPIX_TYPES
   IMPLICIT NONE

   REAL(KIND=DP) :: SCATTER
   REAL(KIND=DP), INTENT(IN) :: AV,LAMBDA

   INTEGER(KIND=I4B) :: I
   REAL(KIND=DP)     :: TAUL,TAUV,XLAMBDA,EXPONENT
   REAL(KIND=DP), DIMENSION(0:5) :: A=(/1.000D0,2.006D0,-1.438D0,0.7364D0,-0.5076D0,-0.0592D0/)
   REAL(KIND=DP), DIMENSION(0:5) :: K=(/0.7514D0,0.8490D0,1.013D0,1.282D0,2.005D0,5.832D0/)

!  Calculate the optical depth at visual wavelength
   TAUV=AV/1.086D0

!  Convert the optical depth to that at the desired wavelength
   TAUL=TAUV*XLAMBDA(LAMBDA)

!  Calculate the attenuation due to scattering by dust (SCATTER)
   SCATTER=0.0D0
   IF(TAUL.LT.1.0D0) THEN
      EXPONENT=K(0)*TAUL
      IF(EXPONENT.LT.100.0D0) THEN
         SCATTER=A(0)*EXP(-EXPONENT)
      ENDIF
   ELSE
      DO I=1,5
         EXPONENT=K(I)*TAUL
         IF(EXPONENT.LT.100.0D0) THEN
            SCATTER=SCATTER+A(I)*EXP(-EXPONENT)
         ENDIF
      ENDDO
   ENDIF

   RETURN
END FUNCTION SCATTER
!=======================================================================

!=======================================================================
!
!  Determine the ratio of the optical depth at a given wavelength to
!  that at visual wavelength (λ=5500Å) using the extinction curve of
!  Savage & Mathis (1979, ARA&A, 17, 73, Table 2)
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  LAMBDA  = wavelength (in Å)
!
!  Program variables:
!  XLAMBDA       = Ratio of tau(λ)/tau(V) at the desired wavelength
!                  by 1D spline interpolation over the grid values
!  XLAMBDA_GRID  = tau(λ)/tau(V) values, determined by dividing the
!                  Aλ/E(B-V) values from Savage & Mathis (1979) by
!                  an assumed reddening of R=AV/E(B-V)=3.1
!  XLAMBDA_DERIV = 2nd derivative of XLAMBDA_GRID values from SPLINE
!  LAMBDA_GRID   = wavelengths (in Å) listed in Table 2
!  NUM_LAMBDA    = number of wavelengths
!  START         = .TRUE. when XLAMBDA is first called
!
!  Functions called:
!  SPLINE = second derivative of the supplied 1D function (in spline.f90)
!  SPLINT = 1-dimensional cubic spline interpolated value (in spline.f90)
!
!-----------------------------------------------------------------------
FUNCTION XLAMBDA(LAMBDA)

   USE HEALPIX_TYPES
   USE SHIELDING_MODULE, ONLY : START,NUM_LAMBDA,LAMBDA_GRID,XLAMBDA_GRID,XLAMBDA_DERIV
   IMPLICIT NONE

   REAL(KIND=DP) :: XLAMBDA
   REAL(KIND=DP), INTENT(IN) :: LAMBDA

   REAL(KIND=DP) :: LAMBDA_VALUE

   IF(START) THEN
      CALL SPLINE(LAMBDA_GRID,XLAMBDA_GRID,NUM_LAMBDA,1.0D30,1.0D30,XLAMBDA_DERIV)
   ENDIF

   LAMBDA_VALUE=LAMBDA
   IF(LAMBDA.LT.LAMBDA_GRID(1)) LAMBDA_VALUE=LAMBDA_GRID(1)
   IF(LAMBDA.GT.LAMBDA_GRID(NUM_LAMBDA)) LAMBDA_VALUE=LAMBDA_GRID(NUM_LAMBDA)

   CALL SPLINT(LAMBDA_GRID,XLAMBDA_GRID,XLAMBDA_DERIV,NUM_LAMBDA,LAMBDA_VALUE,XLAMBDA)
   IF(XLAMBDA.LT.0.0D0) XLAMBDA=0.0D0

   RETURN
END FUNCTION XLAMBDA
!=======================================================================

!=======================================================================
!
!  Calculate the mean wavelength (in Å) of the 33 dissociating bands of
!  CO, weighted by their fractional contribution to the total shielding
!  from van Dishoeck & Black (1988, ApJ, 334, 771, Equation 4)
!
!-----------------------------------------------------------------------
!
!  Input parameters:
!  NCO = CO column density (in cm^-2)
!  NH2 = H2 column density (in cm^-2)
!
!  Program variables:
!  LBAR = mean wavelength (in Å)
!  U    = log10(NCO)
!  W    = log10(NH2)
!
!-----------------------------------------------------------------------
FUNCTION LBAR(NCO,NH2)

   USE HEALPIX_TYPES
   IMPLICIT NONE

   REAL(KIND=DP) :: LBAR
   REAL(KIND=DP), INTENT(IN) :: NCO,NH2

   REAL(KIND=DP) :: U,W

   U=DLOG10(NCO+1.0D0)
   W=DLOG10(NH2+1.0D0)

   LBAR=(5675.0D0 - 200.6D0*W) &
     & - (571.6D0 - 24.09D0*W)*U &
     & + (18.22D0 - 0.7664D0*W)*U**2

!  LBAR cannot be smaller than the wavelength of band 1 (913.6Å)
!  and cannot be larger than the wavelength of band 33 (1076.1Å)
   IF(LBAR.LT.913.6D0)  LBAR=913.6D0
   IF(LBAR.GT.1076.1D0) LBAR=1076.1D0

   RETURN
END FUNCTION LBAR
!=======================================================================

max_stars_repo_name: uclchem/uclpdr

max_issues_repo_name: uclchem/uclpdr

max_forks_repo_name: uclchem/uclpdr


text:
c
c ----------------------------------------------------
c
      subroutine domprep(domflags,lbase,ibase,jbase,kbase)

      use amr_module
      implicit double precision (a-h, o-z)

      integer*1 domflags(0:ibase+1,0:jbase+1,0:kbase+1)

c
c ::::::::::::::::::::::::::: PREPDOM :::::::::::::::::::::
c
c  prepare 3 dimensional array of domain for proper nesting
c
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::

      do 10 k = 0, kbase+1
      do 10 j = 0, jbase+1
      do 10 i = 0, ibase+1
        domflags(i,j,k) = 0
 10   continue

      mptr = lstart(lbase)
 15   continue
      do 20 k = node(ndklo,mptr) + 1, node(ndkhi,mptr) + 1
      do 20 j = node(ndjlo,mptr) + 1, node(ndjhi,mptr) + 1
      do 20 i = node(ndilo,mptr) + 1, node(ndihi,mptr) + 1
        domflags(i,j,k) = 1
 20   continue
      mptr = node(levelptr, mptr)
      if (mptr .ne. 0) go to 15

c
c take care of periodic domains or if border of domain touches a
c  physical boundary then set domain in ghost cell as well
c
      if (xperdom) then
         do 25 k = 0, kbase+1
         do 25 j = 0, jbase+1
          domflags(      0,j,k) = domflags(ibase,j,k)
          domflags(ibase+1,j,k) = domflags(    1,j,k)
 25      continue
       else
         do 65 k = 1, kbase
         do 65 j = 1, jbase
           if (domflags(    1,j,k) .eq. 1) domflags(      0,j,k) = 1
           if (domflags(ibase,j,k) .eq. 1) domflags(ibase+1,j,k) = 1
 65      continue
      endif
      if (yperdom) then
         do 35 k = 0, kbase+1
         do 35 i = 0, ibase+1
           domflags(i,      0,k) = domflags(i,jbase,k)
           domflags(i,jbase+1,k) = domflags(i,    1,k)
 35      continue
       else
         do 55 k = 1, kbase
         do 55 i = 1, ibase
            if (domflags(i,    1,k) .eq. 1) domflags(i,      0,k) = 1
            if (domflags(i,jbase,k) .eq. 1) domflags(i,jbase+1,k) = 1
 55      continue
      endif
      if (zperdom) then
         do 45 j = 0, jbase+1
         do 45 i = 0, ibase+1
           domflags(i,j,      0) = domflags(i,j,kbase)
           domflags(i,j,kbase+1) = domflags(i,j,    1)
 45      continue
      else
         do 75 j = 1, jbase
         do 75 i = 1, ibase
           if (domflags(i,j,    1) .eq. 1) domflags(i,j,      0) = 1
           if (domflags(i,j,kbase) .eq. 1) domflags(i,j,kbase+1) = 1
 75       continue
      endif
c
c     corner stuff
c
      do k = 1, kbase
        if (domflags(0,1,k)+domflags(1,0,k).eq.2)  domflags(0,0,k) = 1
        if (domflags(0,jbase,k)+domflags(1,jbase+1,k).eq.2)
     .                                       domflags(0,jbase+1,k) = 1
        if (domflags(ibase,0,k)+domflags(ibase+1,1,k).eq.2)
     .                                       domflags(ibase+1,0,k) = 1
        if (domflags(ibase,jbase+1,k)+domflags(ibase+1,jbase,k).eq.2)
     .                                 domflags(ibase+1,jbase+1,k) = 1
      end do
      do i = 1, ibase
        if (domflags(i,0,1)+domflags(i,1,0).eq.2)  domflags(i,0,0) = 1
        if (domflags(i,jbase,0)+domflags(i,jbase+1,1).eq.2)
     .                                       domflags(i,jbase+1,0) = 1
        if (domflags(i,0,kbase)+domflags(i,1,kbase+1).eq.2)
     .                                       domflags(i,0,kbase+1) = 1
        if (domflags(i,jbase+1,kbase)+domflags(i,jbase,kbase+1).eq.2)
     .                                 domflags(i,jbase+1,kbase+1) = 1
      end do
      do j = 1, jbase
        if (domflags(0,j,1)+domflags(1,j,0).eq.2)  domflags(0,j,0) = 1
        if (domflags(ibase,j,0)+domflags(ibase+1,j,1).eq.2)
     .                                       domflags(ibase+1,j,0) = 1
        if (domflags(0,j,kbase)+domflags(1,j,kbase+1).eq.2)
     .                                       domflags(0,j,kbase+1) = 1
        if (domflags(ibase+1,j,kbase)+domflags(ibase,j,kbase+1).eq.2)
     .                                 domflags(ibase+1,j,kbase+1) = 1
      end do
c
c     The neigbors of a corner are the three locations each having an
c     address which differs from that of the corner by 1 in just one of
c     the three indexes. If a corner index is 0, the corresponding neighbor
c     index may be 0 or 1. If a corner index is base+1, the corresponding
c     neighbor index may be base+1 or base (where base => [i,j,k]base).
c
      if (( domflags(      0,      0,      1)
     &     +domflags(      0,      1,      0)
     &     +domflags(      1,      0,      0)) .eq. 3)
     .      domflags(      0,      0,      0)   =   1

      if (( domflags(ibase+1,      0,      1)
     &     +domflags(ibase+1,      1,      0)
     &     +domflags(ibase  ,      0,      0)) .eq. 3)
     .      domflags(ibase+1,      0,      0)   =   1

      if (( domflags(      0,jbase+1,      1)
     &     +domflags(      0,jbase  ,      0)
     &     +domflags(      1,jbase+1,      0)) .eq. 3)
     .      domflags(      0,jbase+1,      0)   =   1

      if (( domflags(ibase+1,jbase+1,      1)
     &     +domflags(ibase+1,jbase  ,      0)
     &     +domflags(ibase  ,jbase+1,      0)) .eq. 3)
     .      domflags(ibase+1,jbase+1,      0)   =   1

      if (( domflags(      0,      0,kbase  )
     &     +domflags(      0,      1,kbase+1)
     &     +domflags(      1,      0,kbase+1)) .eq. 3)
     .      domflags(      0,      0,kbase+1)   =   1

      if (( domflags(ibase+1,      0,kbase  )
     &     +domflags(ibase+1,      1,kbase+1)
     &     +domflags(ibase  ,      0,kbase+1)) .eq. 3)
     .      domflags(ibase+1,      0,kbase+1)   =   1

      if (( domflags(      0,jbase+1,kbase  )
     &     +domflags(      0,jbase  ,kbase+1)
     &     +domflags(      1,jbase+1,kbase+1)) .eq. 3)
     .      domflags(      0,jbase+1,kbase+1)   =   1

      if (( domflags(ibase+1,jbase+1,kbase  )
     &     +domflags(ibase+1,jbase  ,kbase+1)
     &     +domflags(ibase  ,jbase+1,kbase+1)) .eq. 3)
     .      domflags(ibase+1,jbase+1,kbase+1)   =   1

      if (dprint) then
         write(outunit,*)" from domprep: domflags at level  ", lbase
         do 40 kk = 1, kbase
         k = kbase + 1 - kk
         write(outunit,*) 'plane k = ',k
         do 40 jj = 1, jbase
         j = jbase + 1 - jj
         write(outunit,100)(domflags(i,j,k),i=1,ibase)
 100     format(80i1)
 40      continue
      endif

      return
      end

ext: f

max_stars_repo_name: ian-r-rose/visualization

max_issues_repo_name: ian-r-rose/visualization

max_forks_repo_name: ian-r-rose/visualization


text:
module ufsLandNoahMPDriverModule

implicit none

contains

subroutine ufsLandNoahMPDriverInit(namelist, static, forcing, noahmp)

  use NamelistRead
  use ufsLandNoahMPType
  use ufsLandStaticModule
  use ufsLandInitialModule
  use ufsLandForcingModule
  use ufsLandNoahMPRestartModule

  implicit none

  type (namelist_type)       :: namelist
  type (noahmp_type)         :: noahmp
  type (static_type)         :: static
  type (initial_type)        :: initial
  type (forcing_type)        :: forcing
  type (noahmp_restart_type) :: restart

  call static%ReadStatic(namelist)

  call noahmp%Init(namelist,namelist%lensub)

  if(namelist%restart_simulation) then
    call restart%ReadRestartNoahMP(namelist, noahmp)
  else
    call initial%ReadInitial(namelist)
    call initial%TransferInitialNoahMP(namelist, noahmp)
  end if

  call static%TransferStaticNoahMP(noahmp)

  call noahmp%TransferNamelist(namelist)

  call forcing%ReadForcingInit(namelist)

end subroutine ufsLandNoahMPDriverInit

subroutine ufsLandNoahMPDriverRun(namelist, static, forcing, noahmp)

use machine , only : kind_phys
use noahmpdrv
use set_soilveg_mod
use funcphys
use namelist_soilveg, only : z0_data
use physcons, only : con_hvap , con_cp, con_jcal, con_eps, con_epsm1,    &
                     con_fvirt, con_rd, con_hfus,                        &
		     tfreeze=> con_t0c, rhoh2o => rhowater

use interpolation_utilities
use time_utilities
use cosine_zenith
use NamelistRead, only         : namelist_type
use ufsLandNoahMPType, only    : noahmp_type
use ufsLandStaticModule, only  : static_type
use ufsLandForcingModule
use ufsLandIOModule
use ufsLandNoahMPRestartModule

type (namelist_type)  :: namelist
type (noahmp_type)    :: noahmp
type (forcing_type)   :: forcing
type (static_type)    :: static
type (output_type)    :: output
type (noahmp_restart_type)    :: restart

integer          :: timestep
double precision :: now_time
character*19     :: now_date  ! format: yyyy-mm-dd hh:nn:ss
integer          :: now_yyyy

real, allocatable, dimension(:) :: rho
real(kind=kind_phys), parameter :: one     = 1.0_kind_phys

associate (                               &
   ps         => forcing%surface_pressure    ,&
   t1         => forcing%temperature         ,&
   q1         => forcing%specific_humidity   ,&
   dlwflx     => forcing%downward_longwave   ,&
   dswsfc     => forcing%downward_shortwave  ,&
   wind       => forcing%wind_speed          ,&
   tprcp      => forcing%precipitation       ,&
   im         => noahmp%static%im        ,&
   km         => noahmp%static%km        ,&
   itime      => noahmp%static%itime     ,&
   delt       => noahmp%static%delt      ,&
   isot       => noahmp%static%isot      ,&
   ivegsrc    => noahmp%static%ivegsrc   ,&
   idveg      => noahmp%static%idveg     ,&
   iopt_crs   => noahmp%static%iopt_crs  ,&
   iopt_btr   => noahmp%static%iopt_btr  ,&
   iopt_run   => noahmp%static%iopt_run  ,&
   iopt_sfc   => noahmp%static%iopt_sfc  ,&
   iopt_frz   => noahmp%static%iopt_frz  ,&
   iopt_inf   => noahmp%static%iopt_inf  ,&
   iopt_rad   => noahmp%static%iopt_rad  ,&
   iopt_alb   => noahmp%static%iopt_alb  ,&
   iopt_snf   => noahmp%static%iopt_snf  ,&
   iopt_tbot  => noahmp%static%iopt_tbot ,&
   iopt_stc   => noahmp%static%iopt_stc  ,&
   iopt_rsf   => noahmp%static%iopt_rsf  ,&
   iopt_gla   => noahmp%static%iopt_gla  ,&
   errmsg     => noahmp%static%errmsg    ,&
   errflg     => noahmp%static%errflg    ,&
   soiltyp    => noahmp%model%soiltyp    ,&
   vegtype    => noahmp%model%vegtype    ,&
   slopetyp   => noahmp%model%slopetyp   ,&
   u1         => noahmp%model%u1         ,&
   v1         => noahmp%model%v1         ,&
   sigmaf     => noahmp%model%sigmaf     ,&
   emiss      => noahmp%model%emiss      ,&
   albdvis    => noahmp%model%albdvis    ,&
   albdnir    => noahmp%model%albdnir    ,&
   albivis    => noahmp%model%albivis    ,&
   albinir    => noahmp%model%albinir    ,&
   snet       => noahmp%model%snet       ,&
   tg3        => noahmp%model%tg3        ,&
   cm         => noahmp%model%cm         ,&
   ch         => noahmp%model%ch         ,&
   prsl1      => noahmp%model%prsl1      ,&
   shdmin     => noahmp%model%shdmin     ,&
   shdmax     => noahmp%model%shdmax     ,&
   snoalb     => noahmp%model%snoalb     ,&
   sfalb      => noahmp%model%sfalb      ,&
   zf         => noahmp%model%zf         ,&
   flag_iter  => noahmp%model%flag_iter  ,&
   flag_guess => noahmp%model%flag_guess ,&
   dry        => noahmp%model%dry        ,&
   prslki     => noahmp%model%prslki     ,&
   weasd      => noahmp%model%weasd      ,&
   snwdph     => noahmp%model%snwdph     ,&
   tskin      => noahmp%model%tskin      ,&
   srflag     => noahmp%model%srflag     ,&
   canopy     => noahmp%model%canopy     ,&
   trans      => noahmp%model%trans      ,&
   tsurf      => noahmp%model%tsurf      ,&
   zorl       => noahmp%model%zorl       ,&
   smc        => noahmp%model%smc        ,&
   stc        => noahmp%model%stc        ,&
   slc        => noahmp%model%slc        ,&
   sncovr1    => noahmp%model%sncovr1    ,&
   qsurf      => noahmp%model%qsurf      ,&
   gflux      => noahmp%model%gflux      ,&
   drain      => noahmp%model%drain      ,&
   evap       => noahmp%model%evap       ,&
   hflx       => noahmp%model%hflx       ,&
   ep         => noahmp%model%ep         ,&
   runoff     => noahmp%model%runoff     ,&
   cmm        => noahmp%model%cmm        ,&
   chh        => noahmp%model%chh        ,&
   evbs       => noahmp%model%evbs       ,&
   evcw       => noahmp%model%evcw       ,&
   sbsno      => noahmp%model%sbsno      ,&
   snowc      => noahmp%model%snowc      ,&
   stm        => noahmp%model%stm        ,&
   snohf      => noahmp%model%snohf      ,&
   smcwlt2    => noahmp%model%smcwlt2    ,&
   smcref2    => noahmp%model%smcref2    ,&
   wet1       => noahmp%model%wet1       ,&
   xlatin     => noahmp%model%xlatin     ,&
   xcoszin    => noahmp%model%xcoszin    ,&
   iyrlen     => noahmp%model%iyrlen     ,&
   julian     => noahmp%model%julian     ,&
   rainn_mp   => noahmp%model%rainn_mp   ,&
   rainc_mp   => noahmp%model%rainc_mp   ,&
   snow_mp    => noahmp%model%snow_mp    ,&
   graupel_mp => noahmp%model%graupel_mp ,&
   ice_mp     => noahmp%model%ice_mp     ,&
   snowxy     => noahmp%model%snowxy     ,&
   tvxy       => noahmp%model%tvxy       ,&
   tgxy       => noahmp%model%tgxy       ,&
   canicexy   => noahmp%model%canicexy   ,&
   canliqxy   => noahmp%model%canliqxy   ,&
   eahxy      => noahmp%model%eahxy      ,&
   tahxy      => noahmp%model%tahxy      ,&
   cmxy       => noahmp%model%cmxy       ,&
   chxy       => noahmp%model%chxy       ,&
   fwetxy     => noahmp%model%fwetxy     ,&
   sneqvoxy   => noahmp%model%sneqvoxy   ,&
   alboldxy   => noahmp%model%alboldxy   ,&
   qsnowxy    => noahmp%model%qsnowxy    ,&
   wslakexy   => noahmp%model%wslakexy   ,&
   zwtxy      => noahmp%model%zwtxy      ,&
   waxy       => noahmp%model%waxy       ,&
   wtxy       => noahmp%model%wtxy       ,&
   tsnoxy     => noahmp%model%tsnoxy     ,&
   zsnsoxy    => noahmp%model%zsnsoxy    ,&
   snicexy    => noahmp%model%snicexy    ,&
   snliqxy    => noahmp%model%snliqxy    ,&
   lfmassxy   => noahmp%model%lfmassxy   ,&
   rtmassxy   => noahmp%model%rtmassxy   ,&
   stmassxy   => noahmp%model%stmassxy   ,&
   woodxy     => noahmp%model%woodxy     ,&
   stblcpxy   => noahmp%model%stblcpxy   ,&
   fastcpxy   => noahmp%model%fastcpxy   ,&
   xlaixy     => noahmp%model%xlaixy     ,&
   xsaixy     => noahmp%model%xsaixy     ,&
   taussxy    => noahmp%model%taussxy    ,&
   smoiseq    => noahmp%model%smoiseq    ,&
   smcwtdxy   => noahmp%model%smcwtdxy   ,&
   deeprechxy => noahmp%model%deeprechxy ,&
   rechxy     => noahmp%model%rechxy     ,&
   t2mmp      => noahmp%model%t2mmp      ,&
   q2mp       => noahmp%model%q2mp        &
   )

allocate(   rho(im))

flag_iter  = .true.
flag_guess = .false.

call set_soilveg(0,isot,ivegsrc,0)
call gpvs()

zorl     = z0_data(vegtype) * 100.0   ! at driver level, roughness length in cm

time_loop : do timestep = 1, namelist%run_timesteps

  now_time = namelist%initial_time + timestep * namelist%timestep_seconds

  call date_from_since("1970-01-01 00:00:00", now_time, now_date)
  read(now_date(1:4),'(i4)') now_yyyy
  iyrlen = 365
  if(mod(now_yyyy,4) == 0) iyrlen = 366

  if(.not.namelist%restart_simulation .and. timestep == 1) &
     call noahmp%InitStates(namelist, now_time)

  call forcing%ReadForcing(namelist, static, now_time)

  call interpolate_monthly(now_time, im, static%gvf_monthly, sigmaf)
  call interpolate_monthly(now_time, im, static%albedo_monthly, sfalb)

  call calc_cosine_zenith(now_time, im, static%latitude, static%longitude, xcoszin, julian)

  u1 = wind
  v1 = 0.0_kind_phys
  snet   = dswsfc * (1.0_kind_phys - sfalb)
  prsl1  = ps * exp(-1.d0*zf/29.25d0/t1)
  srflag = 0.0d0
  where(t1 < tfreeze) srflag = 1.d0
  prslki = (exp(zf/29.25d0/t1))**(2.d0/7.d0)

  rainn_mp = 1000.0 * tprcp / delt
  rainc_mp = 0.0
  snow_mp = 0.0
  graupel_mp = 0.0
  ice_mp = 0.0

      call noahmpdrv_run                                               &
          ( im, km, itime, ps, u1, v1, t1, q1, soiltyp, vegtype,       &  !  ---  inputs
            sigmaf, dlwflx, dswsfc, snet, delt, tg3, cm, ch,           &
            prsl1, prslki, zf, dry, wind, slopetyp,                    &
            shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess,      &
            idveg, iopt_crs, iopt_btr, iopt_run, iopt_sfc, iopt_frz,   &
            iopt_inf, iopt_rad, iopt_alb, iopt_snf, iopt_tbot,         &
            iopt_stc, xlatin, xcoszin, iyrlen, julian,                 &
            rainn_mp, rainc_mp, snow_mp, graupel_mp, ice_mp,           &
            con_hvap, con_cp, con_jcal, rhoh2o, con_eps, con_epsm1,    &
            con_fvirt, con_rd, con_hfus,                               &
            weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc,        &  !  ---  in/outs:
            canopy, trans, tsurf, zorl,                                &
            snowxy, tvxy, tgxy, canicexy, canliqxy, eahxy, tahxy, cmxy,&  ! --- Noah MP specific
            chxy, fwetxy, sneqvoxy, alboldxy, qsnowxy, wslakexy, zwtxy,&
            waxy, wtxy, tsnoxy, zsnsoxy, snicexy, snliqxy, lfmassxy,   &
            rtmassxy, stmassxy, woodxy, stblcpxy, fastcpxy, xlaixy,    &
            xsaixy, taussxy, smoiseq, smcwtdxy, deeprechxy, rechxy,    &
            albdvis, albdnir,  albivis,  albinir,emiss,                &
            sncovr1, qsurf, gflux, drain, evap, hflx, ep, runoff,      &  !  ---  outputs:
            cmm, chh, evbs, evcw, sbsno, snowc, stm, snohf,            &
            smcwlt2, smcref2, wet1, t2mmp, q2mp, errmsg, errflg)

  rho = prsl1 / (con_rd*t1*(one+con_fvirt*q1))
  hflx = hflx * rho * con_cp
  evap = evap * rho * con_hvap

  where(dswsfc>0.0 .and. sfalb<0.0) dswsfc = 0.0

  call output%WriteOutputNoahMP(namelist, noahmp, forcing, now_time)

  if(namelist%restart_timesteps > 0) then
    if(mod(timestep,namelist%restart_timesteps) == 0) then
      call restart%WriteRestartNoahMP(namelist, noahmp, now_time)
    end if
  end if

  if(errflg /= 0) then
    write(*,*) "noahmpdrv_run reporting an error"
    write(*,*) errmsg
    stop
  end if

end do time_loop

end associate

end subroutine ufsLandNoahMPDriverRun

subroutine ufsLandNoahMPDriverFinalize()
end subroutine ufsLandNoahMPDriverFinalize

end module ufsLandNoahMPDriverModule

max_stars_repo_name: JCSDA-internal/UCLDAS

max_issues_repo_name: JCSDA-internal/UCLDAS

max_forks_repo_name: JCSDA-internal/ucldas


text:
program t
  character(len=5)::a='hello'
  character(5)::b
  character*5::c
  b='there'
  c='world'
  print '(a)',a
  print '(a)',b
  print '(a)',c
end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
  subroutine eventLoop(gsmObj, gsmRxn)

! ======================================================================
!
! Main physics loop for spallation events - simulates INC, coalescence,
! preequilibrium, evaporation, and Fermi Break-up emission where
! appropriate
!
!
! Written by CMJ, XCP-3 (03/2019); parsed from original "gsm03.f90" routine
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, int64, real64, error_unit
    use gsm_params, only: zro, hlf, one, two
    use randomNumberGenerator, only: RN_next_particle

    ! Some objects used (mostly for setup):
    use standardDCMData,  only: photonEG
    use evaporationClass, only : setLevelDensity

    ! For initializing things...
    use hist_mod, only: hist_tally

! ======================================================================

    implicit none
    class(GSM),        intent(inout) :: gsmObj
    type(GSMReaction), intent(inout), target :: gsmRxn

    integer(int32) :: ibadren, in, iz, ln
    real(real64)   :: bf0, delu, erotev, pevap, &
         & u, ue, wam
    logical :: inelasticEvent = .TRUE.
    logical :: momentumConserved = .TRUE.
    logical :: mDCMUsed = .FALSE.

    type(GSMResults),    pointer :: results => NULL()
    type(GSMProjectile), pointer :: proj    => NULL()
    type(GSMTarget),     pointer :: targ    => NULL()

! ======================================================================

    results => gsmRxn%results
    proj    => results%initialProj
    targ    => results%initialTarg
    mDCMUsed = (gsmRxn%incFlag == mDCMFlagged)

30  continue   ! for cem reference
35  continue
    inelasticEvent = .TRUE.

    ! Reset simulation results object:
    call results%resetEvent()

    ! ----- for bremss. photons only ----- (sample incident energy)
    !   Sampling the energy of gamma from Schiff's spectra (1/E) and
    !   preparing g + N --> pi + N angular distributions:
    !   [The latter moved to TYPINT  02/25/05]
    if (results%initialProj%particleFlag == bremsProjFlag) then
       ! Sample energy and update energy for sDCM/mDCM simulation
       proj%kinEnergy = proj%brems%sampleEnergy(gsmObj%rang())
       gsmRxn%sDCMProj%kinenergy = proj%kinEnergy
    endif

    ! Sample INC if smooth transitions are desired:
    if ( gsmObj%options%smoothTransition .and. proj%numBaryons <= 1 ) then
       gsmRxn%incFlag = gsmObj%inquireINCModel( proj, targ )
       mDCMUsed = (gsmRxn%incFlag == mDCMFlagged)
    end if

    ! Performing IntraNuclear Cascade (INC) calculation
    select case( gsmRxn%incFlag )
       case( sDCMFlagged )
          call gsmObj%standardDCMInterface (gsmRxn%data, gsmRxn%sDCMProj, results)
       case( mDCMFlagged )
          !>>> CRITICAL REGION
          call gsmObj%modifiedDCMInterface( gsmRxn )
          !>>> END CRITICAL REGION
       case default
          call gsmObj%standardDCMInterface (gsmRxn%data, gsmRxn%sDCMProj, results)
    end select

    ! THESE 'IF' BLOCKS CHECK FOR PHYSICALITY (IF NOT, 'RESETS' THE EVENT
    ! [GO TO 30])

    ! Performing Coalescence, Pre-Equil, Equil, and Fermi Break-Up OR
    ! Evaporation/Fission when appropriate; unallowed reaction re-runs cascade.
    if (results%targRes%kinEnergy <= zro) then
       inelasticEvent = .FALSE.

       ! Negative Exit energy (elastic event); run cascade again (increment num. elastic events with CEM cascade only)
       if ( mDCMUsed ) then
          ! Print warning to user (shouldn't reach here)
          write(gsmObj%io%message, 1000) "An unexpected elastic interaction was flagged."
          call gsmObj%io%print(5, 3, gsmObj%io%message)
       else
          ! CEM cascade allows elastic events in this manner
          if (associated(gsmRxn%outData)) gsmRxn%outData%intel = gsmRxn%outData%intel + 1
       end if
    else
       if (results%numElasticEvents > 0 .and. .not.mDCMUsed) then
          ! Elastic Event; run cascade again (no momentum imparted)
          if (associated(gsmRxn%outData)) gsmRxn%outData%intel = gsmRxn%outData%intel + 1
          inelasticEvent = .FALSE.
       else
          ! Inelastic Event occurred, obtaining data
          ! Rerun cascade if an unphysical particle nucleus was produced via CASCAD:
          if ( ( results%targRes%numBaryons < 4.d0 .or. &
               & results%targRes%numProtons < one  .or. &
               & results%targRes%numBaryons < results%targRes%numProtons) .and. &
               & .not.(mDCMUsed) ) then
             write(gsmObj%io%message,1100) &
                  & "The target size (A<=4) is not allowed when using the Standard DCM INC model in GSM", gsmRxn%eventNum
             call gsmObj%io%print(2, 3, gsmObj%io%message)
             go to 30
          endif

          ! ABOVE PHYSICS IS VALID - CONTINUE ON
          ! Renormalize residuals to ensure physicality
          if ( .not.mDCMUsed ) then
             ! NOTE: renorm() is NOT used with LAQGSM's INC becuase the recalculation of
             !       energy is not performed for the projectile, and in addition would
             !       require modification of the equations used. Mass/charge conservation
             !       is observed some of the time with the use of the LAQGSM INC, however
             !       not often. 'renorm()' has been updated to allow mass/charge conservation
             !       checking with the LAQGSM INC, however should not be used as it will
             !       provide incorrect results (due to A/Z conservation failure in 'CASCAW()')
             !>>> ASSIGN TO AN OPENMP TASK?
             !>>> {
             call gsmObj%renorm (results, ibadren)

             ! Bad Re-Normalization; conservation laws broken
             if ( ibadren > 0) then
                ! call gsmObj%printPartProp (proj, targ, outData%ncas, intel)  ! for debug.
                write(gsmObj%io%message,1100) "A bad renormalization occurred", gsmRxn%eventNum
                call gsmObj%io%print(4, 3, gsmObj%io%message)
                go to 30
             endif
             !>>> }
          endif

          !>>> ASSIGN TO AN OPENMP TASK?
          !>>> {
          ! Form compound nuclei from the emitted INC progeny
          call gsmObj%coalescenceInterface (gsmRxn%models, results, gsmRxn%outData)

          ! Check for momentum conservation:
          call gsmObj%checkMomentum(results, momentumConserved)
          if ( .not. momentumConserved ) then
             write(gsmObj%io%message,1100) "Conservation of momentum failed", gsmRxn%eventNum
             call gsmObj%io%print(4, 3, gsmObj%io%message)
             go to 30
          endif
          !>>> }

          ! Converting to MeV
          results%targRes%kinEnergy = results%targRes%kinEnergy*1000.d0
          results%projRes%kinEnergy = results%projRes%kinEnergy*1000.d0

          ! Checking for allowed residual particle; LAQGSM Checks for this above
          if ( ( results%targRes%numBaryons >= 4.d0 .and. &
               & results%targRes%numProtons >= one  .and. &
               & results%targRes%numBaryons >= results%targRes%numProtons) .or. &
               & (mDCMUsed) ) then

             iz = nint(results%targRes%numProtons)       ! number of protons/charged particles
             in = nint(results%targRes%numBaryons - results%targRes%numProtons)  ! number of neutrons/neutral particles
             pevap = gsmObj%genData%molnixenergies%defineenergy (iz, in, 3)

             ! Momenta (ln), Angular Momenta (amnucl), Energy, and Fission Barrier (bf) Book-keeping
             call gsmRxn%data%preeq%auxl (results%targRes%numBaryons, results%targRes%numProtons, &
                  & results%targRes%angularMom, bf0, ln, erotev, delu)

             u = results%targRes%kinEnergy + delu
             ue = u - pevap - erotev ! rotational energy (ev)

             !>>> CRITICAL REGION (FOR TIME BEING)
             !>>> {
             if(results%tallySim) then
                ! Statistics at Various Stages of the Reaction
                call gsmObj%ststcs (results%targRes%numBaryons, &
                     & results%targRes%numProtons, ue, ln, bf0, 1)

                if (gsmRxn%output%printHist) then
                   ! Tallying Particle Data
                   call hist_tally(results%targRes%linearMom, results%targRes%angularMom, &
                        & results%targRes%kinEnergy, results%targRes%numBaryons, &
                        & results%targRes%numProtons, results%targExc%numNeutrons, &
                        & results%targExc%numProtons, results%targExc%numHoles )
                end if

                ! Valid event occurred, count event and notify user/client
                gsmRxn%outData%ncas = gsmRxn%outData%ncas + 1
                gsmRxn%eventNum = gsmRxn%eventNum + 1
                if ( mod(gsmRxn%eventNum, gsmObj%options%printIncrement) == 0 ) then
                   if ( gsmRxn%eventNum < 99999999 ) then
                      write (gsmObj%io%message, 2000) gsmRxn%eventNum
                   else
                      write (gsmObj%io%message, 2050) real(gsmRxn%eventNum)
                   endif
                   call gsmObj%io%print(2, 5, gsmObj%io%message)
                endif
             end if
             !>>> }

             ! ===================================================
             ! Tally Projectile (where appropriate)
             ! ===================================================
             if ( mDCMUsed ) then
                ! Count scores from residual projectile nucleus where appropriate
                if ( gsmObj%options%tallyResidualProjectile .and. &
                     & results%projRes%kinEnergy > zro .and. &
                     & results%projRes%numBaryons > 0.9 ) then

                   iz = nint(results%projRes%numProtons)
                   in = nint(results%projRes%numBaryons - results%projRes%numProtons)
                   pevap = gsmObj%genData%molnixenergies%defineenergy (iz, in, 3)

                   ! Momenta (ln), Angular Momenta (amnucl), Energy, and Fission Barrier (bf) Book-keeping
                   call gsmRxn%data%preeq%auxl (results%projRes%numBaryons, &
                        & results%projRes%numProtons, &
                        & results%projRes%angularMom, bf0, ln, erotev, delu)
                   u = results%projRes%kinEnergy + delu
                   ue = u - pevap - erotev ! rotational energy (ev)

                   !>>> CRITICAL REGION (OPENMP) FOR NOW
                   !>>> {
                   if(results%tallySim) then
                      ! Statistics at Various Stages of the Reaction
                      call gsmObj%ststcs (results%projRes%numBaryons, &
                           & results%projRes%numProtons, ue, ln, bf0, 1)

                      if (gsmRxn%output%printHist) then
                         ! Tallying Particle Data
                         call hist_tally(results%projRes%linearMom, results%projRes%angularMom, &
                              & results%projRes%kinEnergy, results%projRes%numBaryons, &
                              & results%projRes%numProtons, results%projExc%numNeutrons, &
                              & results%projExc%numProtons, results%projExc%numHoles )
                      end if
                   end if
                   !>>> }

                   ! Check for fermi break-results%projRes%kinEnergy or other emission
                   if ( gsmObj%genModels%fbu%recommendFermiBreakUp(results%projRes%numBaryons) ) then

                      ! Fermi Break-Up Driving Routine
                      call gsmObj%fermiBreakUpInterface (results%projRes, &
                           & gsmRxn)
                   else
                      ! Preequilibrium/Evaporation decay
                      results%info%wf = zro
                      wam = one

                      ! Calculating Pre-Equilibrium and Equilibrium Particle Emssion
                      call gsmObj%simulateDecay (gsmRxn, &
                           & results%projRes, results%projExc, &
                           & proj%afMultiplier, &
                           & proj%czMultiplier, &
                           & wam )

                      if (wam == -13.d0) then
                         wam = one
                         write(gsmObj%io%message,1100) "Projectile pre-equilibrium decay failure", gsmRxn%eventNum
                         call gsmObj%io%print(3, 3, gsmObj%io%message)

                         go to 30
                      endif
                   endif

                   ! Tally the projectile separately
                   results%projRes%kinEnergy = results%projRes%kinEnergy/1000.d0

                   ! if (gsmRxn%eventNum <= nnnp) call gsmObj%printPartProp (proj, &
                   !      & targ, gsmRxn%outData%ncas, intel)
                   if (results%info%fusion > zro .and. &
                        & gsmRxn%output%accountFission > 0) &
                        & gsmRxn%outData%nfis = gsmRxn%outData%nfis + 1
                   gsmRxn%outData%sfu = gsmRxn%outData%sfu + results%info%wf

                elseif ( gsmObj%options%tallyResidualProjectile ) then

                   ! Proj. residual not physical, check for errors
                   ! Check for physical residual projectile particle
                   if ( results%projRes%kinEnergy < zro .or. results%projRes%numProtons < zro ) then

                      ! Unphysical residual produced, rerun the cascade
                      write(gsmObj%io%message,1100) "An unphysical projectile residual was detected", gsmRxn%eventNum
                      call gsmObj%io%print(3, 3, gsmObj%io%message)

                      if(results%tallySim) then
                         ! Decrement counters
                         gsmRxn%outData%ncas = gsmRxn%outData%ncas - 1
                         gsmRxn%eventNum = gsmRxn%eventNum - 1
                      end if
                     if (associated(gsmRxn%outData)) gsmRxn%outData%intel = &
                          & gsmRxn%outData%intel + 1
                      inelasticEvent = .FALSE.

                      go to 35

                   endif
                endif
             endif

             ! ==========================
             ! Tally Target Residual
             ! ==========================
             if ( gsmObj%genModels%fbu%recommendFermiBreakUp(results%targRes%numBaryons) ) then
                ! Fermi Break-Up Driving Routine
                call gsmObj%fermiBreakUpInterface (results%targRes, &
                     & gsmRxn)
             else
                !   Regular decay of target; ap > 12
                wam = one

                !   Determine semiempirical af/an fits for fission cross sections;
                !   gamma energy varies randomly, so need new value for each iteration:
                !   Corrected by NVMokhov, 10/08/03:
                !   NOTE:
                !     afMultiplier = a_f(CEM)/a_f(RAL)
                !     czMultiplier = C(Z)[CEM]/C(Z)[RAL]
                if (proj%particleFlag == bremsProjFlag) then
                   ! Determine use of CEM or LAQGSM values:
                   if( gsmObj%options%useSDCMAfCzMultipliers ) then
                      call setleveldensity( targ%numBaryons, targ%numProtons, &
                           & proj%kinEnergy, targ%afMultiplier, &
                           & targ%czMultiplier )
                   else
                      call setleveldensity( targ%numBaryons, targ%numProtons, &
                           & proj%kinEnergy, targ%afMultiplier, &
                           & targ%czMultiplier, 1 )
                   endif
                endif

                ! Calculating Pre-Equilibrium and Equilibrium Particle Emssion
                call gsmObj%simulateDecay (gsmRxn, &
                     & results%targRes, results%targExc, &
                     & targ%afMultiplier, &
                     & targ%czMultiplier, &
                     & wam )

                if (wam == -13.d0) then
                   wam = one
                   write(gsmObj%io%message,1100) "Target pre-equilibrium decay failure", gsmRxn%eventNum
                   call gsmObj%io%print(3, 3, gsmObj%io%message)
                   go to 30
                endif
             endif

             ! KKG  10/13/04
             if (proj%particleFlag == bremsProjFlag) then
                call proj%brems%incrementTEqv(proj%kinEnergy)
                call proj%brems%incrementSXAbs( &
                     & photonEG%photocrosssection( &
                     &      1000 * proj%kinEnergy, targ%numBaryons))
             endif

          else

             ! Unallowed A, Z values; run Cascade again, unphysical nucleus produced
             write(gsmObj%io%message,1100) "An unphysical target residual was detected", gsmRxn%eventNum
             call gsmObj%io%print(3, 3, gsmObj%io%message)
             go to 30

          endif

          results%targRes%kinEnergy = results%targRes%kinEnergy/1000.d0   ! converting to GeV
          ! if (gsmRxn%eventNum <= nnnp) call gsmObj%printPartProp (proj, targ, intel)
          if (results%info%fusion > zro .and. gsmRxn%output%accountFission > 0) &
               & gsmRxn%outData%nfis = gsmRxn%outData%nfis + 1
          gsmRxn%outData%sfu = gsmRxn%outData%sfu + results%info%wf

          if(results%tallySim) then
             !>>> CRITICAL REGION
             ! Accumulating Events; evaporation and fission model
             call gsmObj%vlobd( gsmRxn )
             !>>> END CRITICAL REGION
          end if
       endif
    endif

    !>>> THE RN_SEED (FROM RN_next_particle) NEEDS TO BE OPENMP PRIVATE
    ! Successfully simulated a single non-elastic event; go to next RN stride
    if ( inelasticEvent .and. gsmObj%useDefaultRNG ) then
       call RN_next_particle( gsmRxn%eventNum+1, 0_int64, gsmRxn%eventNum )
    end if

    return
! ======================================================================
1000 format (A)
1100 format (A, " (event ", i8, ").  Event will be restarted...")
2000 format (2x, 'nc = ', i8)
2050 format (2x, 'nc = ', es11.3)
! ======================================================================
  end subroutine eventLoop

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


submodule (h5fortran:write) writer
    !! This submodule is for writing 0-D..7-D data
    use hdf5, only: h5dwrite_f, H5T_STD_I64LE

    implicit none (type, external)

contains

    module procedure hdf_write_scalar
        integer(HID_T)  :: sid, did
        integer(HSIZE_T), allocatable :: dims(:)
        integer :: ier

        if(.not.self%is_open) error stop 'h5fortran:writer: file handle is not open'

        allocate(dims(0))

        did = 0
        sid = 0
        !! sentinel

        select type (value)
            type is (character(*))
                call h5ltmake_dataset_string_f(self%lid, dname, value, ier)
            type is (real(real64))
                !! NOTE: 0d does not use chunk_size
                call hdf_create(self,dname, H5T_NATIVE_DOUBLE, dims,sid,did, compact=compact)
                call h5dwrite_f(did, H5T_NATIVE_DOUBLE, value, dims, ier)
            type is (real(real32))
                call hdf_create(self,dname, H5T_NATIVE_REAL, dims,sid,did, compact=compact)
                call h5dwrite_f(did, H5T_NATIVE_REAL, value, dims, ier)
            type is (integer(int32))
                call hdf_create(self,dname, H5T_NATIVE_INTEGER, dims,sid,did, compact=compact)
                call h5dwrite_f(did, H5T_NATIVE_INTEGER, value, dims, ier)
            type is (integer(int64))
                call hdf_create(self,dname, H5T_STD_I64LE, dims,sid,did, compact=compact)
                call h5dwrite_f(did, H5T_STD_I64LE, value, dims, ier)
            class default
                error stop 'h5fortran: invalid data type'
        end select
        if (ier/=0) then
            write(stderr,*) 'h5fortran:ERROR: could not write ',dname, ' to ', self%filename
            error stop
        endif

        call hdf_wrapup(did, sid, ier)

        if (present(ierr)) ierr = ier
        if (check(ier, self%filename, dname) .and. .not.present(ierr)) error stop
    end procedure hdf_write_scalar

    module procedure hdf_write_1d
        integer(HID_T) :: did, sid, mem_sid
        integer(HSIZE_T) :: dims(rank(value))
        integer :: ier

        dims = shape(value)
        select type (value)
            type is (real(real64))
                call hdf_create(self,dname, H5T_NATIVE_DOUBLE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (real(real32))
                call hdf_create(self,dname, H5T_NATIVE_REAL, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int32))
                call hdf_create(self,dname, H5T_NATIVE_INTEGER, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int64))
                call hdf_create(self,dname, H5T_STD_I64LE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            class default
                error stop 'h5fortran:write:invalid data type'
        end select

        mem_sid = H5S_ALL_F !< default

        if(present(istart) .and. present(iend)) then
            if(present(stride)) then
                !! necessary to use this present check for Intel and GCC
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend, stride)
            else
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend)
            endif
        endif

        select type (value)
            type is (real(real64))
                call h5dwrite_f(did, H5T_NATIVE_DOUBLE, value, dims, ier, mem_sid, sid)
            type is (real(real32))
                call h5dwrite_f(did, H5T_NATIVE_REAL, value, dims, ier, mem_sid, sid)
            type is (integer(int32))
                call h5dwrite_f(did, H5T_NATIVE_INTEGER, value, dims, ier, mem_sid, sid)
            type is (integer(int64))
                call h5dwrite_f(did, H5T_STD_I64LE, value, dims, ier, mem_sid, sid)
        end select
        if (ier/=0) error stop 'h5fortran:ERROR: could not write ' // dname // ' to ' // self%filename

        call hdf_wrapup(did, sid, ier)

        if (present(ierr)) ierr = ier
        if (check(ier, self%filename, dname) .and. .not.present(ierr)) error stop
    end procedure hdf_write_1d

    module procedure hdf_write_2d
        integer(HID_T) :: did, sid, mem_sid
        integer(HSIZE_T) :: dims(rank(value))
        integer :: ier

        dims = shape(value)
        select type (value)
            type is (real(real64))
                call hdf_create(self,dname, H5T_NATIVE_DOUBLE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (real(real32))
                call hdf_create(self,dname, H5T_NATIVE_REAL, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int32))
                call hdf_create(self,dname, H5T_NATIVE_INTEGER, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int64))
                call hdf_create(self,dname, H5T_STD_I64LE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            class default
                error stop 'h5fortran:write:invalid data type'
        end select

        mem_sid = H5S_ALL_F !< default

        if(present(istart) .and. present(iend)) then
            if(present(stride)) then
                !! necessary to use this present check for Intel and GCC
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend, stride)
            else
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend)
            endif
        endif

        select type (value)
            type is (real(real64))
                call h5dwrite_f(did, H5T_NATIVE_DOUBLE, value, dims, ier, mem_sid, sid)
            type is (real(real32))
                call h5dwrite_f(did, H5T_NATIVE_REAL, value, dims, ier, mem_sid, sid)
            type is (integer(int32))
                call h5dwrite_f(did, H5T_NATIVE_INTEGER, value, dims, ier, mem_sid, sid)
            type is (integer(int64))
                call h5dwrite_f(did, H5T_STD_I64LE, value, dims, ier, mem_sid, sid)
        end select
        if (ier/=0) error stop 'h5fortran:ERROR: could not write ' // dname // ' to ' // self%filename

        call hdf_wrapup(did, sid, ier)

        if (present(ierr)) ierr = ier
        if (check(ier, self%filename, dname) .and. .not.present(ierr)) error stop
    end procedure hdf_write_2d

    module procedure hdf_write_3d
        integer(HID_T) :: did, sid, mem_sid
        integer(HSIZE_T) :: dims(rank(value))
        integer :: ier

        dims = shape(value)
        select type (value)
            type is (real(real64))
                call hdf_create(self,dname, H5T_NATIVE_DOUBLE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (real(real32))
                call hdf_create(self,dname, H5T_NATIVE_REAL, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int32))
                call hdf_create(self,dname, H5T_NATIVE_INTEGER, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int64))
                call hdf_create(self,dname, H5T_STD_I64LE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            class default
                error stop 'h5fortran:write:invalid data type'
        end select

        mem_sid = H5S_ALL_F !< default

        if(present(istart) .and. present(iend)) then
            if(present(stride)) then
                !! necessary to use this present check for Intel and GCC
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend, stride)
            else
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend)
            endif
        endif

        select type (value)
            type is (real(real64))
                call h5dwrite_f(did, H5T_NATIVE_DOUBLE, value, dims, ier, mem_sid, sid)
            type is (real(real32))
                call h5dwrite_f(did, H5T_NATIVE_REAL, value, dims, ier, mem_sid, sid)
            type is (integer(int32))
                call h5dwrite_f(did, H5T_NATIVE_INTEGER, value, dims, ier, mem_sid, sid)
            type is (integer(int64))
                call h5dwrite_f(did, H5T_STD_I64LE, value, dims, ier, mem_sid, sid)
        end select
        if (ier/=0) error stop 'h5fortran:ERROR: could not write ' // dname // ' to ' // self%filename

        call hdf_wrapup(did, sid, ier)

        if (present(ierr)) ierr = ier
        if (check(ier, self%filename, dname) .and. .not.present(ierr)) error stop
    end procedure hdf_write_3d

    module procedure hdf_write_4d
        integer(HID_T) :: did, sid, mem_sid
        integer(HSIZE_T) :: dims(rank(value))
        integer :: ier

        dims = shape(value)
        select type (value)
            type is (real(real64))
                call hdf_create(self,dname, H5T_NATIVE_DOUBLE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (real(real32))
                call hdf_create(self,dname, H5T_NATIVE_REAL, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int32))
                call hdf_create(self,dname, H5T_NATIVE_INTEGER, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int64))
                call hdf_create(self,dname, H5T_STD_I64LE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            class default
              error stop 'h5fortran:write:invalid data type'
        end select

        mem_sid = H5S_ALL_F !< default

        if(present(istart) .and. present(iend)) then
            if(present(stride)) then
                !! necessary to use this present check for Intel and GCC
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend, stride)
            else
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend)
            endif
        endif

        select type (value)
            type is (real(real64))
                call h5dwrite_f(did, H5T_NATIVE_DOUBLE, value, dims, ier, mem_sid, sid)
            type is (real(real32))
                call h5dwrite_f(did, H5T_NATIVE_REAL, value, dims, ier, mem_sid, sid)
            type is (integer(int32))
                call h5dwrite_f(did, H5T_NATIVE_INTEGER, value, dims, ier, mem_sid, sid)
            type is (integer(int64))
                call h5dwrite_f(did, H5T_STD_I64LE, value, dims, ier, mem_sid, sid)
        end select
        if (ier/=0) error stop 'h5fortran:ERROR: could not write ' // dname // ' to ' // self%filename

        call hdf_wrapup(did, sid, ier)

        if (present(ierr)) ierr = ier
        if (check(ier, self%filename, dname) .and. .not.present(ierr)) error stop
    end procedure hdf_write_4d

    module procedure hdf_write_5d
        integer(HID_T) :: did, sid, mem_sid
        integer(HSIZE_T) :: dims(rank(value))
        integer :: ier

        dims = shape(value)
        select type (value)
            type is (real(real64))
                call hdf_create(self,dname, H5T_NATIVE_DOUBLE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (real(real32))
                call hdf_create(self,dname, H5T_NATIVE_REAL, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int32))
                call hdf_create(self,dname, H5T_NATIVE_INTEGER, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int64))
                call hdf_create(self,dname, H5T_STD_I64LE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            class default
                error stop 'h5fortran:write:invalid data type'
        end select

        mem_sid = H5S_ALL_F !< default

        if(present(istart) .and. present(iend)) then
            if(present(stride)) then
                !! necessary to use this present check for Intel and GCC
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend, stride)
            else
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend)
            endif
        endif

        select type (value)
            type is (real(real64))
                call h5dwrite_f(did, H5T_NATIVE_DOUBLE, value, dims, ier, mem_sid, sid)
            type is (real(real32))
                call h5dwrite_f(did, H5T_NATIVE_REAL, value, dims, ier, mem_sid, sid)
            type is (integer(int32))
                call h5dwrite_f(did, H5T_NATIVE_INTEGER, value, dims, ier, mem_sid, sid)
            type is (integer(int64))
                call h5dwrite_f(did, H5T_STD_I64LE, value, dims, ier, mem_sid, sid)
        end select
        if (ier/=0) error stop 'h5fortran:ERROR: could not write ' // dname // ' to ' // self%filename

        call hdf_wrapup(did, sid, ier)

        if (present(ierr)) ierr = ier
        if (check(ier, self%filename, dname) .and. .not.present(ierr)) error stop
    end procedure hdf_write_5d

    module procedure hdf_write_6d
        integer(HID_T) :: did, sid, mem_sid
        integer(HSIZE_T) :: dims(rank(value))
        integer :: ier

        dims = shape(value)
        select type (value)
            type is (real(real64))
                call hdf_create(self,dname, H5T_NATIVE_DOUBLE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (real(real32))
                call hdf_create(self,dname, H5T_NATIVE_REAL, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int32))
                call hdf_create(self,dname, H5T_NATIVE_INTEGER, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int64))
                call hdf_create(self,dname, H5T_STD_I64LE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            class default
                error stop 'h5fortran:write:invalid data type'
        end select

        mem_sid = H5S_ALL_F !< default

        if(present(istart) .and. present(iend)) then
            if(present(stride)) then
                !! necessary to use this present check for Intel and GCC
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend, stride)
            else
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend)
            endif
        endif

        select type (value)
            type is (real(real64))
                call h5dwrite_f(did, H5T_NATIVE_DOUBLE, value, dims, ier, mem_sid, sid)
            type is (real(real32))
                call h5dwrite_f(did, H5T_NATIVE_REAL, value, dims, ier, mem_sid, sid)
            type is (integer(int32))
                call h5dwrite_f(did, H5T_NATIVE_INTEGER, value, dims, ier, mem_sid, sid)
            type is (integer(int64))
                call h5dwrite_f(did, H5T_STD_I64LE, value, dims, ier, mem_sid, sid)
        end select
        if (ier/=0) error stop 'h5fortran:ERROR: could not write ' // dname // ' to ' // self%filename

        call hdf_wrapup(did, sid, ier)

        if (present(ierr)) ierr = ier
        if (check(ier, self%filename, dname) .and. .not.present(ierr)) error stop
    end procedure hdf_write_6d

    module procedure hdf_write_7d
        integer(HID_T) :: did, sid, mem_sid
        integer(HSIZE_T) :: dims(rank(value))
        integer :: ier

        dims = shape(value)
        select type (value)
            type is (real(real64))
                call hdf_create(self,dname, H5T_NATIVE_DOUBLE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (real(real32))
                call hdf_create(self,dname, H5T_NATIVE_REAL, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int32))
                call hdf_create(self,dname, H5T_NATIVE_INTEGER, dims, sid, did, chunk_size, istart, iend, stride, compact)
            type is (integer(int64))
                call hdf_create(self,dname, H5T_STD_I64LE, dims, sid, did, chunk_size, istart, iend, stride, compact)
            class default
                error stop 'h5fortran:write:invalid data type'
        end select

        mem_sid = H5S_ALL_F !< default

        if(present(istart) .and. present(iend)) then
            if(present(stride)) then
                !! necessary to use this present check for Intel and GCC
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend, stride)
            else
                call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend)
            endif
        endif

        select type (value)
            type is (real(real64))
                call h5dwrite_f(did, H5T_NATIVE_DOUBLE, value, dims, ier, mem_sid, sid)
            type is (real(real32))
                call h5dwrite_f(did, H5T_NATIVE_REAL, value, dims, ier, mem_sid, sid)
            type is (integer(int32))
                call h5dwrite_f(did, H5T_NATIVE_INTEGER, value, dims, ier, mem_sid, sid)
            type is (integer(int64))
                call h5dwrite_f(did, H5T_STD_I64LE, value, dims, ier, mem_sid, sid)
        end select
        if (ier/=0) error stop 'h5fortran:ERROR: could not write ' // dname // ' to ' // self%filename

        call hdf_wrapup(did, sid, ier)

        if (present(ierr)) ierr = ier
        if (check(ier, self%filename, dname) .and. .not.present(ierr)) error stop
    end procedure hdf_write_7d

end submodule writer


text:
      subroutine fpfrno(maxtr,up,left,right,info,point,merk,n1,
     * count,ier)
c  subroutine fpfrno collects the free nodes (up field zero) of the
c  triply linked tree the information of which is kept in the arrays
c  up,left,right and info. the maximal length of the branches of the
c  tree is given by n1. if no free nodes are found, the error flag
c  ier is set to 1.
c  ..
c  ..scalar arguments..
      integer maxtr,point,merk,n1,count,ier
c  ..array arguments..
      integer up(maxtr),left(maxtr),right(maxtr),info(maxtr)
c  ..local scalars
      integer i,j,k,l,n,niveau
c  ..
      ier = 1
      if(n1.eq.2) go to 140
      niveau = 1
      count = 2
  10  j = 0
      i = 1
  20  if(j.eq.niveau) go to 30
      k = 0
      l = left(i)
      if(l.eq.0) go to 110
      i = l
      j = j+1
      go to 20
  30  if (i.lt.count) go to 110
      if (i.eq.count) go to 100
      go to 40
  40  if(up(count).eq.0) go to 50
      count = count+1
      go to 30
  50  up(count) = up(i)
      left(count) = left(i)
      right(count) = right(i)
      info(count) = info(i)
      if(merk.eq.i) merk = count
      if(point.eq.i) point = count
      if(k.eq.0) go to 60
      right(k) = count
      go to 70
  60  n = up(i)
      left(n) = count
  70  l = left(i)
  80  if(l.eq.0) go to 90
      up(l) = count
      l = right(l)
      go to 80
  90  up(i) = 0
      i = count
 100  count = count+1
 110  l = right(i)
      k = i
      if(l.eq.0) go to 120
      i = l
      go to 20
 120  l = up(i)
      j = j-1
      if(j.eq.0) go to 130
      i = l
      go to 110
 130  niveau = niveau+1
      if(niveau.le.n1) go to 10
      if(count.gt.maxtr) go to 140
      ier = 0
 140  return
      end

max_stars_repo_name: Ennosigaeon/scipy

max_issues_repo_name: Ennosigaeon/scipy

max_forks_repo_name: Ennosigaeon/scipy


text:
      subroutine gg_h(p,msq)
      implicit none
c----Lowest order matrix element for H production
c----in the heavy quark (mt=Infinity) limit.
C----averaged over initial colours and spins
c     g(-p1)+g(-p2)-->H -->  b(p3)+b(p4))
      include 'constants.f'
      include 'masses.f'
      include 'qcdcouple.f'
      include 'ewcouple.f'
      include 'hdecaymode.f'
      integer j,k
      double precision msq(-nf:nf,-nf:nf),p(mxpart,4),s,s12
      double precision hdecay,gg,Asq,hprod,msqgamgam
c      double precision tn,ftn
      s(j,k)=2d0*(p(j,4)*p(k,4)-p(j,1)*p(k,1)
     .           -p(j,2)*p(k,2)-p(j,3)*p(k,3))

c--- set msq=0 to initialize
      do j=-nf,nf
      do k=-nf,nf
      msq(j,k)=0d0
      enddo
      enddo
      s12=s(1,2)

C   Deal with Higgs decay
      if (hdecaymode == 'tlta') then
          call htautaudecay(p,3,4,hdecay)
      elseif (hdecaymode == 'bqba') then
          call hbbdecay(p,3,4,hdecay)
      elseif (hdecaymode == 'gaga') then
          hdecay=msqgamgam(hmass)
      else
      write(6,*) 'Unimplemented process in gg_hgg_v'
      stop
      endif
      hdecay=hdecay/((s12-hmass**2)**2+(hmass*hwidth)**2)

      hprod=1d0
C--- The following lines may be uncommented in order to implement
c--- the correct dependence on the top quark mass
c      tn = 4d0*(mt/hmass)**2
c      if (tn.lt.(1.0d0)) then
c         ftn = 0.5d0*((dlog((1d0+dsqrt(1d0-tn))
c     &        /(1d0-dsqrt(1d0-tn))))**2-pisq)
c      else
c         ftn = -2d0*(dasin(1.0d0/dsqrt(tn)))**2
c      endif
c      hprod=((3d0*tn/4d0)*(2d0+(tn-1d0)*ftn))**2
C---- end of mass corrections in production
c      write(6,*) 'hprod',hprod

      Asq=(as/(3d0*pi))**2/vevsq
      gg=0.5d0*Asq*V*s12**2
      msq(0,0)=avegg*gg*hdecay*hprod

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
!##############################################################################
!# Tutorial 004f: Datastructures - comparison of linked lists vs. maps
!##############################################################################

module tutorial004f

  ! Include basic Feat-2 modules
  use fsystem
  use genoutput
  use listInt
  use mapInt
  use random
  use statistics

  implicit none
  private

  public :: start_tutorial004f

contains

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

  subroutine start_tutorial004f

    ! Declare some variables
    type(t_random) :: rrandom
    type(t_timer) :: rtimer
    type(t_listInt) :: rlist
    type(it_listInt) :: rlistIterator
    type(t_mapInt) :: rmap
    type(it_mapInt) :: rmapIterator
    integer(I32), dimension(:), allocatable :: Irandom
    integer :: i

    ! Print a message
    call output_lbrk()
    call output_separator (OU_SEP_STAR)
    call output_line ("This is FEAT-2. Tutorial 004f")
    call output_separator (OU_SEP_MINUS)

    ! =================================
    ! Initialise random number generator
    ! =================================

    call rng_initByClock(rrandom)

    ! =================================
    ! Generate 20000 random integer values
    ! =================================

    allocate(Irandom(20000))
    do i=1,size(Irandom)
      call rng_get(rrandom,Irandom(i))
    end do

    ! =================================
    ! Start time measurement for list
    ! =================================

    call stat_clearTimer(rtimer)
    call stat_startTimer(rtimer)

    ! =================================
    ! Create list
    ! =================================

    call list_create(rlist, size(Irandom))

    ! =================================
    ! Insert items at the end
    ! =================================

    do i=1,size(Irandom)
#ifdef USE_LARGEINT
      call list_push_back(rlist,int(Irandom(i), I64))
#else
      call list_push_back(rlist,Irandom(i))
#endif
    end do

    ! =================================
    ! Search for all items in list
    ! =================================

    do i=1,size(Irandom)
#ifdef USE_LARGEINT
      rlistIterator = list_find(rlist, int(Irandom(i), I64))
#else
      rlistIterator = list_find(rlist, Irandom(i))
#endif
    end do

    ! =================================
    ! Release list
    ! =================================

    call list_release(rlist)

    ! =================================
    ! Stop time measurement
    ! =================================

    call stat_stopTimer(rtimer)
    call output_line ("Total time for linked list: "&
        // trim(adjustl(sys_sd(rtimer%delapsedReal,10)))&
        // " seconds")

    ! =================================
    ! Start time measurement for map
    ! =================================

    call stat_clearTimer(rtimer)
    call stat_startTimer(rtimer)

    ! =================================
    ! Create map
    ! =================================

    call map_create(rmap, size(Irandom))

    ! =================================
    ! Insert items
    ! =================================

    do i=1,size(Irandom)
#ifdef USE_LARGEINT
      rmapIterator = map_insert(rmap,int(Irandom(i), I64))
#else
      rmapIterator = map_insert(rmap,Irandom(i))
#endif
    end do

    ! =================================
    ! Search for all items in map
    ! =================================

    do i=1,size(Irandom)
#ifdef USE_LARGEINT
      rmapIterator = map_find(rmap, int(Irandom(i), I64))
#else
      rmapIterator = map_find(rmap, Irandom(i))
#endif
    end do

    ! =================================
    ! Release map
    ! =================================

    call map_release(rmap)

    ! =================================
    ! Stop time measurement
    ! =================================

    call stat_stopTimer(rtimer)
    call output_line ("Total time for map        : "&
        // trim(adjustl(sys_sd(rtimer%delapsedReal,10)))&
        // " seconds")

    ! =================================
    ! Free memory
    ! =================================

    deallocate(Irandom)

  end subroutine

end module