subroutine guardarResultados(saveASmatlab,saveAstxt,sl)
use debugstuff
use resultVars, only : allpoints,nPts
use sourcevars, only: PSV,SH,currentiFte,nfuentes,Po
use modelVars, only: NFREC,NPTSTIME,dfrec,t_vec,Uo
use glovars, only: z0,pi
implicit none
integer :: ip,i,comp
character(LEN=100) :: yax
logical, intent(in) :: saveASmatlab,saveAstxt,sl
real :: t0
complex*16, dimension(NPTSTIME) :: S
complex*16, dimension(:,:,:), allocatable :: ShR
character(LEN=32)  :: name
character(LEN=100) :: titleN
character(LEN=3),dimension(8) :: cname

!    W_to_t  guarda la respuesta en frecuencia y tiempo en los formatos indicados
   allocate(t_vec(NPTSTIME))
   allocate(Uo(NPTSTIME))
   call system("mkdir out")
   call chdir("out")
   CALL getcwd(yax)
      write(6,'(a,/,a)') 'Saving files At:',TRIM(yax)
   call system('rm *.m *.txt')

   if (sl) allocate(ShR(nPts,8,NPTSTIME))

   do currentiFte=1,nFuentes
      call sourceAmplitudeFunction
      if (sl) ShR = 0; cname(1:8) = ''
      ! para la correción del tiempo inicial de la señal
      t0 = Po(currentiFte)%t0
      t_vec(1:NPTSTIME) = z0
      t_vec(1) = exp(cmplx(0.0,-0.01* dfrec * t0*(2*pi),8))
      do i = 2,nfrec+1
        t_vec(i) = exp(cmplx(0.0,-(i-1)* dfrec * t0*(2*pi),8))
      end do!
      do i = NPTSTIME-NFREC+2,NPTSTIME
        t_vec(i) = exp(cmplx(0.0,-(i-NPTSTIME-1)* dfrec * t0*(2*pi),8))
      end do

      if (PSV) then
        do iP = 1,nPts
        !#< r ____________   W   _________________________________________ !#>
          write(yAx,'(a)') '$u_3$ [m]'
!         print*,ip,"---------------------",allpoints(iP)%resp(:,currentiFte)%W; print*," "
          call W_to_t(S,allpoints(iP)%resp(:,currentiFte)%W,'w--',yax,iP,&
                    allpoints(iP)%center%x,&
                    allpoints(iP)%center%z,currentiFte,saveASmatlab,saveAstxt)
          if (sl) shR(iP,3,:) = S; cname(3) = '_w_'
        !#< r ____________   U   _________________________________________ !#>
          write(yAx,'(a)') '$u_1$ [m]'
          call W_to_t(S,allpoints(iP)%resp(:,currentiFte)%U,'u--',yax,iP,&
                    allpoints(iP)%center%x,&
                    allpoints(iP)%center%z,currentiFte,saveASmatlab,saveAstxt)
          if (sl) shR(iP,1,:) = S; cname(1) = '_u_'
        !#< r ____________   sxx   _________________________________________ !#>
          write(yAx,'(a)') '$\sigma_{xx}$ [m]'
          call W_to_t(S,allpoints(iP)%resp(:,currentiFte)%sxx,'sxx',yax,iP,&
                    allpoints(iP)%center%x,&
                    allpoints(iP)%center%z,currentiFte,saveASmatlab,saveAstxt)
          if (sl) shR(iP,4,:) = S; cname(4) = 'sxx'
        !#< r ____________   szx   _________________________________________ !#>
          write(yAx,'(a)') '$\sigma_{zx}$ [m]'
          call W_to_t(S,allpoints(iP)%resp(:,currentiFte)%szx,'szx',yax,iP,&
                    allpoints(iP)%center%x,&
                    allpoints(iP)%center%z,currentiFte,saveASmatlab,saveAstxt)
          if (sl) shR(iP,5,:) = S; cname(5) = 'szx'
        !#< r ____________   szz   _________________________________________ !#>
          write(yAx,'(a)') '$\sigma_{zz}$ [m]'
          call W_to_t(S,allpoints(iP)%resp(:,currentiFte)%szz,'szz',yax,iP,&
                    allpoints(iP)%center%x,&
                    allpoints(iP)%center%z,currentiFte,saveASmatlab,saveAstxt)
          if (sl) shR(iP,6,:) = S; cname(6) = 'szz'
        end do
      end if !psv
      if (SH) then
        do iP = 1,nPts
          !#< r ____________   V   _________________________________________ !#>
          write(yAx,'(a)') '$u_2$ [m]'
          call W_to_t(S,allpoints(iP)%resp(:,currentiFte)%V,'v--',yax,iP,&
                    allpoints(iP)%center%x,&
                    allpoints(iP)%center%z,currentiFte,saveASmatlab,saveAstxt)
          if (sl) shR(iP,2,:) = S; cname(2) = '_v_'
          write(yAx,'(a)') '$\sigma_{xy}$ [m]'
          call W_to_t(S,allpoints(iP)%resp(:,currentiFte)%sxy,'sxy',yax,iP,&
                    allpoints(iP)%center%x,&
                    allpoints(iP)%center%z,currentiFte,saveASmatlab,saveAstxt)
          if (sl) shR(iP,7,:) = S; cname(7) = 'sxy'
          write(yAx,'(a)') '$\sigma_{zy}$ [m]'
          call W_to_t(S,allpoints(iP)%resp(:,currentiFte)%szy,'szy',yax,iP,&
                    allpoints(iP)%center%x,&
                    allpoints(iP)%center%z,currentiFte,saveASmatlab,saveAstxt)
          if (sl) shR(iP,8,:) = S; cname(8) = 'szy'
!          call W_to_t(allpoints(iP)%resp(:,currentiFte)%Ty,'Ty-',yax,iP,&
!                    allpoints(iP)%center%x,&
!                    allpoints(iP)%center%z,saveASmatlab,saveAstxt)
        end do !ip
      end if !sh

      if (sl) then
      do comp=1,8
       if (cname(comp) .eq. '') cycle
       write(titleN,'(a,I0,a,a,a)') 'In_',currentiFte,'_',cname(comp),'.m'
       OPEN(3555,FILE=trim(titleN),FORM="FORMATTED",ACTION='WRITE')
       !write(name,'(a,I0,a,I0)') 'inci_',currentiFte,'_',comp
        write(name,'(a)') 'var'
       call scripToMatlabMNmatrixZ(nPts,NPTSTIME,shR(:,comp,:),name,3555)
       close(3555)
      end do
      end if
    end do !insidencia
end subroutine guardarResultados

subroutine W_to_t(Sout,W,nombre,yax,iP,x_i,z_i,currentiFte,saveASmatlab,saveAstxt)
      use modelVars, only : NFREC, NPTSTIME, OMEI, &
                            t_vec, dt,Uo
      use glovars, only : z0
      use resultvars, only : allpoints
      use sourceVars, only : Po
      use wavelets
      use debugStuff
      implicit none
      complex*16, dimension(Nfrec+1), intent(in) :: W !espectro
      character(LEN=3)   :: nombre ! del componente
      character(LEN=100) :: titleN,yAx, CTIT ! texto para grafica
      integer ,intent(in) :: iP ! indice del receptor
      real*8, intent(in) :: x_i,z_i ! coordenadas del receptor
      integer, intent(in) :: currentiFte ! indice de la incidencia actual
      logical, intent(in) :: saveASmatlab,saveAstxt

      character(LEN=32)  :: name
      complex*16, dimension(NPTSTIME) :: S,Sout
      integer :: i,n_maxtime

      write(CTIT,'(a,F7.2,a,F7.2,a)')'(', x_i,' , ',z_i,')'

      S = z0
      S(1:nfrec+1)= W(1:nfrec+1:+1)
      S(NPTSTIME-NFREC+2:NPTSTIME) = conjg(W(nfrec:2:-1))

      ! (0) espectro sin corregir
      if (saveAstxt .eqv. .true.) then
        write(titleN,'(I0,a,I0,a,a,a,I0,a,I0,a,I0,a,I0,a)') &
               iP,'_in_',currentiFte,'_rawfrec_',nombre,'[', &
               int(x_i),'.',abs(int((x_i-int(x_i))*10)),';', &
               int(z_i),'.',abs(int((z_i-int(z_i))*10)),'].txt'
        write(name,'(a,I0,a,a,I0)') 'rawf_IP',iP,trim(nombre),'_in_',currentiFte
        OPEN(3201,FILE=trim(titleN),FORM="FORMATTED",ACTION='WRITE')
        write(3201,'(a,a,/,a,EN26.9,a,EN26.9,a,/,a)') '%',yAx, &
         '%(', allpoints(iP)%center%x,',',allpoints(iP)%center%z,')','%% '
        call showMNmatrixZ(int(nfrec +1),1,S(1:int(nfrec +1)),name,3201)
        close (3201)
       end if
!        call plotXYcomp(S(1:int(nfrec)),real(DFREC,4), &
!         int(nfrec),titleN, 'frec[hz] ',yAx, CTIT ,1200,800,0.0)
       if (saveASmatlab .eqv. .true.) then
        write(titleN,'(I0,a,I0,a,a,a,I0,a,I0,a,I0,a,I0,a)') &
               iP,'_in_',currentiFte,'_rawfrec_',nombre,'[', &
               int(x_i),'.',abs(int((x_i-int(x_i))*10)),';', &
               int(z_i),'.',abs(int((z_i-int(z_i))*10)),'].m'
        write(name,'(a,I0,a,a,I0)') 'rawf_IP',iP,trim(nombre),'_in_',currentiFte
        OPEN(3202,FILE=trim(titleN),FORM="FORMATTED",ACTION='WRITE')
        write(3202,'(a,a,/,a,EN26.9,a,EN26.9,a,/,a)') '%',yAx, &
         '%(', allpoints(iP)%center%x,',',allpoints(iP)%center%z,')','%% '
        call scripToMatlabMNmatrixZ(int(nfrec +1),1,S(1:int(nfrec +1)),name,3202)
        close (3202)
       end if

      !(1) señal en tiempo

!     if(ip.eq.3) then
!     print*,"*******"
!     DO i=1,NPTSTIME
!     print*,i,S(i),t_vec(i)
!     end do
!     end if

      S = S * t_vec ! t0 tiempo inicial
      S = S * Uo ! conv con fucion de amplitud
      S = FFTW(NPTSTIME,S,+1,1/(NPTSTIME*dt)) !backward (pasar al tiempo)

      ! (2) remover efecto de la frecuencia imaginaria (método DWN)
      !  La frecuencia imaginaria se utiliza como un amortiguamiento artificial
      !  que tiene el efecto de atenuar la amplitud de los arrivos de las
      !  fuentes periódicas en el espacio
      !  (surgen de truncar y discretizar una integral infinita)
      S = S * exp(- OMEI * Dt*((/(i,i=0, NPTSTIME-1)/))) !

      ! tiempo maximo para graficar
         n_maxtime = int(Po(currentiFte)%maxtime/dt)
         if(Po(currentiFte)%maxtime .lt. dt) n_maxtime = 2*nfrec
         if(Po(currentiFte)%maxtime .gt. NPTSTIME * real(dt,4)) n_maxtime = NPTSTIME
         S(n_maxtime+1: NPTSTIME) = z0;
         Sout = S

      !  (3) GRAFICAR en el tiempo
      if (saveAstxt .eqv. .true.) then
        write(titleN,'(I0,a,I0,a,a,a,I0,a,I0,a,I0,a,I0,a)') &
               iP,'_in_',currentiFte,'_sism_',nombre,'[', &
               int(x_i),'.',abs(int((x_i-int(x_i))*10)),';', &
               int(z_i),'.',abs(int((z_i-int(z_i))*10)),'].txt'
        write(name,'(a,I0,a,a,I0)') 'Sism_IP',iP,trim(nombre),'_in_',currentiFte
        OPEN(3203,FILE=trim(titleN),FORM="FORMATTED",ACTION='WRITE')
        write(3203,'(a,a,/,a,EN26.9,a,EN26.9,a,/,a)') '%',yAx, &
         '%(', allpoints(iP)%center%x,',',allpoints(iP)%center%z,')','%% '
        call showMNmatrixZ(NPTSTIME,1,S(1:int(NPTSTIME)),name,3203)
        close (3203)
       end if

       if (saveASmatlab .eqv. .true.) then
        write(titleN,'(I0,a,I0,a,a,a,I0,a,I0,a,I0,a,I0,a)') &
               iP,'_in_',currentiFte,'_sism_',nombre,'[', &
               int(x_i),'.',abs(int((x_i-int(x_i))*10)),';', &
               int(z_i),'.',abs(int((z_i-int(z_i))*10)),'].m'
        write(name,'(a,I0,a,a,I0)') 'Sism_IP',iP,trim(nombre),'_in_',currentiFte
        OPEN(3204,FILE=trim(titleN),FORM="FORMATTED",ACTION='WRITE')
        write(3204,'(a,a,/,a,EN26.9,a,EN26.9,a,/,a)') '%',yAx, &
         '%(', allpoints(iP)%center%x,',',allpoints(iP)%center%z,')','%% '
        call scripToMatlabMNmatrixZ(NPTSTIME,1,S(1:int(NPTSTIME)),name,3204)
        close (3204)
       end if

      !  (4) espectro correcto
!       ! tapper ?

        ! fft
        S = FFTW(NPTSTIME,S(1:int(NPTSTIME)),-1,dt) !forward

        if (saveAstxt .eqv. .true.) then
        write(titleN,'(I0,a,I0,a,a,a,I0,a,I0,a,I0,a,I0,a)') &
               iP,'_in_',currentiFte,'_f_',nombre,'[', &
               int(x_i),'.',abs(int((x_i-int(x_i))*10)),';', &
               int(z_i),'.',abs(int((z_i-int(z_i))*10)),'].txt'
        write(name,'(a,I0,a,a,I0)') 'f_IP',iP,trim(nombre),'_in_',currentiFte
        OPEN(3205,FILE=trim(titleN),FORM="FORMATTED",ACTION='WRITE')
        write(3205,'(a,a,/,a,EN26.9,a,EN26.9,a,/,a)') '%',yAx, &
         '%(', allpoints(iP)%center%x,',',allpoints(iP)%center%z,')','%% '
        call showMNmatrixZ(int(NPTSTIME/2+1),1,S(1:int(NPTSTIME/2+1)),name,3205)
        close (3205)
       end if

       if (saveASmatlab .eqv. .true.) then
        write(titleN,'(I0,a,I0,a,a,a,I0,a,I0,a,I0,a,I0,a)') &
               iP,'_in_',currentiFte,'_f_',nombre,'[', &
               int(x_i),'.',abs(int((x_i-int(x_i))*10)),';', &
               int(z_i),'.',abs(int((z_i-int(z_i))*10)),'].m'
        write(name,'(a,I0,a,a,I0)') 'f_IP',iP,trim(nombre),'_in_',currentiFte
        OPEN(3206,FILE=trim(titleN),FORM="FORMATTED",ACTION='WRITE')
        write(3206,'(a,a,/,a,EN26.9,a,EN26.9,a,/,a)') '%',yAx, &
         '%(', allpoints(iP)%center%x,',',allpoints(iP)%center%z,')','%% '
        call scripToMatlabMNmatrixZ(int(NPTSTIME/2+1),1,S(1:int(NPTSTIME/2+1)),name,3206)
        close (3206)
       end if

