text:
#include "../defines.inc"

subroutine mc_reptation(IT1, IT2, IB1, IB2, rand_stat, forward, super)
!--------------------------------------------------------------*
!
!           Makes Monti Carlo reptation move
!
!    Quinn separated out this file on 8/9/17
!
!---------------------------------------------------------------

! values from wlcsim_data
   use params, only: wlc_VP, wlc_ABP, wlc_V, wlc_AB, wlc_R &
                     , wlc_RP, wlc_U, wlc_UP, wlc_nPointsMoved, wlc_pointsMoved
   use mersenne_twister
   use params, only: dp
   use vector_utils, only: random_perp, cross
   use polydispersity, only: first_bead_of_chain, last_bead_of_chain, get_IB, length_of_chain

   implicit none
   logical, intent(in) :: super
   integer, intent(out) :: IT1   ! Index of test bead 1
   integer, intent(out) :: IT2   ! Index of test bead 2
   integer, intent(out) :: IB1   ! Index of test bead 1
   integer, intent(out) :: IB2   ! Index of test bead 2

   integer IP    ! Test polymer
   integer I  ! Test indices
! Things for random number generator
   type(random_stat), intent(inout) :: rand_stat  ! status of random number generator
   real(dp) urnd(1) ! single random number
   integer irnd(1)
   real(dp) DR(3)    ! Displacement for slide move
   real(dp) Uvec(3) ! parallel component of triad
   real(dp) pDir(3) ! perp component of triad
   real(dp) tDir(3) ! twist component of triad
   real(dp) r_relative(3) ! r in new coordinate system
   real(dp) u_relative(3) ! u in new coordinate system
   real(dp) v_relative(3) ! v in new coordinate system
   logical, intent(out) :: forward

!TOdo saving RP is not actually needed, even in these cases, but Brad's code assumes that we have RP.
   if (WLC_P__RING .OR. WLC_P__INTERP_BEAD_LENNARD_JONES) then
      wlc_RP = wlc_R
      wlc_UP = wlc_U
   endif

! single bead reptation
   call random_index(WLC_P__NP, irnd, rand_stat)
   IP = irnd(1)
   IT1 = first_bead_of_chain(IP)
   IT2 = last_bead_of_chain(IP)
   IB1 = 1
   IB2 = length_of_chain(IP)
! move forward or backward
   call random_number(urnd, rand_stat)
   if (urnd(1) .lt. 0.5_dp) then
      forward = .true.
      dR = wlc_R(:, IT1 + 1) - wlc_R(:, IT1)

      Uvec = wlc_U(:, IT1)
      ! chose coordinate system
      if (WLC_P__LOCAL_TWIST) then
         pDir = wlc_V(:, IT1)
         tDir = cross(Uvec, pDir)
      else
         call random_perp(Uvec, pDir, tDir, rand_stat)
      endif
      ! find next r and u in new coordinate system
      u_relative(1) = dot_product(Uvec, wlc_U(:, IT1 + 1))
      u_relative(2) = dot_product(pDir, wlc_U(:, IT1 + 1))
      u_relative(3) = dot_product(tDir, wlc_U(:, IT1 + 1))
      if (WLC_P__LOCAL_TWIST) then
         v_relative(1) = dot_product(Uvec, wlc_V(:, IT1 + 1))
         v_relative(2) = dot_product(pDir, wlc_V(:, IT1 + 1))
         v_relative(3) = dot_product(tDir, wlc_V(:, IT1 + 1))
      endif
      r_relative(1) = dot_product(Uvec, dR)
      r_relative(2) = dot_product(pDir, dR)
      r_relative(3) = dot_product(tDir, dR)

      ! orient coordinate system with end of chain
      Uvec = wlc_U(:, IT2)
      if (WLC_P__LOCAL_TWIST) then
         pDir = wlc_V(:, IT2)
         tDir = cross(Uvec, pDir)
      else
         call random_perp(Uvec, pDir, tDir, rand_stat)
      endif
      ! update UP and RP
      wlc_UP(:, IT2) = Uvec*u_relative(1) + pDir*u_relative(2) + tDir*u_relative(3)
      wlc_UP(:, IT2) = wlc_UP(:, IT2)/norm2(wlc_UP(:, IT2))
      if (WLC_P__LOCAL_TWIST) then
         wlc_VP(:, IT2) = Uvec*v_relative(1) + pDir*v_relative(2) + tDir*v_relative(3)
         wlc_VP(:, IT2) = wlc_VP(:, IT2)/norm2(wlc_VP(:, IT2))
      endif
      wlc_RP(:, IT2) = wlc_R(:, IT2) + Uvec*r_relative(1) + pDir*r_relative(2) + tDir*r_relative(3)

      do I = IT1, IT2 - 1
         wlc_RP(:, I) = wlc_R(:, I + 1)
         wlc_UP(:, I) = wlc_U(:, I + 1)
         if (WLC_P__LOCAL_TWIST) wlc_VP(:, I) = wlc_V(:, I + 1)
      enddo
      if (super) then
         do I = IT1, IT2 - 1
            wlc_ABP(I) = wlc_AB(I + 1)
         enddo
         wlc_ABP(IT2) = wlc_AB(IT1)
      endif

      ! RperpMag = sqrt(r_relative(2)**2 + r_relative(3)**2)
      ! RparaMag = r_relative(1)
      ! call test_equiv_forward(U,R,UP,RP,WLC_P__NT,IT1,IT2,RparaMag,RperpMag)

   else
      forward = .false.
      dR = wlc_R(:, IT2) - wlc_R(:, IT2 - 1)

      Uvec = wlc_U(:, IT2)
      ! chose coordinate system
      if (WLC_P__LOCAL_TWIST) then
         pDir = wlc_V(:, IT2)
         tDir = cross(Uvec, pDir)
      else
         call random_perp(Uvec, pDir, tDir, rand_stat)
      endif
      ! find next r and u in new coordinate system
      u_relative(1) = dot_product(Uvec, wlc_U(:, IT2 - 1))
      u_relative(2) = dot_product(pDir, wlc_U(:, IT2 - 1))
      u_relative(3) = dot_product(tDir, wlc_U(:, IT2 - 1))
      if (WLC_P__LOCAL_TWIST) then
         v_relative(1) = dot_product(Uvec, wlc_V(:, IT2 - 1))
         v_relative(2) = dot_product(pDir, wlc_V(:, IT2 - 1))
         v_relative(3) = dot_product(tDir, wlc_V(:, IT2 - 1))
      endif
      r_relative(1) = dot_product(Uvec, dR)
      r_relative(2) = dot_product(pDir, dR)
      r_relative(3) = dot_product(tDir, dR)

      ! orient coordinate system with end of chain
      Uvec = wlc_U(:, IT1)
      if (WLC_P__LOCAL_TWIST) then
         pDir = wlc_V(:, IT2)
         tDir = cross(Uvec, pDir)
      else
         call random_perp(Uvec, pDir, tDir, rand_stat)
      endif
      ! update UP and RP
      wlc_UP(:, IT1) = Uvec*u_relative(1) + pDir*u_relative(2) + tDir*u_relative(3)
      wlc_UP(:, IT1) = wlc_UP(:, IT1)/norm2(wlc_UP(:, IT1))
      if (WLC_P__LOCAL_TWIST) then
         wlc_VP(:, IT1) = Uvec*v_relative(1) + pDir*v_relative(2) + tDir*v_relative(3)
         wlc_VP(:, IT1) = wlc_VP(:, IT1)/norm2(wlc_VP(:, IT1))
      endif
      wlc_RP(:, IT1) = wlc_R(:, IT1) - Uvec(:)*r_relative(1) - pDir(:)*r_relative(2) - tDir(:)*r_relative(3)

      do I = IT1 + 1, IT2
         wlc_RP(:, I) = wlc_R(:, I - 1)
         wlc_UP(:, I) = wlc_U(:, I - 1)
         if (WLC_P__LOCAL_TWIST) wlc_VP(:, I) = wlc_V(:, I - 1)
      enddo
      if (super) then
         do I = IT1 + 1, IT2
            wlc_ABP(I) = wlc_AB(I - 1)
         enddo
         wlc_ABP(IT1) = wlc_AB(IT2)
      endif
   endif
   do I = IT1, IT2
      wlc_nPointsMoved = wlc_nPointsMoved + 1
      wlc_pointsMoved(wlc_nPointsMoved) = I
   enddo
end subroutine

max_stars_repo_name: SpakowitzLab/BasicWLC

max_issues_repo_name: riscalab/wlcsim

max_forks_repo_name: riscalab/wlcsim


subroutine general_read_nmmb(grd,filename,mype,g_z,g_ps,g_u,g_v,g_tv,g_tsen,g_q,g_oz)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    general_read_nmmb  adaptation of read_nems_nmmb_guess for general resolutions
!   prgmmr: parrish          org: np22                date: 2003-09-05
!
! abstract: copied from read_nems_nmmb_guess, primarily for reading in NMMB NEMSIO files, for now the
!            input resolution and the grid that variables are reconstructed on MUST be the same
!            as the analysis grid/resolution.
!
! program history log:
!   2011-07-01  carley     - Initial adaptation
!   2015-05-12  wu         - changes to read in multiple guess files for FGAT/4DEnVar
!
!   input argument list:
!     grd      - structure variable containing information about grid
!                    (initialized by general_sub2grid_create_info, located in general_sub2grid_mod.f90)
!     filename - input nemsio file name
!     mype     - mpi task id
!
!   output argument list:
!     g_*      - guess fields
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind,r_single
  use gridmod, only: pdtop_ll,pt_ll,nmmb_verttype,use_gfs_ozone,regional_ozone
  use constants, only: zero,one_tenth,one,fv,r0_01
  use gsi_nemsio_mod, only: gsi_nemsio_open,gsi_nemsio_close,gsi_nemsio_read
  use general_sub2grid_mod, only: sub2grid_info

  implicit none

!   Declare passed variables
    type(sub2grid_info)                   ,intent(in   ) :: grd
    character(70)                         ,intent(in   ) :: filename
    integer(i_kind)                       ,intent(in   ) :: mype
    real(r_kind),dimension(grd%lat2,grd%lon2)     ,intent(  out) :: g_z,g_ps
    real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent(  out) :: g_u,g_v,&
         g_tv,g_tsen,g_q,g_oz

!   Declare local variables
    integer(i_kind) i,j,k,kr,nsig,lat2,lon2,mype_input,ierr
    real(r_kind) pd,psfc_this,pd_to_ps,dumtv
    logical good_o3mr
    real(r_kind),dimension(grd%lat2,grd%lon2) :: g_pd

!     get conversion factor for pd to psfc

  if(nmmb_verttype=='OLD') then
     pd_to_ps=pdtop_ll+pt_ll
  else
     pd_to_ps=pt_ll
  end if

   lat2=grd%lat2
   lon2=grd%lon2
   nsig=grd%nsig

   mype_input=0

     call gsi_nemsio_open(filename,'READ', &
                          'GENERAL_READ_NMMB:  problem with ens input file!',mype,mype_input,ierr)

!                            ! pd

     call gsi_nemsio_read('dpres','hybrid sig lev','H',1,g_pd(:,:),mype,mype_input)
     do i=1,lon2
        do j=1,lat2
!               convert wrf nmm pd variable to psfc in mb, and then to log(psfc) in cb
           pd=r0_01*g_pd(j,i)
           psfc_this=pd+pd_to_ps
           g_ps(j,i)=one_tenth*psfc_this
        end do
     end do

!                          !   fis

     call gsi_nemsio_read('hgt','sfc','H',1,g_z(:,:),mype,mype_input)

