text:
.( load user32.dll dynamic library )
loadlibrary user32.dll ?exit            \ 'loadlibrary' is from Sir yap's demo in KsanaVM forth tutorial.
                                        \ user32.dll at this moment is a string of filename.
                                        \ after loadlibrary will successfully do its job, 'user32.dll' become a new command word.
                                        \ As a command word, 'user32.dll' leaves the pointer to DLL user32.dll.

.( get 'MessageBoxA' function from user32.dll )
user32.dll @                            \ Get 'user32.dll' address
4 0 winapi MessageBoxA ?exit            \ 'winapi' [ DLL-address n -- error ] gets an export function from the given DLL.
                                        \ in this example, 'winapi' gets the exported function 'MessageBoxA' from user32.dll
                                        \ and assign the save name 'MessageBoxA' to the function.
                                        \ Where n is 'MessageBoxA' function's input argument number. In this exampke, it needs
                                        \ 4 arguments as shown below C++ function prototype,
                                        \
                                        \    int WINAPI MessageBoxA(
                                        \      __in_opt  HWND hWnd,
                                        \      __in_opt  LPCTSTR lpText,
                                        \      __in_opt  LPCTSTR lpCaption,
                                        \      __in      UINT uType
                                        \    );
                                        \
                                        \ After 'winapi' done. 'MessageBoxA'

.( Test using it . . . )
: test 0 s" hello, world!" s" WinApi Test" 0 MessageBoxA ;    \ see, there are 4 arguments before calling MessageBoxA
test .ss                                                      \ MessageBoxA returns an error code let's check it out in the stack.

max_stars_repo_name: hcchengithub/weforth

max_issues_repo_name: hcchengithub/weforth

max_forks_repo_name: hcchengithub/weforth


text:
C********************************************************
C
C   TEST PROGRAM FOR FORTRAN AND C CALLABLE SUBROUTINE
C                    PUTSPICE95
C
C  Modified 8/6/99 by Michael Brady
C
C  Added the variable "proj_name".  Previously, the
C  integer "project" was used as a parameter for
C  getlabcon and get proj, causing them to fail.
C
C********************************************************
C
	include 'VICMAIN_FOR'
	subroutine main44

        implicit none

        character*6 proj_name
	integer*4 mode
	character*80 imgfile
	character*10 source
	integer*4 project
	integer*4 system
	character*10 usr_info(10)
	character*5 camera
        character*12 target_name
        character*12 temp
	double precision dbuf(100)
	integer*4 ibuf(200)
	integer*4 ind, count, def
	equivalence (dbuf, ibuf)
	integer*4 scet(6)
	integer*4 data(80)
        integer*4 unit, status, cam, fds

C	Selecting mode to tell putspice95() where to
C	write the data
	print *, 'Input mode (local = 0, remote = 1)'
	read *, mode

	call xvunit(unit, 'INP', 1, status, ' ')
	call xvopen(unit, status, 'OPEN_ACT', 'SA', ' ')
	call getproj(unit, proj_name, cam, fds, ind)
	call getlabcon(unit, proj_name, data, ind)

        call mve(4, 6, data(8), scet, 1, 1)
        call mvlc(data(25), target_name, 12)

        call xvparm('INP', imgfile, count, def, 0)
        call xvparm('SOURCE', source, count, def, 0)
	call xvparm('PROJECT', project, count, def, 0)
        call xvparm('CAMERA', camera, count, def, 0)
        call xvparm('TARGET', temp, count, def, 0)

        if (def .EQ. 0) then
           target_name=temp
        endif

	call xvparm('INSTITUTE', usr_info(1), count, def, 0)
	call xvparm('PURPOSE', usr_info(2), count, def, 0)
	call xvparm('PROG_NAME', usr_info(3), count, def, 0)
	call xvparm('SP_REF', usr_info(4), count, def, 0)
	call xvparm('REQ_NO', usr_info(5), count, def, 0)
	call xvparm('YEAR', usr_info(6), count, def, 0)
	call xvparm('MONTH_DAY', usr_info(7), count, def, 0)
	call xvparm('HOUR_MIN', usr_info(8), count, def, 0)
	call xvparm('FILE_ID', usr_info(9), count, def, 0)
	call xvparm('USR_GRP_ID', usr_info(10), count, def, 0)

	if (project .eq. -77) then
	   system = 2
	else
	   system = 1
	endif

        call xvmessage
     &     ('............................................', ' ')
        call xvmessage(' ',' ')
        call xvmessage('..................', ' ')
        call xvmessage('Calling GetSpice95', ' ')
	call getspice95(mode, project, camera, scet, target_name,
     &        system, source, usr_info, dbuf, ind)
        call xvmessage(' ',' ')

	if (ind .eq. 0) then
	   print*, 'ZGETSPICE95::Test From FORTRAN:FAILURE'
	else
           call xvparm('SOURCE', source, count, def, 0)
	   call mvcl(source, ibuf(11), 4)

        call xvmessage
     &     ('...... Testing PUTSPICE95 In Fortran .......', ' ')
	   ind  = 0
	   call putspice95(project, dbuf, mode, ind)
	   if (ind .eq. 0) then
	      print*, 'PUTSPICE95_ERROR:: Fortran Test FAILED'
	   else
              print*, 'PUTSPICE95:: Fortran Test SUCCESS'
	   endif
       endif

        call xvmessage('           ', ' ')
        call xvmessage
     &     ('...... Testing PUTSPICE95 In C .......', ' ')
        ind  = 0
        call tzputspice95(project, dbuf, mode, ind)
        if (ind .eq. 0) then
           print*, 'PUTSPICE_ERROR:: C Test FAILED'
        else
           print*, 'PUTSPICE95::C Test SUCCESS'
        endif

	return
	end

ext: f

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
module string_converter_module
    use, intrinsic :: iso_fortran_env, only : error_unit, real64, int64
    use :: string_module, only : string
    implicit none
    private

    public :: string_converter

    type :: string_converter
    contains
        procedure, nopass :: toint => toint
        procedure, nopass :: toint64 => toint64
        procedure, nopass :: to_logical => to_logical
        procedure, nopass :: toreal32 => toreal32
        procedure, nopass :: toreal64 => toreal64
        procedure, nopass :: to_string_array => to_string_array
        procedure, nopass :: to_int_array => to_int_array
        procedure, nopass :: to_int64_array => to_int64_array
        procedure, nopass :: fromint => fromint
    end type string_converter

    character, parameter :: default_array_delimiter = ' '
contains
    integer function toint(str)
        type(string), intent(in) :: str

        integer :: error

        read(str%char_array, *, iostat=error) toint
        if ( error /= 0) &
            error stop "string_converter::Error converting string to integer"
    end function toint

    integer(int64) function toint64(str)
        type(string), intent(in) :: str

        integer :: error

        read(str%char_array, *, iostat=error) toint64
        if ( error /= 0) &
            error stop "string_converter::Error converting string to integer64"
    end function toint64

    type(string) function fromint(input)
        integer, intent(in) :: input

        character(len=20) :: dummy

        write(dummy, *) input
        fromint = string(trim(adjustl(dummy)))
    end function fromint

    logical function to_logical(str)
        type(string), intent(in) :: str

        to_logical = .false.
        if ( str%to_lower() == "true" ) to_logical = .true.
        if ( str%to_lower()== "yes" ) to_logical = .true.
    end function to_logical

    function to_string_array(str, delimiter) result(array)
        type(string), intent(in) :: str
        character, intent(in), optional :: delimiter
        type(string), dimension(:), allocatable :: array

        character :: actual_delimiter

        actual_delimiter = default_array_delimiter
        if (present(delimiter)) actual_delimiter = delimiter

        array = str%split_and_strip(actual_delimiter)
    end function to_string_array

    function to_int_array(str, delimiter) result(array)
        type(string), intent(in) :: str
        character, intent(in), optional :: delimiter
        integer, dimension(:), allocatable :: array

        type(string), dimension(:), allocatable :: str_array
        integer :: idx

        str_array = to_string_array(str, delimiter)

        array = [(toint(str_array(idx)), idx = 1, size(str_array))]

        do idx = 1, size(str_array)
            call str_array(idx)%cleanup()
        end do
        deallocate(str_array)
    end function to_int_array

    function to_int64_array(str, delimiter) result(array)
        type(string), intent(in) :: str
        character, intent(in), optional :: delimiter
        integer(int64), dimension(:), allocatable :: array

        type(string), dimension(:), allocatable :: str_array
        integer :: idx

        str_array = to_string_array(str, delimiter)

        array = [(toint64(str_array(idx)), idx = 1, size(str_array))]

        do idx = 1, size(str_array)
            call str_array(idx)%cleanup()
        end do
        deallocate(str_array)
    end function to_int64_array

    real function toreal32(str)
        type(string), intent(in) :: str

        integer :: error

        read(str%char_array, *, iostat=error) toreal32
        if ( error /= 0) &
            error stop "string_converter::Error converting string to real32"
    end function toreal32

    real(real64) function toreal64(str)
        type(string), intent(in) :: str

        integer :: error

        read(str%char_array, *, iostat=error) toreal64
        if ( error /= 0) &
            error stop "string_converter::Error converting string to real64"
    end function toreal64
end module string_converter_module

max_stars_repo_name: cheshyre/ntcl-util

max_issues_repo_name: cheshyre/ntcl-util

max_forks_repo_name: cheshyre/ntcl-util


text:
C***********************************************************************
C    Module:  gridtest.f
C
C    Copyright (C) 1996 Harold Youngren, Mark Drela
C
C    This program is free software; you can redistribute it and/or modify
C    it under the terms of the GNU General Public License as published by
C    the Free Software Foundation; either version 2 of the License, or
C    (at your option) any later version.
C
C    This program is distributed in the hope that it will be useful,
C    but WITHOUT ANY WARRANTY; without even the implied warranty of
C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C    GNU General Public License for more details.
C
C    You should have received a copy of the GNU General Public License
C    along with this program; if not, write to the Free Software
C    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
C
C    Report problems to:    guppy@maine.com
C                        or drela@mit.edu
C***********************************************************************

C Example/Test of grid routine
C  Sets up two plots, each containing a grid with a label and symbol line
C  First  plot comes up as B&W in portrait window
C  Second plot comes up as color in larger landscape window
C
      INCLUDE 'colors.inc'
      CHARACTER INP*10
      DATA LMASK1, LMASK2 / -32640, -30584 /
      DATA LMASK3, LMASK4/  -21846, Z'AAAAAAAA'/
C
C---- number of grid intervals per axis annotation interval
      NGR = 2
C
      XMIN=.1
      YMIN=.1
      XMAX=.9
      YMAX=.6
      XDEL = 0.1
      YDEL = 0.1
      XSF = 1.0
      YSF = 1.0
      NXG = NGR * INT((XMAX-XMIN)/XDEL + 0.001)
      NYG = NGR * INT((YMAX-YMIN)/YDEL + 0.001)
      DXG = XSF*XDEL / FLOAT(NGR)
      DYG = YSF*YDEL / FLOAT(NGR)
C
      CH = 0.03
      SH = 0.2*CH
C
C---What devices to plot to?
      WRITE(*,*) ' '
 1    WRITE(*,*) ' Enter -1 for no PS, 0 for B&W PS, 1 for color PS'
      READ(*,1000,end=2000) INP
      IOUT = -1
      if(INP.NE.' ') READ(INP,*,end=2000,err=2000) IOUT
 2    IF(IOUT.LT.0) IDEV = 1
      IF(IOUT.EQ.0) IDEV = 3
      IF(IOUT.GE.1) IDEV = 5
C
C---plot #1 in B&W portrait mode 0.5 of screen height
      CALL PLINITIALIZE
      CALL PLOPEN(-0.5,0,IDEV)
      CALL PLOT(0.1,0.1,-3)
      CALL NEWPEN(1)
      CALL NEWFACTOR(10.)
      CALL PLCHAR(0.2,0.2,1.2*CH,'TEST FOR GRID',0.0,13)
      CALL PLSYMB(999.,0.2,CH,1,0.0,0)
      CALL PLSYMB(0.2,0.2,CH,1,0.0,-1)
      CALL PLGRID(0.0,0.0, NXG,DXG, NYG,DYG, LMASK4)
      CALL PLFLUSH
      PAUSE
      CALL PLOT(0.,0.,-999)