!     ! grafica logaritmica
!           write(titleN,'(a,a,I0,a,I0,a,I0,a,I0,a,I0,a)') &
!              'fL_',nombre,iP,'[', &
!              int(x_i),'.',abs(int((x_i-int(x_i))*10)),';', &
!              int(z_i),'.',abs(int((z_i-int(z_i))*10)),'].pdf'
!           logflag = 'logx     '
!           write(xAx,'(a)') "frec[Hz]"
!           call plotSpectrum(S(1: NPTSTIME/2),real(DFREC,4), NPTSTIME/2, NPTSTIME/2, &
!               titleN,xAx,yAx,logflag,1200,800,real(DFREC*(NFREC+1),4))

end subroutine W_to_t

subroutine sourceAmplitudeFunction
      use wavelets !las funciones: ricke
      use modelVars, only : Dt,Uo,dt!,maxtime
      use modelVars, only : nfrec, NPTSTIME,OMEI!DFREC
      use gloVars, only : ve => verbose,Ur,Ui,z0,Printnum
      use sourceVars, only:Po,iFte=>currentiFte
      !use ploteo10pesos
      implicit none
      integer  :: i,nval
!      character(LEN=9)          :: logflag
!      character(LEN=100)        :: titleN,xAx,yAx,CTIT

!      CHARACTER(len=32) :: arg
      logical :: inquire, lexist!, argumA
      real*8 :: val1,val2
      complex*16, dimension(NPTSTIME) :: S
      !complex*16 :: FFTW!(NPTSTIME)
      integer :: n_maxtime
!      argumA = .false.
!      CALL get_command_argument(1, arg)
!      IF (LEN_TRIM(arg) .ne. 0) then
!        if (trim(arg) .eq. '-a') then
!        iFte = 1
!        ve = 6
!        argumA = .true.
!        end if
!      end if

       n_maxtime = int(Po(iFte)%maxtime/dt)
       if(Po(iFte)%maxtime .lt. dt) n_maxtime = 2*nfrec
       if(Po(iFte)%maxtime .gt. NPTSTIME * real(dt,4)) n_maxtime = NPTSTIME

      ! grafica del factor de correción
      S = exp(- OMEI * Dt*((/(i,i=0, NPTSTIME-1)/)))
!      write(titleN,'(a)') 'x-LevantonDWN.pdf'
!      write(CTIT,'(a)') 'exp(-omei t)'
!      write(xAx,'(a)') 't[sec]'
!      write(yAx,'(a)') ''
!        call plotXYcomp(S(1:n_maxtime),real(Dt,4), &
!         int(n_maxtime),titleN,xAx,yAx, CTIT ,1200,800,0.0)

!     real :: factor
      ! Amplitude function of incident wave
      ! prepare the signal we will use as incident wave amplitude.
      ! el tamaño del ricker es 2*NFREC porque incluirá el conjugado

! 153
      Uo(:)=z0

      if(Po(iFte)%ampfunction .eq. 1) then
        call ricker(Uo(:),Po(iFte)%Ts,Po(iFte)%Tp) ! Ricker wavelet saved on Uo
!          if (ve .ge. 1) then
!            write(Printnum,'(a)')'   Incident wave amplitude function: Ricker'
!            write(titleN,'(a,I0,a)') 'x-amp',iFte,'-t.pdf'
!            write(CTIT,'(a)') 'WaveAmplitude of Ricker wavelet'
!            xAx = 'time[sec]'
!            write(yAx,'(a)') 'amplitude'
!            call plotXYcomp(Uo(1:n_maxtime),real(Dt,4), &
!                 n_maxtime,titleN,xAx,yAx,CTIT,1200,800,0.0)
!          end if
        ! forward
        Uo(:) = FFTW(NPTSTIME,Uo(:),-1,Dt)
        Uo(1) = 0
!          if (ve .ge. 1) then
!            write(titleN,'(a,I0,a)') 'x-amp',iFte,'-f.pdf'
!            xAx = 'frec[Hz] '
!            write(yAx,'(a)') 'amplitude'
!!           logflag = 'logx     '
!!           logflag = 'none     '
!            call plotXYcomp(Uo(:),real(DFREC,4), &
!                 size(Uo(:)),titleN,xAx,yAx,CTIT,1200,800,0.0)
!!           call plotSpectrum(Uo(:,iFte),real(DFREC,4), size(Uo(:,iFte)),int(size(Uo(:,iFte))/2), &
!!           titleN,xAx,yAx,logflag,1200,800,real(DFREC*(NFREC+1),4))
!          end if
     !********************************************************************************
      elseif(Po(iFte)%ampfunction .eq. 2) then ! Gaussian
        call gaussian(Uo(:),NPTSTIME,Po(iFte)%sigGaus)
!          if (ve .ge. 1) then
!           write(Printnum,'(a)')'   Incident wave amplitude function: Gaussian'
!            write(titleN,'(a,I0,a)') 'x-amp',iFte,'-f.pdf'
!            write(CTIT,'(a)') 'WaveAmplitude of Gaussian wavelet'
!            write(xAx,'(a)') 'frec[Hz] '
!            write(yAx,'(a)') 'amplitude'
!!           call plotXYcomp(Uo(:,iFte),real(DFREC,4), &
!!                n_maxtime,titleN,xAx,yAx,CTIT,1200,800,0.0)
!!           print*,Uo(:,iFte)
!            call plotXYcomp(Uo(:),real(DFREC,4), &
!                 size(Uo(:)),titleN,xAx,yAx,CTIT,1200,800,0.0)
!!           call plotSpectrum(Uo(:,iFte),real(DFREC,4), size(Uo(:,iFte)),int(size(Uo(:,iFte))/2), &
!!           titleN,xAx,yAx,logflag,1200,800,real(DFREC*(NFREC+1),4))
!          end if
      elseif(Po(iFte)%ampfunction .gt. 20) then ! GaussianCúbico
        call gaussian(Uo(:),NPTSTIME,Po(iFte)%sigGaus) !shift =  int((int(sigGaus) - sigGaus) * 1000)
        if (int((int(Po(iFte)%sigGaus) - Po(iFte)%sigGaus) * 1000) .gt. 0) &
        write(Printnum,'(a,I0)')'   One sided Gaussian filter with shift in points:',&
         int((int(Po(iFte)%sigGaus) - Po(iFte)%sigGaus) * 1000)
        Uo(:) = Uo(:) ** (Po(iFte)%ampfunction-20)
!          if (ve .ge. 1) then
!           write(Printnum,'(a,I0)')'   Incident wave amplitude function: Gaussian ',int((Po(iFte)%ampfunction-20))
!            write(titleN,'(a,I0,a)') 'x-amp',iFte,'-f.pdf'
!            write(CTIT,'(a,I0,a)') 'WaveAmplitude of Gaussian',int((Po(iFte)%ampfunction-20)),' wavelet'
!            write(xAx,'(a)') 'frec[Hz] '
!            write(yAx,'(a)') 'amplitude'
!!           call plotXYcomp(Uo(:,iFte),real(DFREC,4), &
!!                n_maxtime,titleN,xAx,yAx,CTIT,1200,800,0.0)
!!           print*,Uo(:,iFte)
!            call plotXYcomp(Uo(:),real(DFREC,4), &
!                 size(Uo(:)),titleN,xAx,yAx,CTIT,1200,800,0.0)
!!           call plotSpectrum(Uo(:,iFte),real(DFREC,4), size(Uo(:,iFte)),int(size(Uo(:,iFte))/2), &
!!           titleN,xAx,yAx,logflag,1200,800,real(DFREC*(NFREC+1),4))
!          end if
      elseif(Po(iFte)%ampfunction .eq. 3) then ! inDispl.txt
     !********************************************************************************
!        CALL chdir("..")
!        CALL chdir("ins")
        inquire(file="inAmplitude.txt",exist=lexist)
        if (lexist) then
        OPEN(77,FILE="inAmplitude.txt",FORM="FORMATTED")
        else
        write(6,'(a)') 'There is a missing input file. '
        stop 'Check "inAmplitude.txt" on this directory'
        end if
        READ(77,*)
        READ(77,*) inquire
        READ(77,*) nval
        if (nval .gt. NPTSTIME) then
        print*, "**** warning inAmplitude trimmed to ", NPTSTIME
        nval = NPTSTIME
        end if
        Uo(:) = 0
        if ( inquire) then
        do i = 1,nval
        READ(77,*) val1
        Uo(i) = UR * val1
        end do
        close(77)
!          if (ve .ge. 1) then
!            CALL chdir("..")
!            call chdir(trim(adjustl(rutaOut)))
!            write(Printnum,'(a)')'   Incident wave amplitude function from file'
!            write(titleN,'(a,I0,a)') 'x-amp',iFte,'-t.pdf'
!            write(CTIT,'(a)') 'WaveAmplitude'
!            xAx = 'time[sec]'
!            write(yAx,'(a)') 'amplitude'
!            call plotXYcomp(Uo(1:NPTSTIME),real(Dt,4), &
!                 n_maxtime,titleN,xAx,yAx,CTIT,1200,800,0.0)
!          end if

        ! forward
        Uo(:) = FFTW(NPTSTIME,Uo(:),-1,Dt)

!          if (ve .ge. 1) then
!            write(titleN,'(a,I0,a)') 'x-amp',iFte,'-f.pdf'
!            xAx = 'frec[Hz] '
!            write(yAx,'(a)') 'amplitude'
!            logflag = 'logx     '
!!           logflag = 'none     '
!            call plotSpectrum(Uo(:),real(DFREC,4), size(Uo(:)),int(size(Uo(:))/2), &
!            titleN,xAx,yAx,logflag,1200,800,real(DFREC*(NFREC+1),4))
!            CALL chdir("..")
!            CALL chdir("ins")
!          end if
        else
        do i = 1,nval
        READ(77,*) val1, val2
        Uo(i) = UR * val1 + UI * val2
        end do
!        if (ve .ge. 1) then
!            CALL chdir("..")
!            call chdir(trim(adjustl(rutaOut)))
!            write(titleN,'(a,I0,a)') 'x-amp',iFte,'-f.pdf'
!            xAx = 'frec[Hz] '
!            write(yAx,'(a)') 'amplitude'
!            logflag = 'logx     '
!!           logflag = 'none     '
!            call plotSpectrum(Uo(:),real(DFREC,4), size(Uo(:)),int(size(Uo(:))/2), &
!            titleN,xAx,yAx,logflag,1200,800,real(DFREC*(NFREC+1),4))
!            CALL chdir("..")
!            CALL chdir("ins")
!          end if
        close(77)
        end if

!        CALL chdir("..")
!        call chdir(trim(adjustl(rutaOut)))

      else  ! DIRAC -----------------------------------------
       if (ve .ge. 1) write(Printnum,'(a)')'   Incident wave amplitude function: Dirac delta'
        Uo(:)=UR
      end if !tipoPulso
      Uo(:) = Uo(:)*Po(iFte)%Escala !escala de la señal
!      write(Printnum,'(a)') &
!       "---------------------------------------------------------------------------------"
!      if (argumA .eqv. .true.) then
!      if (iFte .eq. nfuentes) then
!      stop "argumento -a"
!      else
!      iFte = iFte + 1
!      go to 153
!      end if
!      end if
end subroutine sourceAmplitudeFunction


text:
module mod_monolis
  use mod_monolis_prm
  use mod_monolis_com
  use mod_monolis_mat
  use mod_monolis_solve
  use mod_monolis_eigen_lanczos
  use mod_monolis_eigen_LOBPCG
  use mod_monolis_util
  use mod_monolis_matrix_copy
  use mod_monolis_sparse_util
  use mod_monolis_stdlib
  use mod_monolis_graph
  use mod_monolis_io
  use mod_monolis_io_mtx
  use mod_monolis_c2d3_shape
  use mod_monolis_c2d4_shape
  use mod_monolis_c2d6_shape
  use mod_monolis_c3d10_shape
  use mod_monolis_c3d4_shape
  use mod_monolis_c3d8_shape
  use mod_monolis_shape_util
  use mod_monolis_hash
  use mod_monolis_geom
  use mod_monolis_neighbor_search
  use mod_monolis_wrapper
end module mod_monolis

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis


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

      SUBROUTINE QUAL2 (MXND, MXCORN, MLN, NCORN, LCORN, LNODES, ICOMB,
     &   BNSIZE, ANGLE, LXN, ITEST, LTEST, QUAL, POSBL2, POSBL3, ROWCHN,
     &   ISTART, IEND)
C***********************************************************************

C  SUBROTINE QUAL2 = CHECKS THE QUALITY OF A SEMICIRCLE INTERPRETATION

C***********************************************************************

      DIMENSION LNODES (MLN, MXND), ANGLE (MXND), LCORN (MXCORN)
      DIMENSION BNSIZE (2, MXND)
      DIMENSION ICOMB (MXCORN), ITEST (2), LTEST (2), LXN (4, MXND)

      LOGICAL ERR, POSBL2, POSBL3, ROWCHN, SHRUNK

      REAL NICKS, NICKC

C ... See note below regarding bug...
      ISTEP = 0

C  ASSUME PERFECT QUALITY

      QUAL = 0.
      POSBL2 = .FALSE.
      POSBL3 = .FALSE.
      ROWCHN = .FALSE.

C  FIRST GET THE INTERVAL LENGTHS TO THE CHOSEN CORNERS

      ILEN = 2
      CALL SPACED (MXND, MXCORN, MLN, ILEN, NCORN, LCORN, LNODES, ICOMB,
     &   ITEST, LTEST, ERR)
      IF (ERR) RETURN

C  SEE IF A SEMICIRCLE INTERPRETATION IS POSSIBLE WITH
C  THESE INTERVALS

      IF ( (LTEST(1) .GE. 2) .AND. (LTEST(2) .GE. 2) ) THEN
         POSBL2 = .TRUE.
      ENDIF

C  NOT ADD UP THE NICKS FOR BAD ANGLES

      DO 100 I =1, NCORN
         NODE = LCORN (I)
         IF (ICOMB (I) .EQ. 1) THEN
            QUAL = QUAL + NICKC (ANGLE (NODE), LXN (1, NODE))
         ELSE
            QUAL = QUAL + NICKS (ANGLE (NODE), LXN (1, NODE))
         ENDIF
  100 CONTINUE

C  NOW SEE IF A TRIANGLE INTERPRETATION IS WARRANTED

      IF (LTEST (1) .GT. LTEST (2)) THEN
         I1 = ITEST (1)
         L1 = LTEST (1)
         I2 = ITEST (2)
         L2 = LTEST (2)
      ELSE
         I1 = ITEST (2)
         L1 = LTEST (2)
         I2 = ITEST (1)
         L2 = LTEST (1)
      ENDIF
      LDIF = (L1 - L2) / 2
      IF (LDIF .GT. L1 / 2) LDIF = L1 - LDIF