!                          !   u,v,q,tsen,tv,cwmr,qr,qli
     do kr=1,nsig
        k=nsig+1-kr
        call gsi_nemsio_read('ugrd','mid layer','V',kr,g_u(:,:,k),   mype,mype_input)
        call gsi_nemsio_read('vgrd','mid layer','V',kr,g_v(:,:,k),   mype,mype_input)
        call gsi_nemsio_read('spfh','mid layer','H',kr,g_q(:,:,k),   mype,mype_input)
        call gsi_nemsio_read('tmp' ,'mid layer','H',kr,g_tsen(:,:,k),mype,mype_input)

	do i=1,lon2
           do j=1,lat2
              g_tv(j,i,k) = g_tsen(j,i,k) * (one+fv*g_q(j,i,k))  !To be consistent with read_nems_nmmb_guess
	                                                         ! compute tv prior to enforcing the limit on q.
	      g_q(j,i,k)=max(g_q(j,i,k),1.e-10_r_kind)           !step comes from compute_qvar3d.f90
              dumtv=g_tsen(j,i,k) * (one+fv*g_q(j,i,k))
	      g_tsen(j,i,k)=dumtv/(one+fv*max(zero,g_q(j,i,k)))     !Recompute tsen based on the limit enforced on q.
            end do
        end do

        if(regional_ozone) then
           if(use_gfs_ozone) then
              g_oz(:,:,k)=zero
           else
              good_o3mr=.false.
              call gsi_nemsio_read('o3mr' ,'mid layer','H',kr,g_oz(:,:,k),mype,mype_input,good_o3mr)
              if(.not.good_o3mr) write(6,*)' IN GENERAL_READ_NMMB, O3MR FIELD NOT YET AVAILABLE'
           end if
        end if
     end do

     call gsi_nemsio_close(filename,'GENERAL_READ_NMMB',mype,mype_input)
  return
end subroutine general_read_nmmb


text:
C
C	$Id: skpblk.f,v 1.4 2008-07-27 12:23:44 haley Exp $
C
C			     Copyright (C)  1997
C	     University Corporation for Atmospheric Research
C			     All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C

      SUBROUTINE SKPBLK(IOS, STATUS)
C
C  Position the buffer pointer to point at the next non-blank
C  character in the file.
C
C
C  OUTPUT
C       IOS     --  I/O status flag.  This flag is valid only
C                   if STATUS indicates an error.
C       STATUS  --  The error status as defined in COMMON CAPERR.
C
C
      COMMON /CAPERR/ ALLOK, EOFFL, INTERR, MAXINT, PLBERR, PMBERR,
     1                FABERR, TXBERR, FLTERR, MAXFLT, NOTINT, SIZERR,
     2                UNDASC, DEVERR, DOCERR, TBLERR , STSZER, ENTERR,
     3                TABERR, TABSER, PRSFIL
      INTEGER ALLOK, EOFFL, INTERR, MAXINT, PLBERR, PMBERR,
     1        FABERR, TXBERR, FLTERR, MAXFLT, NOTINT, SIZERR, UNDASC,
     2        DEVERR, DOCERR, STSZER, ENTERR, TABERR, TABSER, TBLERR,
     3        PRSFIL
      COMMON /CAPIOB/ UNIT, IPTR, LSIZE, NFST, NFLG
      COMMON /CAPIO2/ LINE, GRNM
      INTEGER LNMAX
      PARAMETER (LNMAX=80)
      CHARACTER*80 LINE
      CHARACTER*80 GRNM
      INTEGER UNIT, IPTR, LSIZE, NFST, NFLG
C
      INTEGER IOS, STATUS
      CHARACTER*1 BLANK
      INTEGER II
C
      DATA BLANK/' '/
C
C  Scan the current buffer for a non-blank.
C
 10   CONTINUE
      DO 20 II = IPTR,LSIZE
        IF (LINE(II:II) .NE. BLANK) THEN
C
C  Non-blank found.
C
                IPTR = II
                RETURN
        END IF
 20   CONTINUE
C
C  End of current buffer and no non-blank, get another line.
C
      CALL RDLINE (IOS, STATUS)
C
C  Return if EOF or bad read.
C
      IF (STATUS.NE.ALLOK) RETURN
C
C  Go back and scan this line.
C
      GO TO 10
C
      END

ext: f

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


! Code automatically generated by 'PreSARAH'
! Expressions for amplitudes are obtained by FeynArts/FormCalc
! Created at 14:09 on 7.11.2016

Module DecayFFV

Use Control
Use LoopFunctions
Use AddLoopFunctions
Use Settings

Contains

Subroutine TreeAmp_FtoFV(Mex1,Mex2,Mex3,coupT1L, coupT1R,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude(4)
Complex(dp), Intent(in) :: coupT1L, coupT1R

Amplitude = 0._dp

! Now the full amplitude for the different Lorentz Parts
! Coefficient of Mat[<v1|6,ec[3]|v2>]
Amplitude(1)= 1*(0.,1._dp)*coupT1R

! Coefficient of Mat[<v1|7,ec[3]|v2>]
Amplitude(2)= 1*(0.,1._dp)*coupT1L

End Subroutine TreeAmp_FtoFV

! --------------------------------------
! Topology Nr.1
! --------------------------------------

Subroutine Amp_VERTEX_FtoFV_Topology1_SFF(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1L, coup1R, coup2L, coup2R, coup3L, coup3R,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude(4)
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1L, coup1R, coup2L, coup2R, coup3L, coup3R
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

Amplitude = 0._dp

! Necessary Loop Integrals
int1=SA_B0(Mex3**2,ML2**2,ML3**2)
int2=SA_Cget("C0  ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int3=SA_Cget("C00 ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int4=SA_Cget("C1  ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int5=SA_Cget("C12 ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int6=SA_Cget("C2  ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int7=SA_Cget("C22 ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of Mat[<v1|6,ec[3]|v2>]
Amplitude(1)= &
& +(-1*(0.,1._dp)*coup1R*coup2L*coup3R*int1) &
& +(2*(0.,1._dp)*coup1R*coup2L*coup3R*int3) &
& +(1*(0.,1._dp)*coup1R*coup2L*coup3R*int2*Mex1**2) &
& +(1*(0.,1._dp)*coup1R*coup2L*coup3R*int4*Mex1**2) &
& +(1*(0.,1._dp)*coup1R*coup2L*coup3R*int6*Mex1**2) &
& +(-1*(0.,1._dp)*coup1L*coup2R*coup3L*int6*Mex1*Mex2) &
& +(-1*(0.,1._dp)*coup1R*coup2L*coup3R*int4*Mex2**2) &
& +(-1*(0.,1._dp)*coup1R*coup2L*coup3R*int2*ML1**2) &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3R*int2*Mex1*ML2) &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3R*int4*Mex1*ML2) &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3R*int6*Mex1*ML2) &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3L*int2*Mex2*ML2) &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3L*int4*Mex2*ML2) &
& +(-1*(0.,1._dp)*coup1L*coup2L*coup3L*int4*Mex1*ML3) &
& +(-1*(0.,1._dp)*coup1L*coup2L*coup3L*int6*Mex1*ML3) &
& +(-1*(0.,1._dp)*coup1R*coup2R*coup3R*int4*Mex2*ML3) &
& +(1*(0.,1._dp)*coup1R*coup2L*coup3L*int2*ML2*ML3)

! Coefficient of Mat[<v1|7,ec[3]|v2>]
Amplitude(2)= &
& +(-1*(0.,1._dp)*coup1L*coup2R*coup3L*int1) &
& +(2*(0.,1._dp)*coup1L*coup2R*coup3L*int3) &
& +(1*(0.,1._dp)*coup1L*coup2R*coup3L*int2*Mex1**2) &
& +(1*(0.,1._dp)*coup1L*coup2R*coup3L*int4*Mex1**2) &
& +(1*(0.,1._dp)*coup1L*coup2R*coup3L*int6*Mex1**2) &
& +(-1*(0.,1._dp)*coup1R*coup2L*coup3R*int6*Mex1*Mex2) &
& +(-1*(0.,1._dp)*coup1L*coup2R*coup3L*int4*Mex2**2) &
& +(-1*(0.,1._dp)*coup1L*coup2R*coup3L*int2*ML1**2) &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3L*int2*Mex1*ML2) &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3L*int4*Mex1*ML2) &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3L*int6*Mex1*ML2) &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3R*int2*Mex2*ML2) &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3R*int4*Mex2*ML2) &
& +(-1*(0.,1._dp)*coup1R*coup2R*coup3R*int4*Mex1*ML3) &
& +(-1*(0.,1._dp)*coup1R*coup2R*coup3R*int6*Mex1*ML3) &
& +(-1*(0.,1._dp)*coup1L*coup2L*coup3L*int4*Mex2*ML3) &
& +(1*(0.,1._dp)*coup1L*coup2R*coup3R*int2*ML2*ML3)

! Coefficient of Mat[<v1|6|v2>] Pair[ec[3], k[1]]
Amplitude(3)= &
& +(-2*(0.,1._dp)*coup1L*coup2R*coup3L*int5*Mex1) &
& +(-2*(0.,1._dp)*coup1L*coup2R*coup3L*int6*Mex1) &
& +(-2*(0.,1._dp)*coup1L*coup2R*coup3L*int7*Mex1) &
& +(2*(0.,1._dp)*coup1R*coup2L*coup3R*int5*Mex2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3L*int2*ML2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3L*int4*ML2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3L*int6*ML2) &
& +(-2*(0.,1._dp)*coup1R*coup2R*coup3R*int4*ML3)

! Coefficient of Mat[<v1|7|v2>] Pair[ec[3], k[1]]
Amplitude(4)= &
& +(-2*(0.,1._dp)*coup1R*coup2L*coup3R*int5*Mex1) &
& +(-2*(0.,1._dp)*coup1R*coup2L*coup3R*int6*Mex1) &
& +(-2*(0.,1._dp)*coup1R*coup2L*coup3R*int7*Mex1) &
& +(2*(0.,1._dp)*coup1L*coup2R*coup3L*int5*Mex2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3R*int2*ML2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3R*int4*ML2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3R*int6*ML2) &
& +(-2*(0.,1._dp)*coup1L*coup2L*coup3L*int4*ML3)

If ((Abs(Amplitude(1))).ne.(Abs(Amplitude(1)))) Then
Amplitude=0._dp
End if

!Field Order
! FieldOrder[1,{S, F, F},1]={External[1][gt1], AntiField[Internal[2]][i2], Internal[1][i1]};
! FieldOrder[1,{S, F, F},2]={Internal[3][i3], External[2][gt2], AntiField[Internal[1]][i1]};
! FieldOrder[1,{S, F, F},3]={Internal[2][i2], AntiField[Internal[3]][i3], External[3][gt3]};
End Subroutine Amp_VERTEX_FtoFV_Topology1_SFF

Subroutine Amp_VERTEX_FtoFV_Topology1_FSS(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1L, coup1R, coup2L, coup2R, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude(4)
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1L, coup1R, coup2L, coup2R, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

Amplitude = 0._dp

! Necessary Loop Integrals
int1=SA_Cget("C0  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int2=SA_Cget("C00 ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int3=SA_Cget("C1  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int4=SA_Cget("C11 ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int5=SA_Cget("C12 ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int6=SA_Cget("C2  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int7=SA_Cget("C22 ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of Mat[<v1|6,ec[3]|v2>]
Amplitude(1)= 2*(0.,1._dp)*coup1R*coup2L*coup3*int2

! Coefficient of Mat[<v1|7,ec[3]|v2>]
Amplitude(2)= 2*(0.,1._dp)*coup1L*coup2R*coup3*int2

! Coefficient of Mat[<v1|6|v2>] Pair[ec[3], k[1]]
Amplitude(3)= &
& +(-2*(0.,1._dp)*coup1L*coup2R*coup3*int5*Mex1) &
& +(-2*(0.,1._dp)*coup1L*coup2R*coup3*int6*Mex1) &
& +(-2*(0.,1._dp)*coup1L*coup2R*coup3*int7*Mex1) &
& +(-2*(0.,1._dp)*coup1R*coup2L*coup3*int3*Mex2) &
& +(-2*(0.,1._dp)*coup1R*coup2L*coup3*int4*Mex2) &
& +(-2*(0.,1._dp)*coup1R*coup2L*coup3*int5*Mex2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3*int1*ML1) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3*int3*ML1) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3*int6*ML1)

! Coefficient of Mat[<v1|7|v2>] Pair[ec[3], k[1]]
Amplitude(4)= &
& +(-2*(0.,1._dp)*coup1R*coup2L*coup3*int5*Mex1) &
& +(-2*(0.,1._dp)*coup1R*coup2L*coup3*int6*Mex1) &
& +(-2*(0.,1._dp)*coup1R*coup2L*coup3*int7*Mex1) &
& +(-2*(0.,1._dp)*coup1L*coup2R*coup3*int3*Mex2) &
& +(-2*(0.,1._dp)*coup1L*coup2R*coup3*int4*Mex2) &
& +(-2*(0.,1._dp)*coup1L*coup2R*coup3*int5*Mex2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3*int1*ML1) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3*int3*ML1) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3*int6*ML1)

If ((Abs(Amplitude(1))).ne.(Abs(Amplitude(1)))) Then
Amplitude=0._dp
End if

!Field Order
! FieldOrder[1,{F, S, S},1]={External[1][gt1], Internal[1][i1], AntiField[Internal[2]][i2]};
! FieldOrder[1,{F, S, S},2]={AntiField[Internal[1]][i1], External[2][gt2], Internal[3][i3]};
! FieldOrder[1,{F, S, S},3]={Internal[2][i2], AntiField[Internal[3]][i3], External[3][gt3]};
End Subroutine Amp_VERTEX_FtoFV_Topology1_FSS

Subroutine Amp_VERTEX_FtoFV_Topology1_VFF(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1L, coup1R, coup2L, coup2R, coup3L, coup3R,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude(4)
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1L, coup1R, coup2L, coup2R, coup3L, coup3R
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

Amplitude = 0._dp

! Necessary Loop Integrals
int1=SA_B0(Mex3**2,ML2**2,ML3**2)
int2=SA_Cget("C0  ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int3=SA_Cget("C00 ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int4=SA_Cget("C1  ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int5=SA_Cget("C12 ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int6=SA_Cget("C2  ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)
int7=SA_Cget("C22 ",Mex3**2,Mex2**2,Mex1**2,ML2**2,ML3**2,ML1**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of Mat[<v1|6,ec[3]|v2>]
Amplitude(1)= &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3L*rMS) &
& +(-2*(0.,1._dp)*coup1L*coup2L*coup3L*int1) &
& +(4*(0.,1._dp)*coup1L*coup2L*coup3L*int3) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3L*int2*Mex1**2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3L*int4*Mex1**2) &
& +(4*(0.,1._dp)*coup1L*coup2L*coup3L*int6*Mex1**2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3R*int6*Mex1*Mex2) &
& +(-2*(0.,1._dp)*coup1L*coup2L*coup3L*int4*Mex2**2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3L*int6*Mex2**2) &
& +(-2*(0.,1._dp)*coup1L*coup2L*coup3L*int6*Mex3**2) &
& +(-2*(0.,1._dp)*coup1L*coup2L*coup3L*int2*ML1**2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3R*int2*ML2*ML3)

! Coefficient of Mat[<v1|7,ec[3]|v2>]
Amplitude(2)= &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3R*rMS) &
& +(-2*(0.,1._dp)*coup1R*coup2R*coup3R*int1) &
& +(4*(0.,1._dp)*coup1R*coup2R*coup3R*int3) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3R*int2*Mex1**2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3R*int4*Mex1**2) &
& +(4*(0.,1._dp)*coup1R*coup2R*coup3R*int6*Mex1**2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3L*int6*Mex1*Mex2) &
& +(-2*(0.,1._dp)*coup1R*coup2R*coup3R*int4*Mex2**2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3R*int6*Mex2**2) &
& +(-2*(0.,1._dp)*coup1R*coup2R*coup3R*int6*Mex3**2) &
& +(-2*(0.,1._dp)*coup1R*coup2R*coup3R*int2*ML1**2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3L*int2*ML2*ML3)

! Coefficient of Mat[<v1|6|v2>] Pair[ec[3], k[1]]
Amplitude(3)= &
& +(-4*(0.,1._dp)*coup1R*coup2R*coup3R*int5*Mex1) &
& +(-4*(0.,1._dp)*coup1R*coup2R*coup3R*int7*Mex1) &
& +(4*(0.,1._dp)*coup1L*coup2L*coup3L*int5*Mex2) &
& +(4*(0.,1._dp)*coup1L*coup2L*coup3L*int6*Mex2) &
& +(-4*(0.,1._dp)*coup1L*coup2R*coup3R*int6*ML2) &
& +(-4*(0.,1._dp)*coup1L*coup2R*coup3L*int6*ML3)

! Coefficient of Mat[<v1|7|v2>] Pair[ec[3], k[1]]
Amplitude(4)= &
& +(-4*(0.,1._dp)*coup1L*coup2L*coup3L*int5*Mex1) &
& +(-4*(0.,1._dp)*coup1L*coup2L*coup3L*int7*Mex1) &
& +(4*(0.,1._dp)*coup1R*coup2R*coup3R*int5*Mex2) &
& +(4*(0.,1._dp)*coup1R*coup2R*coup3R*int6*Mex2) &
& +(-4*(0.,1._dp)*coup1R*coup2L*coup3L*int6*ML2) &
& +(-4*(0.,1._dp)*coup1R*coup2L*coup3R*int6*ML3)

If ((Abs(Amplitude(1))).ne.(Abs(Amplitude(1)))) Then
Amplitude=0._dp
End if

!Field Order
! FieldOrder[1,{V, F, F},1]={External[1][gt1], AntiField[Internal[2]][i2], Internal[1][i1]};
! FieldOrder[1,{V, F, F},2]={Internal[3][i3], External[2][gt2], AntiField[Internal[1]][i1]};
! FieldOrder[1,{V, F, F},3]={Internal[2][i2], AntiField[Internal[3]][i3], External[3][gt3]};
End Subroutine Amp_VERTEX_FtoFV_Topology1_VFF

Subroutine Amp_VERTEX_FtoFV_Topology1_FSV(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1L, coup1R, coup2L, coup2R, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude(4)
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1L, coup1R, coup2L, coup2R, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

Amplitude = 0._dp

! Necessary Loop Integrals
int1=SA_Cget("C0  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int2=SA_Cget("C1  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int3=SA_Cget("C2  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of Mat[<v1|6,ec[3]|v2>]
Amplitude(1)= &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3*int3*Mex1) &
& +(-1*(0.,1._dp)*coup1R*coup2R*coup3*int2*Mex2) &
& +(-1*(0.,1._dp)*coup1R*coup2L*coup3*int1*ML1)

! Coefficient of Mat[<v1|7,ec[3]|v2>]
Amplitude(2)= &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3*int3*Mex1) &
& +(-1*(0.,1._dp)*coup1L*coup2L*coup3*int2*Mex2) &
& +(-1*(0.,1._dp)*coup1L*coup2R*coup3*int1*ML1)

! Coefficient of Mat[<v1|6|v2>] Pair[ec[3], k[1]]
Amplitude(3)= -2*(0.,1._dp)*coup1R*coup2R*coup3*int2

! Coefficient of Mat[<v1|7|v2>] Pair[ec[3], k[1]]
Amplitude(4)= -2*(0.,1._dp)*coup1L*coup2L*coup3*int2

If ((Abs(Amplitude(1))).ne.(Abs(Amplitude(1)))) Then
Amplitude=0._dp
End if
!Field Order
! FieldOrder[1,{F, S, V},1]={External[1][gt1], Internal[1][i1], AntiField[Internal[2]][i2]};
! FieldOrder[1,{F, S, V},2]={AntiField[Internal[1]][i1], External[2][gt2], Internal[3][i3]};
! FieldOrder[1,{F, S, V},3]={Internal[2][i2], External[3][gt3], AntiField[Internal[3]][i3]};
End Subroutine Amp_VERTEX_FtoFV_Topology1_FSV

Subroutine Amp_VERTEX_FtoFV_Topology1_FVS(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1L, coup1R, coup2L, coup2R, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude(4)
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1L, coup1R, coup2L, coup2R, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

Amplitude = 0._dp

! Necessary Loop Integrals
int1=SA_Cget("C0  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int2=SA_Cget("C1  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int3=SA_Cget("C2  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of Mat[<v1|6,ec[3]|v2>]
Amplitude(1)= &
& +(-1*(0.,1._dp)*coup1R*coup2L*coup3*int3*Mex1) &
& +(1*(0.,1._dp)*coup1L*coup2R*coup3*int2*Mex2) &
& +(-1*(0.,1._dp)*coup1L*coup2L*coup3*int1*ML1)

! Coefficient of Mat[<v1|7,ec[3]|v2>]
Amplitude(2)= &
& +(-1*(0.,1._dp)*coup1L*coup2R*coup3*int3*Mex1) &
& +(1*(0.,1._dp)*coup1R*coup2L*coup3*int2*Mex2) &
& +(-1*(0.,1._dp)*coup1R*coup2R*coup3*int1*ML1)

! Coefficient of Mat[<v1|6|v2>] Pair[ec[3], k[1]]
Amplitude(3)= -2*(0.,1._dp)*coup1L*coup2R*coup3*int3

! Coefficient of Mat[<v1|7|v2>] Pair[ec[3], k[1]]
Amplitude(4)= -2*(0.,1._dp)*coup1R*coup2L*coup3*int3

If ((Abs(Amplitude(1))).ne.(Abs(Amplitude(1)))) Then
Amplitude=0._dp
End if
!Field Order
! FieldOrder[1,{F, V, S},1]={External[1][gt1], Internal[1][i1], AntiField[Internal[2]][i2]};
! FieldOrder[1,{F, V, S},2]={AntiField[Internal[1]][i1], External[2][gt2], Internal[3][i3]};
! FieldOrder[1,{F, V, S},3]={AntiField[Internal[3]][i3], External[3][gt3], Internal[2][i2]};
End Subroutine Amp_VERTEX_FtoFV_Topology1_FVS

Subroutine Amp_VERTEX_FtoFV_Topology1_FVV(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1L, coup1R, coup2L, coup2R, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude(4)
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1L, coup1R, coup2L, coup2R, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

Amplitude = 0._dp

! Necessary Loop Integrals
int1=SA_B0(Mex3**2,ML2**2,ML3**2)
int2=SA_Cget("C0  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int3=SA_Cget("C00 ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int4=SA_Cget("C1  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int5=SA_Cget("C11 ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int6=SA_Cget("C12 ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int7=SA_Cget("C2  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int8=SA_Cget("C22 ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of Mat[<v1|6,ec[3]|v2>]
Amplitude(1)= &
& +(-1*(0.,1._dp)*coup1L*coup2L*coup3*rMS) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3*int1) &
& +(4*(0.,1._dp)*coup1L*coup2L*coup3*int3) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3*int4*Mex1**2) &
& +(3*(0.,1._dp)*coup1L*coup2L*coup3*int7*Mex1**2) &
& +(3*(0.,1._dp)*coup1R*coup2R*coup3*int4*Mex1*Mex2) &
& +(3*(0.,1._dp)*coup1R*coup2R*coup3*int7*Mex1*Mex2) &
& +(3*(0.,1._dp)*coup1L*coup2L*coup3*int4*Mex2**2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3*int7*Mex2**2) &
& +(-2*(0.,1._dp)*coup1L*coup2L*coup3*int4*Mex3**2) &
& +(-2*(0.,1._dp)*coup1L*coup2L*coup3*int7*Mex3**2) &
& +(3*(0.,1._dp)*coup1R*coup2L*coup3*int2*Mex1*ML1) &
& +(3*(0.,1._dp)*coup1L*coup2R*coup3*int2*Mex2*ML1) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3*int2*ML1**2)

! Coefficient of Mat[<v1|7,ec[3]|v2>]
Amplitude(2)= &
& +(-1*(0.,1._dp)*coup1R*coup2R*coup3*rMS) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3*int1) &
& +(4*(0.,1._dp)*coup1R*coup2R*coup3*int3) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3*int4*Mex1**2) &
& +(3*(0.,1._dp)*coup1R*coup2R*coup3*int7*Mex1**2) &
& +(3*(0.,1._dp)*coup1L*coup2L*coup3*int4*Mex1*Mex2) &
& +(3*(0.,1._dp)*coup1L*coup2L*coup3*int7*Mex1*Mex2) &
& +(3*(0.,1._dp)*coup1R*coup2R*coup3*int4*Mex2**2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3*int7*Mex2**2) &
& +(-2*(0.,1._dp)*coup1R*coup2R*coup3*int4*Mex3**2) &
& +(-2*(0.,1._dp)*coup1R*coup2R*coup3*int7*Mex3**2) &
& +(3*(0.,1._dp)*coup1L*coup2R*coup3*int2*Mex1*ML1) &
& +(3*(0.,1._dp)*coup1R*coup2L*coup3*int2*Mex2*ML1) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3*int2*ML1**2)

! Coefficient of Mat[<v1|6|v2>] Pair[ec[3], k[1]]
Amplitude(3)= &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3*int4*Mex1) &
& +(-4*(0.,1._dp)*coup1R*coup2R*coup3*int6*Mex1) &
& +(-4*(0.,1._dp)*coup1R*coup2R*coup3*int8*Mex1) &
& +(-4*(0.,1._dp)*coup1L*coup2L*coup3*int5*Mex2) &
& +(-4*(0.,1._dp)*coup1L*coup2L*coup3*int6*Mex2) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3*int7*Mex2) &
& +(-6*(0.,1._dp)*coup1L*coup2R*coup3*int4*ML1) &
& +(-6*(0.,1._dp)*coup1L*coup2R*coup3*int7*ML1)

! Coefficient of Mat[<v1|7|v2>] Pair[ec[3], k[1]]
Amplitude(4)= &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3*int4*Mex1) &
& +(-4*(0.,1._dp)*coup1L*coup2L*coup3*int6*Mex1) &
& +(-4*(0.,1._dp)*coup1L*coup2L*coup3*int8*Mex1) &
& +(-4*(0.,1._dp)*coup1R*coup2R*coup3*int5*Mex2) &
& +(-4*(0.,1._dp)*coup1R*coup2R*coup3*int6*Mex2) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3*int7*Mex2) &
& +(-6*(0.,1._dp)*coup1R*coup2L*coup3*int4*ML1) &
& +(-6*(0.,1._dp)*coup1R*coup2L*coup3*int7*ML1)

If ((Abs(Amplitude(1))).ne.(Abs(Amplitude(1)))) Then
Amplitude=0._dp
End if
!Field Order
! FieldOrder[1,{F, V, V},1]={External[1][gt1], Internal[1][i1], AntiField[Internal[2]][i2]};
! FieldOrder[1,{F, V, V},2]={AntiField[Internal[1]][i1], External[2][gt2], Internal[3][i3]};
! FieldOrder[1,{F, V, V},3]={External[3][gt3], Internal[2][i2], AntiField[Internal[3]][i3]};
End Subroutine Amp_VERTEX_FtoFV_Topology1_FVV

Subroutine SquareAmp_FtoFV(Mex1,Mex2,Mex3,Amplitude,AmplitudeC,AmpSquared)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(in) :: Amplitude(4), AmplitudeC(4)
Real(dp), Intent(out) :: AmpSquared
Real(dp) :: LorentzFactor, et

AmpSquared = 0._dp

If (Mex3.gt.1._dp) Then

! Mat[<v1|6,ec[3]|v2>] x Mat[<v1|6,ec[3]|v2>]
LorentzFactor = (Mex1**4 + Mex2**4 + Mex2**2*Mex3**2 - 2*Mex3**4 + Mex1**2*(-2*Mex2**2 + Mex3**2))/(4._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(1)*Conjg(AmplitudeC(1))

! Mat[<v1|6,ec[3]|v2>] x Mat[<v1|7,ec[3]|v2>]
LorentzFactor = (-3*Mex1*Mex2)/2._dp
AmpSquared = AmpSquared + LorentzFactor*Amplitude(1)*Conjg(AmplitudeC(2))

! Mat[<v1|6,ec[3]|v2>] x Mat[<v1|6|v2>] Pair[ec[3], k[1]]
LorentzFactor = -(Mex2*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(1)*Conjg(AmplitudeC(3))

! Mat[<v1|6,ec[3]|v2>] x Mat[<v1|7|v2>] Pair[ec[3], k[1]]
LorentzFactor = -(Mex1*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(1)*Conjg(AmplitudeC(4))

! Mat[<v1|7,ec[3]|v2>] x Mat[<v1|6,ec[3]|v2>]
LorentzFactor = (-3*Mex1*Mex2)/2._dp
AmpSquared = AmpSquared + LorentzFactor*Amplitude(2)*Conjg(AmplitudeC(1))

! Mat[<v1|7,ec[3]|v2>] x Mat[<v1|7,ec[3]|v2>]
LorentzFactor = (Mex1**4 + Mex2**4 + Mex2**2*Mex3**2 - 2*Mex3**4 + Mex1**2*(-2*Mex2**2 + Mex3**2))/(4._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(2)*Conjg(AmplitudeC(2))

! Mat[<v1|7,ec[3]|v2>] x Mat[<v1|6|v2>] Pair[ec[3], k[1]]
LorentzFactor = -(Mex1*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(2)*Conjg(AmplitudeC(3))

! Mat[<v1|7,ec[3]|v2>] x Mat[<v1|7|v2>] Pair[ec[3], k[1]]
LorentzFactor = -(Mex2*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(2)*Conjg(AmplitudeC(4))

! Mat[<v1|6|v2>] Pair[ec[3], k[1]] x Mat[<v1|6,ec[3]|v2>]
LorentzFactor = -(Mex2*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(3)*Conjg(AmplitudeC(1))

! Mat[<v1|6|v2>] Pair[ec[3], k[1]] x Mat[<v1|7,ec[3]|v2>]
LorentzFactor = -(Mex1*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(3)*Conjg(AmplitudeC(2))

! Mat[<v1|6|v2>] Pair[ec[3], k[1]] x Mat[<v1|6|v2>] Pair[ec[3], k[1]]
LorentzFactor = ((Mex1**2 + Mex2**2 - Mex3**2)*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(16._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(3)*Conjg(AmplitudeC(3))

! Mat[<v1|6|v2>] Pair[ec[3], k[1]] x Mat[<v1|7|v2>] Pair[ec[3], k[1]]
LorentzFactor = (Mex1*Mex2*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(3)*Conjg(AmplitudeC(4))

! Mat[<v1|7|v2>] Pair[ec[3], k[1]] x Mat[<v1|6,ec[3]|v2>]
LorentzFactor = -(Mex1*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(4)*Conjg(AmplitudeC(1))

! Mat[<v1|7|v2>] Pair[ec[3], k[1]] x Mat[<v1|7,ec[3]|v2>]
LorentzFactor = -(Mex2*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(4)*Conjg(AmplitudeC(2))

! Mat[<v1|7|v2>] Pair[ec[3], k[1]] x Mat[<v1|6|v2>] Pair[ec[3], k[1]]
LorentzFactor = (Mex1*Mex2*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(8._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(4)*Conjg(AmplitudeC(3))

! Mat[<v1|7|v2>] Pair[ec[3], k[1]] x Mat[<v1|7|v2>] Pair[ec[3], k[1]]
LorentzFactor = ((Mex1**2 + Mex2**2 - Mex3**2)*(Mex1**4 + (Mex2**2 - Mex3**2)**2 - 2*Mex1**2*(Mex2**2 + Mex3**2)))/(16._dp*Mex3**2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(4)*Conjg(AmplitudeC(4))

Else

et = Mex1-Mex2

! Mat[<v1|6,ec[3]|v2>] x Mat[<v1|6,ec[3]|v2>]
LorentzFactor = (Mex1**2 + Mex2**2)/2._dp
AmpSquared = AmpSquared + LorentzFactor*Amplitude(1)*Conjg(AmplitudeC(1))

! Mat[<v1|6,ec[3]|v2>] x Mat[<v1|7,ec[3]|v2>]
LorentzFactor = (-Mex1*Mex2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(1)*Conjg(AmplitudeC(2))

! Mat[<v1|6,ec[3]|v2>] x Mat[<v1|6|v2>] Pair[ec[3], k[1]]
LorentzFactor = -Mex2*et
AmpSquared = AmpSquared + LorentzFactor*Amplitude(1)*Conjg(AmplitudeC(3))

! Mat[<v1|6,ec[3]|v2>] x Mat[<v1|7|v2>] Pair[ec[3], k[1]]
LorentzFactor = -Mex1*et
AmpSquared = AmpSquared + LorentzFactor*Amplitude(1)*Conjg(AmplitudeC(4))

! Mat[<v1|7,ec[3]|v2>] x Mat[<v1|6,ec[3]|v2>]
LorentzFactor = (-Mex1*Mex2)
AmpSquared = AmpSquared + LorentzFactor*Amplitude(2)*Conjg(AmplitudeC(1))

! Mat[<v1|7,ec[3]|v2>] x Mat[<v1|7,ec[3]|v2>]
LorentzFactor = (Mex1**2 + Mex2**2)/2._dp
AmpSquared = AmpSquared + LorentzFactor*Amplitude(2)*Conjg(AmplitudeC(2))

! Mat[<v1|7,ec[3]|v2>] x Mat[<v1|6|v2>] Pair[ec[3], k[1]]
LorentzFactor = -Mex1*et
AmpSquared = AmpSquared + LorentzFactor*Amplitude(2)*Conjg(AmplitudeC(3))

! Mat[<v1|7,ec[3]|v2>] x Mat[<v1|7|v2>] Pair[ec[3], k[1]]
LorentzFactor = -Mex2*et
AmpSquared = AmpSquared + LorentzFactor*Amplitude(2)*Conjg(AmplitudeC(4))

! Mat[<v1|6|v2>] Pair[ec[3], k[1]] x Mat[<v1|6,ec[3]|v2>]
LorentzFactor = -Mex2*et
AmpSquared = AmpSquared + LorentzFactor*Amplitude(3)*Conjg(AmplitudeC(1))

! Mat[<v1|6|v2>] Pair[ec[3], k[1]] x Mat[<v1|7,ec[3]|v2>]
LorentzFactor = -Mex1*et
AmpSquared = AmpSquared + LorentzFactor*Amplitude(3)*Conjg(AmplitudeC(2))

! Mat[<v1|6|v2>] Pair[ec[3], k[1]] x Mat[<v1|6|v2>] Pair[ec[3], k[1]]
LorentzFactor = (Mex1**2 + Mex2**2)/2._dp
AmpSquared = AmpSquared + LorentzFactor*Amplitude(3)*Conjg(AmplitudeC(3))

! Mat[<v1|6|v2>] Pair[ec[3], k[1]] x Mat[<v1|7|v2>] Pair[ec[3], k[1]]
LorentzFactor = Mex1*Mex2
AmpSquared = AmpSquared + LorentzFactor*Amplitude(3)*Conjg(AmplitudeC(4))

! Mat[<v1|7|v2>] Pair[ec[3], k[1]] x Mat[<v1|6,ec[3]|v2>]
LorentzFactor = -Mex1*et
AmpSquared = AmpSquared + LorentzFactor*Amplitude(4)*Conjg(AmplitudeC(1))

! Mat[<v1|7|v2>] Pair[ec[3], k[1]] x Mat[<v1|7,ec[3]|v2>]
LorentzFactor = -Mex2*et
AmpSquared = AmpSquared + LorentzFactor*Amplitude(4)*Conjg(AmplitudeC(2))

! Mat[<v1|7|v2>] Pair[ec[3], k[1]] x Mat[<v1|6|v2>] Pair[ec[3], k[1]]
LorentzFactor = Mex1*Mex2
AmpSquared = AmpSquared + LorentzFactor*Amplitude(4)*Conjg(AmplitudeC(3))

! Mat[<v1|7|v2>] Pair[ec[3], k[1]] x Mat[<v1|7|v2>] Pair[ec[3], k[1]]
LorentzFactor = (Mex1**2 + Mex2**2)/2._dp
AmpSquared = AmpSquared + LorentzFactor*Amplitude(4)*Conjg(AmplitudeC(4))

End if

End Subroutine SquareAmp_FtoFV
End Module DecayFFV


text:
// Generated Files
${PROJ_DIR}/lpif/rtl/lpif_txrx/lpif_txrx_x8_q2_master_top.sv
${PROJ_DIR}/lpif/rtl/lpif_txrx/lpif_txrx_x8_q2_master_concat.sv
${PROJ_DIR}/lpif/rtl/lpif_txrx/lpif_txrx_x8_q2_master_name.sv

// Logic Link files
-f ${PROJ_DIR}/llink/rtl/llink.f

// Common Files
-f ${PROJ_DIR}/common/rtl/common.f

max_stars_repo_name: chipsalliance/aib-protocols

max_issues_repo_name: chipsalliance/aib-protocols

max_forks_repo_name: chipsalliance/aib-protocols


text:
! Copyright (C) 2011
! Free Software Foundation, Inc.

! This file is part of the gtk-fortran GTK+ Fortran Interface library.

! This is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 3, or (at your option)
! any later version.

! This software is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.

! Under Section 7 of GPL version 3, you are granted additional
! permissions described in the GCC Runtime Library Exception, version
! 3.1, as published by the Free Software Foundation.

! You should have received a copy of the GNU General Public License along with
! this program; see the files COPYING3 and COPYING.RUNTIME respectively.
! If not, see <http://www.gnu.org/licenses/>.
!
! Contributed by James Tappin
! Last modification: 11-21-2011

! --------------------------------------------------------
! gtk-hl-combobox.f90
! Generated: Tue Oct 29 17:02:03 2013 GMT
! Please do not edit this file directly,
! Edit gtk-hl-combobox-tmpl.f90, and use ./mk_gtk_hl.pl to regenerate.
! Generated for GTK+ version: 2.24.0.
! Generated for GLIB version: 2.38.0.
! --------------------------------------------------------

!*
! ComboBox
module gtk_hl_combobox
  ! This interface implements the GtkComboBoxText widget for making a chooser.
  ! While this has more limited capabilities than the full GtkComboBox, it
  ! is adequate for the vast majority of uses.
  !/

  use gtk_sup
  use iso_c_binding
  ! Autogenerated use's
  use gtk, only: gtk_combo_box_get_active, gtk_combo_box_new,&
       & gtk_combo_box_set_active, &
       & gtk_widget_set_sensitive, gtk_widget_set_tooltip_text,&
!!$GTK< 2.24!       & gtk_combo_box_get_active_text, &
!!$GTK< 2.24!       & gtk_combo_box_new_with_entry, &
!!$GTK< 2.24!       & gtk_combo_box_new_text, gtk_combo_box_entry_new_text, &
!!$GTK< 2.24!       & gtk_combo_box_append_text,&
!!$GTK< 2.24!       & gtk_combo_box_insert_text,&
!!$GTK< 2.24!       & gtk_combo_box_prepend_text,&
!!$GTK< 2.24!       & gtk_combo_box_remove_text,&
       & gtk_combo_box_text_get_active_text, gtk_combo_box_text_new,&
       & gtk_combo_box_text_new_with_entry, &
       & gtk_combo_box_text_remove, &
       & gtk_combo_box_text_append_text,&
       & gtk_combo_box_text_insert_text,&
       & gtk_combo_box_text_prepend_text,&
!!$GTK>=3.0!       & gtk_combo_box_text_remove_all, &
       & gtk_combo_box_get_model, gtk_tree_model_iter_n_children, &
       & TRUE, FALSE, g_signal_connect

  implicit none

  interface hl_gtk_combo_box_delete
     module procedure hl_gtk_combo_box_delete_single
     module procedure hl_gtk_combo_box_delete_multi
     module procedure hl_gtk_combo_box_delete_all
  end interface hl_gtk_combo_box_delete

contains
  !+
  function hl_gtk_combo_box_new(has_entry, changed, data, initial_choices, &
       & sensitive, tooltip, active) result(cbox)

    type(c_ptr) :: cbox
    integer(kind=c_int), intent(in), optional :: has_entry
    type(c_funptr), optional :: changed
    type(c_ptr), intent(in), optional :: data
    character(len=*), dimension(:), intent(in), optional :: initial_choices
    integer(kind=c_int), intent(in), optional :: sensitive
    character(kind=c_char), dimension(*), optional, intent(in) :: tooltip
    integer(kind=c_int), optional, intent(in) :: active

    ! Creator for the combobox.
    !
    ! HAS_ENTRY: boolean: optional: Set to TRUE to add an entry field.
    ! CHANGED: c_funptr: optional: Callback routine for the "changed" signal.
    ! DATA: c_ptr: optional: User data for the changed callback.
    ! INITIAL_CHOICES: string(): optional: Initial list of choices.
    ! SENSITIVE: boolean: optional: Set to FALSE to make the widget start in an
    ! 		insensitive state.
    ! TOOLTIP: string: optional: A tooltip to display when the pointer is
    ! 		held over the widget.
    ! ACTIVE: c_int: optional: The initial active selection.
    !-

    integer(kind=c_int) :: ientry
    integer(kind=c_int) :: i

    if (present(has_entry)) then
       ientry = has_entry
    else
       ientry = FALSE
    end if

    if (ientry == TRUE) then
!GTK3
       cbox = gtk_combo_box_text_new_with_entry()
!GTK2
!!$GTK< 2.24!       cbox = gtk_combo_box_entry_new_text()
    else
!GTK3
       cbox = gtk_combo_box_text_new()
!GTK2
!!$GTK< 2.24!       cbox =  gtk_combo_box_new_text()
    end if

    if (present(initial_choices)) then
       do i=1,size(initial_choices)
!GTK3
          call gtk_combo_box_text_append_text(cbox, &
               & trim(initial_choices(i))//C_NULL_CHAR)
!GTK2
!!$GTK< 2.24!          call gtk_combo_box_append_text(cbox, &
!!$GTK< 2.24!               & trim(initial_choices(i))//C_NULL_CHAR)
       end do
    end if

    if (present(active)) call gtk_combo_box_set_active(cbox, active)

    if (present(changed)) then
       if (present(data)) then
          call g_signal_connect(cbox, "changed"//C_NULL_CHAR, changed, data)
       else
          call g_signal_connect(cbox, "changed"//C_NULL_CHAR, changed)
       end if
    end if

    if (present(sensitive)) call gtk_widget_set_sensitive(cbox, sensitive)
    if (present(tooltip)) call gtk_widget_set_tooltip_text(cbox, tooltip)
  end function hl_gtk_combo_box_new

  !+
  subroutine hl_gtk_combo_box_add_text(cbox, text, index, at_start)

    type(c_ptr), intent(in) :: cbox
    character(kind=c_char), dimension(*), optional :: text
    integer(kind=c_int), intent(in), optional :: index
    integer(kind=c_int), intent(in), optional :: at_start

    ! Add a new choice to a combo box.
    !
    ! CBOX: c_ptr: required: The combo box to modify.
    ! TEXT: string: required: The text to add.
    ! INDEX: c_int: optional: The location at which to add the text.
    ! AT_START: boolean: optional: If set to TRUE and INDEX is not given
    ! 		then add the text at the start of the list.
    !
    ! If neither INDEX nor AT_START is present the text is appended.
    !-

    integer(kind=c_int) :: prepend

    if (present(index)) then
!GTK3
       call gtk_combo_box_text_insert_text(cbox, index, text)
!GTK2
!!$GTK< 2.24!       call gtk_combo_box_insert_text(cbox, index, text)
    else
       if (present(at_start)) then
          prepend = at_start
       else
          prepend = FALSE
       end if
       if (prepend == TRUE) then
!GTK3
          call gtk_combo_box_text_prepend_text(cbox, text)
!GTK2
!!$GTK< 2.24!          call gtk_combo_box_prepend_text(cbox, text)
       else
!GTK3
          call gtk_combo_box_text_append_text(cbox, text)
!GTK2
!!$GTK< 2.24!          call gtk_combo_box_append_text(cbox, text)
       end if
    end if
  end subroutine hl_gtk_combo_box_add_text

  !+
  subroutine hl_gtk_combo_box_delete_single(cbox, index)

    type(c_ptr), intent(in) :: cbox
    integer(kind=c_int), intent(in) :: index

    ! Delete a line from a combo box
    !
    ! CBOX: c_ptr: required: The combo box to update
    ! INDEX: c_int: required: The index of the choce to remove
    !
    ! Usually called via the generic hl_gtk_combo_box_delete interface.
    !-

!GTK3
    call gtk_combo_box_text_remove(cbox, index)
!GTK2
!!$GTK< 2.24!    call gtk_combo_box_remove_text(cbox, index)

  end subroutine hl_gtk_combo_box_delete_single

  !+
  subroutine hl_gtk_combo_box_delete_multi(cbox, index)

    type(c_ptr), intent(in) :: cbox
    integer(kind=c_int), dimension(:), intent(in) :: index

    ! Delete lines from a combo box
    !
    ! CBOX: c_ptr: required: The combo box to update
    ! INDEX: c_int(): required: The index of the choce to remove
    !
    ! Usually called via the generic hl_gtk_combo_box_delete interface.
    !-

    integer(kind=c_int) :: i

    do i = size(index),1,-1
!GTK3
       call gtk_combo_box_text_remove(cbox, index(i))
!GTK2
!!$GTK< 2.24!       call gtk_combo_box_remove_text(cbox, index(i))
    end do
  end subroutine hl_gtk_combo_box_delete_multi

  !+
  subroutine hl_gtk_combo_box_delete_all(cbox)

    type(c_ptr), intent(in) :: cbox

    ! Delete all lines from a combo box
    !
    ! CBOX: c_ptr: required: The combo box to update
    !
    ! Usually called via the generic hl_gtk_combo_box_delete interface.
    !-

    integer(kind=c_int) :: n, i

!GTK3
!!$GTK>=3.0!    call gtk_combo_box_text_remove_all(cbox)
!GTK2
! Note that by using the hl delete we avoid 3-level conditional compilation.
    n = hl_gtk_combo_box_n_entries(cbox)
    do i = 1, n
       call hl_gtk_combo_box_delete(cbox, 0_c_int)
    end do

  end subroutine hl_gtk_combo_box_delete_all
  !+
  function hl_gtk_combo_box_get_active(cbox, text, ftext) result(index)

    integer(kind=c_int) :: index
    type(c_ptr), intent(in) :: cbox
    type(c_ptr), intent(out), optional :: text
    character(len=*), intent(out), optional :: ftext

    ! Get the selection from a combo box
    !
    ! CBOX: c_ptr: required: The combo box to query.
    ! TEXT: c_ptr: optional: C pointer to the text.
    ! FTEXT: fstring: optional: The string as a Fortran string.
    !-

    type(c_ptr), target :: ctext

    index = gtk_combo_box_get_active(cbox)

    if (present(text) .or. present(ftext)) then

!GTK3
      ctext = gtk_combo_box_text_get_active_text(cbox)
!GTK2
!!$GTK< 2.24!       ctext = gtk_combo_box_get_active_text(cbox)

       ! This is a bit ugly
       if (present(ftext)) &
            & call convert_c_string(ctext, ftext)

       if (present(text)) text=ctext
    end if
  end function hl_gtk_combo_box_get_active

  !+
  function hl_gtk_combo_box_n_entries(cbox) result(count)
    integer(kind=c_int) :: count
    type(c_ptr), intent(in) :: cbox

    ! Find the number of entries in a combo box.
    !
    ! CBOX: c_ptr: required: The combo box to query.
    !-

    type(c_ptr) :: model

    model = gtk_combo_box_get_model(cbox)
    count = gtk_tree_model_iter_n_children (model, C_NULL_PTR)

  end function hl_gtk_combo_box_n_entries

end module gtk_hl_combobox

max_stars_repo_name: jshahbazi/fortran-gui-demo

max_issues_repo_name: jshahbazi/fortran-gui-demo

max_forks_repo_name: jshahbazi/fortran-gui-demo


text:
      program calcSAL

!    note that the documentation for spherepack programs shaec.f, shsec.f
!          along with the example program on page 33 of the spherepack2.0
!          manual helped guide the writing of this program

!     set grid size with parameter statements

      PARAMETER ( NI=4500, NJ=NI/2+1)  !0.08 degrees
      PARAMETER (nnlat=NJ,nnlon=NI)

!     set saved and unsaved work space lengths in terms of nnlat, nnlon

      parameter (nn15=nnlon+15)
!      parameter (llsave=nnlat*(nnlat+1)+3*((nnlat-2)*(nnlat-1)+nn15))
!      parameter (llwork=(nnlat+1)*(nnlon+3*nnlat)+nnlat*(2*nnlat+1))
      parameter (l1use=nnlat,l2=(nnlat+1)/2)
      parameter (llsave=2*nnlat*l2+3*((l1use-2)*(nnlat+nnlat-l1use-1))
     &        /2+nnlon+15)
      parameter (ldwork=nnlat*(nnlon+nnlon))
      parameter (lshaec=llsave)
      parameter (lshsec=llsave)

!     dimension arrays

      real*8 LOADNUMBERS(10000)

      real*4 amp(ni,nj)
      real*4 phase(ni,nj)
      real*4 loadtide(ni,nj)

      REAL*8 TIDEPAD(nnlon,nnlat)
      real*8 tidepad2(nnlat,nnlon)

      real*8 loadtide2(nnlat,nnlon)

      real*8 dwork(ldwork)
      real*8 wshaec(llsave)
      real*8 wshsec(llsave)
      real*8 a(nnlat,nnlat)
      real*8 b(nnlat,nnlat)

      INTEGER ISTART(4), ICOUNT(4)

      character*60,parameter :: infile = "loadnumbers_binary"

!     define pi

      pi=4.0*atan(1.0)

!     read in loadnumbers

      write(6,'(a,a)') 'reading ',trim(infile)
      call flush(6)
      open(10, file=infile, form='unformatted')
      read(10) loadnumbers
      close (10)

!     read in amplitude and phase

      write(6,'(a,a)') 'reading tide_AmPh.A'
      call flush(6)
      INQUIRE( IOLENGTH=MRECL) amp
      open(10, file='tide_AmPh.A', form='unformatted', access="direct", RECL=MRECL)
      read(10, rec=1) amp
      read(10, rec=2) phase
      close (10)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!     calculate loadtide for omegat=0 degrees:
      write(6,'(a)') 'calculate loadtide for omegat=0 degrees:'
      call flush(6)

         do j=1,NJ
           do i=1,NI
              if (amp(i,j) .gt. 1.0d10) then
                tidepad(i,j)=0.0
             else
                tidepad(i,j)=amp(i,j)*cos((pi/180.0)*(0.0-phase(i,j)))
             endif
           enddo
         enddo

!     translate tidepad from geophysical coordinates (which use longitude
!       and latitude, the latter measured from south to north pole) to
!       tidepad2 in mathematical spherical coordinates (which use colatitude
!       and longitude, the former measured from north to south pole)

      do j=1,nnlon
         do i=1,nnlat
            tidepad2(i,j)=tidepad(j,nnlat+1-i)
         enddo
      enddo

!
! initialize spherepack routines
!
      call shaeci(nnlat,nnlon,wshaec,lshaec,dwork,ldwork,ierror)
      if (ierror .ne. 0) print *, 'shaeci ierror = ',ierror
      call shseci(nnlat,nnlon,wshsec,lshsec,dwork,ldwork,ierror)
      if (ierror .ne. 0) print *, 'shseci ierror = ',ierror

!  calculate spherical harmonics

      call shaec(nnlat,nnlon,0,1,tidepad2,nnlat,nnlon,a,b,nnlat,
     &                 nnlat,wshaec,lshaec,dwork,ldwork,ierror)
      if (ierror .ne. 0) print *, 'shaec ierror = ',ierror

!  multiply the harmonics by the load numbers; note start at j=2 since
!        j=1, i = 1 denotes P00; note set a(1,1)=b(1,1)=0 to eliminate
!        possibility of nonzero mean

      a(1,1)=0.0
      b(1,1)=0.0

      do j = 2,nnlat
         do i = 1,j
             a(i,j) = a(i,j)*loadnumbers(j-1)
             b(i,j) = b(i,j)*loadnumbers(j-1)
         enddo
      enddo

!  transform back to physical coordinates

      call shsec(nnlat,nnlon,0,1,loadtide2,nnlat,nnlon,a,b,nnlat,
     &                 nnlat,wshsec,lshsec,dwork,ldwork,ierror)
      if (ierror .ne. 0) print *, 'shsec ierror = ',ierror

!  translate loadtide2 from mathematical spherical coordinates to
!      loadtide in geophysical coordinates.
!  also flip the sign so that this is the SAL *correction* to the tides
!   which is approximately 180 degrees out of phase with the tides.

      do j=1,nnlon
         do i=1,nnlat
            loadtide(j,i)=-loadtide2(nnlat+1-i,j)
         enddo
      enddo

!  sal_ReIm follows the OSU (TPXO) convention for the Imaginary
!   component: phase = atan2(-Im,Re) and amp = sqrt( Re**2 + Im**2 ).

      write(6,'(a)') 'writing Re to sal_ReIm.A'
      call flush(6)
      open( 80, file='sal_ReIm.A', form='unformatted', access="direct", RECL=MRECL)
      write(80, rec=1)loadtide

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!     calculate loadtide for omegat=90 degrees:
      write(6,'(a)') 'calculate loadtide for omegat=90 degrees:'
      call flush(6)

         do j=1,NJ
           do i=1,NI
              if (amp(i,j) .gt. 1.0d10) then
                tidepad(i,j)=0.0
             else
                tidepad(i,j)=amp(i,j)*cos((pi/180.0)*(90.0-phase(i,j)))
             endif
           enddo
         enddo

!     translate tidepad from geophysical coordinates (which use longitude
!       and latitude, the latter measured from south to north pole) to
!       tidepad2 in mathematical spherical coordinates (which use colatitude
!       and longitude, the former measured from north to south pole)

      do j=1,nnlon
         do i=1,nnlat
            tidepad2(i,j)=tidepad(j,nnlat+1-i)
         enddo
      enddo

!
! initialize spherepack routines
!
      call shaeci(nnlat,nnlon,wshaec,lshaec,dwork,ldwork,ierror)
      if (ierror .ne. 0) print *, 'shaeci ierror = ',ierror
      call shseci(nnlat,nnlon,wshsec,lshsec,dwork,ldwork,ierror)
      if (ierror .ne. 0) print *, 'shseci ierror = ',ierror

!  calculate spherical harmonics

      call shaec(nnlat,nnlon,0,1,tidepad2,nnlat,nnlon,a,b,nnlat,
     &                 nnlat,wshaec,lshaec,dwork,ldwork,ierror)
      if (ierror .ne. 0) print *, 'shaec ierror = ',ierror

!  multiply the harmonics by the load numbers; note start at j=2 since
!        j=1, i = 1 denotes P00; note set a(1,1)=b(1,1)=0 to eliminate
!        possibility of nonzero mean

      a(1,1)=0.0
      b(1,1)=0.0

      do j = 2,nnlat
         do i = 1,j
             a(i,j) = a(i,j)*loadnumbers(j-1)
             b(i,j) = b(i,j)*loadnumbers(j-1)
         enddo
      enddo

!  transform back to physical coordinates

      call shsec(nnlat,nnlon,0,1,loadtide2,nnlat,nnlon,a,b,nnlat,
     &                 nnlat,wshsec,lshsec,dwork,ldwork,ierror)
      if (ierror .ne. 0) print *, 'shsec ierror = ',ierror

!  translate loadtide2 from mathematical spherical coordinates to
!      loadtide in geophysical coordinates.
!  sal_ReIm follows the OSU (TPXO) convention for the Imaginary
!   component: phase = atan2(-Im,Re).
!  this implies that the sign of loadtide should be flipped, but we
!  then flip the sign again so that this is the SAL *correction* to the
!  tides which is approximately 180 degrees out of phase with the tides.
!  overall, leave the sign of loadtide2 "as is".

      do j=1,nnlon
         do i=1,nnlat
            loadtide(j,i)=loadtide2(nnlat+1-i,j)
         enddo
      enddo

      write(6,'(a)') 'writing Im to sal_ReIm.A'
      call flush(6)
      write(80, rec=2) loadtide
      close(80)
      endprogram calcSAL

max_stars_repo_name: HYCOM/HYCOM-examples

max_issues_repo_name: HYCOM/HYCOM-examples

max_forks_repo_name: HYCOM/HYCOM-examples


text:
!*******************************************************************
!                                                                  *
      INTEGER FUNCTION JTHN(K,N,I)
!                                                                  *
!   Written by G. Gaigalas,                                        *
!   Vanderbilt University,  Nashville               October  1996  *
!   Transform to fortran 90/95 by G. Gaigalas       December 2012  *
!   The last modification made by G. Gaigalas       October  2017  *
!                                                                  *
!*******************************************************************
!
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: I, N, K
!-----------------------------------------------
      IF(N == 1) THEN
        JTHN=MOD(K,I)
      ELSEIF(N == 2) THEN
        JTHN=MOD(K/I,I)
      ELSEIF(N == 3) THEN
        JTHN=MOD(K/(I*I),I)
      ELSEIF(N == 4) THEN
        JTHN=MOD(K/(I*I*I),I)
      ELSEIF(N == 5) THEN
        JTHN=MOD(K/(I*I*I*I),I)
      ELSE
        WRITE(6,'(A)') ' ERROR IN JTHN '
        STOP
      ENDIF
      RETURN
      END FUNCTION JTHN

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


module eqempirical

character(len=512) :: relation
double precision :: magnitude

end module

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

program main
!----
! Print fault characteristics (dimensions, slip) from published empirical relations
!----

use io, only: stdout
use eq, only: empirical
use eqempirical, only: relation, magnitude

implicit none

! Local variables
double precision :: len, wid

! Parse command line inputs
call gcmdln()
if (relation.eq.'') then
    call usage('eqempirical: empirical relation not defined')
endif
if (magnitude.lt.-10.0d0 .and. relation.ne.'list') then
    call usage('eqempirical: magnitude not defined')
endif

if (relation.eq.'WC') then
    call empirical(magnitude,wid,len,'WC:print','')
elseif (relation.eq.'MB') then
    call empirical(magnitude,wid,len,'MB:print','')
elseif (relation.eq.'B') then
    call empirical(magnitude,wid,len,'B:print','')
elseif (relation.eq.'YM') then
    call empirical(magnitude,wid,len,'YM:print','')
elseif (relation.eq.'AH') then
    call empirical(magnitude,wid,len,'AH:print','')
elseif (relation.eq.'list') then
    write(stdout,*) 'Empirical Relations:'
    write(stdout,*) '  WC: Wells and Coppersmith, BSSA (1994)'
    write(stdout,*) '  MB: Mai and Beroza, BSSA (2000)'
    write(stdout,*) '  B: Blaser et al., BSSA (2010)'
    write(stdout,*) '  YM: Yen and Ma, BSSA (2011)'
    write(stdout,*) '  AH: Allen and Hayes, BSSA (2017)'
else
    call usage('eqempirical: no empirical relationship named '//trim(relation))
endif

end

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

subroutine gcmdln()

use eqempirical, only: relation, &
                       magnitude
implicit none

! Local variables
integer :: i, ios, narg
character(len=512) :: tag, arg

! Initialize control parameters
ios = 0
relation = ''
magnitude = -1d10

! Number of arguments
narg = command_argument_count()
if (narg.eq.0) then
  call usage('')
endif

! Read arguments
i = 1
do while (i.le.narg)

    call get_command_argument(i,tag)

    if (tag.eq.'-model') then
        i = i + 1
        call get_command_argument(i,arg,status=ios)
        read (arg,*,err=9032,iostat=ios) relation
    elseif (tag.eq.'-mag') then
        i = i + 1
        call get_command_argument(i,arg,status=ios)
        read (arg,*,err=9032,iostat=ios) magnitude
    elseif (tag.eq.'-list') then
        relation = 'list'
    else
        call usage('')
    endif

    9032 if (ios.ne.0) then
        call usage('eqempirical: error parsing "'//trim(tag)//'" flag arguments')
    endif

    i = i + 1
enddo

return
end subroutine

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

subroutine usage(str)

use io, only: stderr

implicit none

character(len=*) :: str

if (str.ne.'') then
    write(stderr,*) trim(str)
    write(stderr,*)
endif

write(stderr,*) 'Usage: eqempirical -model MODEL -mag MAGNITUDE'
write(stderr,*)
write(stderr,*) '-model MODEL      Empirical relation name'
write(stderr,*) '-mag MAGNITUDE    Earthquake magnitude'
write(stderr,*) '-list             List available models'
write(stderr,*)

call error_exit(1)
end subroutine


text:
! should fail 'implicit-none' test: 1 case detected
! should fail 'implicit-none-all' test: 2 cases detected
program example1
  implicit none ! ERROR: IMPLICIT NONE must follow all USE statements

  use file1

  != unit(m) :: dist2
  real :: dist2 = 20

  != unit(m) :: dist3
  real :: dist3

  != unit(s) :: time2
  real :: time2

  ! Should be consistent
  dist2 = add(dist1, dist2)

  ! Should be inconsistent
  dist2 = add(time1, dist2)

  ! Should be inconsistent
  dist3 = add(time2, dist3)

contains
  real function add(a,b)
    real :: a, b
    add = a + b
  end function add
end program example1

max_stars_repo_name: raehik/camfort

max_issues_repo_name: apthorpe/camfort

max_forks_repo_name: raehik/camfort


text:
!
! Copyright (c) 2020 Bálint Aradi, Universität Bremen. All rights reserved.
! Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
!

!> Wrapper for node.
module saydx_node
  use saydx_ciface, only : c_node_dereference, c_node_get_name, c_node_has_parent
  use saydx_cinterop, only : c_ptr, c_null_ptr, c_associated, c_f_string_pointer
  use saydx_common, only : fatal_error
  use saydx_commontypes, only : c_ptr_movable_wrapper_t
  implicit none
  private

  public :: node_t

  type, extends(c_ptr_movable_wrapper_t) :: node_t
  contains
    procedure :: get_name => node_get_name
    procedure :: has_parent => node_has_parent
    final :: node_final
  end type node_t

contains

  function node_get_name(this) result(pstr)
    class(node_t), intent(in) :: this
    character(:), pointer :: pstr

    call c_f_string_pointer(c_node_get_name(this%cptr), pstr)

  end function node_get_name

  function node_has_parent(this) result(hasparent)
    class(node_t), intent(in) :: this
    logical :: hasparent

    hasparent = logical(c_node_has_parent(this%cptr))

  end function node_has_parent

  subroutine node_final(this)
    type(node_t), intent(inout) :: this

    if (.not. c_associated(this%cptr)) then
      return
    end if
    call c_node_dereference(this%cptr)
    this%cptr = c_null_ptr

  end subroutine node_final

end module saydx_node

max_stars_repo_name: saydx/libsaydx

max_issues_repo_name: saydx/libsaydx

max_forks_repo_name: saydx/libsaydx


text: James Tommy Kunze served on the City Council from 4/20/1956 to 9/20/1958.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine tstrdx(cdrs,iodb,iopt,jflag,jsflag,narn1,narn2,
     $ ndrs,ndrsmx,ndrsr,nodbmx,noptmx,noutpt,nrdxsp,nstmax,
     $ qredox,uspec)
c
c     This subroutine determines if the chemical model to be
c     computed has a redox aspect. This will be determined to be
c     so if a species in the model has an associated reaction
c     that is a redox reaction and the species is not in the
c     active basis set.
c
c     An auxiliary basis species (say Oxalate-) that is in the
c     model but is included in the active basis set is effectively.
c     treated as detached from a corresponding strict basis species
c     (e.g., the concentration/activity of Oxalate- is not determined
c     by the assumption of equilibrium for its associated reaction,
c     which would link it to HCO3-). In effect, Oxalate- is treated
c     as being composed of a pseudo-element, and its presence in the
c     model does not require a redox variable.
c
c     This subroutine is called by:
c
c       EQ6/eq6.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c
c     Principal output:
c
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ndrsmx,nodbmx,noptmx,nstmax
c
      integer noutpt
c
      integer iodb(nodbmx),iopt(noptmx),jflag(nstmax),jsflag(nstmax),
     $ ndrs(ndrsmx),ndrsr(2,nstmax)
c
      integer narn1,narn2,nrdxsp
c
      logical qredox
c
      character(len=48) uspec(nstmax)
c
      real(8) cdrs(ndrsmx)
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer jlen,ns
c
      character(len=56) uspn56
c
      real(8) cx
c
      real(8) coefdr
c
c-----------------------------------------------------------------------
c
      qredox = .false.
c
      if (iopt(15) .le. 0) then
        do ns = narn1,narn2
c
          if (jflag(ns) .eq. 30) then
c
c           The species is not in the active basis set. It is a
c           "dependent" species whose concentration/activity is
c           computed assuming its associated reaction is in a state
c           of equilibrium.
c
            if (jsflag(ns) .lt. 2) then
c
c             The species not hard suppressed.
c
c             Calling sequence substitutions:
c               nrdxsp for nse
c
              cx = coefdr(cdrs,ndrs,ndrsmx,ndrsr,nrdxsp,ns,nstmax)
              if (cx .ne. 0.) then
                qredox = .true.
                if (iodb(1) .ge. 1) then
c
c                 Calling sequence substitutions:
c                   uspec(ns) for unam48
c
                  call fmspnm(jlen,uspec(ns),uspn56)
                  write (noutpt,1000) uspn56(1:jlen)
 1000             format(/' * Note - (EQ6/tstrdx) The reaction',
     $            ' associated with the species',/7x,a,' is a redox',
     $            ' reaction.')
                  go to 999
                endif
              endif
            endif
          endif
        enddo
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
  999 continue
      end


text:
C Copyright(C) 1999-2021 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE SIORPT(MODULE,MESS,DISP)
      IMPLICIT INTEGER (A-Z)
      CHARACTER*(*) MODULE
      CHARACTER*(*) MESS
      CHARACTER*256 LOCLIN
      CHARACTER*20 ERRORT
      CHARACTER*120 ERRMSG
      CHARACTER*8 ERRMOD

      ERRMOD = MODULE
      ERRMSG = MESS
      CALL CPUDAC(ERRORT)
      CALL CHRTRM(ERRORT,LT)
      CALL CHRTRM(ERRMOD,L)
      IF (L.EQ.0) THEN
         L = 1
      END IF

      IF (DISP.EQ.1) THEN

      ELSE IF (DISP.EQ.2) THEN
         LOCLIN = '#PLT error (warning) in module '//ERRMOD(1:L)//
     *            ' at '//ERRORT(1:LT)
         WRITE (6,10) LOCLIN

         LOCLIN = '#Error message: '//ERRMSG
         WRITE (6,10) LOCLIN

      ELSE IF (DISP.EQ.3) THEN
         LOCLIN = '#PLT error (traceback) in module '//ERRMOD(1:L)//
     *            ' at '//ERRORT(1:LT)
         WRITE (6,10) LOCLIN
         LOCLIN = '#Error message: '//ERRMSG
         WRITE (6,10) LOCLIN

      ELSE IF (DISP.EQ.4) THEN
         LOCLIN = '#PLT error (fatal) in module '//ERRMOD(1:L)//' at '//
     *            ERRORT(1:LT)
         WRITE (6,10) LOCLIN
         LOCLIN = '#Error message: '//ERRMSG
         WRITE (6,10) LOCLIN

      ELSE
         LOCLIN = '#PLT error (fatal) in module '//ERRMOD(1:L)//' at '//
     *            ERRORT(1:LT)
         WRITE (6,10) LOCLIN
         LOCLIN = '#Error message: '//ERRMSG
         WRITE (6,10) LOCLIN
      END IF

 10   FORMAT (1X,A)

      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
      program test
      real*16 sum,sum2,rmsd,energy,ave
      character*16 fname

      ave=0.0
      rmsd=0.0
      sum=0.0
      sum2=0.0
      read(5,*)timestep,nline
      read(5,*)fname

      open(10,file=fname)

      read(10,*)
      do i=1,nline
         read(10,20)energy
         sum=sum+energy
         sum2=sum2+energy**2
      enddo
      close(10)
 20   format(5x,f12.6)

c      write(6,*)sum**2/nline,sum2
      ave=sum/nline
      rmsd=(sum2-sum*sum/nline)/nline
      rmsd=sqrt(rmsd)

      write(6,10)timestep,rmsd,ave
 10   format(3f16.7)
      end

max_issues_repo_name: IBM/BlueMatter

max_forks_repo_name: IBM/BlueMatter


module Constants
  implicit none

  real, parameter :: pi = 4 * atan(1.0)
  real, parameter :: E_e = 510998.91013

  ! we now have enumerators in F2003/8, for the sake of interop with C
  enum, bind(c) ! unnamed 1
    enumerator :: red =1, blue, black =5
    enumerator yellow
    enumerator gold, silver, bronze
    enumerator :: purple
    enumerator :: pink, lavender
  end enum

  enum ! unnamed 2
    enumerator :: a, b, c
  end enum

  enum :: Named1
    enumerator :: x1, y1, z1
  end enum

  enum Named2
    enumerator :: x2, y2, z2
  end enum

  enum(8) Named3
    enumerator :: x3, y3, z3
  end enum

  enum*8 Named4
    enumerator :: x4, y4, z4
  end enum

  enum(8) :: Named5
    enumerator :: x5, y5, z5
  end enum

  enum*8 :: Named6
    enumerator :: x6, y6, z6
  end enum

  enum, bind(c) :: Named7
    enumerator :: x7, y7, z7
  end enum

  real, parameter :: hc = 12398.4193

  public

end module Constants


text:
!*robodoc*u* tests/testABN
!  NAME
!    testABN
!  SYNOPSIS

!$Id: testABN.f90 380 2017-03-22 11:03:09Z mexas $

program testABN

!  PURPOSE
!    Testing solidification with MPI/IO, cgca_m2mpiio/cgca_pswci
!  DESCRIPTION
!    Solidification with Cray reduction, CO_SUM (cgca_sld3).
!    Timing output of MPI/IO (cgca_pswci) against
!    the serial version (cgca_swci).
!  NOTES
!    The program must be called with 2 command line arguments,
!    both positive integers. These are codimensions along 1 and 2.
!    The number of images must be such that
!    codimension3 = num_images()/( codimension1 * codimension3 )
!    is a positive integer. Example:
!      ./testABN.x 2 2
!    which will make the third codimension equal to 16/(2*2)=4.
!  AUTHOR
!    Anton Shterenlikht
!  COPYRIGHT
!    See LICENSE
!  USES
!    cgca testaux
!  USED BY
!    Part of CGPACK test suite
!  SOURCE

use testaux

implicit none

real,parameter :: gigabyte=real(2**30), resolution=1.0e-5, &
 loge2 = log(real(2))
logical(kind=ldef),parameter :: yesdebug = .true., nodebug = .false.

integer(kind=idef) :: l1,u1,l2,u2,l3,u3,col1,cou1,col2,cou2,col3,cou3, &
 nuc,    & ! number of nuclei in the model
 nimages, codim(3)[*], p
integer(kind=iarr),allocatable :: space(:,:,:,:)[:,:,:]
integer(kind=ilrg) :: icells, mcells, img

real :: time1, time2, fsizeb, fsizeg, tdiff

logical(kind=ldef) :: solid

!*********************************************************************72
! first executable statement

img = this_image()
nimages = num_images()

! do a check on image 1
if ( img .eq. 1 ) then
 call getcodim(nimages,codim)
 ! print a banner
 call banner("ABN")
 ! print the parameter values
 call cgca_pdmp
 write (*,'(a,i0,a)') "running on ", nimages, " images in a 3D grid"
 ! dump the value of p
 write (*,"(a,i0)") "p=",p
end if

! all images read codim from image 1
sync all
codim(:) = codim(:)[1]

l1 = 1
l2 = l1
l3 = l1

! The array size is only controlled by this value in this program.
u1 = 2**7 ! 128
u2 = u1
u3 = u1

col1 = 1
col2 = col1
col3 = col1

cou1 = codim(1)-col1+1
cou2 = codim(2)-col2+1
cou3 = codim(3)-col3+1

! total number of cells in a coarray
icells = int(u1-l1+1,kind=ilrg) * int(u2-l2+1,kind=ilrg) *            &
 int(u3-l3+1,kind=ilrg)

! total number of cells in the model
mcells = icells * int(codim(1),kind=ilrg) * int(codim(2),kind=ilrg) * &
  int(codim(3),kind=ilrg)

! total number of nuclei
nuc = int( resolution * real( mcells ) )

if ( img .eq. 1 ) then
  write (*,'(a,2(i0,":",i0,","),i0,":",i0,")")') &
    "bounds: (",l1,u1,l2,u2,l3,u3
  write (*,'(a,2(i0,":",i0,","),i0,":",i0,")")') &
    "cobounds: (",col1,cou1,col2,cou2,col3,cou3

! Total output file size, in B and in GB.
fsizeb = real( mcells*storage_size(space,kind=ilrg)/8_ilrg )
fsizeg = fsizeb / gigabyte

  write (*,'(a,i0,a)') "Each image has ",icells, " cells"
  write (*,'(a,i0,a)') "The model has ", mcells, " cells"
  write (*,'(a,i0,a)') "The model has ", nuc, " nuclei"
  write (*,'(2(a,es9.2),a)') "The output file size is ", fsizeb, &
   " B, or ", fsizeg, "GB."
end if

! calculate file size in MB
! allocate coarray
call cgca_as(l1,u1,l2,u2,l3,u3,col1,cou1,col2,cou2,col3,1,space)

! initialise coarray to liquid
space = cgca_liquid_state

! populate nuclei
call cgca_nr(space,nuc,nodebug)

! solidify
!subroutine cgca_sld3(coarray,iter,heartbeat,solid)
call cgca_sld3(space,0,1,solid)

! dump the model
call cpu_time(time1)
call cgca_pswci(space, cgca_state_type_grain, 'mpiio.raw')
call cpu_time(time2)
tdiff = time2-time1
if (img .eq. 1) write (*,*) "MPI/IO: ", tdiff, "s, rate: ", fsizeg/tdiff, "GB/s."

call cpu_time(time1)
call cgca_swci (space, cgca_state_type_grain, 10, 'serial.raw')
call cpu_time(time2)
tdiff = time2-time1
if (img .eq. 1) write (*,*) "Serial IO: ", tdiff, "s, rate: ", fsizeg/tdiff, "GB/s."

! deallocate all arrays
call cgca_ds(space)

end program testABN

!*roboend*

max_stars_repo_name: lcebaman/casup

max_issues_repo_name: lcebaman/casup

max_forks_repo_name: lcebaman/casup


text:
      subroutine nodmas(y,z,nodm,xmass,yz,nodd,xdamp)
c     implicit double precision (a-h,o-z)                                    dp
c
c     read in concentrated nodal masses and dampers
c
      common/bk03/numdc,imassn,idampn,irller,penstf
      common/bk14/lfna(15),lfnt(6)
      dimension y(*),z(*),nodm(*),xmass(2,*),yz(2,*),nodd(*),xdamp(2,*)
      character*80 txts,mssg
c
c     n o d a l   m a s s e s
c
      if (imassn.eq.0) go to 20
c
      call header
      write(lfnt(2),50)
      do 10 in=1,imassn
      call gttxsg (txts,lcount)
      read(unit=txts,fmt=40,err=34) n,(xmass(i,in),i=1,2)
      yz(1,in)=y(n)
      yz(2,in)=z(n)
      nodm(in)=n
   10 write(lfnt(2),60) n,(xmass(i,in),i=1,2)
c
   20 if (idampn.eq.0) return
c
      call header
      write(lfnt(2),70)
      do 30 in=1,idampn
      call gttxsg (txts,lcount)
      read(unit=txts,fmt=40,err=36) n,(xdamp(i,in),i=1,2)
      nodd(in)=n
   30 write(lfnt(2),60) n,(xdamp(i,in),i=1,2)
c
      return
c
   34 write (unit=mssg,fmt=80) in
      call termin (txts,mssg,lcount,1)
   36 write (unit=mssg,fmt=90) in
      call termin (txts,mssg,lcount,1)
c
   40 format(i5,2e10.1)
   50 format(///' n o d a l   m a s s   d a t a   '//
     1          '   node     y-mass     z-mass')
   60 format(i7,6e11.3)
   70 format(///' n o d a l   d a m p e r   d a t a   '//
     1          '   node   y-damper   z-damper')
   80 format (' error reading nodal mass card #',i5)
   90 format (' error reading discrete damper card #',i5)
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
      subroutine gradctr(iflg)
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C***********************************************************************
CD1
CD1  PURPOSE
CD1
CD1  This subroutine manages initization and boundary conditions
CD1  that involve gradients
CD1
C***********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Initial implementation: 19-APR-02, Programmer: George Zyvoloski
CD2
CD2 $Log:   /pvcs.config/fehm90/src/gradctr.f_a  $
CD2
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:12   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:07:20   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
CD2
C***********************************************************************
CD3
CD3  REQUIREMENTS TRACEABILITY
CD3
CD3 2.6 Provide Input/Output Data Files
CD3 3.0 INPUT AND OUTPUT REQUIREMENTS
CD3
C***********************************************************************
CD4
CD4  SPECIAL COMMENTS AND REFERENCES
CD4
CD4 Requirements from SDN: 10086-RD-2.20-00
CD4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD4   FEHM Application Version 2.20
CD4
**********************************************************************

      use comai
      use combi
      use comco2
      use comdi
      use comdti
      use comfi
      use commeth
      implicit none

      integer iflg,icode, izone, inode, idir, iroot,neq_total
      integer mi,neqp1,i,i1,i2,j,jj,ja,ngrad_max,i_2nd
      integer open_file,igradm, imodel, jk, node_second_0
c gaz 052019
      integer, allocatable ::  izonef_tmp(:)
      character*80 gradm_name
      character*5 dgradm
      character*80 gradmod_root
      logical grad_dum
      character*3 grad_all
      parameter(ngrad_max= 50)
      real*8 dist,var_inode

c================================================================
      if(igrad.eq.0) return
c================================================================
      if(iflg.eq.0) then
c
c read input
c cordg : reference coordinate
c idirg : gradient coordinate direction (1 = x,2 = y, 3 = z)
c igradf=1 : initial pressure is a variable
c igradf=2 : initial temperature is a variable
c igradf=3 : initial saturation  is a variable
c igradf=4 : specified pressure (or head) is a variable
c igradf=5 : specified inflow temperture (or enthalpy) is a variable
c xg1 is x coordinate at reference point
c yg1 is y coordinate at reference point
c zg1 is z coordinate at reference point
c var0 is variable at cordg
c grad1 is the linear gradient
c
         if(.not.allocated(izone_grad)) then
            allocate(izone_grad(ngrad_max))
            allocate(igradf(ngrad_max))
            allocate(cordg(ngrad_max))
            allocate(var0(ngrad_max))
            allocate(grad1(ngrad_max))
            allocate(idirg(ngrad_max))
            allocate(gradmod_filename(ngrad_max))
            allocate(igradmodnamlen(ngrad_max))
            allocate(igradmodelfile(ngrad_max))
            igradmd = 0
            gradmod_filename = ''
         end if
         iroot = 8
         gradmod_root(1:8)='gradtemp'
         read(inpt,*) ngrad
        do i = 1,ngrad
         grad_dum = .false.
         read(inpt,'(a80)') wdd1
          do jj = 1,78
           if(wdd1(jj:jj+2).eq.'all') grad_dum = .true.
          enddo
         backspace inpt
         if(grad_dum) then
          read(inpt,*)
     &        grad_all,cordg(i),idirg(i),
     &        igradf(i),var0(i),grad1(i)
          izone_grad(i) = -1
         else
          read(inpt,*)
     &        izone_grad(i),cordg(i),idirg(i),
     &        igradf(i),var0(i),grad1(i)
          if(izone_grad(i).ne.-1) then
           izone_grad(i) = abs(izone_grad(i))
          endif
         endif
        enddo
c
c     save zone file and parameters
c
         gradm_name = ''
         igradmd = igradmd + 1
         idgradmc = 10000+igradmd
         write(dgradm,'(i5)')idgradmc
         gradm_name(1:iroot) = gradmod_root(1:iroot)
         gradm_name(iroot+1:iroot+1) ='.'
         gradm_name(iroot+2:iroot+5) = dgradm(2:5)
         gradm_name(iroot+6:iroot+11) = '.gradf'
         gradmod_filename(igradmd)(1:iroot+11)
     &   = gradm_name(1:iroot+11)
         igradmodnamlen(igradmd)= iroot+11
c
c complete name here
c
         igradmodelfile(igradmd) = open_file(gradm_name, 'unknown')
         j = igradmodelfile(igradmd)
         write(j,'(a4, 1x, i9)') 'grad',ngrad
         do i = 1, ngrad
            write(j,*)
     &           izone_grad(i),cordg(i),idirg(i),
     &           igradf(i),var0(i),grad1(i)
         enddo
         write(j,'(a4)') 'end '
         write(j,'(8(1x,i9))') (izonef(i),i=1,n0)
         close (j)
      else if(iflg.eq.1) then
c
c modify initial values and BC's
c
        do jj = 1,igradmd
c
c read gradient info from auxillary files
c and zone list for that request
c
         j = igradmodelfile(jj)
         i1=igradmodnamlen(jj)
         j = open_file (gradmod_filename(jj), 'old')
c         open(j,file=gradmod_filename(jj)(1:i1),
c     &    status ='unknown')

         read(j,*)dgradm(1:4) ,ngrad
         do i = 1, ngrad
            read(j,*)
     &           izone_grad(i),cordg(i),idirg(i),
     &           igradf(i),var0(i),grad1(i)
         enddo
c gaz 052019
         if(.not.allocated(izonef_tmp)) allocate(izonef_tmp(n0))
         read(j,'(a4)') dgradm(1:4)
         read(j,*) (izonef_tmp(i),i=1,n0)
c
         if(iread.le.0) then
c code with no restart file is present
          if(gdpm_flag.ge.3.and.gdpm_flag.le.6) then
           neq_total = neq_primary
          else
           neq_total = neq
          endif
            do izone=1,ngrad
               do inode=1,neq_total
                  if(izonef_tmp(inode).eq.izone_grad(izone)
     &                  .or.izone_grad(izone).eq.-1) then
                     idir = idirg(izone)
                     dist = cord(inode,idir)-cordg(izone)
                     if(igradf(izone).eq.1) then
                        pho(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.2) then
                        to(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.3) then
                        so(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.4) then
                        pflow(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.5) then
                        esk(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.-5) then
                        esk(inode)=-(var0(izone) + grad1(izone)*dist)
                     else if(igradf(izone).eq.6) then
                        phometh(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.7) then
                        pflowmeth(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.8) then
                        qflux(inode)=var0(izone) + grad1(izone)*dist
c RJP 04/10/07 added the following part for CO2
                     else if(igradf(izone).eq.9) then
                        phoco2(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.10) then
                        pflowco2(inode)=var0(izone) + grad1(izone)*dist
                     endif
                  endif
               enddo
c  calculate variables from gradients for gdpm and gdkm
             if(gdpm_flag.ge.3.and.gdpm_flag.le.6) then
               do inode=1,neq_primary
                 imodel = igdpm(inode)
                 if(ngdpm_layers(imodel).ne.0) then
                  if(izonef_tmp(inode).eq.izone_grad(izone))then
c notice that the gdpm nodes are 1D
                   node_second_0 = nelm(nelm(inode+1))
                   do jk = 1, ngdpm_layers(imodel)
                     i_2nd = node_second_0 + (jk-1)
                     idir = idirg(izone)
                     dist =  cord(inode,1) - cord(i_2nd,1)
                     if(igradf(izone).eq.11) then
                        var_inode = pho(inode)
                        pho(i_2nd) = var_inode + grad1(izone)*dist
                     else if(igradf(izone).eq.12) then
                        var_inode = to(inode)
                        to(i_2nd)= var_inode + grad1(izone)*dist
                     endif
                   enddo
                   endif
                  endif
               enddo
             endif
            enddo
         else
            do izone=1,ngrad
               do inode=1,n0
                  if(izonef_tmp(inode).eq.izone_grad(izone)) then
                     idir = idirg(izone)
                     dist = cord(inode,idir)-cordg(izone)
                     if(igradf(izone).eq.4) then
                        pflow(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.5) then
                        esk(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.-5) then
                        esk(inode)=-(var0(izone) + grad1(izone)*dist)
                     else if(igradf(izone).eq.7) then
                        pflowmeth(inode)=var0(izone) + grad1(izone)*dist
                     else if(igradf(izone).eq.8) then
                        qflux(inode)=var0(izone) + grad1(izone)*dist
c RJP 04/10/07 added following for CO2
                     else if(igradf(izone).eq.10) then
                        pflowco2(inode)=var0(izone) + grad1(izone)*dist
                     endif
                  endif
               enddo
            enddo
         endif
c        rewind j
c        write(j,*) 'gradctr file ', jj, ' read'
        close(j, status = 'delete')
       enddo
         if(allocated(izone_grad)) then
            deallocate(izone_grad)
            deallocate(izonef_tmp)
            deallocate(igradf)
            deallocate(var0)
            deallocate(grad1)
            deallocate(cordg)
            deallocate(idirg)
            deallocate(gradmod_filename)
            deallocate(igradmodnamlen)
            deallocate(igradmodelfile)
         end if

      else if(iflg.eq.2) then
c

      else if(iflg.eq.3) then
c

      endif
c
      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
! wrapfscope_ns1.f
! This file is generated by Shroud nowrite-version. Do not edit.
! Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
! other Shroud Project Developers.
! See the top-level COPYRIGHT file for details.
!
! SPDX-License-Identifier: (BSD-3-Clause)
!
!>
!! \file wrapfscope_ns1.f
!! \brief Shroud generated wrapper for ns1 namespace
!<
! splicer begin namespace.ns1.file_top
! splicer end namespace.ns1.file_top
module scope_ns1_mod
    use iso_c_binding, only : C_INT
    ! splicer begin namespace.ns1.module_use
    ! splicer end namespace.ns1.module_use
    implicit none

    ! splicer begin namespace.ns1.module_top
    ! splicer end namespace.ns1.module_top

    !  enum ns1::Color
    integer(C_INT), parameter :: red = 20
    integer(C_INT), parameter :: blue = 21
    integer(C_INT), parameter :: white = 22

    interface

        ! splicer begin namespace.ns1.additional_interfaces
        ! splicer end namespace.ns1.additional_interfaces
    end interface

contains

    ! splicer begin namespace.ns1.additional_functions
    ! splicer end namespace.ns1.additional_functions

end module scope_ns1_mod

max_stars_repo_name: ExternalRepositories/shroud

max_issues_repo_name: ExternalRepositories/shroud

max_forks_repo_name: ExternalRepositories/shroud