C
C---plot #2 in landscape mode 0.7 of screen width
C   (green grid with red lettering with blue symbol line)
      CALL PLOPEN(0.7,0,IDEV)
      CALL PLOT(0.1,0.1,-3)
      CALL NEWFACTOR(10.)
      CALL NEWPEN(1)
      NXG = NGR * INT((XMAX-XMIN)/XDEL + 0.001)
      NYG = NGR * INT((YMAX-YMIN)/YDEL + 0.001)
      DXG = XSF*XDEL / FLOAT(NGR)
      DYG = YSF*YDEL / FLOAT(NGR)

      call NEWCOLORNAME('green')
      CALL PLGRID(0.0,0.0, NXG,DXG, NYG,DYG, LMASK2)

      CALL NEWPEN(1)
      CALL NEWCOLORRGB(255,0,0)
      CALL PLCHAR(0.2,0.2,1.2*CH,'TEST FOR GRID',0.0,13)
      call NEWCOLOR(BLUE)
      call NEWPEN(1)
      CALL PLSYMB(999.,0.2,CH,1,0.0,0)
      CALL PLSYMB(0.2,0.2,CH,1,0.0,-1)
      CALL PLFLUSH
      PAUSE
      CALL PLCLOSE
C
 1000 FORMAT(A)
 2000 STOP
      END


text:
! Bug report from Jean (ANL) nameListDecl_2.f90
module m_test2011_70

   implicit none
   save

   logical  ::  &
      iage_on, dye_on

   namelist /passive_tracers_on_nml/  &
      iage_on, dye_on
end module

program p
  use m_test2011_70
  print *,'OK' ! parser test for namelist only
end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


! Code automatically generated by 'PreSARAH'
! Expressions for amplitudes are obtained by FeynArts/FormCalc
! Created at 11:55 on 9.11.2016

Module DecaySSS

Use Control
Use LoopFunctions
Use AddLoopFunctions
Use Settings

Contains

Subroutine TreeAmp_StoSS(Mex1,Mex2,Mex3,coupT1,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Complex(dp), Intent(in) :: coupT1
! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= 1*(0.,1._dp)*coupT1

End Subroutine TreeAmp_StoSS

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

Subroutine Amp_VERTEX_StoSS_Topology1_FFF(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
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

! 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("C1  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int4=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 {}[[1]]
Amplitude= &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3L*int1*ML1) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3R*int1*ML1) &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3L*int2*Mex1**2*ML1) &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3R*int2*Mex1**2*ML1) &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3L*int3*Mex1**2*ML1) &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3R*int3*Mex1**2*ML1) &
& +(3*(0.,1._dp)*coup1R*coup2R*coup3L*int4*Mex1**2*ML1) &
& +(3*(0.,1._dp)*coup1L*coup2L*coup3R*int4*Mex1**2*ML1) &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3L*int2*Mex2**2*ML1) &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3R*int2*Mex2**2*ML1) &
& +(3*(0.,1._dp)*coup1R*coup2R*coup3L*int3*Mex2**2*ML1) &
& +(3*(0.,1._dp)*coup1L*coup2L*coup3R*int3*Mex2**2*ML1) &
& +(1*(0.,1._dp)*coup1R*coup2R*coup3L*int4*Mex2**2*ML1) &
& +(1*(0.,1._dp)*coup1L*coup2L*coup3R*int4*Mex2**2*ML1) &
& +(-1*(0.,1._dp)*coup1R*coup2R*coup3L*int2*Mex3**2*ML1) &
& +(-1*(0.,1._dp)*coup1L*coup2L*coup3R*int2*Mex3**2*ML1) &
& +(-1*(0.,1._dp)*coup1R*coup2R*coup3L*int3*Mex3**2*ML1) &
& +(-1*(0.,1._dp)*coup1L*coup2L*coup3R*int3*Mex3**2*ML1) &
& +(-1*(0.,1._dp)*coup1R*coup2R*coup3L*int4*Mex3**2*ML1) &
& +(-1*(0.,1._dp)*coup1L*coup2L*coup3R*int4*Mex3**2*ML1) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3L*int2*ML1**3) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3R*int2*ML1**3) &
& +(2*(0.,1._dp)*coup1L*coup2R*coup3L*int1*ML2) &
& +(2*(0.,1._dp)*coup1R*coup2L*coup3R*int1*ML2) &
& +(1*(0.,1._dp)*coup1L*coup2R*coup3L*int4*Mex1**2*ML2) &
& +(1*(0.,1._dp)*coup1R*coup2L*coup3R*int4*Mex1**2*ML2) &
& +(2*(0.,1._dp)*coup1L*coup2R*coup3L*int3*Mex2**2*ML2) &
& +(2*(0.,1._dp)*coup1R*coup2L*coup3R*int3*Mex2**2*ML2) &
& +(1*(0.,1._dp)*coup1L*coup2R*coup3L*int4*Mex2**2*ML2) &
& +(1*(0.,1._dp)*coup1R*coup2L*coup3R*int4*Mex2**2*ML2) &
& +(-1*(0.,1._dp)*coup1L*coup2R*coup3L*int4*Mex3**2*ML2) &
& +(-1*(0.,1._dp)*coup1R*coup2L*coup3R*int4*Mex3**2*ML2) &
& +(2*(0.,1._dp)*coup1L*coup2R*coup3L*int2*ML1**2*ML2) &
& +(2*(0.,1._dp)*coup1R*coup2L*coup3R*int2*ML1**2*ML2) &
& +(2*(0.,1._dp)*coup1R*coup2L*coup3L*int1*ML3) &
& +(2*(0.,1._dp)*coup1L*coup2R*coup3R*int1*ML3) &
& +(1*(0.,1._dp)*coup1R*coup2L*coup3L*int3*Mex1**2*ML3) &
& +(1*(0.,1._dp)*coup1L*coup2R*coup3R*int3*Mex1**2*ML3) &
& +(2*(0.,1._dp)*coup1R*coup2L*coup3L*int4*Mex1**2*ML3) &
& +(2*(0.,1._dp)*coup1L*coup2R*coup3R*int4*Mex1**2*ML3) &
& +(1*(0.,1._dp)*coup1R*coup2L*coup3L*int3*Mex2**2*ML3) &
& +(1*(0.,1._dp)*coup1L*coup2R*coup3R*int3*Mex2**2*ML3) &
& +(-1*(0.,1._dp)*coup1R*coup2L*coup3L*int3*Mex3**2*ML3) &
& +(-1*(0.,1._dp)*coup1L*coup2R*coup3R*int3*Mex3**2*ML3) &
& +(2*(0.,1._dp)*coup1R*coup2L*coup3L*int2*ML1**2*ML3) &
& +(2*(0.,1._dp)*coup1L*coup2R*coup3R*int2*ML1**2*ML3) &
& +(2*(0.,1._dp)*coup1L*coup2L*coup3L*int2*ML1*ML2*ML3) &
& +(2*(0.,1._dp)*coup1R*coup2R*coup3R*int2*ML1*ML2*ML3)

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

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