C  THIS TESTS THE FORCED TRIANGLE - THE NEW ROW MUST BE
C  ENDED AT A CURRENT SIDE NODE

      IF (L1 .EQ. L2) THEN
         NCHG1 = JUMPLP (MXND, MLN, LNODES, I1, LDIF)
         NCHG2 = JUMPLP (MXND, MLN, LNODES, I2, LDIF)
         QUAL13 = QUAL + NICKC (ANGLE (NCHG1), LXN (1, NCHG1))
         QUAL23 = QUAL + NICKC (ANGLE (NCHG2), LXN (1, NCHG2))
         IF ( (SHRUNK (BNSIZE (2, NCHG1), LNODES(8, NCHG1))) .AND.
     &      (SHRUNK (BNSIZE (2, NCHG2), LNODES(8, NCHG2))) ) THEN
            POSBL3 = .TRUE.
            POSBL2 = .FALSE.
            IF (QUAL13 .LT. QUAL23) THEN
               IEND = NCHG1
               ISTART = I1
            ELSE
               IEND = NCHG2
               ISTART = I2
            ENDIF
         ELSEIF (SHRUNK (BNSIZE (2, NCHG1), LNODES(8, NCHG1))) THEN
            POSBL3 = .TRUE.
            POSBL2 = .FALSE.
            IEND = NCHG1
            ISTART = I1
         ELSEIF (SHRUNK (BNSIZE (2, NCHG2), LNODES(8, NCHG2))) THEN
            POSBL3 = .TRUE.
            POSBL2 = .FALSE.
            IEND = NCHG2
            ISTART = I2
         ELSE
            POSBL3 = .FALSE.
         ENDIF

C  SEE IF THE ROW NEEDS ADJUSTED SO THAT A RECTANGLE REMAINS POSSIBLE
C  WITH A SIGNIFICANTLY REDUCED ELEMENT SIZE ON THE LONG SIDE

      ELSE
C ... There is a bug here since ISTEP is not defined
C     Since it has been 'kindof' working for several years,
C     we assume that ISTEP=0 will work the first time through
C     and we add that initialization to the beginning of this routine
         NCHG1 = JUMPLP (MXND, MLN, LNODES, I1, ISTEP)
         NCHG2 = JUMPLP (MXND, MLN, LNODES, I1, L1 - ISTEP)
         IF ( (SHRUNK (BNSIZE (2, NCHG1), LNODES(8, NCHG1))) .AND.
     &      (SHRUNK (BNSIZE (2, NCHG2), LNODES(8, NCHG2))) ) THEN
            ROWCHN = .TRUE.
            ISTART = I2
            IEND = I1

C  CHECK THE SIZE REDUCTIONS AND TRIANGLE INTERPRETATION

         ELSE
            DO 110 ISTEP = LDIF + 1, L1 / 2 - 1
               NCHG1 = JUMPLP (MXND, MLN, LNODES, I1, ISTEP)
               NCHG2 = JUMPLP (MXND, MLN, LNODES, I1, L1 - ISTEP)
               IF ( (SHRUNK (BNSIZE (2, NCHG1), LNODES(8, NCHG1))) .AND.
     &            (SHRUNK (BNSIZE (2, NCHG2), LNODES(8, NCHG2))) ) THEN
                  ROWCHN = .TRUE.
                  ISTART = I2
                  IEND = I1
                  GOTO 120
               ENDIF

  110       CONTINUE
  120       CONTINUE
         ENDIF
      ENDIF

      RETURN

      END

ext: f

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
      subroutine psplint(ncat,x,y)
      include "interp.h"
      integer ncat,klo,khi,k
      real*8 h,a,b,y,x
c
c  this started out as numerical recipies splint, uses
c  derivitives and points passed in splinearray through interp.h
c
      if(x .lt. splinearray(ncat,1,1) .or.
     $     x .gt. splinearray(ncat,1,npoints)) then
         print *,'interpolation out of range, for category ',ncat,
     $        ' asking for: ',x,'range is: ',splinearray(ncat,1,1),
     $        splinearray(ncat,1,npoints)
         stop
      endif

      klo=1
      khi=npoints
 1    if (khi-klo.gt.1) then
         k=(khi+klo)/2
         if(splinearray(ncat,1,k).gt.x)then
            khi=k
         else
            klo=k
         endif
         goto 1
      endif
      h=splinearray(ncat,1,khi)-splinearray(ncat,1,klo)
c     if (h.eq.0.) pause 'bad splinearray input.'
      a=(splinearray(ncat,1,khi)-x)/h
      b=(x-splinearray(ncat,1,klo))/h
      y=a*splinearray(ncat,2,klo)+b*splinearray(ncat,2,khi)+
     *     ((a**3-a)*splinearray(ncat,3,klo)+(b**3-b)
     $     *splinearray(ncat,3,khi))*(h**2)/6.
      return
      end

max_stars_repo_name: phaustin/empm

max_issues_repo_name: phaustin/empm

max_forks_repo_name: phaustin/empm


text:
Note: University facilities are subject to certain public access regulations...

 Open to all
  Haring Hall Haring 1301 (near southwest doors)  At the back of this UCD Bathrooms mens restroom is a single small and dirty shower. It has a shower curtain.
  Hickey Gym  There are showers you can access from the locker rooms (one in the womens locker room and one in the mens). The showers are also accessible from the Hickey Pool pool as well. Not really restricted since nobody checks your ID when you enter so...
  Plant and Environmental Sciences Building Plant and Environmental Sciences Building 1244B  this is a lactation station, with a small shower for women to rinse themselves (not a normal, fullsize shower)
  Swimming Swimming Pools
  Veihmeyer Hall  at the east end, in the restroom under the stairwell (no room number) is a single shower.

 Restricted access
  ARC  new, private, comfortable showers with plenty of traffic.
  King Hall  The basement bathrooms have single showers.  Rarely used and clean.
  Rec Pool  heavily used and not so nice.
  Sciences Laboratory Building  Off of the main hallway on the second floor one can find an individual shower. Its relatively new and very clean. It seems to be used mainly for breastfeeding. A sign on the door says This is not a public bathroom.
  Veterinary Medical Teaching Hospital  second floor, female locker room
  Davis Athletic Club  its got tons of showers, but you have to be a member

20110426 17:05:09 nbsp It seems bizarre that the most bicyclefriendly city in the country would not have even a single public shower for bicycle commuters. Users/SteveDavison

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!> Reference: https://refactoring.guru/design-patterns/visitor/go/example
module visitor_pattern

    implicit none
    private

    public :: square, circle, rectangle, area_calculator, middle_coordinates

    !> Two abstract classes

    type, abstract :: shape
    contains
        procedure(get_type_procedure), deferred :: get_type
        procedure(accept_procedure), deferred :: accept
    end type shape

    type, abstract :: visitor
    contains
        procedure(visit_procedure), deferred :: visit
    end type visitor

    abstract interface
        function get_type_procedure(self) result(result)
            import shape
            class(shape), intent(inout) :: self
            character(:), allocatable :: result
        end function get_type_procedure
        subroutine accept_procedure(self, v)
            import shape, visitor
            class(shape), intent(inout) :: self
            class(visitor), intent(inout) :: v
        end subroutine accept_procedure
        subroutine visit_procedure(self, s)
            import visitor, shape
            class(visitor), intent(inout) :: self
            class(shape), intent(inout) :: s
        end subroutine visit_procedure
    end interface

    !> Specific shapes

    type, extends(shape) :: square
        integer :: side
    contains
        procedure :: get_type => square_get_type
        procedure :: accept => square_accept
    end type square

    type, extends(shape) :: circle
        integer :: radius
    contains
        procedure :: get_type => circle_get_type
        procedure :: accept => circle_accept
    end type circle

    type, extends(shape) :: rectangle
        integer :: l
        integer :: b
    contains
        procedure :: get_type => rectangle_get_type
        procedure :: accept => rectangle_accept
    end type rectangle

    !> Specific visitors

    type, extends(visitor) :: area_calculator
        integer :: area
    contains
        procedure :: visit => area_calculator_visit
    end type area_calculator

    type, extends(visitor) :: middle_coordinates
        integer :: x, y
    contains
        procedure :: visit => middle_coordinates_visit
    end type middle_coordinates

contains

    function square_get_type(self) result(result)
        class(square), intent(inout) :: self
        character(:), allocatable :: result
        result = "Square"
    end function square_get_type

    function circle_get_type(self) result(result)
        class(circle), intent(inout) :: self
        character(:), allocatable :: result
        result = "Circle"
    end function circle_get_type

    function rectangle_get_type(self) result(result)
        class(rectangle), intent(inout) :: self
        character(:), allocatable :: result
        result = "Rectangle"
    end function rectangle_get_type

    subroutine square_accept(self, v)
        class(square), intent(inout) :: self
        class(visitor), intent(inout) :: v
        call v%visit(self)
    end subroutine square_accept

    subroutine circle_accept(self, v)
        class(circle), intent(inout) :: self
        class(visitor), intent(inout) :: v
        call v%visit(self)
    end subroutine circle_accept

    subroutine rectangle_accept(self, v)
        class(rectangle), intent(inout) :: self
        class(visitor), intent(inout) :: v
        call v%visit(self)
    end subroutine rectangle_accept

    subroutine area_calculator_visit(self, s)
        class(area_calculator), intent(inout) :: self
        class(shape), intent(inout) :: s
        select type (s)
        type is (square)
            print *, "Calculating area for square.🔥"
        type is (circle)
            print *, "Calculating area for circle.🔥"
        type is (rectangle)
            print *, "Calculating area for rectangle.🔥"
        end select
    end subroutine area_calculator_visit

    subroutine middle_coordinates_visit(self, s)
        class(middle_coordinates), intent(inout) :: self
        class(shape), intent(inout) :: s
        select type (s)
        type is (square)
            print *, "Calculating middle point coordinates for square.💠"
        type is (circle)
            print *, "Calculating middle point coordinates for circle.💠"
        type is (rectangle)
            print *, "Calculating middle point coordinates for rectangle.💠"
        end select
    end subroutine middle_coordinates_visit

end module visitor_pattern

max_stars_repo_name: zoziha/Fortran-Design-Patterns

max_issues_repo_name: zoziha/Fortran-Design-Patterns

max_forks_repo_name: zoziha/Fortran-Design-Patterns


text:
module fe_mesh

! The quad is enumerated as follows:
!
! 4       3
!  +-----+
!  |     |
!  |     |
!  +-----+
! 1       2
!
! The hexahedron is enumerated as follows
!
! The bottom face:             The top face:
! 4       3                      8       7
!  +-----+                        +-----+
!  |     |                        |     |
!  |     |                        |     |
!  +-----+                        +-----+
! 1       2                      5       6
!

use types
use utils, only: stop_error, newunit, zeros
use quadrature, only: gauss_pts, gauss_wts, lobatto_wts, lobatto_pts
use solvers, only: solve_sym
use feutils, only: phih
use isolve, only: solve_cg
use utils, only: assert
implicit none
private
public cartesian_mesh_2d, cartesian_mesh_3d, define_connect_tensor_2d, &
    define_connect_tensor_3d, c2fullc_2d, c2fullc_3d, fe2quad_2d, fe2quad_3d, &
    vtk_save, get_element_3d, fe_eval_xyz, line_save, fe2quad_3d_lobatto, &
    quad2fe_3d, quad2fe_3d_lobatto, cartesian_mesh_3d_mask

contains

subroutine cartesian_mesh_2d(nx, ny, p1, p2, nodes, elems)
! Returns a uniform 2D cartesian mesh of nx x ny elements
!
! The lower left corner is p1(2), the upper right corner is p2(2).
integer, intent(in) :: nx, ny ! the number of elements in x and y-direction
real(dp), intent(in) :: p1(:), p2(:)
real(dp), allocatable, intent(out) :: nodes(:, :)
integer, allocatable, intent(out) :: elems(:, :)
integer :: Ne, Nn, i, j, idx
Ne = nx*ny
Nn = (nx+1)*(ny+1)
allocate(nodes(2, Nn), elems(4, Ne))
! Loop over nodes:
idx = 1
do i = 1, nx + 1
    do j = 1, ny + 1
        nodes(:, idx) = (p2-p1) * [i-1, j-1] / [nx, ny] + p1
        idx = idx + 1
    end do
end do
! Loop over elements:
idx = 1
elems = 0
do i = 1, nx
    do j = 1, ny
        elems(1, idx) = p(i  , j  )
        elems(2, idx) = p(i+1, j  )
        elems(3, idx) = p(i+1, j+1)
        elems(4, idx) = p(i  , j+1)
        idx = idx + 1
    end do
end do

contains

    integer pure function p(i, j)
    integer, intent(in) :: i, j
    p = (i-1)*(ny+1) + j
    end function

end subroutine

subroutine cartesian_mesh_3d(nx, ny, nz, p1, p2, nodes, elems)
! Returns a uniform 3D cartesian mesh of nx x ny x nz elements
!
! The lower left front corner is p1(3), the upper right back corner is p2(3).
integer, intent(in) :: nx, ny, nz ! number of elements in x, y and z-directions
real(dp), intent(in) :: p1(:), p2(:)
real(dp), allocatable, intent(out) :: nodes(:, :)
integer, allocatable, intent(out) :: elems(:, :)
integer :: Ne, Nn, i, j, k, idx
Ne = nx*ny*nz
Nn = (nx+1)*(ny+1)*(nz+1)
allocate(nodes(3, Nn), elems(8, Ne))
! Loop over nodes:
idx = 1
do i = 1, nx + 1
    do j = 1, ny + 1
        do k = 1, nz + 1
            nodes(:, idx) = (p2-p1) * [i-1, j-1, k-1] / [nx, ny, nz] + p1
            idx = idx + 1
        end do
    end do
end do
! Loop over elements:
idx = 1
elems = 0
do i = 1, nx
    do j = 1, ny
        do k = 1, nz
            elems(1, idx) = p(i  , j  , k  )
            elems(2, idx) = p(i+1, j  , k  )
            elems(3, idx) = p(i+1, j+1, k  )
            elems(4, idx) = p(i  , j+1, k  )
            elems(5, idx) = p(i  , j  , k+1)
            elems(6, idx) = p(i+1, j  , k+1)
            elems(7, idx) = p(i+1, j+1, k+1)
            elems(8, idx) = p(i  , j+1, k+1)
            idx = idx + 1
        end do
    end do
end do

contains

    integer pure function p(i, j, k)
    integer, intent(in) :: i, j, k
    p = (i-1)*(ny+1)*(nz+1) + (j-1)*(nz+1) + k
    end function

end subroutine

subroutine define_connect_tensor_2d(nex, ney, p, ibc, gn)
! 2D connectivity table for tensor-product order p elements
integer, intent(in) :: nex, ney ! Number of elements in x and y directions
integer, intent(in) :: p ! Polynomial order of elements
! Boundary condition: 1 = Neumann, 2 = Dirichlet, 3 = periodic
integer, intent(in) :: ibc
! gn(i, j, e) is the global node number of the local (i,j) node in e-th element
integer, allocatable, intent(out) :: gn(:, :, :)
integer :: nodes(nex*p+1, ney*p+1)
integer :: inode, ix, iy, iel, iex, iey
! Construct array of global nodes
! 0 = no associated basis function as in e.g. Dirichlet BCs
nodes = 0
inode = 0
do iy = 1, ney*p+1
    do ix = 1, nex*p+1
        if (ibc >= 2 .and. (ix == 1 .or. ix == nex*p+1 .or. &
                            iy == 1 .or. iy == ney*p+1)) cycle
        inode = inode + 1
        nodes(ix, iy) = inode
    end do
end do
if (ibc == 3) then
    ! Now we need to connect the basis functions on the opposite sites of the
    ! boundary.
    ! top-bottom sides
    do ix = 2, nex*p
        inode = inode + 1
        nodes(ix, 1)       = inode
        nodes(ix, ney*p+1) = inode
    end do
    ! left-right sides
    do iy = 2, ney*p
        inode = inode + 1
        nodes(1, iy)       = inode
        nodes(nex*p+1, iy) = inode
    end do
    ! Corners
    inode = inode + 1
    nodes(1, 1)             = inode
    nodes(1, ney*p+1)       = inode
    nodes(nex*p+1, 1)       = inode
    nodes(nex*p+1, ney*p+1) = inode
end if
! Construct connectivity table of global nodes in each element
allocate(gn(p+1, p+1, nex*ney))
iel = 0
do iex = 1, nex
    do iey = 1, ney
        iel = iel + 1 ! Element number
        ix = (iex-1)*p+1 ! Lower left corner
        iy = (iey-1)*p+1
        ! Get global node numbers in element
        gn(:, :, iel) = nodes(ix:ix+p, iy:iy+p)
    end do
end do
end subroutine

subroutine define_connect_tensor_3d(nex, ney, nez, p, ibc, gn)
! 3D connectivity table for tensor-product order p elements
integer, intent(in) :: nex, ney, nez ! Number of elements in x, y, z directions
integer, intent(in) :: p ! Polynomial order of elements
! Boundary condition: 1 = Neumann, 2 = Dirichlet, 3 = periodic
integer, intent(in) :: ibc
! gn(i, j, k, e) is the global node number of the local (i,j,k) node
! in the e-th element:
integer, allocatable, intent(out) :: gn(:, :, :, :)
integer :: nodes(nex*p+1, ney*p+1, nez*p+1)
integer :: inode, ix, iy, iz, iel, iex, iey, iez
! Construct array of global nodes
! 0 = no associated basis function as in e.g. Dirichlet BCs
nodes = 0
inode = 0
do iz = 1, nez*p+1
    do iy = 1, ney*p+1
        do ix = 1, nex*p+1
            if (ibc >= 2 .and. (ix == 1 .or. ix == nex*p+1 .or. &
                iy == 1 .or. iy == ney*p+1 .or. &
                iz == 1 .or. iz == nez*p+1)) cycle
            inode = inode + 1
            nodes(ix, iy, iz) = inode
        end do
    end do
end do
if (ibc == 3) then
    ! Now we need to connect the basis functions on the opposite sites of the
    ! boundary.
    ! top-bottom faces (middle part)
    do ix = 2, nex*p
        do iy = 2, ney*p
            inode = inode + 1
            nodes(ix, iy, 1)       = inode
            nodes(ix, iy, nez*p+1) = inode
        end do
    end do
    ! top-bottom faces (4 edges connecting them)
    do iz = 2, nez*p
        inode = inode + 1
        nodes(1, 1, iz)             = inode
        nodes(1, ney*p+1, iz)       = inode
        nodes(nex*p+1, 1, iz)       = inode
        nodes(nex*p+1, ney*p+1, iz) = inode
    end do
    ! left-right faces (middle part)
    do iy = 2, ney*p
        do iz = 2, nez*p
            inode = inode + 1
            nodes(1, iy, iz)       = inode
            nodes(nex*p+1, iy, iz) = inode
        end do
    end do
    ! left-right faces (4 edges connecting them)
    do ix = 2, nex*p
        inode = inode + 1
        nodes(ix, 1, 1)             = inode
        nodes(ix, 1, nez*p+1)       = inode
        nodes(ix, ney*p+1, 1)       = inode
        nodes(ix, ney*p+1, nez*p+1) = inode
    end do
    ! front-back faces (middle part)
    do ix = 2, nex*p
        do iz = 2, nez*p
            inode = inode + 1
            nodes(ix, 1, iz)       = inode
            nodes(ix, ney*p+1, iz) = inode
        end do
    end do
    ! front-back faces (4 edges connecting them)
    do iy = 2, ney*p
        inode = inode + 1
        nodes(1, iy, 1)             = inode
        nodes(1, iy, nez*p+1)       = inode
        nodes(nex*p+1, iy, 1)       = inode
        nodes(nex*p+1, iy, nez*p+1) = inode
    end do
    ! Corners (vertices)
    inode = inode + 1
    nodes(1, 1, 1)                   = inode
    nodes(1, ney*p+1, 1)             = inode
    nodes(nex*p+1, 1, 1)             = inode
    nodes(nex*p+1, ney*p+1, 1)       = inode
    nodes(1, 1, nez*p+1)             = inode
    nodes(1, ney*p+1, nez*p+1)       = inode
    nodes(nex*p+1, 1, nez*p+1)       = inode
    nodes(nex*p+1, ney*p+1, nez*p+1) = inode
end if
! Construct connectivity table of global nodes in each element
allocate(gn(p+1, p+1, p+1, Nex*Ney*Nez))
iel = 0
do iex = 1, Nex
    do iey = 1, Ney
        do iez = 1, Nez
            iel = iel + 1 ! Element number
            ix = (iex-1)*p+1 ! Lower left front corner
            iy = (iey-1)*p+1
            iz = (iez-1)*p+1
            ! Get global node numbers in element
            gn(:, :, :, iel) = nodes(ix:ix+p, iy:iy+p, iz:iz+p)
        end do
    end do
end do
end subroutine

subroutine c2fullc_2d(in, ib, c, fullc)
! Converts FE coefficient vector to full coefficient vector
! It puts 0 for Dirichlet boundary conditions (ib==0), otherwise it just copies
! the coefficients.
integer, intent(in) :: in(:, :, :)
integer, intent(in) :: ib(:, :, :)
real(dp), intent(in) :: c(:) ! coefficient vector with regards to ib
real(dp), intent(out) :: fullc(:) ! full coefficients vector with regards to in
integer :: e, i, j
do e = 1, size(in, 3)
    do i = 1, size(in, 1)
        do j = 1, size(in, 2)
            if (ib(i, j, e) == 0) then
                fullc(in(i, j, e)) = 0 ! Dirichlet
            else
                fullc(in(i, j, e)) = c(ib(i, j, e))
            end if
        end do
    end do
end do
end subroutine

subroutine c2fullc_3d(in, ib, c, fullc)
! Converts FE coefficient vector to full coefficient vector
! It puts 0 for Dirichlet boundary conditions (ib==0), otherwise it just copies
! the coefficients.
integer, intent(in) :: in(:, :, :, :)
integer, intent(in) :: ib(:, :, :, :)
real(dp), intent(in) :: c(:) ! coefficient vector with regards to ib
real(dp), intent(out) :: fullc(:) ! full coefficients vector with regards to in
integer :: e, i, j, k
do e = 1, size(in, 4)
    do i = 1, size(in, 1)
        do j = 1, size(in, 2)
            do k = 1, size(in, 3)
                if (ib(i, j, k, e) == 0) then
                    fullc(in(i, j, k, e)) = 0 ! Dirichlet
                else
                    fullc(in(i, j, k, e)) = c(ib(i, j, k, e))
                end if
            end do
        end do
    end do
end do
end subroutine

subroutine fe2quad_2d(elems, xin, xiq, phihq, in, fullu, uq)
! Transforms fullu from FE-coefficient to quadrature-grid representation.
! fullu is a full FE coefficient vector, having values for all nodes in the
! mesh, including domain-boundary nodes.
integer, intent(in) :: elems(:, :)
real(dp), intent(in) :: xin(:)
real(dp), intent(in) :: xiq(:)
real(dp), intent(in) :: phihq(:, :)
integer, intent(in) :: in(:, :, :)
real(dp), intent(in) :: fullu(:)
real(dp), intent(out) :: uq(:, :, :)
integer :: ie, ilnx, ilny, iqx, iqy

! evaluate at quad points in each element
do ie = 1, size(elems, 2)
    uq(:, :, ie) = 0
    do ilnx = 1, size(xin)
    do ilny = 1, size(xin)
        do iqx = 1, size(xiq)
        do iqy = 1, size(xiq)
            uq(iqx, iqy, ie) = uq(iqx, iqy, ie) + &
                fullu(in(ilnx, ilny, ie)) * phihq(iqx, ilnx) * phihq(iqy, ilny)
        end do
        end do
    end do
    end do
end do
end subroutine

subroutine fe2quad_3d(elems, xin, xiq, phihq, in, fullu, uq)
! Transforms fullu from FE-coefficient to quadrature-grid representation.
! fullu is a full FE coefficient vector, having values for all nodes in the
! mesh, including domain-boundary nodes.
integer, intent(in) :: elems(:, :)
real(dp), intent(in) :: xin(:)
real(dp), intent(in) :: xiq(:)
real(dp), intent(in) :: phihq(:, :)
integer, intent(in) :: in(:, :, :, :)
real(dp), intent(in) :: fullu(:)
real(dp), intent(out) :: uq(:, :, :, :)
integer :: ie, ilnx, ilny, ilnz, iqx, iqy, iqz

! evaluate at quad points in each element
uq = 0
do ie = 1, size(elems, 2)
    do ilnz = 1, size(xin)
    do ilny = 1, size(xin)
    do ilnx = 1, size(xin)
        do iqz = 1, size(xiq)
        do iqy = 1, size(xiq)
        do iqx = 1, size(xiq)
            uq(iqx, iqy, iqz, ie) = uq(iqx, iqy, iqz, ie) + &
                fullu(in(ilnx, ilny, ilnz, ie)) * &
                phihq(iqx, ilnx) * phihq(iqy, ilny) * phihq(iqz, ilnz)
        end do
        end do
        end do
    end do
    end do
    end do
end do
end subroutine

subroutine quad2fe_3d(Ne, Nb, p, jac_det, wtq3, Sp, Sj, Sx, phi_v, in, ib, &
        uq, fullu)
! Transforms quadrature-grid representation to FE-coefficient (fullu).
! fullu is a full FE coefficient vector, having values for all nodes in the
! mesh, including domain-boundary nodes.
! Assumes the same 'jac_det' for each finite element.
use poisson3d_assembly, only: assemble_3d_coo_rhs
integer, intent(in) :: Ne, Nb, p
real(dp), intent(in) :: jac_det ! |J| for each element (must be constant)
real(dp), intent(in) :: wtq3(:, :, :)
integer, intent(in) :: Sp(:), Sj(:)
real(dp), intent(in) :: Sx(:)
real(dp), intent(in) :: phi_v(:, :, :, :, :, :)
real(dp), intent(in) :: uq(:, :, :, :)
integer, intent(in) :: in(:, :, :, :), ib(:, :, :, :)
real(dp), intent(out) :: fullu(:)
real(dp) :: rhs(Nb), sol(Nb)
call assemble_3d_coo_rhs(Ne, p, uq, jac_det, wtq3, ib, phi_v, rhs)
sol = solve_cg(Sp, Sj, Sx, rhs, zeros(size(rhs)), 1e-12_dp, 400)
call c2fullc_3d(in, ib, sol, fullu)
end subroutine

subroutine quad2fe_3d_lobatto(Ne, p, in, uq, fullu)
! Transforms quadrature-grid representation to FE-coefficient (fullu).
! fullu is a full FE coefficient vector, having values for all nodes in the
! mesh, including domain-boundary nodes.
! Assumes the same 'jac_det' for each finite element and spectral elements.
integer, intent(in) :: Ne, p
real(dp), intent(in) :: uq(:, :, :, :)
integer, intent(in) :: in(:, :, :, :)
real(dp), intent(out) :: fullu(:)
integer :: ie, iqx, iqy, iqz
! The quad points are the coefficients for spectral elements
do ie = 1, Ne
    do iqz = 1, p+1
    do iqy = 1, p+1
    do iqx = 1, p+1
        fullu(in(iqx, iqy, iqz, ie)) = uq(iqx, iqy, iqz, ie)
    end do
    end do
    end do
end do
end subroutine

subroutine fe2quad_3d_lobatto(elems, xiq, in, fullu, uq)
integer, intent(in) :: elems(:, :)
real(dp), intent(in) :: xiq(:)
integer, intent(in) :: in(:, :, :, :)
real(dp), intent(in) :: fullu(:)
real(dp), intent(out) :: uq(:, :, :, :)
integer :: ie, iqx, iqy, iqz
! evaluate at quad points in each element
do ie = 1, size(elems, 2)
    do iqz = 1, size(xiq)
    do iqy = 1, size(xiq)
    do iqx = 1, size(xiq)
        uq(iqx, iqy, iqz, ie) = fullu(in(iqx, iqy, iqz, ie))
    end do
    end do
    end do
end do
end subroutine

subroutine vtk_save(filename, nx, ny, nz, nodes, elems, xiq, fq)
character(len=*), intent(in) :: filename
integer, intent(in) :: nx, ny, nz
real(dp), intent(in):: nodes(:, :)
integer, intent(in) :: elems(:, :)
real(dp), intent(in) :: xiq(:)
real(dp), intent(in) :: fq(:, :, :, :)
integer :: Ne, e, iqx, iqy, iqz, i, j, k, Nq
real(dp), dimension(size(xiq)) :: x, y, z, xp, yp, zp
real(dp) :: qdata(nx*size(xiq), ny*size(xiq), nz*size(xiq))
real(dp) :: lx, ly, lz
real(dp) :: jacx, jacy, jacz, jac_det
integer :: u

Nq = size(xiq)

open(newunit(u), file=filename, status="replace")
write(u, "(a)") "# vtk DataFile Version 2.0"
write(u, "(a)") "Function at quadrature points"
write(u, "(a)") "ASCII"
write(u, "(a)") "DATASET RECTILINEAR_GRID"
write(u, "(a, i8, i8, i8)") "DIMENSIONS", nx*Nq, ny*Nq, nz*Nq

Ne = size(elems, 2)
lx = nodes(1, elems(7, 1)) - nodes(1, elems(1, 1)) ! Element sizes
ly = nodes(2, elems(7, 1)) - nodes(2, elems(1, 1))
lz = nodes(3, elems(7, 1)) - nodes(3, elems(1, 1))
jacx = lx/2
jacy = ly/2
jacz = lz/2
jac_det = abs(jacx*jacy*jacz)
xp = (xiq + 1) * jacx
yp = (xiq + 1) * jacy
zp = (xiq + 1) * jacz
! Note: this depends on the exact order 'elems' is constructed
write(u, "(a, ' ', i8, ' ', a)") "X_COORDINATES", nx*Nq, "float"
do e = 1, nx*ny*nz, ny*nz
    x = xp + nodes(1, elems(1, e))
    do iqz = 1, Nq
        write(u, *) x(iqz)
    end do
end do
write(u, "(a, ' ', i8, ' ', a)") "Y_COORDINATES", ny*Nq, "float"
do e = 1, ny*nz, nz
    y = yp + nodes(2, elems(1, e))
    do iqz = 1, Nq
        write(u, *) y(iqz)
    end do
end do
write(u, "(a, ' ', i8, ' ', a)") "Z_COORDINATES", nz*Nq, "float"
do e = 1, nz
    z = zp + nodes(3, elems(1, e))
    do iqz = 1, Nq
        write(u, *) z(iqz)
    end do
end do

e = 1
do i = 1, nx
do j = 1, ny
do k = 1, nz
    do iqx = 1, Nq
    do iqy = 1, Nq
    do iqz = 1, Nq
        qdata((i-1)*Nq + iqx, (j-1)*Nq + iqy, (k-1)*Nq + iqz) = &
            fq(iqx, iqy, iqz, e)
    end do
    end do
    end do
    e = e + 1
end do
end do
end do

write(u, "(a, ' ', i10)") "POINT_DATA", size(qdata)
write(u, "(a)") "SCALARS MyFunction float"
write(u, "(a)") "LOOKUP_TABLE default"
write(u, *) qdata

close(u)
end subroutine

subroutine line_save(filename, xin, nodes, elems, in, fullu, x1, x2, n)
character(len=*), intent(in) :: filename
real(dp), intent(in):: xin(:), nodes(:, :), fullu(:), x1(:), x2(:)
integer, intent(in):: elems(:, :), in(:, :, :, :), n
integer :: i, u
real(dp) :: x(3), val
open(newunit(u), file=filename, status="replace")
do i = 1, n
    x = x1 + (i-1) * (x2-x1) / (n-1)
    val = fe_eval_xyz(xin, nodes, elems, in, fullu, x)
    write(u, *) x, val
end do
close(u)
end subroutine

real(dp) function fe_eval_xyz(xin, nodes, elems, in, fullu, x) result(val)
real(dp), intent(in):: xin(:), nodes(:, :), x(:), fullu(:)
integer, intent(in):: elems(:, :), in(:, :, :, :)
real(dp), dimension(3) :: l, jac, xi, xa, xb
integer :: nx, ny, nz, p, e

e = get_element_3d(nodes, elems, x)

! xa, xb are two opposite corners of the hexahedron
xa = nodes(:, elems(1, e))
xb = nodes(:, elems(7, e))
! l is the diagonal vector:
l = xb - xa
! Assume hexahedral shape:
jac = l / 2
! xi is in the reference domain [-1, 1]^3
xi = (x-xa) / jac - 1

p = size(xin) - 1
val = 0
do nz = 1, p+1
do ny = 1, p+1
do nx = 1, p+1
    val = val + fullu(in(nx, ny, nz, e)) * &
        phih(xin, nx, xi(1)) * phih(xin, ny, xi(2)) * phih(xin, nz, xi(3))
end do
end do
end do
end function

integer function get_element_3d(nodes, elems, x) result(e)
! Returns the element that containst the point 'x'
real(dp), intent(in) :: nodes(:, :), x(:)
integer, intent(in) :: elems(:, :)
real(dp) :: xa(3), xb(3)
integer :: i, Ne
Ne = size(elems, 2)
do i = 1, Ne
    ! xa, xb are two opposite corners of the hexahedron
    xa = nodes(:, elems(1, i))
    xb = nodes(:, elems(7, i))
    if (point_in_hex(xa, xb, x)) then
        e = i
        return
    end if
end do
e = 0
call stop_error("The point 'x' is not inside the mesh.")
end function

logical pure function point_in_hex(p1, p2, x) result(r)
! Returns .true. if the point 'x' is in the hexahedron specified by the two
! opposite corners p1 (lower, left, front) and p2 (upper, right, back).
real(dp), intent(in) :: p1(:), p2(:), x(:)
r = all(p1 <= x) .and. all(x <= p2)
end function

subroutine cartesian_mesh_3d_mask(myid, nx, ny, nz, nsubx, nsuby, nsubz, &
        mask_elems)
! Returns a mask for a subdomain "myid" within a uniform 3D cartesian mesh of
! (nx, ny, nz) elements, where subdomains have (nsubx, nsuby, nsubz) elements.
! Total number of elements in the domain in each direction
integer, intent(in) :: nx, ny, nz
! number of subdomains elements in x, y and z-directions
integer, intent(in) :: nsubx, nsuby, nsubz
! Processor id = 0, 1, 2, ...
integer, intent(in) :: myid
integer, allocatable, intent(out) :: mask_elems(:)
integer :: Ne, i, j, k, idx
integer :: nesubx, nesuby, nesubz
integer :: iminmax(6)
call assert(modulo(nx, nsubx) == 0)
call assert(modulo(ny, nsuby) == 0)
call assert(modulo(nz, nsubz) == 0)
Ne = nx*ny*nz
nesubx = nx / nsubx
nesuby = ny / nsuby
nesubz = nz / nsubz
allocate(mask_elems(Ne))

iminmax = get_minmax(myid+1, nsubx, nsuby, nsubz, nesubx, nesuby, nesubz)

mask_elems = 0
idx = 1
do i = 1, nx
    do j = 1, ny
        do k = 1, nz
            if     (i >= iminmax(1) .and. i <= iminmax(2) .and. &
                    j >= iminmax(3) .and. j <= iminmax(4) .and. &
                    k >= iminmax(5) .and. k <= iminmax(6)) then
                mask_elems(idx) = 1
            end if
            idx = idx + 1
        end do
    end do
end do
end subroutine

function get_minmax(id, nsubx, nsuby, nsubz, nesubx, nesuby, nesubz) result(r)
! Calculates element indices belonging to the subdomain 'id'
integer, intent(in) :: id ! The subdomain ID = 1, 2, 3, ...
! number of subdomains in each direction:
integer, intent(in) :: nsubx, nsuby, nsubz
! number of elements in a subdomain in each direction:
integer, intent(in) :: nesubx, nesuby, nesubz
integer :: r(6) ! r = [iminx, imaxx, iminy, imaxy, iminz, imaxz]
integer :: i, j, k, idx
idx = 1
do i = 1, nsubx
    do j = 1, nsuby
        do k = 1, nsubz
            if (idx == id) then
                r = [(i-1)*nesubx + 1, i*nesubx, &
                     (j-1)*nesuby + 1, j*nesuby, &
                     (k-1)*nesubz + 1, k*nesubz]
                return
            end if
            idx = idx + 1
        end do
    end do
end do
call stop_error("get_minmax(): 'id' not found")
end function

end module

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
!-------------------------------------------------------------
!--  Modul for the Saving Control Values
!-------------------------------------------------------------
 MODULE Control_Mod
   USE Kind_Mod
!
!-----------------------------------------------------------------
!---  Scenario
!-----------------------------------------------------------------
!
!--- Identifier for scenario
      CHARACTER(80) :: Bsp        = ''      ! Identifier for scenario

!--- Files
      CHARACTER(80) :: RunFile    = ''    & ! Simulation data file
&                    , MetFile    = ''    & ! Meteorology file
&                    , SysFile    = ''    & ! Chemical mechanism
&                    , ChemFile   = ''    & ! Chemical mechenism (spaccim input)
&                    , InitFile   = ''    & ! Initial concentrations
&                    , DataFile   = ''    & ! Gas and Aqueous DATA
&                    , MWFile     = ''    & ! molecular weights of species
&                    , NetcdfFile = ''    & ! NetCDF output file
&                    , ODEsolver  = ''    & ! Method for Rosenbrock Integration
&                    , TargetFile = ''      ! file for reductions analysis (target species)

      CHARACTER(7)  :: OutputPath   = 'OUTPUT/'  ! path to output folder
      CHARACTER(80) :: FluxMetaFile = ''         ! meta data for unformatted flux data
      CHARACTER(80) :: FluxFile     = ''         ! flux data (unformatted)

      REAL(dp) :: Red_TStart, Red_TEnd     ! time interval for redcution procedure
      REAL(dp) :: eps_red                  ! threshold for reduction procedure
!
!--- Unit Numbers
      INTEGER, PARAMETER :: RunUnit      = 101  &
&                         , MetUnit      = 102  &
&                         , SysUnit      = 103  &
&                         , ChemUnit     = 104  &
&                         , MWUnit       = 105  &
&                         , InitUnit     = 106  &
&                         , DataUnit     = 107  &
&                         , FluxMetaUnit = 109  &
&                         , FluxUnit     = 110  &
&                         , TikZUnit     = 111
!
!-- Set Levels and Parameters for Processes
      REAL(dp) :: LWCLevelmin    & ! Lower level for LWC
&               , LWCLevelmax    & ! Upper level for LWC
&               , Temperature0   & ! Initial temperature
&               , Pressure0        ! Initial pressure for combustion mechanisms

!-- print matrices of the system (alpha, beta, miter, lu_miter, permutvector,..)
      LOGICAL :: MatrixPrint     & ! print certain matrices to file
&              , DebugPrint      & ! debugging rosenbrock steps
&              , NetCdfPrint     & ! print out concs and other stuff to netcdf file
&              , constLWC        & ! true if lwc value is fixed
&              , Teq             & ! if Teq=.TRUE. simulate combustion mechanism
&              , pHSet           & ! Initial pH by charge balance (1=on, 0=off)
&              , WaitBar         & ! ladebalken im terminal bei simulation (=1, default=0)
&              , FluxDataPrint   & ! writing flux data and analyse after simulaiton -> print new reaction file
&              , Simulation      & ! calculation of species concentration
&              , Reduction         ! reduction of chemical species and reactions

      LOGICAL :: KPP_Conversion

      INTEGER :: Error_Est         ! error estimation 1 = inf norm  , 2 = euklid norm

!--- Output of reaction fluxes
      REAL(dp)                    :: StpFlux
      INTEGER                     :: iStpFlux ! number of written flux steps

      INTEGER,          PARAMETER :: newReac_nr    = 1977
      CHARACTER(LEN=*), PARAMETER :: newReac_name  = 'BigRates.dat'
      INTEGER, ALLOCATABLE        :: newReac_List(:)
!-----------------------------------------------------------------
!---  Times
!-----------------------------------------------------------------
!-----------------------------------------------------------------
!--- Times:  < 0 in seconds, > 0 in hours, = 0  meteorology
      REAL(dp) :: tBegin              & ! Model start time
&               , tEnd              & ! Model end time
&               , StpNetcdf           ! Time step for Netcdf output
      INTEGER  :: nOutP

!--- NetCDF globals
      INTEGER, ALLOCATABLE :: iNcdfGas(:), iNcdfAqua(:), iNcdfSolid(:), iNcdfParti(:)
      INTEGER, ALLOCATABLE :: iNCout_G(:), iNCout_A_l(:), iNCout_A_m3(:), iNCout_S(:), iNCout_P(:)
      INTEGER  :: nNcdfGas=0, nNcdfAqua=0, nNcdfSolid=0, nNcdfParti=0 ! number of output spc for each phase
!

!---  Photolysis
      INTEGER :: iDate               ! Current date
      REAL(dp) :: rlat, rlon          ! Latitude, Longitude

!---  Liquid water content parameter

      REAL(dp) :: LWCconst=2.0d-8       ! [l/m3] no cloud
      REAL(dp) :: NCC=1000.0d0

!--  dust factor for damping of photolysis rates, measured JNO2
      REAL(dp) :: Dust = 1.0d0
      REAL(dp) :: minStp = 1.0d-25
      REAL(dp) :: maxStp = 100.0d0

!-- initialize timers
      REAL(dp) :: Timer_Start=0.0d0, Timer_Finish=0.d0
      REAL(dp) :: Tspan(2)
      REAL(dp) :: Tspan_tot(2)

      REAL(dp) :: Time_Read=0.0d0
      REAL(dp) :: TimeRates=0.0d0
      REAL(dp) :: TimeSymbolic=0.0d0
      REAL(dp) :: TimeFac=0.0d0
      REAL(dp) :: TimeSolve=0.0d0
      REAL(dp) :: TimeJac=0.0d0
      REAL(dp) :: TimeNetCDF=0.0d0
      REAL(dp) :: TimeErrCalc=0.0d0
      REAL(dp) :: TimeFluxWrite=0.0d0
      REAL(dp) :: TimeRhsCalc=0.0d0
      REAL(dp) :: TimeReduction=0.0d0

      REAL(dp) :: TimeIntegration=0.0d0
      REAL(dp) :: TimeRateA=0.0d0
      REAL(dp) :: TimeRateE=0.0d0
      REAL(dp) :: TimeJacobianA=0.0d0
      REAL(dp) :: TimeJacobianE=0.0d0
      REAL(dp) :: TimeNetCDFA=0.0d0
!

!--- type for some statistics
      TYPE Output_T
        REAL(dp), ALLOCATABLE :: y(:)    ! y-vector at Tend
        INTEGER :: nsteps     = 0              ! # succ. steps
        INTEGER :: nfailed    = 0              ! # failed steps
        INTEGER :: nRateEvals = 0              ! # Rate evaluation
        INTEGER :: npds       = 0              ! # Jacobian evaluation
        INTEGER :: ndecomps   = 0              ! # LU factorisation
        INTEGER :: nsolves    = 0              ! # solved lin algebra
        REAL(dp) :: Ttimestep = 0.0d0  ! mean Time for one ROW step
      END TYPE Output_T

      TYPE(Output_T) :: Out
!-----------------------------------------------------------------
!---  Numerics
!-----------------------------------------------------------------
!
!
!--- Tolerances for ROW Scheme
      REAL(dp) :: RtolROW        & ! Relative tolerance for ROW method
&               , AtolGas        & ! Absolute tolerance for gas phase
&               , AtolAqua       & ! Absolute tolerance for liquid phase
&               , AtolTemp         ! Absolute tolerance for Temperatur

      INTEGER  :: Ordering     ! default = 8
      INTEGER  :: ParOrdering  ! < 0 serial ordering , 0,1,2 parallel

      REAL(dp), ALLOCATABLE :: ATolAll(:)

!--- Logical variable for PI setsize controler
      LOGICAL :: PI_StepSize

!--- Control Parameter for solving linear algebra
      ! 'cl' for classic jacobian calc or
      ! 'ex' for extended matrix without calc of jac
      CHARACTER(2) :: LinAlg = '??'

!
!-----------------------------------------------------------------
!---  Set Constants and Unit Conversion
!-----------------------------------------------------------------
!
!---  constants
    REAL(dp), PARAMETER :: HOUR       = 3600.0             &
&                        , hourday    = 86400.0            &
&                        , secday     = 4.32d04            &
&                        , mONE       = -1.d0              &
&                        , Pi         = 4.0d0*ATAN(1.0d0)  &
&                        , DR         = Pi / 180.d0        &
&                        , PiHalf     = 2.0d0*ATAN(1.0d0)  &
&                        , Pi34       = 3.0d0/4.0d0/Pi     &
&                        , eps        = EPSILON(1.0d0)     &  ! such that 1+eps>1 with working precision
&                        , small      = TINY(1.0d0)        &  ! smallest pos. real value
&                        , big        = HUGE(1.0d0)        &  ! largest pos. real value
&                        , epsY       = 1.0d-7             &
&                        , SI_am      = 1.66053892173d-27  &  ! Atomic mass unit  [kg]
&                        , SI_na      = 6.0221412927d+23   &  ! Avogadro's number [1/mol]
&                        , SI_kB      = 1.380648813d-23    &  ! Bolzmann constant [J/K]
&                        , SI_Gas     = SI_na * SI_kB         ! Gas constant      [J/mol/K]
!
!--- Real number constants
  REAL(dp), PARAMETER ::   ZERO    =     0.0d0   &
&                      ,  ONE     =     1.0d0   &
&                      ,  TWO     =     2.0d0  ,   rTWO     =   ONE/tWO     &
&                      ,  THREE   =     3.0d0  ,   rTHREE   =   ONE/THREE   &
&                      ,  FOUR    =     4.0d0  ,   rFOUR    =   ONE/FOUR    &
&                      ,  FIVE    =     5.0d0  ,   rFIVE    =   ONE/FIVE    &
&                      ,  SIX     =     6.0d0  ,   rSIX     =   ONE/SIX     &
&                      ,  SEVEN   =     7.0d0  ,   rSEVEN   =   ONE/SEVEN   &
&                      ,  EIGHT   =     8.0d0  ,   rEIGHT   =   ONE/EIGHT   &
&                      ,  NINE    =     9.0d0  ,   rNINE    =   ONE/NINE    &
&                      ,  TEN     =    10.0d0  ,   rTEN     =   ONE/TEN     &
&                      ,  ELEVN   =    11.0d0  ,   rELEVN   =   ONE/ELEVN   &
&                      ,  TWELV   =    12.0d0  ,   rTWELV   =   ONE/TWELV   &
&                      , TWENTY   =    20.0d0  ,   rTWENTY  =   ONE/TWENTY  &
&                      , mTHIRTY  =   -30.0d0  &
&                      , rm300    = mONE/300.d0,   r300     =   ONE/300.d0
!
!--- Orders of magnitude
  REAL(dp), PARAMETER ::  nano    =     1.0d-09    &
&                      , micro    =     1.0d-06    &
&                      , milli    =     1.0d-03    &
&                      , kilo     =     1.0d+03    &
&                      , mega     =     1.0d+06    &
&                      , giga     =     1.0d+09
!
!--- Natural logarithms
  REAL(dp), PARAMETER ::   ln10   =     LOG(TEN)    &
&                      ,  rln10   = ONE/LOG(TEN)
!
!--- minimum values if there is no sun
  REAL(dp), PARAMETER ::    EyChiZmin  =  9.357d-14
!
!
!--- Unit Conversion constants
!
!     Pressure
  REAL(dp), PARAMETER :: bar_to_dyncm2 = 1.0d06
  REAL(dp), PARAMETER :: dyncm2_to_Pa  = 1.0d-01
  REAL(dp), PARAMETER :: Pa_to_dyncm2  = 1.0d+01
  REAL(dp), PARAMETER :: bar_to_Pa     = 1.0d05
  REAL(dp), PARAMETER :: atm_to_Pa     = 101325.0d0
!
!     Energy
  REAL(dp), PARAMETER :: cal_to_joule  = 4.184d0
  REAL(dp), PARAMETER :: joule_to_cal  = ONE / cal_to_joule
  REAL(dp), PARAMETER :: joule_to_kcal = milli * joule_to_cal
  REAL(dp), PARAMETER :: kcal_to_joule = kilo * cal_to_joule
  REAL(dp), PARAMETER :: joule_to_erg  = TEN * mega
  REAL(dp), PARAMETER :: erg_to_joule  = rTEN * micro
!
!--- Physical constants ********************************************

!     Universal gas constant [J / mol / K]
  REAL(dp), PARAMETER :: R         = 8.31446210000000D0
  REAL(dp), PARAMETER :: rR        = ONE/R

!
!     Universal gas constant, calorie units [cal / mol / K]
  REAL(dp), PARAMETER :: Rcal      = R * joule_to_cal
  REAL(dp), PARAMETER :: rRcal     = ONE/Rcal
!
!     Universal gas constant, CGS units [erg / mol / K]
  REAL(dp), PARAMETER :: Rerg      = R * joule_to_erg
  REAL(dp), PARAMETER :: rRerg      = ONE/Rerg
!
!     Standard pressure [Pa]
  REAL(dp), PARAMETER :: Patm      = atm_to_Pa
  REAL(dp), PARAMETER :: rPatm     = ONE/Patm

!--- Unit Conversion
    INTEGER :: GasUnit, AquaUnit, GasRateUnit
    REAL(dp) :: GasFac

    INTEGER :: UnitGas=0
    INTEGER :: UnitAqua=0

!-----------------------------------------------------------------
!---  Output control
!-----------------------------------------------------------------

!
!-- input chemical reaction mechnism
    INTEGER, PARAMETER :: LenLine=400
    INTEGER, PARAMETER :: LenName=100
    INTEGER, PARAMETER :: LenType=20

!-- logicals for classic or extended matrix case
    LOGICAL :: CLASSIC  = .FALSE.
    LOGICAL :: EXTENDED = .FALSE.

    LOGICAL :: useMUMPS = .FALSE.
    LOGICAL :: useSparseLU = .FALSE.

    LOGICAL :: ChemKin = .FALSE.

    REAL(dp), ALLOCATABLE :: integrated_rates(:)
    REAL(dp), ALLOCATABLE :: mixing_ratios_spc(:,:)

!-- Type declaration
    TYPE List
      INTEGER,   ALLOCATABLE  :: List(:)
      Real(dp),  ALLOCATABLE  :: ListE(:)
      Real(dp),  ALLOCATABLE  :: ListP(:)
      INTEGER                :: len
    END TYPE List

    TYPE Chain
      CHARACTER(100), ALLOCATABLE :: sName(:,:)
      INTEGER,        ALLOCATABLE :: sIdx(:,:)
      Type(List),     ALLOCATABLE :: rIdx(:)
    END TYPE Chain

    TYPE Families_T
      CHARACTER(LenName), ALLOCATABLE :: Name(:)
      INTEGER,            ALLOCATABLE :: Index(:)
    END TYPE Families_T

    INTEGER, ALLOCATABLE :: maxErrorCounter(:)

 END MODULE Control_Mod

max_stars_repo_name: KarlJohnsonnn/AtCSol

max_issues_repo_name: KarlJohnsonnn/AtCSol

max_forks_repo_name: KarlJohnsonnn/AtCSol


text:
!! R624 image-selector
!    is lbracket cosubscript-list rbracket
!
!! R625 cosubscript
!    is scalar-int-expr
!
a[3] = 1
b = a[3]

! NOTE 6.17
A(:)[1,4] = 1
b = A(:)[1,4]

END

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
C
C***********************************************************************
C                                                                      *
C            NACHOS II - A Finite Element Computer Program             *
C                        for Incompressible Flow Problems              *
C                                                                      *
C     Copyright (c) 1986,2019   National Technology & Engineering      *
C                               Solutions of Sandia, LLC (NTESS)       *
C                                                                      *
C                            All rights reserved.                      *
C                                                                      *
C     This software is distributed under the BSD 3-Clause License.     *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE DEFORM (KIND,X,Y,U,V,SCNDII)
C
C     ******************************************************************
C
C     SUBROUTINE TO EVALUATE THE INVARIANTS OF THE RATE OF DEFORMATION
C     TENSOR AT THE ELEMENT INTEGRATION POINTS
C
C     ******************************************************************
C
      COMMON /PRBDAT/ IAXSYM,ITMDEP,IFORCE,IFREE,IVAR1,IVAR2,IPFUNC,
     1                IPNLTY ,PNLTY
      COMMON /GAUSS1/ GSPT1(4),GSPT2(4)
      COMMON /HAMMR1/ HAMPT1(3),HAMPT2(3)
      COMMON /TRI6/   F6Q(6,7),F6L(3,7),DF6QDS(6,7),DF6QDT(6,7)
      COMMON /QUAD8/  F8Q(8,9),F8L(4,9),DF8QDS(8,9),DF8QDT(8,9)
      COMMON /QUAD9/  F9Q(9,9),F9L(4,9),DF9QDS(9,9),DF9QDT(9,9)
C
      DIMENSION X(9), Y(9), U(9), V(9), SCNDII(4)
C
      DATA NGAUSS/4/,NHAM/3/
C
C     STATEMENT FUNCTIONS FOR THE JACOBIAN OF THE PARAMETRIC MAPPING
C     AND OTHER ELEMENT FUNCTIONS
C
      F11Q(S,T)=.25*(D1+S*D2+T*D3+S*T*D4+S**2*D5+S**2*T*D6)
      F12Q(S,T)=-.25*(B1+T*B2+S*B3+S*T*B4+T**2*B5+S*T**2*C6)
      F21Q(S,T)=-.25*(C1+S*C2+T*C3+S*T*C4+S**2*C5+S**2*T*B6)
      F22Q(S,T)=.25*(A1+T*A2+S*A3+S*T*A4+T**2*A5+S*T**2*A6)
      DETJQ(S,T)=F11Q(S,T)*F22Q(S,T)-F12Q(S,T)*F21Q(S,T)
C
      F11T(S,T)=S*A1+T*A2+A3
      F12T(S,T)=-(S*C1+T*C2+C3)
      F21T(S,T)=-(S*B1+T*B2+B3)
      F22T(S,T)=S*D1+T*D2+D3
      DETJT(S,T)=F11T(S,T)*F22T(S,T)-F12T(S,T)*F21T(S,T)
C
      U6Q(L)=F6Q(1,L)*U(1)+F6Q(2,L)*U(2)+F6Q(3,L)*U(3)+F6Q(4,L)*U(4)+
     1F6Q(5,L)*U(5)+F6Q(6,L)*U(6)
      U8Q(L)=F8Q(1,L)*U(1)+F8Q(2,L)*U(2)+F8Q(3,L)*U(3)+F8Q(4,L)*U(4)+
     1F8Q(5,L)*U(5)+F8Q(6,L)*U(6)+F8Q(7,L)*U(7)+F8Q(8,L)*U(8)
      U9Q(L)=F9Q(1,L)*U(1)+F9Q(2,L)*U(2)+F9Q(3,L)*U(3)+F9Q(4,L)*U(4)+
     1F9Q(5,L)*U(5)+F9Q(6,L)*U(6)+F9Q(7,L)*U(7)+F9Q(8,L)*U(8)+
     2F9Q(9,L)*U(9)
C
      RAD6L(L)=F6L(1,L)*X(1)+F6L(2,L)*X(2)+F6L(3,L)*X(3)
      RAD8L(L)=F8L(1,L)*X(1)+F8L(2,L)*X(2)+F8L(3,L)*X(3)+F8L(4,L)*X(4)
      RAD9L(L)=F9L(1,L)*X(1)+F9L(2,L)*X(2)+F9L(3,L)*X(3)+F9L(4,L)*X(4)
C
      RAD6Q(L)=F6Q(1,L)*X(1)+F6Q(2,L)*X(2)+F6Q(3,L)*X(3)+F6Q(4,L)*X(4)+
     1F6Q(5,L)*X(5)+F6Q(6,L)*X(6)
      RAD8Q(L)=F8Q(1,L)*X(1)+F8Q(2,L)*X(2)+F8Q(3,L)*X(3)+F8Q(4,L)*X(4)+
     1F8Q(5,L)*X(5)+F8Q(6,L)*X(6)+F8Q(7,L)*X(7)+F8Q(8,L)*X(8)
      RAD9Q(L)=F9Q(1,L)*X(1)+F9Q(2,L)*X(2)+F9Q(3,L)*X(3)+F9Q(4,L)*X(4)+
     1F9Q(5,L)*X(5)+F9Q(6,L)*X(6)+F9Q(7,L)*X(7)+F9Q(8,L)*X(8)+
     2F9Q(9,L)*X(9)
C
C     ******************************************************************
C
C     SELECT ELEMENT TYPE
C
      GO TO (110, 120, 10, 20, 10, 30), KIND
C
C     QUADRILATERAL ELEMENTS
C
   10 CONTINUE
      A1=-X(1)+X(2)+X(3)-X(4)
      A2=X(1)-X(2)+X(3)-X(4)
      A3=0.
      A4=0.
      A5=0.
      A6=0.
      B1=-Y(1)+Y(2)+Y(3)-Y(4)
      B2=Y(1)-Y(2)+Y(3)-Y(4)
      B3=0.
      B4=0.
      B5=0.
      B6=0.
      C1=-X(1)-X(2)+X(3)+X(4)
      C2=X(1)-X(2)+X(3)-X(4)
      C3=0.
      C4=0.
      C5=0.
      C6=0.
      D1=-Y(1)-Y(2)+Y(3)+Y(4)
      D2=Y(1)-Y(2)+Y(3)-Y(4)
      D3=0.
      D4=0.
      D5=0.
      D6=0.
      GO TO 40
   20 CONTINUE
      A1=2.*(X(6)-X(8))
      A2=X(1)-X(2)+X(3)-X(4)
      A3=2.*(X(1)+X(2)+X(3)+X(4)-2.*X(5)-2.*X(7))
      A4=2.*(-X(1)-X(2)+X(3)+X(4)+2.*X(5)-2.*X(7))
      A5=-X(1)+X(2)+X(3)-X(4)-2.*X(6)+2.*X(8)
      A6=0.
      B1=2.*(Y(6)-Y(8))
      B2=Y(1)-Y(2)+Y(3)-Y(4)
      B3=2.*(Y(1)+Y(2)+Y(3)+Y(4)-2.*Y(5)-2.*Y(7))
      B4=2.*(-Y(1)-Y(2)+Y(3)+Y(4)+2.*Y(5)-2.*Y(7))
      B5=-Y(1)+Y(2)+Y(3)-Y(4)-2.*Y(6)+2.*Y(8)
      B6=0.
      C1=2.*(X(7)-X(5))
      C2=X(1)-X(2)+X(3)-X(4)
      C3=2.*(X(1)+X(2)+X(3)+X(4)-2.*X(6)-2.*X(8))
      C4=2.*(-X(1)+X(2)+X(3)-X(4)-2.*X(6)+2.*X(8))
      C5=-X(1)-X(2)+X(3)+X(4)+2.*X(5)-2.*X(7)
      C6=0.
      D1=2.*(Y(7)-Y(5))
      D2=Y(1)-Y(2)+Y(3)-Y(4)
      D3=2.*(Y(1)+Y(2)+Y(3)+Y(4)-2.*Y(6)-2.*Y(8))
      D4=2.*(-Y(1)+Y(2)+Y(3)-Y(4)-2.*Y(6)+2.*Y(8))
      D5=-Y(1)-Y(2)+Y(3)+Y(4)+2.*Y(5)-2.*Y(7)
      D6=0.
      GO TO 40
   30 CONTINUE
      A1=2.*(X(6)-X(8))
      A2=X(1)-X(2)+X(3)-X(4)
      A3=4.*(X(6)+X(8)-2.*X(9))
      A4=2.*(-X(1)-X(2)+X(3)+X(4)+2.*X(5)-2.*X(7))
      A5=-X(1)+X(2)+X(3)-X(4)-2.*X(6)+2.*X(8)
      A6=2.*(X(1)+X(2)+X(3)+X(4)-2.*(X(5)+X(6)+X(7)+X(8))+4.*X(9))
      B1=2.*(Y(6)-Y(8))
      B2=Y(1)-Y(2)+Y(3)-Y(4)
      B3=4.*(Y(6)+Y(8)-2.*Y(9))
      B4=2.*(-Y(1)-Y(2)+Y(3)+Y(4)+2.*Y(5)-2.*Y(7))
      B5=-Y(1)+Y(2)+Y(3)-Y(4)-2.*Y(6)+2.*Y(8)
      B6=2.*(Y(1)+Y(2)+Y(3)+Y(4)-2.*(Y(5)+Y(6)+Y(7)+Y(8))+4.*Y(9))
      C1=2.*(X(7)-X(5))
      C2=X(1)-X(2)+X(3)-X(4)
      C3=4.*(X(5)+X(7)-2.*X(9))
      C4=2.*(-X(1)+X(2)+X(3)-X(4)-2.*X(6)+2.*X(8))
      C5=-X(1)-X(2)+X(3)+X(4)+2.*X(5)-2.*X(7)
      C6=2.*(X(1)+X(2)+X(3)+X(4)-2.*(X(5)+X(6)+X(7)+X(8))+4.*X(9))
      D1=2.*(Y(7)-Y(5))
      D2=Y(1)-Y(2)+Y(3)-Y(4)
      D3=4.*(Y(5)+Y(7)-2.*Y(9))
      D4=2.*(-Y(1)+Y(2)+Y(3)-Y(4)-2.*Y(6)+2.*Y(8))
      D5=-Y(1)-Y(2)+Y(3)+Y(4)+2.*Y(5)-2.*Y(7)
      D6=2.*(Y(1)+Y(2)+Y(3)+Y(4)-2.*(Y(5)+Y(6)+Y(7)+Y(8))+4.*Y(9))
   40 CONTINUE
C
C     COMPUTE STRAIN RATES AT GAUSS POINTS
C
      DO 100 I=1,NGAUSS
      EXX=0.
      EYY=0.
      EXY=0.
      A=F11Q(GSPT1(I),GSPT2(I))
      B=F21Q(GSPT1(I),GSPT2(I))
      C=F12Q(GSPT1(I),GSPT2(I))
      D=F22Q(GSPT1(I),GSPT2(I))
      EI=1./DETJQ(GSPT1(I),GSPT2(I))
      GO TO (50, 50, 50, 50, 70, 70), KIND
   50 CONTINUE
      RG=RAD8L(I)
      IF (KIND.EQ.4) RG=RAD8Q(I)
      UG=U8Q(I)
      DO 60 J=1,8
      AA=DF8QDS(J,I)*A+DF8QDT(J,I)*C
      BB=DF8QDS(J,I)*B+DF8QDT(J,I)*D
      EXX=EXX+AA*U(J)*EI
      EYY=EYY+BB*V(J)*EI
      EXY=EXY+AA*V(J)*EI+BB*U(J)*EI
   60 CONTINUE
      GO TO 90
   70 CONTINUE
      RG=RAD9L(I)
      IF (KIND.EQ.6) RG=RAD9Q(I)
      UG=U9Q(I)
      DO 80 J=1,9
      AA=DF9QDS(J,I)*A+DF9QDT(J,I)*C
      BB=DF9QDS(J,I)*B+DF9QDT(J,I)*D
      EXX=EXX+AA*U(J)*EI
      EYY=EYY+BB*V(J)*EI
      EXY=EXY+AA*V(J)*EI+BB*U(J)*EI
   80 CONTINUE
   90 CONTINUE
      SCNDII(I)=2.*(EXX**2+EYY**2)+EXY**2+IAXSYM*2.*(UG/RG)**2
  100 CONTINUE
C
C     EXTRAPOLATE GAUSS POINT VALUES TO NODAL POINTS
C
      CALL EXTRAP (KIND,SCNDII)
      RETURN
C
C     TRIANGULAR ELEMENTS
C
  110 CONTINUE
      A1=0.
      A2=0.
      A3=Y(2)-Y(3)
      B1=0.
      B2=0.
      B3=-X(2)+X(3)
      C1=0.
      C2=0.
      C3=-Y(1)+Y(3)
      D1=0.
      D2=0.
      D3=X(1)-X(3)
      GO TO 130
  120 CONTINUE
      A1=4.*(Y(3)+Y(4)-Y(5)-Y(6))
      A2=4.*(Y(2)+Y(3)-2.*Y(5))
      A3=4.*Y(5)-3.*Y(3)-Y(2)
      B1=4.*(X(3)+X(4)-X(5)-X(6))
      B2=4.*(X(2)+X(3)-2.*X(5))
      B3=4.*X(5)-3.*X(3)-X(2)
      C1=4.*(Y(1)+Y(3)-2.*Y(6))
      C2=4.*(Y(3)+Y(4)-Y(5)-Y(6))
      C3=4.*Y(6)-3.*Y(3)-Y(1)
      D1=4.*(X(1)+X(3)-2.*X(6))
      D2=4.*(X(3)+X(4)-X(5)-X(6))
      D3=4.*X(6)-3.*X(3)-X(1)
  130 CONTINUE
C
C     COMPUTE STRAIN RATES AT INTEGRATION POINTS
C
      DO 150 I=1,NHAM
      EXX=0.
      EYY=0.
      EXY=0.
      A=F11T(HAMPT1(I),HAMPT2(I))
      B=F12T(HAMPT1(I),HAMPT2(I))
      C=F21T(HAMPT1(I),HAMPT2(I))
      D=F22T(HAMPT1(I),HAMPT2(I))
      EI=1./DETJT(HAMPT1(I),HAMPT2(I))
      RG=RAD6L(I)
      IF (KIND.EQ.2) RG=RAD6Q(I)
      UG=U6Q(I)
      DO 140 J=1,6
      AA=DF6QDS(J,I)*A+DF6QDT(J,I)*B
      BB=DF6QDS(J,I)*C+DF6QDT(J,I)*D
      EXX=EXX+AA*U(J)*EI
      EYY=EYY+BB*V(J)*EI
      EXY=EXY+AA*V(J)*EI+BB*U(J)*EI
  140 CONTINUE
      SCNDII(I)=2.*(EXX**2+EYY**2)+EXY**2+IAXSYM*2.*(UG/RG)**2
  150 CONTINUE
C
C     EXTRAPOLATE INTEGRATION POINT VALUES TO NODAL POINTS
C
      CALL EXTRAP (KIND,SCNDII)
C
      RETURN
      END

max_stars_repo_name: nachos-ii/nachos-ii

max_issues_repo_name: nachos-ii/nachos-ii

max_forks_repo_name: nachos-ii/nachos-ii


text:
C***********************************************************************
C***********************************************************************
C
C   Version:        0.3
C   Last modified:  December 27, 1994
C   Authors:        Esmond G. Ng and Barry W. Peyton
C
C   Mathematical Sciences Section, Oak Ridge National Laboratory
C
C***********************************************************************
C***********************************************************************
C*************     MMPY2  .... MATRIX-MATRIX MULTIPLY     **************
C***********************************************************************
C***********************************************************************
C
C   PURPOSE -
C       THIS ROUTINE PERFORMS A MATRIX-MATRIX MULTIPLY, Y = Y + XA,
C       ASSUMING DATA STRUCTURES USED IN SOME OF OUR SPARSE CHOLESKY
C       CODES.
C
C       LOOP UNROLLING: LEVEL 2
C
C   INPUT PARAMETERS -
C       M               -   NUMBER OF ROWS IN X AND IN Y.
C       N               -   NUMBER OF COLUMNS IN X AND NUMBER OF ROWS
C                           IN A.
C       Q               -   NUMBER OF COLUMNS IN A AND Y.
C       XPNT(*)         -   XPNT(J+1) POINTS ONE LOCATION BEYOND THE
C                           END OF THE J-TH COLUMN OF X.  XPNT IS ALSO
C                           USED TO ACCESS THE ROWS OF A.
C       X(*)            -   CONTAINS THE COLUMNS OF X AND THE ROWS OF A.
C       LDY             -   LENGTH OF FIRST COLUMN OF Y.
C
C   UPDATED PARAMETERS -
C       Y(*)            -   ON OUTPUT, Y = Y + AX.
C
C***********************************************************************
C
      SUBROUTINE  MMPY2  (  M     , N     , Q     , XPNT  , X     ,
     &                      Y     , LDY                             )
C
C***********************************************************************
C
C       -----------
C       PARAMETERS.
C       -----------
C
        INTEGER             LDY   , M     , N     , Q
        INTEGER             XPNT(*)
        DOUBLE PRECISION    X(*)          , Y(*)
C
C       ----------------
C       LOCAL VARIABLES.
C       ----------------
C
        INTEGER             I1    , I2
        INTEGER             J1    , J2
        INTEGER             IY    , IYLAST, IYSTRT, IYSTOP, LENY  ,
     &                      MM    , REMAIN, XCOL  , YCOL
        DOUBLE PRECISION    A1    , A2
C
        INTEGER             LEVEL
        PARAMETER       (   LEVEL = 2 )
C
C***********************************************************************
C
C       -----------------------------------------------------------
C       INITIAL OFFSETS, COLUMN LENGTHS, AND INDEX RANGE VARIABLES.
C       -----------------------------------------------------------
        REMAIN = MOD ( N, LEVEL ) + 1
        MM = M
        IYLAST = 0
        LENY = LDY
C
C       ------------------------------------
C       TO COMPUTE EACH COLUMN YCOL OF Y ...
C       ------------------------------------
C
        DO  500  YCOL = 1, Q
C
            IYSTRT = IYLAST + 1
            IYSTOP = IYSTRT + MM - 1
            IYLAST = IYLAST + LENY
C
C           --------------------------------------------------
C           ... PERFORM THE APPROPRATE MATRIX VECTOR MULTIPLY:
C               X * A(*,YCOL) WITH LEVEL 2 LOOP-UNROLLING.
C           --------------------------------------------------
            GO TO ( 200, 100 ), REMAIN
C
  100           CONTINUE
                J1 = XPNT(1)
                I1 = XPNT(1+1) - MM
                A1  = - X(J1)*X(I1)
                DO  150  IY = IYSTRT, IYSTOP
                    Y(IY) = Y(IY) + A1*X(I1)
                    I1 = I1 + 1
  150           CONTINUE
                GO TO 200
C
  200           CONTINUE
                DO  400  XCOL = REMAIN, N, LEVEL
                    J1 = XPNT(XCOL)
                    J2 = XPNT(XCOL+1)
                    I1 = J2           - MM
                    I2 = XPNT(XCOL+2) - MM
                    A1  = - X(J1)*X(I1)
                    A2  = - X(J2)*X(I2)
                    DO  300  IY = IYSTRT, IYSTOP
                        Y(IY) = ( (Y(IY))
     &                          + A1*X(I1)) + A2*X(I2)
                        I1 = I1 + 1
                        I2 = I2 + 1
  300               CONTINUE
  400           CONTINUE
C
            MM = MM - 1
            LENY = LENY - 1
C
  500   CONTINUE
C
        RETURN
        END

max_stars_repo_name: Pangqiyuangh/SeIInv

max_issues_repo_name: Pangqiyuangh/SeIInv

max_forks_repo_name: Pangqiyuangh/SeIInv


text:
Quik is a pets kept cats feline who lives with his Users/TiffanySnell mom, Users/SharonZimmerman Sharon and Users/AmyZimmerman Amy, the infamous Zimmerman twins, Users/ScottBiggart Scott of the legendary DavisCoeds DavisCoeds.com and Fish at the Sycamore Lane Apartments.  Originally from LA, he moved to Davis with his previous Dad.   At some point, probably in 2003, he developed an extremely painful abcess in 2 teeth and spent his time in a depressed haze.  Frustrated by the lack of veterinary care he was receiving, he decide to search for a new mommy.  Eventually, he chose Users/TiffanySnell for her School of Veterinary Medicine veterinary aspirations.  She is a proud member of the Class of 2010.  After a tooth extraction surgery in March 2005 at Midtown Animal Clinic, a unique and princely personality was revealed.  Quik soon became a different kitty.  He generally spends his days chasing after bugs, eating, and sleeping.

Quik used to live on A St with Users/JonBafus Jon, the drummer extrodinare for Sholi, who coined the nickname Puff Pants.

Unfortunately, on September 6, 2008, Quik lost his life after being hit by a car :(


text:
      subroutine singleatoponshell(q1,q2,q7,p,iswitch,m)
      implicit none
************************************************************************
*     Author: R.K. Ellis, January 2012                                 *
*                                                                      *
*     q(-p1)+g(-p2) --> t~(-->e^-(p3)+nu~(p4)+b~(p5))+b(p6)+q'(p7)     *
*                                                                      *
*     keeping polarization information for t~,pb and gluon             *
*     iswitch= 0 for no gluon emission                                 *
*     iswitch=-1 for gluon emission in atop decay                      *
************************************************************************
      include 'constants.f'
      include 'masses.f'
      include 'zprods_decl.f'
      double precision p(mxpart,4),q(mxpart,4),
     & ala,alb,dot,s17,mt2,twopaDp2,twopbDp2
      double complex m(2,2,2),iza,izb,cprop
      integer p1,p2,a,p7,b,e,si,aa,bb,iswitch,q1,q2,q7
      parameter(p1=1,p2=2,a=3,e=4,b=5,p7=6)
C-----matrix element for p1+p2(g) -> t+c where both t and c are on shell
C-----t rendered massless wrt eb, and c rendered massless wrt p2
c--- statement functions
      iza(aa,bb)=cone/za(aa,bb)
      izb(aa,bb)=cone/zb(aa,bb)
c--- end statement functions
c      write(6,*) 'iswitch',iswitch
c      call writeout(p)
c      pause
C---zero all arrays
      m(:,:,:)=czip
      do si=1,4
      q(p1,si)=p(q1,si)
      q(p2,si)=p(q2,si)
      if (iswitch .eq. 0) then
      q(a,si)=p(3,si)+p(4,si)+p(5,si)
      elseif (iswitch .eq. -1) then
      q(a,si)=p(3,si)+p(4,si)+p(5,si)+p(8,si)
      endif
      q(e,si)=p(3,si)
      q(b,si)=p(6,si)
      q(p7,si)=p(q7,si)
      enddo
      mt2=mt**2
C---- now render "a" massless wrt to vector e
C---- now render "b" massless wrt to vector p2
      ala=mt2/(2d0*dot(q,a,e))
      alb=mb**2/(2d0*dot(q,b,p2))
      s17=2d0*dot(q,p1,p7)
      cprop=dcmplx(s17-wmass**2,zip)
      twopaDp2=2d0*dot(q,a,p2)
      twopbDp2=2d0*dot(q,b,p2)
      do si=1,4
      q(a,si)=q(a,si)-ala*q(e,si)
      q(b,si)=q(b,si)-alb*q(p2,si)
      enddo
      call spinoru(6,q,za,zb)

C----order of indices is polb,polg,polatop
      m(1,1,2)= + cprop**(-1)*twopbDp2**(-1) * (  - za(b,p2)*za(p2,p7)*
     &    zb(a,p1) - za(b,p2)*za(p2,p7)*zb(a,p1)*alb )
      m(1,1,2) = m(1,1,2) + cprop**(-1)*twopbDp2**(-1)*mb**2 * (  - za(
     &    p2,p7)*zb(a,p1)*izb(b,p2) )
      m(1,1,2) = m(1,1,2) + cprop**(-1)*twopaDp2**(-1) * ( za(e,p2)*za(
     &    b,p7)*zb(e,p1)*zb(b,a)*izb(b,p2)*ala + za(b,p7)*za(a,p2)*zb(b
     &    ,a)*zb(a,p1)*izb(b,p2) )
      m(1,1,2) = m(1,1,2) + cprop**(-1)*mt2*twopaDp2**(-1) * ( za(e,p2)
     &    *za(b,p7)*zb(b,p1)*iza(e,a)*izb(b,p2) )

      m(1,2,2)= + cprop**(-1)*twopaDp2**(-1) * ( za(e,b)*za(b,p7)*zb(e,
     &    p1)*zb(a,p2)*iza(b,p2)*ala - za(b,a)*za(b,p7)*zb(a,p1)*zb(a,
     &    p2)*iza(b,p2) + za(b,p7)*zb(a,p2)*zb(p1,p2) )
      m(1,2,2) = m(1,2,2) + cprop**(-1)*mt2*twopaDp2**(-1) * ( za(e,b)*
     &    za(b,p7)*zb(p1,p2)*iza(e,a)*iza(b,p2) )

      m(1,1,1)= + cprop**(-1)*mt*twopbDp2**(-1) * (  - za(b,p2)*za(p2,
     &    p7)*zb(e,p1)*izb(e,a) - za(b,p2)*za(p2,p7)*zb(e,p1)*izb(e,a)*
     &    alb )
      m(1,1,1) = m(1,1,1) + cprop**(-1)*mt*twopbDp2**(-1)*mb**2 * (  -
     &    za(p2,p7)*zb(e,p1)*izb(e,a)*izb(b,p2) )
      m(1,1,1) = m(1,1,1) + cprop**(-1)*mt*twopaDp2**(-1) * (  - za(e,
     &    p2)*za(b,p7)*zb(e,b)*zb(e,p1)*izb(e,a)*izb(b,p2)*ala - za(b,
     &    p7)*za(a,p2)*zb(e,b)*zb(a,p1)*izb(e,a)*izb(b,p2) + za(b,p7)*
     &    za(a,p2)*zb(b,p1)*izb(b,p2) )

      m(1,2,1)= + cprop**(-1)*mt*twopaDp2**(-1) * ( za(e,b)*za(b,p7)*
     &    zb(e,p1)*zb(e,p2)*iza(b,p2)*izb(e,a)*ala - za(b,a)*za(b,p7)*
     &    zb(e,p2)*zb(a,p1)*iza(b,p2)*izb(e,a) - za(b,a)*za(b,p7)*zb(p1
     &    ,p2)*iza(b,p2) + za(b,p7)*zb(e,p2)*zb(p1,p2)*izb(e,a) )

      m(2,1,2)= + cprop**(-1)*twopaDp2**(-1)*mb * (  - za(e,p2)*za(p2,
     &    p7)*zb(e,p1)*zb(b,a)*iza(b,p2)*izb(b,p2)*ala - za(a,p2)*za(p2
     &    ,p7)*zb(b,a)*zb(a,p1)*iza(b,p2)*izb(b,p2) )
      m(2,1,2) = m(2,1,2) + cprop**(-1)*mt2*twopaDp2**(-1)*mb * (  -
     &    za(e,p2)*za(p2,p7)*zb(b,p1)*iza(e,a)*iza(b,p2)*izb(b,p2) )

      m(2,2,2)= + cprop**(-1)*twopaDp2**(-1)*mb * (  - za(e,b)*za(p2,p7
     &    )*zb(e,p1)*zb(a,p2)*iza(b,p2)**2*ala + za(b,a)*za(p2,p7)*zb(a
     &    ,p1)*zb(a,p2)*iza(b,p2)**2 - za(p2,p7)*zb(a,p2)*zb(p1,p2)*
     &    iza(b,p2) )
      m(2,2,2) = m(2,2,2) + cprop**(-1)*mt2*twopaDp2**(-1)*mb * (  -
     &    za(e,b)*za(p2,p7)*zb(p1,p2)*iza(e,a)*iza(b,p2)**2 )

      m(2,1,1)= + cprop**(-1)*mt*twopaDp2**(-1)*mb * ( za(e,p2)*za(p2,
     &    p7)*zb(e,b)*zb(e,p1)*iza(b,p2)*izb(e,a)*izb(b,p2)*ala + za(a,
     &    p2)*za(p2,p7)*zb(e,b)*zb(a,p1)*iza(b,p2)*izb(e,a)*izb(b,p2)
     &     - za(a,p2)*za(p2,p7)*zb(b,p1)*iza(b,p2)*izb(b,p2) )

      m(2,2,1)= + cprop**(-1)*mt*twopaDp2**(-1)*mb * (  - za(e,b)*za(p2
     &    ,p7)*zb(e,p1)*zb(e,p2)*iza(b,p2)**2*izb(e,a)*ala + za(b,a)*
     &    za(p2,p7)*zb(e,p2)*zb(a,p1)*iza(b,p2)**2*izb(e,a) + za(b,a)*
     &    za(p2,p7)*zb(p1,p2)*iza(b,p2)**2 - za(p2,p7)*zb(e,p2)*zb(p1,
     &    p2)*iza(b,p2)*izb(e,a) )

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
C    %W% %G%
      subroutine int3fc (a)
C
C  THIS SUBROUTINE CALCULATES INITIAL STATE VECTORS AND PAST VALUES
C  FOR EXCITER TYPE FC
C
      include 'tspinc/param.inc'
      include 'tspinc/prt.inc'
      include 'tspinc/params.inc'
      include 'tspinc/gentbla.inc'
      include 'tspinc/pointr.inc'
      include 'tspinc/gentblb.inc'
      include 'tspinc/int3.inc'
      dimension array(40),a(40)

       equivalence         (array( 1),   vrmax), (array( 2),  vrmin),
     1 (array( 3),  vref), (array( 4),   esatx), (array( 5),  csatx),
     2 (array( 6),   hbr), (array( 7),     cke), (array( 8),    cka),
     3 (array( 9),    dr), (array(10),     cfd), (array(11),    ckf),
     4 (array(12),    df), (array(13),     akf), (array(14),    hbf),
     5 (array(15),   hcb), (array(16),     ckc), (array(17),     db),
     6 (array(18),    dc), (array(19),     hbe), (array(20),     de),
     7 (array(21),  rcex), (array(22),      vr), (array(23),   xcex),
     8 (array(24),    va), (array(25),      da), (array(26),    hba),
     9 (array(27),   vco), (array(28), cfldpst), (array(29),    efd),
     * (array(30),   ckd), (array(31),   x7pst), (array(32),  vcpst),
     * (array(33),   vfe), (array(34),      ve), (array(35),     vf),
     * (array(36),    a4), (array(37),      a5), (array(38),     a6)

C TRANSFER CITER DATA TO ARRAY
      l6 = 6
      do itr = 1,40
        array(itr) = a(itr)
      end do

      dtsk = 0.1

C CHECK FOR LIMIT VIOLATION
c     if (efdo .ge. efdmax .or. efdo .le. efdmin) then
c       write(errbuf,2000) name,base,id,efdo
c2000   format('0 INITIAL FIELD VOLTS VIOLATES LIMITS FOR ',
c    1 a8,2x,f5.1,2x,a1,5x,' INITIAL VOLTS = ',f6.3)
c       call prterr('E',1)
c       iabort = 1
c       go to 2300
c     end if

      efd = efdo

C CALCULATE VOLTAGE TRANSDUCER CIRCUIT
      if(xt .eq. 0.0)  then
        vtr = vhd
        vti = vhq
      else
        vtr = vhd + rt*agend - xt*agenq
        vti = vhq + rt*agenq + xt*agend
      endif
      crt = agend
      cit = agenq
      zcr = rcex*crt - xcex*cit
      zci = rcex*cit + xcex*crt
      vcr = zcr + vtr
      vci = zci + vti
      vc = sqrt(vcr**2 + vci**2)
      vco = vc
      vcpst = vco
      vc1 = vc

C CALCULATE R-BLOCK
      dr = 2.0*dr/dtsk + 1.0
      hbr = (dr - 2.0)*vc1 + vc

C CALCULATE FEX-BLOCK ( INPUT: EFDO,CFD,CKC;
C                      OUTPUT: FEX,VE)
C Note: Initial field current CFD equals initial field EFDO.
      cfd = efdo
      cfldo = cfd
      cfldpst = cfldo

      call getfex (fex,efdo,ve,cfd,ckc)

C CALCULATE SATURATION BLOCK
      ckep = csatx * exp(esatx*abs(efdo))
      ck1 = ckep + ckep*efdo*esatx
      ck2 = efdo**2 * esatx*ckep
      vfe = (cke+ck1)*ve + ckd*cfd - ck2

C CALCULATE A-BLOCK (REGULATOR)
      vr = vfe
      x3 = vr
      va = vr/cka
      da = 2.0*da/dtsk + 1.0
      hba = (da - 2.0)*vr + cka*va

C CHECK REGULATOR LIMITS
      if (vr .ge. vrmax .or. vr .le. vrmin) then
        write (errbuf, 2036) name, base, id, vr
        call prterr ('E',1)
 2036   format('0 INITIAL VR VIOLATES LIMITS FOR ', a8,
     1           2x, f5.1, 2x, a1, 5x, 'INITIAL VR = ', f6.3)
        iabort = 1
        go to 2300
      endif

C CALCULATE E-BLOCK
      de = 2.0*de/dtsk
      hbe = de*ve + vr - vfe
      veo = ve
      vfeo = vfe

C CALCULATE F-BLOCK
      vf = 0.0
      df = 2.0*df/dtsk + 1.0
      akf = 2.0*ckf/dtsk
      hbf = (df - 2.0)*vf - akf*vfe

C CALCULATE B-C BLOCK
      vai = va
      db = 2.0*db/dtsk + 1.0
      dc = 2.0*dc/dtsk + 1.0
      hcb = (db - 2.0)*va - (dc - 2.0)*vai
      verr = vai
      vref = verr + vc1
      ckep = 1.0
      ckap = cka

C Store coefficients for field solution:
C   VE(t) = A4*VERR(t) + A4*X7(t) + A5*CFD(t) + A6
      if (ve .le. 0.0) then
        ckep = 0.0
      else
        ckep = 1.0
      endif
      a1 = ckep*cka*dc/(da*db)
      a2 = ckep + a1*akf/df
      a3 = 1.0/(de + a2*(cke+ck1))
      a4 = a1*a3
      a5 = -a2*a3*ckd
      a6 = hbe*a3 + ckep*hba*a3/da + a1*hcb*a3/dc - a1*hbf*a3/df
     1       + a2*a3*ck2

c line to abort to
 2300 continue

C TRANSFER DATA TO CITER AND RETURN TO INITAL3
      do itr = 1,40
        a(itr) = array(itr)
      end do

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
module blas_c2
  use blas_kinds, only : ik, sp
  implicit none
  private :: ik, sp

  interface
    subroutine cgbmv(trans, m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: trans
      integer(ik) :: m
      integer(ik) :: n
      integer(ik) :: kl
      integer(ik) :: ku
      complex(sp) :: alpha
      complex(sp) :: a(lda, *)
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: beta
      complex(sp) :: y(*)
      integer(ik) :: incy
    end subroutine cgbmv
  end interface
  interface
    subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: trans
      integer(ik) :: m
      integer(ik) :: n
      complex(sp) :: alpha
      complex(sp) :: a(lda, *)
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: beta
      complex(sp) :: y(*)
      integer(ik) :: incy
    end subroutine cgemv
  end interface
  interface
    subroutine cgerc(m, n, alpha, x, incx, y, incy, a, lda)
      import :: ik, sp
      integer(ik) :: lda
      integer(ik) :: m
      integer(ik) :: n
      complex(sp) :: alpha
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: y(*)
      integer(ik) :: incy
      complex(sp) :: a(lda, *)
    end subroutine cgerc
  end interface
  interface
    subroutine cgeru(m, n, alpha, x, incx, y, incy, a, lda)
      import :: ik, sp
      integer(ik) :: lda
      integer(ik) :: m
      integer(ik) :: n
      complex(sp) :: alpha
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: y(*)
      integer(ik) :: incy
      complex(sp) :: a(lda, *)
    end subroutine cgeru
  end interface
  interface
    subroutine chbmv(uplo, n, k, alpha, a, lda, x, incx, beta, y, incy)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: uplo
      integer(ik) :: n
      integer(ik) :: k
      complex(sp) :: alpha
      complex(sp) :: a(lda, *)
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: beta
      complex(sp) :: y(*)
      integer(ik) :: incy
    end subroutine chbmv
  end interface
  interface
    subroutine chemv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: uplo
      integer(ik) :: n
      complex(sp) :: alpha
      complex(sp) :: a(lda, *)
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: beta
      complex(sp) :: y(*)
      integer(ik) :: incy
    end subroutine chemv
  end interface
  interface
    subroutine cher2(uplo, n, alpha, x, incx, y, incy, a, lda)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: uplo
      integer(ik) :: n
      complex(sp) :: alpha
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: y(*)
      integer(ik) :: incy
      complex(sp) :: a(lda, *)
    end subroutine cher2
  end interface
  interface
    subroutine chpmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
      import :: ik, sp
      character(len=1) :: uplo
      integer(ik) :: n
      complex(sp) :: alpha
      complex(sp) :: ap(*)
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: beta
      complex(sp) :: y(*)
      integer(ik) :: incy
    end subroutine chpmv
  end interface
  interface
    subroutine chpr2(uplo, n, alpha, x, incx, y, incy, ap)
      import :: ik, sp
      character(len=1) :: uplo
      integer(ik) :: n
      complex(sp) :: alpha
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: y(*)
      integer(ik) :: incy
      complex(sp) :: ap(*)
    end subroutine chpr2
  end interface
  interface
    subroutine chpr(uplo, n, alpha, x, incx, ap)
      import :: ik, sp
      character(len=1) :: uplo
      integer(ik) :: n
      real(sp) :: alpha
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: ap(*)
    end subroutine chpr
  end interface
  interface
    subroutine cher(uplo, n, alpha, x, incx, a, lda)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: uplo
      integer(ik) :: n
      real(sp) :: alpha
      complex(sp) :: x(*)
      integer(ik) :: incx
      complex(sp) :: a(lda, *)
    end subroutine cher
  end interface
  interface
    subroutine ctbmv(uplo, trans, diag, n, k, a, lda, x, incx)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: uplo
      character(len=1) :: trans
      character(len=1) :: diag
      integer(ik) :: n
      integer(ik) :: k
      complex(sp) :: a(lda, *)
      complex(sp) :: x(*)
      integer(ik) :: incx
    end subroutine ctbmv
  end interface
  interface
    subroutine ctbsv(uplo, trans, diag, n, k, a, lda, x, incx)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: uplo
      character(len=1) :: trans
      character(len=1) :: diag
      integer(ik) :: n
      integer(ik) :: k
      complex(sp) :: a(lda, *)
      complex(sp) :: x(*)
      integer(ik) :: incx
    end subroutine ctbsv
  end interface
  interface
    subroutine ctpmv(uplo, trans, diag, n, ap, x, incx)
      import :: ik, sp
      character(len=1) :: uplo
      character(len=1) :: trans
      character(len=1) :: diag
      integer(ik) :: n
      complex(sp) :: ap(*)
      complex(sp) :: x(*)
      integer(ik) :: incx
    end subroutine ctpmv
  end interface
  interface
    subroutine ctpsv(uplo, trans, diag, n, ap, x, incx)
      import :: ik, sp
      character(len=1) :: uplo
      character(len=1) :: trans
      character(len=1) :: diag
      integer(ik) :: n
      complex(sp) :: ap(*)
      complex(sp) :: x(*)
      integer(ik) :: incx
    end subroutine ctpsv
  end interface
  interface
    subroutine ctrmv(uplo, trans, diag, n, a, lda, x, incx)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: uplo
      character(len=1) :: trans
      character(len=1) :: diag
      integer(ik) :: n
      complex(sp) :: a(lda, *)
      complex(sp) :: x(*)
      integer(ik) :: incx
    end subroutine ctrmv
  end interface
  interface
    subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
      import :: ik, sp
      integer(ik) :: ldb
      integer(ik) :: lda
      character(len=1) :: side
      character(len=1) :: uplo
      character(len=1) :: transa
      character(len=1) :: diag
      integer(ik) :: m
      integer(ik) :: n
      complex(sp) :: alpha
      complex(sp) :: a(lda, *)
      complex(sp) :: b(ldb, *)
    end subroutine ctrsm
  end interface
  interface
    subroutine ctrsv(uplo, trans, diag, n, a, lda, x, incx)
      import :: ik, sp
      integer(ik) :: lda
      character(len=1) :: uplo
      character(len=1) :: trans
      character(len=1) :: diag
      integer(ik) :: n
      complex(sp) :: a(lda, *)
      complex(sp) :: x(*)
      integer(ik) :: incx
    end subroutine ctrsv
  end interface

end module blas_c2

max_stars_repo_name: awvwgk/blas-interface

max_issues_repo_name: awvwgk/blas-interface

max_forks_repo_name: awvwgk/blas-interface