Subroutine Amp_VERTEX_StoSS_Topology1_SSS(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1, coup2, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1, coup2, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

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

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= -1*(0.,1._dp)*coup1*coup2*coup3*int1

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

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

Subroutine Amp_VERTEX_StoSS_Topology1_UUU(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1, coup2, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1, coup2, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

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

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= 1*(0.,1._dp)*coup1*coup2*coup3*int1

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

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

Subroutine Amp_VERTEX_StoSS_Topology1_VVV(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1, coup2, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1, coup2, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

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

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= 4*(0.,1._dp)*coup1*coup2*coup3*int1

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

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

Subroutine Amp_VERTEX_StoSS_Topology1_SSV(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1, coup2, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1, coup2, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! 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("C1  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int4=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 {}[[1]]
Amplitude= &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int1) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int2*Mex1**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int3*Mex1**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int4*Mex1**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int3*Mex2**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int4*Mex2**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int2*Mex3**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int3*Mex3**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int4*Mex3**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int2*ML1**2)

Amplitude=-Amplitude

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

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

Subroutine Amp_VERTEX_StoSS_Topology1_SVS(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1, coup2, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1, coup2, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! Necessary Loop Integrals
int1=SA_B0(Mex3**2,ML3**2,ML2**2)
int2=SA_Cget("C0  ",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("C2  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int1) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int3*Mex1**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int4*Mex1**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int2*Mex2**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int3*Mex2**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int4*Mex2**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int2*Mex3**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int3*Mex3**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int4*Mex3**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int2*ML1**2)

Amplitude=-Amplitude

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

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

Subroutine Amp_VERTEX_StoSS_Topology1_VSS(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1, coup2, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1, coup2, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! 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("C1  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int4=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 {}[[1]]
Amplitude= &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int1) &
& +(-2*(0.,1._dp)*coup1*coup2*coup3*int2*Mex1**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int3*Mex1**2) &
& +(-3*(0.,1._dp)*coup1*coup2*coup3*int4*Mex1**2) &
& +(-2*(0.,1._dp)*coup1*coup2*coup3*int2*Mex2**2) &
& +(-3*(0.,1._dp)*coup1*coup2*coup3*int3*Mex2**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int4*Mex2**2) &
& +(2*(0.,1._dp)*coup1*coup2*coup3*int2*Mex3**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int3*Mex3**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int4*Mex3**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int2*ML1**2)

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

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

Subroutine Amp_VERTEX_StoSS_Topology1_SVV(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1, coup2, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1, coup2, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! 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("C1  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int4=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 {}[[1]]
Amplitude= &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int1) &
& +(1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int2*Mex1**2) &
& +(-1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int3*Mex1**2) &
& +(-3/2._dp*(0.,1._dp)*coup1*coup2*coup3*int4*Mex1**2) &
& +(1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int2*Mex2**2) &
& +(-3/2._dp*(0.,1._dp)*coup1*coup2*coup3*int3*Mex2**2) &
& +(-1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int4*Mex2**2) &
& +(-1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int2*Mex3**2) &
& +(1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int3*Mex3**2) &
& +(1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int4*Mex3**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int2*ML1**2)

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

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

Subroutine Amp_VERTEX_StoSS_Topology1_VSV(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1, coup2, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1, coup2, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! 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("C1  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)
int4=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 {}[[1]]
Amplitude= &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int1) &
& +(3*(0.,1._dp)*coup1*coup2*coup3*int2*Mex1**2) &
& +(2*(0.,1._dp)*coup1*coup2*coup3*int3*Mex1**2) &
& +(7/2._dp*(0.,1._dp)*coup1*coup2*coup3*int4*Mex1**2) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int2*Mex2**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int3*Mex2**2) &
& +(-1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int4*Mex2**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int2*Mex3**2) &
& +(-2*(0.,1._dp)*coup1*coup2*coup3*int3*Mex3**2) &
& +(1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int4*Mex3**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int2*ML1**2)

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

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

Subroutine Amp_VERTEX_StoSS_Topology1_VVS(Mex1,Mex2,Mex3,ML1, ML2, ML3,coup1, coup2, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2, ML3
Complex(dp), Intent(in) :: coup1, coup2, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! Necessary Loop Integrals
int1=SA_B0(Mex3**2,ML3**2,ML2**2)
int2=SA_Cget("C0  ",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("C2  ",Mex2**2,Mex3**2,Mex1**2,ML1**2,ML3**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int1) &
& +(-1*(0.,1._dp)*coup1*coup2*coup3*int2*Mex1**2) &
& +(-1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int3*Mex1**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int4*Mex1**2) &
& +(3*(0.,1._dp)*coup1*coup2*coup3*int2*Mex2**2) &
& +(7/2._dp*(0.,1._dp)*coup1*coup2*coup3*int3*Mex2**2) &
& +(2*(0.,1._dp)*coup1*coup2*coup3*int4*Mex2**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int2*Mex3**2) &
& +(1/2._dp*(0.,1._dp)*coup1*coup2*coup3*int3*Mex3**2) &
& +(-2*(0.,1._dp)*coup1*coup2*coup3*int4*Mex3**2) &
& +(1*(0.,1._dp)*coup1*coup2*coup3*int2*ML1**2)

Amplitude=-Amplitude

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

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

! --------------------------------------
! Topology Nr.2
! --------------------------------------

Subroutine Amp_VERTEX_StoSS_Topology2_SS(Mex1,Mex2,Mex3,ML1, ML2,coup1, coup2,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2
Complex(dp), Intent(in) :: coup1, coup2
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! Necessary Loop Integrals
int1=SA_B0(Mex1**2,ML1**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= 1/2._dp*(0.,1._dp)*coup1*coup2*int1

Amplitude=2._dp*Amplitude

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

!Field Order
! FieldOrder[2,{S, S},1]=Null[[1]];
! FieldOrder[2,{S, S},2]=Null[[2]];
! FieldOrder[2,{S, S},3]=Null[[3]];
End Subroutine Amp_VERTEX_StoSS_Topology2_SS

Subroutine Amp_VERTEX_StoSS_Topology2_VV(Mex1,Mex2,Mex3,ML1, ML2,coup1, coup2,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2
Complex(dp), Intent(in) :: coup1, coup2
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! Necessary Loop Integrals
int1=SA_B0(Mex1**2,ML1**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= &
& +(-1*(0.,1._dp)*coup1*coup2*rMS) &
& +(2*(0.,1._dp)*coup1*coup2*int1)

Amplitude=2._dp*Amplitude

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

!Field Order
! FieldOrder[2,{V, V},1]=Null[[1]];
! FieldOrder[2,{V, V},2]=Null[[2]];
! FieldOrder[2,{V, V},3]=Null[[3]];
End Subroutine Amp_VERTEX_StoSS_Topology2_VV

! --------------------------------------
! Topology Nr.3
! --------------------------------------

Subroutine Amp_VERTEX_StoSS_Topology3_SS(Mex1,Mex2,Mex3,ML1, ML2,coup1, coup2,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2
Complex(dp), Intent(in) :: coup1, coup2
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! Necessary Loop Integrals
int1=SA_B0(Mex2**2,ML1**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= 1/2._dp*(0.,1._dp)*coup1*coup2*int1

Amplitude=2._dp*Amplitude

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

!Field Order
! FieldOrder[3,{S, S},1]=Null[[1]];
! FieldOrder[3,{S, S},2]=Null[[2]];
! FieldOrder[3,{S, S},3]=Null[[3]];
End Subroutine Amp_VERTEX_StoSS_Topology3_SS

Subroutine Amp_VERTEX_StoSS_Topology3_VV(Mex1,Mex2,Mex3,ML1, ML2,coup1, coup2,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2
Complex(dp), Intent(in) :: coup1, coup2
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! Necessary Loop Integrals
int1=SA_B0(Mex2**2,ML1**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= &
& +(-1*(0.,1._dp)*coup1*coup2*rMS) &
& +(2*(0.,1._dp)*coup1*coup2*int1)

Amplitude=2._dp*Amplitude

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

!Field Order
! FieldOrder[3,{V, V},1]=Null[[1]];
! FieldOrder[3,{V, V},2]=Null[[2]];
! FieldOrder[3,{V, V},3]=Null[[3]];
End Subroutine Amp_VERTEX_StoSS_Topology3_VV

! --------------------------------------
! Topology Nr.4
! --------------------------------------

Subroutine Amp_VERTEX_StoSS_Topology4_SS(Mex1,Mex2,Mex3,ML1, ML2,coup1, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2
Complex(dp), Intent(in) :: coup1, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! Necessary Loop Integrals
int1=SA_B0(Mex3**2,ML1**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= 1/2._dp*(0.,1._dp)*coup1*coup3*int1
Amplitude=2._dp*Amplitude

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

!Field Order
! FieldOrder[4,{S, S},1]=Null[[1]];
! FieldOrder[4,{S, S},2]=Null[[2]];
! FieldOrder[4,{S, S},3]=Null[[3]];
End Subroutine Amp_VERTEX_StoSS_Topology4_SS

Subroutine Amp_VERTEX_StoSS_Topology4_VV(Mex1,Mex2,Mex3,ML1, ML2,coup1, coup3,Amplitude)
Implicit None
Real(dp), Intent(in) :: Mex1,Mex2,Mex3
Complex(dp), Intent(out) :: Amplitude
Real(dp), Intent(in) :: ML1, ML2
Complex(dp), Intent(in) :: coup1, coup3
Complex(dp) :: int1,int2,int3,int4,int5,int6,int7,int8,int9,int10

! Necessary Loop Integrals
int1=SA_B0(Mex3**2,ML1**2,ML2**2)

! Now the full amplitude for the different Lorentz Parts
! Coefficient of {}[[1]]
Amplitude= &
& +(-1*(0.,1._dp)*coup1*coup3*rMS) &
& +(2*(0.,1._dp)*coup1*coup3*int1)

Amplitude=2*Amplitude
! Amplitude=0

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

!Field Order
! FieldOrder[4,{V, V},1]=Null[[1]];
! FieldOrder[4,{V, V},2]=Null[[2]];
! FieldOrder[4,{V, V},3]=Null[[3]];
End Subroutine Amp_VERTEX_StoSS_Topology4_VV

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

AmpSquared = 0._dp
AmpSquared = Amplitude*Conjg(AmplitudeC)
End Subroutine SquareAmp_StoSS
End Module DecaySSS


text:
MODULE hwcart_mod
  use iso_c_binding
  use mpi

  implicit none

  enum, bind(C)
     enumerator :: HWCART_MD_HWTHREAD = 0
     enumerator :: HWCART_MD_CORE
     enumerator :: HWCART_MD_L1CACHE
     enumerator :: HWCART_MD_L2CACHE
     enumerator :: HWCART_MD_L3CACHE
     enumerator :: HWCART_MD_SOCKET
     enumerator :: HWCART_MD_NUMA
     enumerator :: HWCART_MD_NODE
  end enum

  enum, bind(C)
     enumerator :: HWCartOrderXYZ = 0
     enumerator :: HWCartOrderXZY
     enumerator :: HWCartOrderZYX
     enumerator :: HWCartOrderYZX
     enumerator :: HWCartOrderZXY
     enumerator :: HWCartOrderYXZ
  end enum

  type, bind(c) :: hwcart_topo_t
     type(c_ptr) :: ptr = c_null_ptr
  end type hwcart_topo_t

  interface
     function hwcart_init(hwcart_topo) bind(C)
       use iso_c_binding
       import hwcart_topo_t
       integer :: hwcart_init
       type(hwcart_topo_t), intent(out) :: hwcart_topo
     end function hwcart_init

     subroutine hwcart_topo_free(hwcart_topo) bind(C)
       use iso_c_binding
       import hwcart_topo_t
       type(hwcart_topo_t) :: hwcart_topo
     end subroutine hwcart_topo_free

     subroutine hwcart_comm_free(hwcart_comm) bind(C, name="hwcart_comm_free_f")
       use iso_c_binding
       integer :: hwcart_comm
     end subroutine hwcart_comm_free

     function hwcart_create_f(hwcart_topo, mpi_comm, nlevels, domain, topo, cart_order, hwcart_comm_out) bind(C)
       use iso_c_binding
       import hwcart_topo_t
       type(hwcart_topo_t), value :: hwcart_topo
       integer, value :: mpi_comm
       integer, value :: nlevels
       type(c_ptr), value :: domain
       type(c_ptr), value :: topo
       integer, value :: cart_order
       integer :: hwcart_comm_out
       integer :: hwcart_create_f
     end function hwcart_create_f

     function hwcart2mpicart_f(hwcart_comm, nlevels, topo, periodic, cart_order, mpicart_comm_out) bind(C)
       use iso_c_binding
       integer, value :: hwcart_comm
       integer, value :: nlevels
       type(c_ptr), value :: topo
       type(c_ptr), value :: periodic
       integer, value :: cart_order
       integer :: mpicart_comm_out
       integer :: hwcart2mpicart_f
     end function hwcart2mpicart_f

     function hwcart_sub_f(hwcart_comm, dims, rank, cart_order, belongs, hwcart_comm_out) bind(C)
       use iso_c_binding
       import hwcart_topo_t
       integer, value :: hwcart_comm
       type(c_ptr), value :: dims
       integer, value :: rank
       integer, value :: cart_order
       type(c_ptr), value :: belongs
       integer :: hwcart_comm_out
       integer :: hwcart_sub_f
     end function hwcart_sub_f

     function hwcart_print_rank_topology_f(hwcart_topo, mpi_comm, nlevels, domain, topo, cart_order) bind(C)
       use iso_c_binding
       import hwcart_topo_t
       type(hwcart_topo_t), value :: hwcart_topo
       integer, value :: mpi_comm
       integer, value :: nlevels
       type(c_ptr), value :: domain
       type(c_ptr), value :: topo
       integer, value :: cart_order
       integer :: hwcart_print_rank_topology_f
     end function hwcart_print_rank_topology_f

     function hwcart_rank2coord_f(hwcart_comm, dims, rank, cart_order, coord_out) bind(C)
       use iso_c_binding
       integer, value :: hwcart_comm
       type(c_ptr), value :: dims
       integer, value :: rank
       integer, value :: cart_order
       type(c_ptr), value :: coord_out
       integer :: hwcart_rank2coord_f
     end function hwcart_rank2coord_f

     function hwcart_coord2rank_f(hwcart_comm, dims, periodic, coord, cart_order, rank_out) bind(C)
       use iso_c_binding
       integer, value :: hwcart_comm
       type(c_ptr), value :: dims
       type(c_ptr), value :: periodic
       type(c_ptr), value :: coord
       integer, value :: cart_order
       integer :: rank_out
       integer :: hwcart_coord2rank_f
     end function hwcart_coord2rank_f

  end interface
CONTAINS

  function hwcart_create(hwcart_topo, mpi_comm, domain, topo, cart_order, hwcart_comm_out)
    type(hwcart_topo_t), value :: hwcart_topo
    integer, value :: mpi_comm
    integer, dimension(:), target :: domain
    integer, dimension(:,:), target :: topo
    integer, value :: cart_order
    integer :: hwcart_comm_out
    integer :: hwcart_create
    integer :: nlevels

    nlevels = size(domain)
    hwcart_create = &
      hwcart_create_f(hwcart_topo, mpi_comm, nlevels, c_loc(domain), c_loc(topo), cart_order, hwcart_comm_out)
  end function hwcart_create

  function hwcart2mpicart(hwcart_comm, topo, periodic, cart_order, mpicart_comm_out)
    integer, value :: hwcart_comm
    integer, dimension(:,:), target :: topo
    logical, dimension(:), target :: periodic(3)
    integer, value :: cart_order
    integer :: nlevels
    integer :: mpicart_comm_out
    integer :: hwcart2mpicart

    nlevels = size(topo)/3
    hwcart2mpicart = &
      hwcart2mpicart_f(hwcart_comm, nlevels, c_loc(topo), c_loc(periodic), cart_order, mpicart_comm_out)
  end function hwcart2mpicart

  function hwcart_sub(hwcart_comm, dims, rank, cart_order, belongs, hwcart_comm_out)
    integer, value :: hwcart_comm
    integer, dimension(:), target :: dims(3)
    integer, value :: rank
    integer, value :: cart_order
    logical, dimension(:), target :: belongs(3)
    integer :: hwcart_comm_out
    integer :: hwcart_sub

    hwcart_sub = &
      hwcart_sub_f(hwcart_comm, c_loc(dims), rank, cart_order, c_loc(belongs), hwcart_comm_out)
  end function hwcart_sub

  function hwcart_print_rank_topology(hwcart_topo, mpi_comm, domain, topo, cart_order)
    type(hwcart_topo_t), value :: hwcart_topo
    integer, value :: mpi_comm
    integer, dimension(:), target :: domain
    integer, dimension(:,:), target :: topo
    integer, value :: cart_order
    integer :: hwcart_comm_out
    integer :: hwcart_print_rank_topology
    integer :: nlevels

    nlevels = size(domain)
    hwcart_print_rank_topology = &
      hwcart_print_rank_topology_f(hwcart_topo, mpi_comm, nlevels, c_loc(domain), c_loc(topo), cart_order)
  end function hwcart_print_rank_topology

  function hwcart_rank2coord(hwcart_comm, dims, rank, cart_order, coord_out)
    use iso_c_binding
    integer, value :: hwcart_comm
    integer, dimension(:), target :: dims(3)
    integer, value :: rank
    integer, value :: cart_order
    integer, dimension(:), target :: coord_out(3)
    integer :: hwcart_rank2coord

    hwcart_rank2coord = &
      hwcart_rank2coord_f(hwcart_comm, c_loc(dims), rank, cart_order, c_loc(coord_out))
  end function hwcart_rank2coord

  function hwcart_coord2rank(hwcart_comm, dims, periodic, coord, cart_order, rank_out)
    use iso_c_binding
    integer, value :: hwcart_comm
    integer, dimension(:), target :: dims(3)
    logical, dimension(:), target :: periodic(3)
    integer, dimension(:), target :: coord(3)
    integer, value :: cart_order
    integer :: rank_out
    integer :: hwcart_coord2rank

    hwcart_coord2rank = &
      hwcart_coord2rank_f(hwcart_comm, c_loc(dims), c_loc(periodic), c_loc(coord), cart_order, rank_out)
  end function hwcart_coord2rank

END MODULE hwcart_mod

max_stars_repo_name: ghex-org/hwcart

max_issues_repo_name: ghex-org/hwcart

max_forks_repo_name: ghex-org/hwcart


text:
!===============================================================================
! One of Andy Nowacki's Fortran utility modules for dealing with seismic
! anisotropy and other problems.
!
! Andy Nowacki <andy.nowacki@bristol.ac.uk>
!
! See the file LICENCE for licence details.
!===============================================================================
!
!  spherical splines contains a few functions for fitting and evaluating
!  spherical splines, as given by
!     Wang & Dahlen (1995)  Spherical-spline parameterization of three-dimensional
!     Earth models.  GRL (22) 3099-3102.
!
!  This module uses LAPACK, and so must be compiled on a system with a version
!  available.  On Mac OS X, one can use the optimised 'vecLib' framework
!  by compiling with the options '-framework vecLib -llapack -lblas
!
!-------------------------------------------------------------------------------
!
!  Andy Nowacki,
!  University of Bristol
!
!-------------------------------------------------------------------------------
!  HISTORY:
!     2011-01-30: Inception date.  First working version (of sorts).
!     2012-02-24: Added alternative spline basis function, the Abel-Poisson
!                 kernel (sph_slines_func_AP), and subroutine to calculate the
!                 most appropriate value for its parameter, h, for a desired
!                 spline width.  Change splines_set type to keep track of which
!                 kernel we're using.
!
!-------------------------------------------------------------------------------
!  TODO:
!     + Allow for strongly irregular knot spacing by calculating f using a
!       local value for dbar, the mean knot spacing (in that region).
!===============================================================================
module spherical_splines
!===============================================================================

   implicit none

!  ** size constants
      integer, parameter, private :: i4 = selected_int_kind(9)       ! long int
      integer, parameter, private :: r4 = selected_real_kind(6,37)   ! SP
      integer, parameter, private :: r8 = selected_real_kind(15,307) ! DP

!  ** precision selector
      integer, parameter, private :: rs = r8

!  ** maths constants and other useful things
      real(rs), parameter, private :: pi = 3.141592653589793238462643_rs ;
      real(rs), parameter, private :: pi2 = pi/2._rs
      real(rs), parameter, private :: ZERO = 0._rs, ONE = 1._rs, TWO = 2._rs, &
         HALF = 0.5_rs, FOURTH = 0.25_rs, THIRD = 1._rs/3._rs

!  ** logical units
      integer, parameter, private :: lu_stdin  = 5
      integer, parameter, private :: lu_stdout = 6
      integer, parameter, private :: lu_stderr = 0
      integer, parameter, private :: lu_file = 10  ! LU to be used for reading and writing

!  ** Names and number of different kernels available
      integer, parameter :: NUM_KERNELS = 2  ! Number of kernel types available
      integer, parameter :: KERNEL_NAME_LEN = 2, &
                            KERNEL_LONG_NAME_LEN = 30
      character(len=KERNEL_NAME_LEN),dimension(NUM_KERNELS) :: &
                                          KERNEL_NAME = (/ 'B ','AP' /)
      character(len=KERNEL_LONG_NAME_LEN),dimension(NUM_KERNELS) :: &
        KERNEL_LONG_NAME = (/ 'B-spline                      ', &
                              'Abel-Poisson                  ' /)

!  ** Hide the helper functions
      private :: sph_splines_delta
      private :: sph_splines_fatal
      private :: sph_splines_info

!  ** derived types used by the module, containing knot/data points, spline
!     coefficients and so on
      type sph_splines_set
         real(rs),dimension(:),allocatable :: phi,theta, & ! Lon & colat
                                              c, &  ! Value at control points
                                              a     ! Spline coefficient
         real(rs) :: dbar = -1._rs, & ! Mean knot spacing globally (for B-splines) &
                     h = -1._rs       !  Abel-Poisson coefficient (-1 if unset)
         integer :: n                 ! Number of control knots
         logical :: degrees = .true.  ! Knot units; defaults to degrees
         character(len=KERNEL_NAME_LEN) :: kernel ! Type of kernel being used
      end type

!-------------------------------------------------------------------------------
CONTAINS

!===============================================================================
subroutine sph_splines_fit(s,quiet)
!===============================================================================
!  Calculate the coefficients of the spherical spline representation of a surface
!  field, with knot points as given in the data.
!
!  INPUT: sph_splines_set structure.  The following are unchanged
!         phi  \  Column vectors of implicit length which contain coordinates
!         theta } of the knot points and the surface values there.
!         c    /
!         degress Logical variable specifying whether input coordinates
!                 are in degrees or not.  (Default .true.)
!
!  OUTPUT: sph_splines_set structure.  The following are changed
!         a      Column vector of coefficients which can be used to calculate
!                the fitted function.  (Re-)Allocated within the subroutine.
!         dbar   Mean distance between knots points globally.
!
!  Constructs the set of knot coefficients by inverting the linear equations
!     c(theta_i,phi_i) = f_j(Delta_ij) * a(i)
!  ->                a = F^-1 * c
!  For this, I use LAPACK

   implicit none

   type(sph_splines_set) :: s
   logical,optional,intent(in) :: quiet
   logical :: silent
   integer :: n  ! number of knots; length of phi,theta,r
   integer :: i,j,ier,nd,npoints
   integer,allocatable :: ipiv(:)  ! Pivot rows used by LAPACK
   real(rs),allocatable :: F(:,:)   ! Matrix of distances between point i and j
   real(rs) :: dmin

!  Get the level of verbosity.  Default to noisy
   silent = .false.
   if (present(quiet)) silent = quiet

!  Check for the right things being present
   if (.not.(allocated(s%phi) .and. allocated(s%theta) .and. allocated(s%c))) &
      call sph_splines_fatal('sph_splines_fit: location or values of knot points not allocated')

!  Get length of input arrays and make sure they're the same length
   n = size(s%phi)
   if (s%n /= n .or. size(s%theta) /= n .or. size(s%c) /= n) &
      call sph_splines_fatal('sph_splines_fit: structure arrays must all be same length.')

!  (Re-)Allocate s%a if already allocated to a length != n
   if (allocated(s%a)) then
      if (size(s%a) /= n) then
         deallocate(s%a);  allocate(s%a(n))
         if (.not.silent) call sph_splines_info( &
               'sph_splines_fit: allocating memory for spline coefficients.')
      endif
   else
      allocate(s%a(n))
      if (.not.silent) call sph_splines_info( &
            'sph_splines_fit: allocating memory for spline coefficients.')
   endif

!  Construct matrix of distances between knots, F_ij = Delta_ij
   if (.not.silent) call sph_splines_info( &
      'sph_splines_fit: Calculating inter-knot spacings.')
   allocate(F(n,n))
   do i=1,n
      do j=i,n
         if (i==j) then
            F(i,j) = 0._rs
         else
            F(i,j) = sph_splines_delta(s%phi(i),s%theta(i),s%phi(j),s%theta(j),degrees=s%degrees)
         endif
         F(j,i) = F(i,j)
      enddo
   enddo

!  Work out average distance between adjacent points globally, dbar
   if (.not.silent) call sph_splines_info( &
      'sph_splines_fit: Calculating mean knot spacing.')
   s%dbar = 0._rs
   nd = 0
   do i=1,n
!  Take the nearest 6 or so points and take the average distance of these
      dmin = 0._rs
      npoints = 0
      do while (npoints <= 6)
         dmin = minval(F(i,:), F(i,:) > dmin)
         s%dbar = s%dbar + sum(F(i,:), F(i,:) == dmin)
         npoints = npoints + count(F(i,:) == dmin)
         nd = nd + count(F(i,:) == dmin)
      enddo
   enddo
   s%dbar = s%dbar/nd

!  Update F to contain the value of the B-spline function, not just the distance
!  F_ij = f_j(Delta_ij)
   if (.not.silent) call sph_splines_info('sph_splines_fit: Constructing F matrix.')
   do i=1,n
      do j=i,n
         F(i,j) = sph_splines_func(F(i,j),s%dbar)
!         F(j,i) = F(i,j)  ! Only need the upper half in passing to ?posv
      enddo
   enddo

!  Solve the set of equations F*a = c  (In LAPACK terms, A*X = B).
!  The knot coefficients, a, are placed in B, so we copy c into a for this.
!  F is replaced by the L and U factors from the factorisation.
   s%a = s%c

!  Create array for pivot rows
   allocate(ipiv(n))
!  Call the appropriate LAPACK routine for single or double precision
!  The matrix F is real, symmetric, and the diagonal is 0.  All eigenvalues are
!  positive, so it is positive definite.  (Please tell me if this is wrong!)
!  Hence solving by Cholesky decomposition should be faster than calling the more
!  general ?gesv, which uses LU decomposition.
!  If the input values do not give a positive definite matrix, LAPACK will report
!  ier > 0 and the users has probably entered the wrong knot locations or values.
   if (.not.silent) write(lu_stderr,'(a,i0,a,i0,a)') &
      'spherical_splines: sph_splines_fit: Finding spline coefficients by inverting ',&
      n,' x ',n,' matrix.'
   if (rs == r8) then
!      call dgesv(n,1,F,n,ipiv,s%a,n,ier)  ! Works for all non-singular matrices
      call dposv('U',n,1,F,n,s%a,n,ier)   ! Works for positive definite matrices
   else if (rs == r4) then
!      call sgesv(n,1,F,n,ipiv,s%a,n,ier)
      call sposv('U',n,1,F,n,s%a,n,ier)
   else
      call sph_splines_fatal('module compiled with neither real or double precision')
   endif

!  Check status
!   if (ier > 0) call sph_splines_fatal( &
!      'sph_splines_fit: trying to solve singular matrix.  Check knots and values.')
   if (ier > 0) call sph_splines_fatal( &
      'sph_splines_fit: F matrix is not positive definite.  Check knots and values.')
   if (ier < 0) call sph_splines_fatal( &
      'sph_splines_fit: illegal value somewhere in call to LAPACK.')
   if (ier == 0 .and..not.silent) call sph_splines_info( &
      'sph_splines_fit: Successfully solved for knot coefficients.')

!  Deallocate things we won't use outside the subroutine
   deallocate(F,ipiv)

!  Mark down which kernel we've fitted
   s%kernel = 'B '

end subroutine sph_splines_fit
!-------------------------------------------------------------------------------

!===============================================================================
subroutine sph_splines_fit_AP(s,half_width,quiet)
!===============================================================================
!  Calculate the spline coefficients for the Abel-Poisson kernel representation
!  h is the desired scaling parameter.  Using sph_splines_AP_width2h can give you
!  a value of h for the desired half-width in degrees of the kernel.
   implicit none
   type(sph_splines_set),intent(inout) :: s
   real(rs),intent(in),optional :: half_width  ! Degrees
   logical,intent(in),optional :: quiet
   logical :: silent
   integer :: n  ! number of knots; length of phi,theta,r
   integer :: i,j,ier
   real(rs),allocatable :: F(:,:)   ! Matrix of distances between point i and j
!   real(rs) :: d,dmin

!  Get the level of verbosity.  Default to noisy
   silent = .false.
   if (present(quiet)) silent = quiet

!  Check for conflict between supplied h or requested half_width
   if (s%h > 0._rs .and. present(half_width) .and. .not.silent) &
      call sph_splines_info('sph_splines_fit_AP: overwriting value of h with that ' &
                         // 'generated by call to sph_splines_AP_width2h')

!  Calculate the desired value of h if we don't know it already
   if (present(half_width)) then
      s%h = sph_splines_AP_width2h(half_width,degrees=.true.)
      if (.not.silent) write(lu_stderr,'(2a,f0.3)') 'spherical_splines: ', &
            'sph_splines_fit_AP: Using h-value of ',s%h
   endif

!  Check for the right things being present
   if (.not.(allocated(s%phi) .and. allocated(s%theta) .and. allocated(s%c))) &
      call sph_splines_fatal('sph_splines_fit_AP: location or values of knot points not allocated')

!  Get length of input arrays and make sure they're the same length
   n = size(s%phi)
   if (s%n /= n .or. size(s%theta) /= n .or. size(s%c) /= n) &
      call sph_splines_fatal('sph_splines_fit_AP: structure arrays must all be same length.')

!  (Re-)Allocate s%a if already allocated to a length != n
   if (allocated(s%a)) then
      if (size(s%a) /= n) then
         deallocate(s%a);  allocate(s%a(n))
         if (.not.silent) call sph_splines_info( &
               'sph_splines_fit_AP: allocating memory for spline coefficients.')
      endif
   else
      allocate(s%a(n))
      if (.not.silent) call sph_splines_info( &
            'sph_splines_fit_AP: allocating memory for spline coefficients.')
   endif

!  Construct matrix of distances between knots, F_ij = Delta_ij
   if (.not.silent) call sph_splines_info( &
      'sph_splines_fit_AP: Calculating inter-knot spacings.')
   allocate(F(n,n))
   do i=1,n
      do j=i,n
         if (i==j) then
            F(i,j) = 0._rs
         else
            F(i,j) = sph_splines_delta(s%phi(i),s%theta(i),s%phi(j),s%theta(j),degrees=s%degrees)
         endif
         F(j,i) = F(i,j)
      enddo
   enddo

!  Work out average distance between adjacent points globally, dbar
!   if (.not.silent) call sph_splines_info( &
!      'sph_splines_fit: Calculating mean knot spacing.')
!
!   s%dbar = 0._rs
!   nd = 0
!   do i=1,n
!  Take the nearest 6 or so points and take the average distance of these
!      dmin = 0._rs
!      npoints = 0
!      do while (npoints <= 6)
!         dmin = minval(F(i,:), F(i,:) > dmin)
!         s%dbar = s%dbar + sum(F(i,:), F(i,:) == dmin)
!         npoints = npoints + count(F(i,:) == dmin)
!         nd = nd + count(F(i,:) == dmin)
!      enddo
!   enddo
!   s%dbar = s%dbar/nd

!  If necessary, convert distances into radians, as assumed by sph_splines_func_AP
   if (s%degrees) F = F*pi/180_rs

!  Update F to contain the value of the AP-spline function, not just the distance
!  F_ij = f_j(Delta_ij)
   if (.not.silent) call sph_splines_info('sph_splines_fit_AP: Constructing F matrix.')
   do i=1,n
      do j=i,n
         F(i,j) = sph_splines_func_AP(F(i,j),s%h)
!         F(j,i) = F(i,j)  ! Only need the upper half in passing to ?posv
      enddo
   enddo

!  Solve the set of equations F*a = c  (In LAPACK terms, A*X = B).
!  The knot coefficients, a, are placed in B, so we copy c into a for this.
!  F is replaced by the L and U factors from the factorisation.
   s%a = s%c

!  Call the appropriate LAPACK routine for single or double precision
!  The matrix F is real, symmetric, and the diagonal is 0.  All eigenvalues are
!  positive, so it is positive definite.  (Please tell me if this is wrong!)
!  Hence solving by Cholesky decomposition should be faster than calling the more
!  general ?gesv, which uses LU decomposition.
!  If the input values do not give a positive definite matrix, LAPACK will report
!  ier > 0 and the users has probably entered the wrong knot locations or values.
   if (.not.silent) write(lu_stderr,'(a,i0,a,i0,a)') &
      'spherical_splines: sph_splines_fit_AP: Finding spline coefficients by inverting ',&
      n,' x ',n,' matrix.'
   if (rs == r8) then
      call dposv('U',n,1,F,n,s%a,n,ier)   ! Works for positive definite matrices
   else if (rs == r4) then
      call sposv('U',n,1,F,n,s%a,n,ier)
   else
      call sph_splines_fatal('module compiled with neither real or double precision')
   endif

!  Check status
!   if (ier > 0) call sph_splines_fatal( &
!      'sph_splines_fit: trying to solve singular matrix.  Check knots and values.')
   if (ier > 0) call sph_splines_fatal( &
      'sph_splines_fit_AP: F matrix is not positive definite.  Check knots and values.')
   if (ier < 0) call sph_splines_fatal( &
      'sph_splines_fit_AP: illegal value somewhere in call to LAPACK.')
   if (ier == 0 .and..not.silent) call sph_splines_info( &
      'sph_splines_fit_AP: Successfully solved for knot coefficients.')

!  Mark down which kernel we've fitted
   s%kernel = 'AP'

end subroutine sph_splines_fit_AP
!-------------------------------------------------------------------------------

!===============================================================================
subroutine sph_splines_eval(phi,theta,s,c,degrees)
!===============================================================================
!  Evaluate a set of splines with coefficients at a set of arbitrary points
!  on a sphere
   implicit none
   type(sph_splines_set),intent(in) :: s
   real(rs),intent(in),dimension(:) :: phi,theta
   real(rs),dimension(:),intent(out) :: c
   logical,optional,intent(in) :: degrees
   real(rs) :: conversion,d,torad
   integer :: n  ! length of phi and theta
   integer :: i,j

!  Get size of arrays and check a few things
   n = size(phi)
   if (size(theta) /= n .or. size(c) /= n) call sph_splines_fatal( &
      'sph_splines_eval: evaluation coordinates and values must be in arrays of same length')
   if (.not.(allocated(s%phi) .and. allocated(s%theta) .and. allocated(s%a) &
       .and. allocated(s%c))) call sph_splines_fatal( &
      'sph_splines_eval: one or more arrays in the supplied spline set not allocated.')
   if (.not.any(KERNEL_NAME == s%kernel)) call sph_splines_fatal( &
      'sph_splines_eval: kernel type "' // s%kernel // '" not supported.')

!  If necessary, convert input phi and theta to the same units as the
!  splines are given in.  Assume they are the same
   conversion = 1._rs
   if (present(degrees)) then
      if (degrees .eqv. s%degrees) then
         conversion = 1._rs
      else if (degrees) then      ! Must mean s%degrees = .false.
         conversion = pi/180._rs
      else if (.not.degrees) then ! Must mean s%degrees = .true.
         conversion = 180._rs/pi
      else
         call sph_splines_fatal('sph_splines_eval: Andy''s coding is wrong.')
      endif
   else  ! Assume input is in degrees
      if (s%degrees) conversion = 1._rs
      if (.not. s%degrees) conversion = pi/180._rs
   endif

!  We need to pass radians to the AP spline function, so one final conversion.
!  The B-sline function merely needs d and dbar in the same units, so no need to
!  change units for this function.
   torad = 1._rs
   if (s%degrees .and. s%kernel=="AP") torad = pi/180._rs

!  Evaluate spline at each point
   do i=1,n
!      c(i) = sum(s%a*sph_splines_func( &
!               sph_splines_delta(phi(i),theta(i),s%phi,s%theta,degrees=s%degrees),&
!               s%dbar)))
      c(i) = 0._rs
      do j=1,s%n
         d = torad * sph_splines_delta(conversion*phi(i),conversion*theta(i),&
                               s%phi(j),s%theta(j),degrees=s%degrees)
         if (s%kernel == 'B ') then
            c(i) = c(i) + sph_splines_func(d,s%dbar)*s%a(j)
         else if (s%kernel == 'AP') then
            c(i) = c(i) + sph_splines_func_AP(d,s%h)*s%a(j)
         else
            call sph_splines_fatal('sph_splines_fit: kernel "' // s%kernel // &
               '" is not recognised.')
         endif
      enddo
   enddo

end subroutine sph_splines_eval
!-------------------------------------------------------------------------------

!===============================================================================
function sph_splines_func(d,dbar) result(f)
!===============================================================================
!  Compute the spherical B-splines function:
!         { (3/4)d^3/(dbar^3) - (6/4)d^2/(dbar^2) + 1,         d <= dbar
!     f = {-(1/4)d1^3 + (3/4)d1^2 - (3/4)d1 + (1/4)  , dbar <= d <= 2dbar
!         { 0                                        ,         d >  2dbar
!  where
!     d1 = (d - dbar)/dbar.

   implicit none
   real(rs),intent(in)  :: d    ! Distance between this point and knot point
   real(rs),intent(in)  :: dbar ! Mean knot spacing
   real(rs) :: f
   real(rs) :: d1

!  Where the distance is zero, the value of the function is 1.
   if (d < 0._rs) then
      call sph_splines_fatal('sph_splines_func: distance must be > 0.')
   else if (d == 0.) then
      f = 1._rs
   else if (d <= dbar) then
      f = (3._rs/4._rs)*d**3/dbar**3 - (6._rs/4._rs)*d**2/dbar**2 + 1._rs
   else if (d > dbar .and. d <= 2._rs*dbar) then
      d1 = (d - dbar)/dbar
      f = -(1._rs/4._rs)*d1**3 + (3._rs/4._rs)*d1**2 - &
          (3._rs/4._rs)*d1 + (1._rs/4._rs)
   else  ! Outside 2dbar, the value and its derivatives are zero.
      f = 0._rs
   endif

end function sph_splines_func
!-------------------------------------------------------------------------------

!===============================================================================
function sph_splines_func_AP(d,h) result(K)
!===============================================================================
!  Compute the Abel-Poisson spherical splines kernel:
!     Kh(x,y) = (1/4*pi)*(1 - h**2)/(1 + h**2 - 2*h*dot_product(x,y))**(3/2),
!  where
!     x and y are points on the unit sphere.  (Hence x.y is the cosine of the
!  angular distance (d) between the two points.)
!     h is the scaling factor, 0<=h<=1.  Values of h=0.96 correspond to a
!  half-width of maximum intensity at ~1-h (radians)

   implicit none
   real(rs),intent(in) :: d,h  ! Angle between points (radians); and scaling
   real(rs) :: K

   if (d < 0._rs) &
      call sph_splines_fatal('sph_splines_func_AP: distance must be > 0.')
   K = (1._rs/(4._rs*pi))*(1._rs - h**2)/(1._rs + h**2 - 2._rs*h*cos(d))**(1.5_rs)

end function sph_splines_func_AP
!-------------------------------------------------------------------------------

!===============================================================================
subroutine sph_splines_write_coefs(s,outfile)
!===============================================================================
!  Write out the coefficients for a spherical splines set.
!
!  Format is (ascii):
!     [1]      n kernel g     ! n = num coefficients
!                             ! kernel = type of kernel
!                             ! g is dbar or h, depending on kernel
!     [2..n+1] lon lat c      ! lon, colat (degrees), value and coefficient

   implicit none
   type(sph_splines_set),intent(in) :: s
   character(len=*), intent(in) :: outfile
   real(rs) :: temp_val,conversion
   logical :: opened
   integer :: unit,i,iostatus

   call sph_splines_fatal('sph_splines_write_ceofs: Subroutine not tested yet!')

!  Look for available logical unit
   unit = 10
   opened = .true.
   do while(opened)
      inquire(unit,opened=opened)
      if (.not.opened) exit
      if (unit == 100) then
         call sph_splines_fatal( &
            'sph_splines_write_coefs: Unable to find available unit for writing.')
      endif
      unit = unit + 1
   enddo

!  Open file
   open(unit,file=trim(outfile),iostat=iostatus)
   if (iostatus /= 0) call sph_splines_fatal( &
      'sph_splines_write_coefs: Problem opening file "' // trim(outfile) // '"')

!  Write header line
   if (s%kernel == 'B') then
      temp_val = s%dbar
      if (.not.s%degrees) temp_val = s%dbar * 180._rs/pi
   elseif (s%kernel == 'AP') then
      temp_val = s%h
   else
      call sph_splines_fatal( &
         'sph_splines_write_coefs: Unknown kernel "' // s%kernel // '"')
   endif

!  If necessary, convert phi and theta to degrees for output
   conversion = 1._rs
   if (.not.s%degrees) conversion = 180._rs/pi
!  Write header
   write(unit,*) s%n, s%kernel, temp_val
!  Write knot points and coefficients
   do i=1,s%n
      write(unit,*) conversion*s%phi(i), conversion*s%theta(i), s%c(i), s%a(i)
   enddo

   close(unit)

end subroutine sph_splines_write_coefs
!-------------------------------------------------------------------------------

!===============================================================================
function sph_splines_read_coefs(infile) result(s)
!===============================================================================
!  Read the coefficients for a spherical splines set
!  For format of files, see sph_splines_write_coefs

   implicit none
   type(sph_splines_set) :: s
   character(len=*),intent(in) :: infile
   real(rs) :: temp_val
   integer :: i,iostatus,unit
   logical :: opened

   call sph_splines_fatal('sph_splines_read_coefs: Subroutine not tested yet!')

!  Look for available logical unit
   unit = 10
   opened = .true.
   do while(opened)
      inquire(unit,opened=opened)
      if (.not.opened) exit
      if (unit == 100) then
         call sph_splines_fatal( &
            'sph_splines_read_coefs: Unable to find available unit for reading.')
      endif
      unit = unit + 1
   enddo

!  Open file
   open(unit,file=trim(infile),iostat=iostatus)
   if (iostatus /= 0) call sph_splines_fatal( &
      'sph_splines_read_coefs: Problem opening file "' // trim(infile) // '"')

!  Read header
   read(unit,*) s%n, s%kernel, temp_val
   if (s%kernel == 'B') then
      s%dbar = temp_val
   elseif (s%kernel == 'AP') then
      s%h = temp_val
   else
      call sph_splines_fatal( &
         'sph_splines_read_coefs: Unknown kernel "' // trim(s%kernel) // '"')
   endif

!  Allocate space for coordinates and coefficients
   call sph_splines_realloc(s%n,s)

!  Read in data
   do i=1,s%n
      read(unit,*) s%phi(i), s%theta(i), s%c(i), s%a(i)
   enddo

   close(unit)

!  Set the units of the coefficients to be degrees
   s%degrees = .true.

end function sph_splines_read_coefs
!-------------------------------------------------------------------------------

!===============================================================================
function sph_splines_AP_width2h(x,degrees) result(h)
!===============================================================================
!  Using an empirical power law, return the value of h which produces a
!  Abel-Poisson kernel which has the requested half-width.
!  Optionally, specify desired half-width in degrees; default is radians.
   implicit none
   real(rs),intent(in) :: x
   logical,intent(in),optional :: degrees
   real(rs) :: h,r
   real(rs),parameter :: a = -1.31743_rs, b = 0.89017_rs, c = -0.315546_rs, &
                         d = 0.0468784_rs

   r = x  ! Input in radians
   if (present(degrees)) then
      if (degrees) r = pi*x/180._rs  ! Input in degrees
   endif

   h = 1._rs + a*r + b*r**2 + c*r**3 + d*r**4

end function sph_splines_AP_width2h
!-------------------------------------------------------------------------------

!===============================================================================
function sph_splines_delta(phi1,theta1,phi2,theta2,degrees) result(delta)
!===============================================================================
!  delta returns the angular distance between two points on a sphere given
!  the colat and lon of each using the Haversine formula

   implicit none
   real(rs),intent(in) :: phi1,theta1,phi2,theta2
   real(rs) :: delta
   real(rs) :: lon1,lon2,lat1,lat2,conversion
   logical,optional,intent(in) :: degrees

!   real(rs) :: cos_lat1,cos_lat2,sin_lat1,sin_lat2,&
!               cos_lon1,cos_lon2,sin_lon1,sin_lon2

!  Default to degrees input
   conversion = pi/180._rs
   if (present(degrees)) then
      if (.not.degrees) conversion = 1._rs
   endif

!  Convert to radians and from colat to lat for calculation
   lon1 = phi1*conversion         ;  lon2 = phi2*conversion
   lat1 = pi2 - theta1*conversion ;  lat2 = pi2 - theta2*conversion

   delta = atan2( sqrt( (cos(lat2)*sin(lon2-lon1))**2 + (cos(lat1)*sin(lat2) - &
           sin(lat1)*cos(lat2)*cos(lon2-lon1))**2) , &
           sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon2-lon1))

!  Convert back to desired units
   delta = delta/conversion

end function sph_splines_delta
!-------------------------------------------------------------------------------

!===============================================================================
subroutine sph_splines_fatal(string)
!===============================================================================
!  Utility to flag up a fatal error which stops the program.  Input should be:
!  <name of subroutine>: <error message>
   implicit none
   character(len=*) :: string

   write(lu_stderr,'(2a)') "spherical_splines: ",trim(string)
!  Assuming compiler compatibility, exit with non-0 status
   stop 1
end subroutine sph_splines_fatal
!-------------------------------------------------------------------------------

!===============================================================================
subroutine sph_splines_info(string)
!===============================================================================
!  Utility to write out some information to stderr.  Input should be:
!  <name of subroutine>: <info message>
   implicit none
   character(len=*) :: string

   write(lu_stderr,'(2a)') "spherical_splines: ",trim(string)
end subroutine sph_splines_info
!-------------------------------------------------------------------------------

!===============================================================================
subroutine sph_splines_realloc(n,s)
!===============================================================================
!  (Re)allocate the right amount of space for the arrays in the sph_splines_set
!  structure.
!  Defaults to verbose, so supply quiet=.true. to hush the subroutine up.

   implicit none
   type(sph_splines_set),intent(inout) :: s
   integer,intent(in) :: n

   if (allocated(s%phi)) then
      if (size(s%phi) /= n) then
         deallocate(s%phi);  allocate(s%phi(n))
      endif
   else
      allocate(s%phi(n))
   endif

   if (allocated(s%theta)) then
      if (size(s%theta) /= n) then
         deallocate(s%theta);  allocate(s%theta(n))
      endif
   else
      allocate(s%theta(n))
   endif

   if (allocated(s%a)) then
      if (size(s%a) /= n) then
         deallocate(s%a);  allocate(s%a(n))
      endif
   else
      allocate(s%a(n))
   endif

   if (allocated(s%c)) then
      if (size(s%c) /= n) then
         deallocate(s%c);  allocate(s%c(n))
      endif
   else
      allocate(s%c(n))
   endif

end subroutine sph_splines_realloc
!-------------------------------------------------------------------------------

!_______________________________________________________________________________
end module spherical_splines

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork


text:
! mem1.f90 --
!     Simply include the source text
!

public :: printline
private

contains

! printline --
!     Subroutine to print a line of text
!
subroutine printline( string )
    character(len=*) :: string

    write(*,*) string

end subroutine printline

max_stars_repo_name: dorrabbkhan/atlnts-infrastructure

max_issues_repo_name: TerribleDev/Fortran-docker-mvc

max_forks_repo_name: TerribleDev/Fortran-docker-mvc


text:
MODULE ghex_structured_mod
  use iso_c_binding
  use ghex_defs
  use ghex_comm_mod

  implicit none

  interface
     ! callback type
     subroutine f_cart_rank_neighbor (id, offset_x, offset_y, offset_z, nbid_out, nbrank_out)
       use iso_c_binding
       integer(c_int), value, intent(in) :: id, offset_x, offset_y, offset_z
       integer(c_int), intent(out) :: nbid_out, nbrank_out
     end subroutine f_cart_rank_neighbor
  end interface

  ! ---------------------
  ! --- module types
  ! ---------------------
  ! field descriptor - always associated with a domain
  type, bind(c) :: ghex_struct_field
     type(c_ptr)    :: data = c_null_ptr
     integer(c_int) ::  offset(3) = -1         ! by default - values from the halo
     integer(c_int) :: extents(3) = -1         ! by default - size of the local extents + halos
     integer(c_int) ::    halo(6) = -1         ! halo to be used for this field
     logical(c_bool) :: periodic(3) = .false.
     integer(c_int) :: n_components = 1        ! number of field components
     integer(c_int) ::     layout = GhexLayoutFieldLast
  end type ghex_struct_field

  ! computational domain: defines the iteration space, and the fields with their halos
  type, bind(c) :: ghex_struct_domain
     type(c_ptr)    :: fields = c_null_ptr  ! computational field data, opaque field not to be accessed by the user
     integer(c_int) :: id = -1
     integer(c_int) :: device_id = GhexDeviceUnknown
     integer(c_int) :: first(3)             ! indices of the first LOCAL grid point, in global index space
     integer(c_int) :: last(3)              ! indices of the last LOCAL grid point, in global index space
     integer(c_int) :: gfirst(3) = [1,1,1]  ! indices of the first GLOBAL grid point, (1,1,1) by default
     integer(c_int) :: glast(3)             ! indices of the last GLOBAL grid point (model dimensions)

     ! optional, used by the staged communication object
     type(c_funptr) :: cart_nbor
  end type ghex_struct_domain

  ! structured grid communication object
  type, bind(c) :: ghex_struct_communication_object
     type(c_ptr) :: ptr = c_null_ptr
  end type ghex_struct_communication_object

  ! definition of the exchange, including physical fields and the communication pattern
  type, bind(c) :: ghex_struct_exchange_descriptor
     type(c_ptr) :: ptr = c_null_ptr
  end type ghex_struct_exchange_descriptor

  ! a handle to track a particular communication instance, supports wait()
  type, bind(c) :: ghex_struct_exchange_handle
     type(c_ptr) :: ptr = c_null_ptr
  end type ghex_struct_exchange_handle

  ! ---------------------
  ! --- module C interfaces
  ! ---------------------
  interface
     type(ghex_struct_exchange_descriptor) function ghex_struct_exchange_desc_new_wrapped(domains_desc, n_domains) &
       bind(c, name="ghex_struct_exchange_desc_new")
       use iso_c_binding
       import ghex_struct_domain, ghex_struct_exchange_descriptor
       type(c_ptr), value :: domains_desc
       integer(c_int), value :: n_domains
     end function ghex_struct_exchange_desc_new_wrapped
  end interface

  ! ---------------------
  ! --- generic ghex interfaces
  ! ---------------------
  interface ghex_free
     subroutine ghex_struct_exchange_desc_free(exchange_desc) bind(c, name="ghex_obj_free")
       use iso_c_binding
       import ghex_struct_exchange_descriptor
       type(ghex_struct_exchange_descriptor) :: exchange_desc
     end subroutine ghex_struct_exchange_desc_free

     subroutine ghex_struct_exchange_handle_free(exchange_handle) bind(c, name="ghex_obj_free")
       use iso_c_binding
       import ghex_struct_exchange_handle
       type(ghex_struct_exchange_handle) :: exchange_handle
     end subroutine ghex_struct_exchange_handle_free

     subroutine ghex_struct_domain_free(domains_desc) bind(c)
       use iso_c_binding
       import ghex_struct_domain
       type(ghex_struct_domain) :: domains_desc
     end subroutine ghex_struct_domain_free

     subroutine ghex_struct_co_free(co) bind(c, name="ghex_obj_free")
       use iso_c_binding
       import ghex_struct_communication_object
       ! reference, not a value - fortran variable is reset to null
       type(ghex_struct_communication_object) :: co
     end subroutine ghex_struct_co_free

     procedure :: ghex_struct_field_free
  end interface ghex_free

  interface ghex_domain_init
     procedure :: ghex_struct_domain_init
  end interface ghex_domain_init

  interface ghex_field_init
     procedure :: ghex_struct_field_init
     procedure :: ghex_struct_field_init_comp
  end interface ghex_field_init

  interface ghex_co_init
     subroutine ghex_struct_co_init(co, comm) bind(c)
       import ghex_struct_communication_object, ghex_communicator
       type(ghex_struct_communication_object) :: co
       type(ghex_communicator), value :: comm
     end subroutine ghex_struct_co_init
  end interface ghex_co_init

  interface ghex_domain_add_field
     subroutine ghex_struct_domain_add_field(domain_desc, field_desc) bind(c)
       use iso_c_binding
       import ghex_struct_domain, ghex_struct_field
       type(ghex_struct_domain) :: domain_desc
       type(ghex_struct_field)  :: field_desc
     end subroutine ghex_struct_domain_add_field
  end interface ghex_domain_add_field

  interface ghex_exchange_desc_new
     procedure :: ghex_struct_exchange_desc_array_new
     procedure :: ghex_struct_exchange_desc_new
  end interface ghex_exchange_desc_new

  interface ghex_exchange
     type(ghex_struct_exchange_handle) function ghex_struct_exchange(co, exchange_desc) bind(c)
       import ghex_struct_exchange_handle, ghex_struct_communication_object, ghex_struct_exchange_descriptor
       type(ghex_struct_communication_object), value :: co
       type(ghex_struct_exchange_descriptor), value :: exchange_desc
     end function ghex_struct_exchange
  end interface ghex_exchange

  interface ghex_wait
     subroutine ghex_struct_exchange_handle_wait(exchange_handle) bind(c)
       use iso_c_binding
       import ghex_struct_exchange_handle
       type(ghex_struct_exchange_handle) :: exchange_handle
     end subroutine ghex_struct_exchange_handle_wait
  end interface ghex_wait

CONTAINS

  subroutine ghex_struct_domain_init(domain_desc, id, first, last, gfirst, glast, cart_nbor, device_id)
    use iso_c_binding
    type(ghex_struct_domain) :: domain_desc
    integer :: id
    integer :: first(3)
    integer :: last(3)
    integer :: gfirst(3)
    integer :: glast(3)
    procedure(f_cart_rank_neighbor), optional :: cart_nbor
    integer, optional :: device_id
    type(c_funptr) :: ptest

    if (present(cart_nbor)) then
      domain_desc%cart_nbor = c_funloc(cart_nbor)
    end if

    if (present(device_id)) then
      domain_desc%device_id = device_id
    else
      domain_desc%device_id = GhexDeviceCPU
    end if

    domain_desc%id = id
    domain_desc%first = first
    domain_desc%last = last
    domain_desc%gfirst = gfirst
    domain_desc%glast = glast
  end subroutine ghex_struct_domain_init

  subroutine ghex_struct_field_init(field_desc, data, halo, offset, periodic)
    type(ghex_struct_field) :: field_desc
    real(ghex_fp_kind), dimension(:,:,:), target :: data
    integer :: halo(6)
    integer, optional :: offset(3)
    logical, optional :: periodic(3)

    field_desc%data = c_loc(data)
    field_desc%halo = halo
    field_desc%extents = shape(data, 4)
    field_desc%n_components = 1
    field_desc%layout = GhexLayoutFieldLast

    if (present(offset)) then
      field_desc%offset = offset
    else
      field_desc%offset = [halo(1), halo(3), halo(5)]
    endif

    if (present(periodic)) then
      field_desc%periodic = periodic
    endif
  end subroutine ghex_struct_field_init

  subroutine ghex_struct_field_init_comp(field_desc, data, halo, offset, periodic, layout)
    type(ghex_struct_field) :: field_desc
    real(ghex_fp_kind), dimension(:,:,:,:), target :: data
    integer :: halo(6)
    integer, optional :: offset(3)
    logical, optional :: periodic(3)
    integer, optional :: layout
    integer(4) :: extents(4)

    field_desc%data = c_loc(data)
    field_desc%halo = halo

    if (present(offset)) then
      field_desc%offset = offset
    else
      field_desc%offset = [halo(1), halo(3), halo(5)]
    endif

    if (present(periodic)) then
      field_desc%periodic = periodic
    endif

    if (present(layout)) then
      field_desc%layout = layout
    else
      field_desc%layout = GhexLayoutFieldLast
    endif

    extents = shape(data, 4)
    if (field_desc%layout == GhexLayoutFieldLast) then
      field_desc%extents = extents(1:3)
      field_desc%n_components = size(data, 4)
    else
      field_desc%extents = extents(2:4)
      field_desc%n_components = size(data, 1)
    end if
  end subroutine ghex_struct_field_init_comp

  subroutine ghex_struct_field_free(field_desc)
    type(ghex_struct_field) :: field_desc
    field_desc%data         = c_null_ptr
    field_desc%offset(:)    = -1
    field_desc%extents(:)   = -1
    field_desc%halo(:)      = -1
    field_desc%periodic(:)  = .false.
    field_desc%layout       = GhexLayoutFieldLast
  end subroutine ghex_struct_field_free

  type(ghex_struct_exchange_descriptor) function ghex_struct_exchange_desc_array_new(domains_desc)
    type(ghex_struct_domain), dimension(:), target :: domains_desc
    ghex_struct_exchange_desc_array_new = ghex_struct_exchange_desc_new_wrapped(c_loc(domains_desc), size(domains_desc, 1));
  end function ghex_struct_exchange_desc_array_new

  type(ghex_struct_exchange_descriptor) function ghex_struct_exchange_desc_new(domains_desc)
    type(ghex_struct_domain), target :: domains_desc
    ghex_struct_exchange_desc_new = ghex_struct_exchange_desc_new_wrapped(c_loc(domains_desc), 1);
  end function ghex_struct_exchange_desc_new

END MODULE ghex_structured_mod

max_stars_repo_name: tehrengruber/GHEX

max_issues_repo_name: tehrengruber/GHEX

max_forks_repo_name: tehrengruber/GHEX


text:
C
C	$Id: gtxtat.f,v 1.6 2008-07-27 00:21:03 haley Exp $
C
C                Copyright (C)  2000
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 GTXTAT(IOS,STATUS)
C
C  Set the TEXT attributes.
C
C  Valid attributes:
C    TEXT INDEX BUNDLE
C    TEXT FONT INDEX
C    TEXT PRECISION
C    CHARACTER EXPANSION FACTOR
C    CHARACTER SPACING
C    TEXT COLOR
C    CHARACTER HEIGHT
C    CHARACTER ORIENTATION
C    CHARACTER PATH
C    TEXT ALIGNMENT
C
      include 'trstat.h'
      include 'trpars.h'
      include 'trinst.h'
      include 'trcode.h'
      include 'gkscom.h'
C
      INTEGER IOS, STATUS
      INTEGER TEMP
C
      IF (OPID .EQ. ATELTI) THEN
C
C  TEXT BUNDLE INDEX
C
        CALL GOPDEC(TXTIDX,MIXCPR,1,IOS,STATUS)
        CALL GSTXI(TXTIDX)
      ELSE IF (OPID .EQ. ATELTF) THEN
C
C  FONT INDEX
C
        CALL GOPDEC(TEMP,MIXCPR,1,IOS,STATUS)
        CALL GENNEG(TEMP,FINDEX)
        IF (FINDEX .LE. 0) FINDEX = 1
C
C  Set the font index to point to the appropriate internal font table.
C
        CALL GQTXFP(IER,IFNT,IPREC)
        CALL GSTXFP(FINDEX,IPREC)
C
      ELSE IF (OPID .EQ. ATELTP) THEN
C
C  TEXT PRECISION
C
        CALL GOPDEC(TXTPRE,MENCPR,1,IOS,STATUS)
        CALL GQTXFP(IER,IFNT,IPREC)
        CALL GSTXFP(IFNT,TXTPRE)
      ELSE IF (OPID .EQ. ATELCE) THEN
C
C  CHARACTER EXPANSION FACTOR
C
        CALL GTFLT(CEXPN,MFLCPR,IOS,STATUS)
        CALL GSCHXP(CEXPN)
      ELSE IF (OPID .EQ. ATELCS) THEN
C
C  CHARACTER SPACING
C
        CALL GTFLT(CSPACE,MFLCPR,IOS,STATUS)
        CALL GSCHSP(CSPACE)
      ELSE IF (OPID .EQ. ATELTC) THEN
C
C  Set the text color.
C
        CALL GOPDEC(TXTCOL,MCICPR,1,IOS,STATUS)
        CALL GSTXCI(TXTCOL)
      ELSE IF (OPID .EQ. ATELCH) THEN
C
C  CHARACTER HEIGHT
C
        CALL GOPDEC(CHIGHT,MOPLEN,1,IOS,STATUS)
        XHIGHT = REAL(CHIGHT)/32767.
C
C  Scale the character height by the original scale factor in the
C  Y direction.  Take a vactor (0,XHIGHT), transform it with the
C  current transformation matrix, then find the length of the
C  transformed vector.
C
        YSCALE = SQRT(CURTM(1,2)*CURTM(1,2)+CURTM(2,2)*CURTM(2,2))
        HN = XHIGHT*YSCALE
        IF (HN .LT. 0.00004) HN = 0.00004
        CALL GSCHH(HN)
      ELSE IF (OPID .EQ. ATELCO) THEN
C
C  CHARACTER ORIENTATION (given as up and base vectors).
C
        CALL GOPDEC(TEMP,MOPLEN,1,IOS,STATUS)
        CALL GENNEG(TEMP,XU)
        IF (STATUS .NE. 0) RETURN
        CALL GOPDEC(TEMP,MOPLEN,1,IOS,STATUS)
        CALL GENNEG(TEMP,YU)
        IF (STATUS .NE. 0) RETURN
        CALL GOPDEC(TEMP,MOPLEN,1,IOS,STATUS)
        CALL GENNEG(TEMP,XB)
        IF (STATUS .NE. 0) RETURN
        CALL GOPDEC(TEMP,MOPLEN,1,IOS,STATUS)
        CALL GENNEG(TEMP,YB)
C
C  Modify the up vector by the current segment transformation.
C
        RXUT = CURTM(1,1)*REAL(XU)+CURTM(1,2)*REAL(YU)
        RYUT = CURTM(2,1)*REAL(XU)+CURTM(2,2)*REAL(YU)
C
C  If the modified values are both zero, then the transformation
C  scaling is such as to make the character orientation of no
C  observable consequence, so set it to horizontal.
C
        IF (RXUT.EQ.0. .AND. RYUT.EQ.0.) THEN
          RXUT = 0.
          RYUT = 1.
        ENDIF
        CALL GSCHUP(RXUT,RYUT)
      ELSE IF (OPID.EQ.ATELTH) THEN
C
C  TEXT PATH
C
        CALL GOPDEC(PATH,MENCPR,1,IOS,STATUS)
        CALL GSTXP(PATH)
      ELSE IF (OPID .EQ. ATELTA) THEN
C
C  TEXT ALIGNMENT
C
        CALL GOPDEC(HORIZ,MENCPR,1,IOS,STATUS)
        IF (STATUS .NE. 0) RETURN
        CALL GOPDEC(VERT,MENCPR,1,IOS,STATUS)
        CALL GSTXAL(HORIZ,VERT)
C
C  Skip over the rest of the alignment information.
C
        CALL GSKPOP(8,LEN-(MENCPR*2/8),IOS,STATUS)
      END IF
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


! { dg-do compile }
!
! PR fortran/45170
! PR fortran/52158

module test
 implicit none
 type t
   procedure(deferred_len), pointer, nopass :: ppt
 end type t
contains
 function deferred_len()
   character(len=:), allocatable :: deferred_len
   deferred_len = 'abc'
 end function deferred_len
 subroutine doIt()
   type(t) :: x
   character(:), allocatable :: temp
   x%ppt => deferred_len
   temp = deferred_len()
   if ("abc" /= temp) STOP 1
 end subroutine doIt
end module test

use test
call doIt ()
end


text:
SUBROUTINE VARZEROInt

  ! *** THIS SUBROUTINE ZERO'S ALL OF THE ARRAYS AFTER ALLOCATION

  !----------------------------------------------------------------------!
  ! CHANGE RECORD
  ! DATE MODIFIED     BY               DESCRIPTION
  !----------------------------------------------------------------------!
  !    2015-06       PAUL M. CRAIG     IMPLEMENTED SIGMA-Z (SGZ) IN EE7.3

  USE GLOBAL

  WRITE(*,'(A)')'ZEROING INTEGER ARRAYS'

  ! *** INTEGER ARRAYS
  IACTLR = 0
  IAIJ = 0
  IAKL = 0
  IAP = 0
  ICALJP = 0
  ICBE = 0
  ICBN = 0
  ICBS = 0
  ICBW = 0
  ICCDA = 0
  ICDA = 0
  ICFLMP = 0
  ICSMTS = 0
  IGWSER = 0
  IJCT = 0
  IJCTLT = 0
  IL = 0
  ILLSHA = 0
  ILLT = 0
  ILTMSR = 0
  !IMASKDRY = 0
  IMDCHH = 0
  IMDCHU = 0
  IMDCHV = 0
  INPNS = 0
  INTPSER = 0
  IOUTJP = 0
  IPART = 0
  IPBE = 0
  IPBN = 0
  IPBS = 0
  IPBW = 0
  IQCAX = 0
  IQCTLD = 0
  IQCTLU = 0
  IQJP = 0
  IQS = 0
  IQWRD = 0
  IQWRU = 0
  ISCDRY = 0
  ISDJP = 0
  ISENT = 0
  ISPBE = 0
  ISPBN = 0
  ISPBS = 0
  ISPBW = 0
  ISPNS = 0
  ISSBCP = 0
  ISTJP = 0
  ISUDPC = 0
  ITPCDA = 0
  IUPCJP = 0
  IWGG = 0
  JCBE = 0
  JCBN = 0
  JCBS = 0
  JCBW = 0
  JCCDA = 0
  JCDA = 0
  JL = 0
  JLLSHA = 0
  JLLT = 0
  JLTMSR = 0
  JMDCHH = 0
  JMDCHU = 0
  JMDCHV = 0
  JNPNS = 0
  JPBE = 0
  JPBN = 0
  JPBS = 0
  JPBW = 0
  JQCAX = 0
  JQCTLD = 0
  JQCTLU = 0
  JQJP = 0
  JQS = 0
  JQWRD = 0
  JQWRU = 0
  JSPNS = 0
  JUPCJP = 0
  JWGG = 0
  JWQCBE = 0
  JWQCBN = 0
  JWQCBS = 0
  JWQCBW = 0
  KBT = 0
  KCEFDC = 0
  KEFFJP = 0
  KFEFDC = 0
  KPS = 0
  KQJP = 0
  KQWRD = 0
  KQWRU = 0
  KUPCJP = 0
  KUPW = 0
  LBERC = 0
  LBNRC = 0
  LBSRC = 0
  LBWRC = 0
  LCBE = 0
  LCBN = 0
  LCBS = 0
  LCBW = 0
  LCDA = 0
  LCONSOL = 0
  LCT = 0
  LCTLT = 0
  LIJ = 0
  LIJLT = 0
  LJUNX = 0
  LJUNY = 0
  LLBC = 0
  LLRC = 0
  LLSHA = 0
  LMDCHH = 0
  LMDCHU = 0
  LMDCHV = 0
  LNC = 0
  LNCLT = 0
  LNEC = 0
  LNWC = 0
  LOBCS = 0
  LOBCS2 = 0
  LBCS = 0
  LORDER = 0
  LPBE = 0
  LPBN = 0
  LPBS = 0
  LPBW = 0
  LQCTLD = 0
  LQCTLU = 0
  LQS = 0
  LQWRD = 0
  LQWRU = 0
  LSBLBCD = 0
  LSBLBCU = 0
  LSC = 0
  LSCLT = 0
  LSEC = 0
  LSHAB = 0
  LSHAP = 0
  LSHAU = 0
  LSHAUE = 0
  LSMTS = 0
  LSWC = 0
  LUPU = 0
  LUPV = 0
  LWGG = 0
  LWVCELL = 0
  MCNTLR = 0
  MCSER = 0
  MCTLAST = 0
  MFDCHZ = 0
  MDCHTYP = 0
  MGWSER = 0
  MGWTLAST = 0
  MPSER = 0
  MPTLAST = 0
  MQCTL = 0
  MQSER = 0
  MQTLAST = 0
  MQWRSR = 0
  MQWRTLST = 0
  MTMSRA = 0
  MTMSRC = 0
  MTMSRP = 0
  MTMSRQ = 0
  MTMSRQE = 0
  MTMSRU = 0
  MTMSRUE = 0
  MTMSRUT = 0
  MTSCUR = 0
  MTSSTSP = 0
  MVEGL = 0
  MVEGSER = 0
  MVEGTLAST = 0
  NATDRY = 0
  NCSER = 0
  NCSERA = 0
  NCSERE = 0
  NCSERJP = 0
  NCSERN = 0
  NCSERQ = 0
  NCSERS = 0
  NCSERW = 0
  NDRY = 0
  NGWSL = 0
  NJEL = 0
  NJPMX = 0
  NLOE = 0
  NLON = 0
  NLOS = 0
  NLOW = 0
  NLRPDRT = 0
  NPORTJP = 0
  NPSERE = 0
  NPSERN = 0
  NPSERS = 0
  NPSERW = 0
  NPSERE1 = 0
  NPSERN1 = 0
  NPSERS1 = 0
  NPSERW1 = 0
  NQCMINS = 0
  NQCMUL = 0
  NQCTLQ = 0
  NQCTYP = 0
  NQSERJP = 0
  NQSERQ = 0
  NQSMF = 0
  NQSMUL = 0
  NQWRMFD = 0
  NQWRMFU = 0
  NQWRSERJP = 0
  NQWRSERQ = 0
  NTSCRE = 0
  NTSCRN = 0
  NTSCRS = 0
  NTSCRW = 0
  NTSSSS = 0
  NTVSFP = 0
  NUDJP = 0
  NUDJPC = 0
  IF( ISVEG > 0 )NVEGSERV = 0
  NXYSDAT = 0
  NZPRJP = 0

  ! ** TOXIC TRANSPORT VARIABLES
  ISDIFBW = 0
  ISPMXZ = 0
  ISTOC = 0
  ISTOXR = 0
  ITOXBU = 0
  ITOXKIN = 0
  ITOXWU = 0
  ITXBDUT = 0
  ITXINT = 0
  ITXPARB = 0
  ITXPARBC = 0
  ITXPARW = 0
  ITXPARWC = 0
  QCTLGRP = 0

  IF( ISTRAN(5) >= 1 )THEN
    NSP2 = 0
  ENDIF   ! *** END OF TOXIC VARIABLE DECLARATIONS

  ! ** SED VARIABLES
  IF( ISTRAN(6) >= 1 )THEN
  ENDIF   ! *** END OF SED VARIABLE DECLARATIONS

  ! ** SND VARIABLES
  IF( ISTRAN(7) >= 1 )THEN
  ENDIF   ! *** END OF SND VARIABLE DECLARATIONS

  ! *** EITHER SED OR SND
  IBLTAUC = 0
  ISEDWU = 0
  IROUSE = 0
  ISBDLD = 0
  ISEDSCOR = 0
  ISLTAUC = 0
  ISNDEQ = 0
  ISNDM1 = 0
  ISNDM2 = 0
  ISPROBDEP = 0
  IWRSPB = 0

  IF( ISTRAN(6) >= 1 .OR. ISTRAN(7) >= 1 )THEN
    ISEDBU = 0
    ISNDBU = 0
    ISNDWU = 0

    ISDBLDIR = 0

    JSS3DSED = 0
    JSS3DSND = 0
    JSS3DTOX = 0

    ! *** BANK EROSION
    IBANKBE = 0
    JBANKBE = 0
    ICHANBE = 0
    JCHANBE = 0
    NBESERN = 0
    MBESER = 0
    MBETLAST = 0
    LPMXZ = 0
  ENDIF

  ! ** WATER QUALITY VARIABLES
  IF( ISTRAN(8) >= 1 )THEN
    IBENMAP = 0
    ICPSL = 0
    ICWQTS = 0
    IMWQZT = 0
    IMWQZT1 = 0
    IMWQZT2 = 0
    ISMT = 0
    ISMZMAP = 0
    JCPSL = 0
    KCPSL = 0
    MVPSL = 0
    MWQCTLT = 0
    MWQPSR = 0
    MWQPTLT = 0

    IWQCBE = 0
    IWQCBN = 0
    IWQCBS = 0
    IWQCBW = 0
    IWQKA = 0
    IWQOBE = 0
    IWQOBN = 0
    IWQOBS = 0
    IWQOBW = 0
    IWQPSC = 0
    IWQPSV = 0
    IWQT = 0
    IWQZMAP = 0
  ENDIF   ! *** END OF WQ VARIABLE DECLARATIONS

  IF( ISWQFLUX == 1 )THEN
    ISPV = 0
    IVPV = 0
    JSPV = 0
    JVPV = 0
  ENDIF

  IF( ISWASP > 0 )THEN
    LCEFDC = 0
    LCHNC = 0
    NCHNC = 0
  ENDIF

  ! Begin MHK variables SCJ
  IJLTURB = 0

  ! *** ALLOCATE MEMORY FOR VARIABLE TO STORE CONCENTRATIONS AT OPEN BOUNDARIES

  ! *** NEW VARIABLES FOR QCTL NQCTYP=3 & 4
  NLOWCHORD = 0

  ! *** SIGMA-ZED - SGZ
  KSZ = 1
  KSZU = 1
  KSZV = 1
  IF( IGRIDV > 0 )THEN
    KSZE = 1
    KSZW = 1
    KSZN = 1
    KSZS = 1
  ENDIF

  LKSZ = .FALSE.

  ! *** WET/DRY BYPASS VARIABLES
  LWET = 0
  LDRY = 0

  ! *** WIND
  IF( NWSER >= 1 )THEN
    ISWDINT = 0
  ENDIF

  ! *** ICE
  IF( ISICE > 0 )THEN
    MITLAST = 0
  ENDIF

  ! *** LOGICALS
  LMASKDRY = .FALSE.
  OLDMASK = .FALSE.
  LWVMASK = .FALSE.
  LOPENBCDRY = .FALSE.

  ! *** CHARACTER
  CLSL = ' '
  CLTMSR = ' '
  SYMBOL = ' '

END

max_stars_repo_name: dsi-llc/EFDCPlus

max_forks_repo_name: dsi-llc/EFDCPlus


text:
	SUBROUTINE GR_INTP  ( inttyp, gx, gy, npts, kx, ky, grid,
     +			      sdint,  iret )
C************************************************************************
C* GR_INTP								*
C*									*
C* This subroutine interpolates data from a grid to a set of points	*
C* defined in GX, GY.  Bilinear interpolation is the only interpolation	*
C* type implemented.							*
C*									*
C* GR_INTP  ( INTTYP, GX, GY, NPTS, KX, KY, GRID, SDINT, IRET )		*
C*									*
C* Input parameters:							*
C*	INTTYP		INTEGER		Interpolation type		*
C*	GX (NPTS)	REAL		Grid x coordinates		*
C*	GY (NPTS)	REAL		Grid y coordinates		*
C*	NPTS		INTEGER		Number of coordinates		*
C*	KX		INTEGER		Number of x grid points		*
C*	KY		INTEGER		Number of y grid points		*
C*	GRID (KX,KY)	REAL		Grid data			*
C*									*
C* Output parameters:							*
C*	SDINT (npts)	REAL		Interpolated data		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	11/88	GEMPAK 4.1 GDPINT			*
C* G. Huffman/USRA	 5/89	Add npts, inttyp			*
C* T. Piper/GSC		11/98	Updated prolog				*
C************************************************************************
        INCLUDE         'GEMPRM.PRM'
C*
	REAL		gx (*), gy (*), grid (KX,KY), sdint (*)
C*
        INCLUDE         'ERMISS.FNC'
C------------------------------------------------------------------------
	iret  = 0
C*
	fkx = FLOAT ( kx )
	fky = FLOAT ( ky )
	DO  i = 1, npts
C
C*	    Define values for interpolation.  The truncation to integers
C*	    works because the locations are always non-negative.
C*	    If the point is exactly on the right or top boundary,
C*	    reset pointers so that valid grid points are used.
C
	    icol   = gx ( i )
	    irow   = gy ( i )
	    IF  ( gx ( i ) .eq. fkx )  icol = icol - 1
	    IF  ( gy ( i ) .eq. fky )  irow = irow - 1
	    icolp1 = icol + 1
	    irowp1 = irow + 1
	    c      = gx ( i ) - FLOAT ( icol )
	    r      = gy ( i ) - FLOAT ( irow )
	    omc    = 1. - c
	    omr    = 1. - r
C
C*	    Check that surrounding grid points exist.
C
	    IF  ( ( ERMISS ( grid ( icol,   irow   ) ) ) .or.
     +	          ( ERMISS ( grid ( icolp1, irow   ) ) ) .or.
     +	          ( ERMISS ( grid ( icol,   irowp1 ) ) ) .or.
     +	          ( ERMISS ( grid ( icolp1, irowp1 ) ) ) .or.
     +		  ( gy (i) .gt. fky ) .or. ( gx (i) .gt. fkx) )  THEN
		sdint ( i ) = RMISSD
	      ELSE
	        sdint ( i ) = ( grid ( icol, irow ) * omc +
     +		                grid ( icolp1, irow ) * c ) * omr +
     +		              ( grid ( icol, irowp1 ) * omc +
     +		                grid ( icolp1, irowp1 ) * c ) * r
	    END IF
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak