text:
      subroutine setrun
c----------------------------------------------------------------------
c     Sets the run parameters reading them from the run_card.dat
c
c 1. PDF set
c 2. Collider parameters
c 3. cuts
c----------------------------------------------------------------------
      implicit none
c
c     include
c
      include 'genps.inc'
      include 'PDF/pdf.inc'
      include 'run.inc'
      include 'alfas.inc'
c
c     local
c
      integer npara
      character*20 param(maxpara),value(maxpara)
      character*20 ctemp
      integer k,i,l1,l2
      character*132 buff
      real*8 sf1,sf2
      integer lp1,lp2
      real*8 eb1,eb2
      real*8 pb1,pb2
C
C     input cuts
C
      include 'cuts.inc'
c
c----------
c     start
c----------
c
c     read the run_card.dat
c
      call load_para(npara,param,value)

c*********************************************************************
c     Jet measure cuts                                               *
c*********************************************************************

      call get_real   (npara,param,value," xqcut ",xqcut,0d0)

c************************************************************************
c    Collider energy and type                                           *
c************************************************************************
c     lpp  = -1 (antiproton), 0 (no pdf), 1 (proton)
c     lpp  =  2 (proton emitting a photon without breaking)
c     lpp  =  3 (electron emitting a photon)
c     ebeam= energy of each beam in GeV

      call get_integer(npara,param,value," lpp1 "   ,lp1,1  )
      call get_integer(npara,param,value," lpp2 "   ,lp2,1  )
      call get_real   (npara,param,value," ebeam1 " ,eb1,7d3)
      call get_real   (npara,param,value," ebeam2 " ,eb2,7d3)

      lpp(1)=lp1
      lpp(2)=lp2
      ebeam(1)=eb1
      ebeam(2)=eb2

c************************************************************************
c    Collider pdf                                                       *
c************************************************************************

      call get_string (npara,param,value," pdlabel ",pdlabel,'cteq6l1')
c
c     if lhapdf is used the following number identifies the set
c
      if(pdlabel.eq.'''lhapdf''')
     $     call get_integer(npara,param,value," lhaid  ",lhaid,10042)

      call pdfwrap

      return
 99   write(*,*) 'error in reading'
      return
      end


text:
module Mod_nsm_FreeSurfaceMatrices2nd
   use Mod_nsm_BaseElmope
   use Mod_CutMesh
   use typre
   implicit none
   private
   public :: SetPointersFreeSurfaceMatrices2nd

   !External Procedures
   procedure() :: NULLSUB

   !Coupling with levelset
   real(rp)                 :: gplev(1)
   integer(ip)              :: elemStatus(1),ngauss_minus,ngauss_plus,ngaus_total
   real(rp), allocatable    :: weigp(:)
   real(rp), pointer        :: xloc(:,:) => NULL()
   integer(ip), allocatable :: kfl_IsSet
   !To define the pressure boundary condition over free surface
   !Nodal coordinates
   real(rp), pointer       :: coord(:) => NULL()

contains

   !For setting the pointers
   subroutine SetPointersFreeSurfaceMatrices2nd(itask)
      implicit none
      integer(ip) :: itask
      select case (itask)

      case(0)
         allocate(kfl_IsSet)
         call a%Memor%allocObj(0,'kfl_IsSet','InitProcedurePointer',1)
         kfl_IsSet = -1

      case(1)
         if (kfl_IsSet == -1) then
            !Pointers are already set
            kfl_IsSet = 1

            if (a%kfl_colev == 1) then
               !Level Set Coupling
               call ConcatenateProcedures(ProcHook_Initializations,AllocLev)
               call ConcatenateProcedures(ProcHook_Finalizations,DeallocLev)

               if(a%kfl_fsurf==1)then
                  call ConcatenateProcedures(ProcHook_PreDirichlet,FreeSurfMats)

                  if(a%kfl_fsurfLapla==1)then
                     call PrependProcedure(ProcHook_PreDirichlet,ElmatsToLapla)
                  end if

                  !Always after FreeSurfMatsTF
                  call ConcatenateProcedures(ProcHook_PreDirichlet,PressureZeroFurf_2nd)
               end if

            endif

         endif

      case(100)
         deallocate(kfl_IsSet)
         call a%Memor%deallocObj(0,'kfl_IsSet','InitProcedurePointer',1)

      end select
   end subroutine

   subroutine AllocLev
      implicit none
      integer(ip) :: auxdim,auxdim1

      auxdim=(e%ndime+1)
      auxdim1=(e%ndime-1)*(e%ndime-1)+2

      call a%Memor%alloc(auxdim*auxdim1,weigp,'weigp','nsm_elmope')

   end subroutine

   subroutine DeallocLev
      implicit none
      integer(ip) :: auxdim,auxdim1

      auxdim=(e%ndime+1)
      auxdim1=(e%ndime-1)*(e%ndime-1)+2

      call a%Memor%dealloc(auxdim*auxdim1,weigp,'weigp','nsm_elmope')

   end subroutine

   subroutine FreeSurfMats
      implicit none
      integer(ip)  :: elemStatus,inode,poinStatus,ipoin

      call a%CutMesh%GetElementType(ielem,elemStatus)

      if(elemStatus==-1)then
         do inode=1,e%pnode
            ipoin=e%lnods(inode)
            call a%CutMesh%GetPointType(ipoin,poinStatus)
            if(poinStatus==-1)then

               elmat(:,inode,:,1:e%pnode) = 0.0_rp
               elrhs(:,inode)=0.0_rp

            end if
         end do
      end if
   end subroutine

   subroutine ElmatsToLapla
      implicit none
      integer(ip)  :: elemStatus,inode,poinStatus,ipoin,idime

      call a%CutMesh%GetElementType(ielem,elemStatus)

      if(elemStatus==-1)then

         elmat(:,:,:,:) = 0.0_rp
         elrhs(:,:)=0.0_rp
         wrmat1=0.0_rp

         ! Viscosity terms : we only consider mu*(grad v, grad u)
         call elmvis(e,dvolt0,acvis,wrmat1)

         elmat(1,1:e%pnode,1,1:e%pnode) = elmat(1,1:e%pnode,1,1:e%pnode) + wrmat1(1:e%pnode,1:e%pnode)

         do inode = 1,e%pnode
            !Rhs
            elrhs(1,inode) = elrhs(1,inode) + acvis*dot_product(e%cartd(:,inode),acden*a%grnor*a%gravi(1:e%ndime))*dvolt0
         enddo

      end if

   end subroutine

   subroutine PressureZeroFurf_2nd
      implicit none
      integer(ip) :: elemStatus,poinStatus,JpoinStatus,ndone,ntodo
      integer(ip) :: nipoin,ipoin,jpoin,jnode
      real(rp)    :: pdsurf,toleralce
      !real(rp)    :: boundaryVector(3),hipo,tmp
      !real(rp)    :: UnitBoundaryVector(3),,BoundaryCenter(3),checkvector(3)
      !3d case
      !real(rp)    :: boundaryVector2(3),hipo2,localx3(3),hipo3,projCheck
      !integer(ip) :: ialone,idime,jdime,inode !node alone
      real(rp)    :: Residuo(3),alfa,gppre(1),beta,grpreN
      integer(ip) :: idime,inode
      real(rp) :: tmp, UnitOrthogonalVector(3)

      call a%CutMesh%GetElementType(ielem,elemStatus)

      toleralce=1e-6

      if(elemStatus==0)then

         !To integrate over the surface or line intersection
         call a%CutMesh%GetLocalIntersectionPoints(ielem,xloc)
         call a%CutMesh%GetSurfaceIntersection(ielem,e,pdsurf)
         call a%CutMesh%GetNinters(ielem,nipoin)
         call a%CutMesh%GetExternalNormal(ielem,e,UnitOrthogonalVector)

         !Compute shape functions at the intersection points
         do jpoin = 1, nipoin
            weigp(jpoin) = 1.0_rp/nipoin
         end do

         !The rutine give the needed shape functions associated
         call e%SetParticularGaussPoints(a%Memor,nipoin,xloc,weigp)

         !-----------------------------------------------------------------------------------------------
!             !Approximate imposition of boundary conditions in immersed boundary methods

            !Ramon-Joan Method
         if(pdsurf>toleralce)then

            do inode=1,e%pnode
               ipoin=e%lnods(inode)
               call a%CutMesh%GetPointType(ipoin,poinStatus)
               !if(poinStatus==-1)then
                  !elmat(:,inode,:,1:e%pnode) = 0.0_rp
                  !elrhs(:,inode) = 0.0_rp
                  do jnode = 1,e%pnode
                     do igaus = 1,e%pgaus
                        e%igaus = igaus
                        tmp=1.0_rp/(a%MatProp(1)%densi*LHSdtinv)

                        elmat(1,inode,1,jnode) = elmat(1,inode,1,jnode) &
                        + 100*e%shape(inode,e%igaus)*e%shape(jnode,e%igaus)*(tmp/chale(1))*weigp(e%igaus)*pdsurf

                     enddo
                  enddo

               !end if
            end do

            !terms from the integrations by parts
            do inode=1,e%pnode
               ipoin=e%lnods(inode)
               call a%CutMesh%GetPointType(ipoin,poinStatus)
               if(poinStatus==1)then
                  do jnode = 1,e%pnode
                     do idime=1,e%ndime
                        do igaus = 1,e%pgaus
                           e%igaus = igaus
                           tmp=1.0_rp/(a%MatProp(1)%densi*LHSdtinv)

                           elmat(1,inode,1,jnode) = elmat(1,inode,1,jnode) &
                              - e%shape(inode,e%igaus)*e%cartd(idime,jnode)*UnitOrthogonalVector(idime)*(tmp)*(weigp(e%igaus)*pdsurf)
                        end do
                     enddo
                  enddo

               end if
            end do

            !terms from the integrations by parts
            do inode=1,e%pnode
               ipoin=e%lnods(inode)
               call a%CutMesh%GetPointType(ipoin,poinStatus)
               if(poinStatus==1)then
                  do idime=1,e%ndime
                     do igaus = 1,e%pgaus
                        e%igaus = igaus
                        tmp=1.0_rp/(a%MatProp(1)%densi*LHSdtinv)

                        elrhs(1,inode) = elrhs(1,inode) &
                           - e%shape(inode,e%igaus)*grpre(1,idime)*tmp*UnitOrthogonalVector(idime)*(weigp(e%igaus)*pdsurf)

                     end do
                  enddo

               end if
            end do

         elseif(pdsurf<=toleralce)then

         end if

      end if

   end subroutine

end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
! File:  psb_srwextd.f90
! Subroutine:
! Arguments:
!
! We have a problem here: 1. How to handle well all the formats?
!                         2. What should we do with rowscale? Does it only
!                            apply when a%fmt()='COO' ??????
!
!
subroutine psb_srwextd(nr,a,info,b,rowscale)
  use psb_mat_mod
  use psb_serial_mod, psb_protect_name => psb_srwextd
  implicit none

  ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
  integer(psb_ipk_), intent(in)                 :: nr
  type(psb_sspmat_type), intent(inout)        :: a
  integer(psb_ipk_),intent(out)                 :: info
  type(psb_sspmat_type), intent(in), optional :: b
  logical,intent(in), optional                 :: rowscale

  integer(psb_ipk_) :: i,j,ja,jb,err_act,nza,nzb
  character(len=20)                 :: name, ch_err
  type(psb_s_coo_sparse_mat)        :: actmp
  logical  rowscale_

  name='psb_srwextd'
  info  = psb_success_
  call psb_erractionsave(err_act)

  if (nr > a%get_nrows()) then
    select type(aa=> a%a)
    type is (psb_s_csr_sparse_mat)
      if (present(b)) then
        call psb_rwextd(nr,aa,info,b%a,rowscale)
      else
        call psb_rwextd(nr,aa,info,rowscale=rowscale)
      end if
    type is (psb_s_coo_sparse_mat)
      if (present(b)) then
        call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale)
      else
        call psb_rwextd(nr,aa,info,rowscale=rowscale)
      end if
    class default
      call aa%mv_to_coo(actmp,info)
      if (info == psb_success_) then
        if (present(b)) then
          call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale)
        else
          call psb_rwextd(nr,actmp,info,rowscale=rowscale)
        end if
      end if
      if (info == psb_success_) call aa%mv_from_coo(actmp,info)
    end select
  end if
  if (info /= psb_success_) goto 9999

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return

end subroutine psb_srwextd
subroutine psb_sbase_rwextd(nr,a,info,b,rowscale)
  use psb_mat_mod
  use psb_serial_mod, psb_protect_name => psb_sbase_rwextd
  implicit none

  ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
  integer(psb_ipk_), intent(in)                        :: nr
  class(psb_s_base_sparse_mat), intent(inout)        :: a
  integer(psb_ipk_),intent(out)                        :: info
  class(psb_s_base_sparse_mat), intent(in), optional :: b
  logical,intent(in), optional                       :: rowscale

  integer(psb_ipk_) :: i,j,ja,jb,err_act,nza,nzb, ma, mb, na, nb
  character(len=20)                 :: name, ch_err
  logical  rowscale_

  name='psb_sbase_rwextd'
  info  = psb_success_
  call psb_erractionsave(err_act)

  if (present(rowscale)) then
    rowscale_ = rowscale
  else
    rowscale_ = .true.
  end if

  ma  = a%get_nrows()
  na  = a%get_ncols()
  nza = a%get_nzeros()

  select type(a)
  type is (psb_s_csr_sparse_mat)

    call psb_ensure_size(nr+1,a%irp,info)

    if (present(b)) then
      mb = b%get_nrows()
      nb = b%get_ncols()
      nzb = b%get_nzeros()

      select type (b)
      type is (psb_s_csr_sparse_mat)
        call psb_ensure_size(nza+nzb,a%ja,info)
        call psb_ensure_size(nza+nzb,a%val,info)
        a%ja(nza+1:nza+nzb)  = b%ja(1:nzb)
        a%val(nza+1:nza+nzb) = b%val(1:nzb)
        do i=1, min(nr-ma,mb)
          a%irp(ma+i+1) =  a%irp(ma+i) + b%irp(i+1) - b%irp(i)
        end do
        do j=i,nr-ma
          a%irp(ma+j+1) = a%irp(ma+j)
        end do
      class default

        write(psb_err_unit,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
      end select
      call a%set_ncols(max(na,nb))

    else

      do j=ma+2,nr+1
        a%irp(j) = a%irp(j-1)
      end do

    end if

    call a%set_nrows(nr)

  type is (psb_s_coo_sparse_mat)
    nza = a%get_nzeros()

    if (present(b)) then
      mb = b%get_nrows()
      nb = b%get_ncols()
      nzb = b%get_nzeros()
      call a%ensure_size(nza+nzb)

      select type(b)
      type is (psb_s_coo_sparse_mat)

        if (rowscale_) then
          do j=1,nzb
            if ((ma + b%ia(j)) <= nr) then
              nza = nza + 1
              a%ia(nza)  = ma + b%ia(j)
              a%ja(nza)  = b%ja(j)
              a%val(nza) = b%val(j)
            end if
          enddo
        else
          do j=1,nzb
            if ((ma + b%ia(j)) <= nr) then
              nza = nza + 1
              a%ia(nza)  = b%ia(j)
              a%ja(nza)  = b%ja(j)
              a%val(nza) = b%val(j)
            end if
          enddo
        endif
        call a%set_nzeros(nza)

      type is (psb_s_csr_sparse_mat)

        do i=1, min(nr-ma,mb)
          do jb = b%irp(i), b%irp(i+1)-1
            nza = nza + 1
            a%val(nza) = b%val(jb)
            a%ia(nza)  = ma + i
            a%ja(nza)  = b%ja(jb)
          end do
        end do
        call a%set_nzeros(nza)

      class default
        write(psb_err_unit,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'

      end select

      call a%set_ncols(max(na,nb))
    endif

    call a%set_nrows(nr)

  class default
    info = psb_err_unsupported_format_
    ch_err=a%get_fmt()
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999

  end select

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return

end subroutine psb_sbase_rwextd

subroutine psb_lsrwextd(nr,a,info,b,rowscale)
  use psb_mat_mod
  use psb_serial_mod, psb_protect_name => psb_lsrwextd
  implicit none

  ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
  integer(psb_lpk_), intent(in)                  :: nr
  type(psb_lsspmat_type), intent(inout)        :: a
  integer(psb_ipk_),intent(out)                  :: info
  type(psb_lsspmat_type), intent(in), optional :: b
  logical,intent(in), optional                 :: rowscale

  integer(psb_lpk_) :: i,j,ja,jb,nza,nzb
  integer(psb_ipk_) :: err_act
  character(len=20)                 :: name, ch_err
  type(psb_ls_coo_sparse_mat)     :: actmp
  logical  rowscale_

  name='psb_lsrwextd'
  info  = psb_success_
  call psb_erractionsave(err_act)

  if (nr > a%get_nrows()) then
    select type(aa=> a%a)
    type is (psb_ls_csr_sparse_mat)
      if (present(b)) then
        call psb_rwextd(nr,aa,info,b%a,rowscale)
      else
        call psb_rwextd(nr,aa,info,rowscale=rowscale)
      end if
    type is (psb_ls_coo_sparse_mat)
      if (present(b)) then
        call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale)
      else
        call psb_rwextd(nr,aa,info,rowscale=rowscale)
      end if
    class default
      call aa%mv_to_coo(actmp,info)
      if (info == psb_success_) then
        if (present(b)) then
          call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale)
        else
          call psb_rwextd(nr,actmp,info,rowscale=rowscale)
        end if
      end if
      if (info == psb_success_) call aa%mv_from_coo(actmp,info)
    end select
  end if
  if (info /= psb_success_) goto 9999

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return

end subroutine psb_lsrwextd
subroutine psb_lsbase_rwextd(nr,a,info,b,rowscale)
  use psb_mat_mod
  use psb_serial_mod, psb_protect_name => psb_lsbase_rwextd
  implicit none

  ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
  integer(psb_lpk_), intent(in)                         :: nr
  class(psb_ls_base_sparse_mat), intent(inout)        :: a
  integer(psb_ipk_),intent(out)                         :: info
  class(psb_ls_base_sparse_mat), intent(in), optional :: b
  logical,intent(in), optional                       :: rowscale

  integer(psb_lpk_) :: i,j,ja,jb,nza,nzb, ma, mb, na, nb
  integer(psb_ipk_) :: err_act
  character(len=20) :: name, ch_err
  logical  rowscale_

  name='psb_lsbase_rwextd'
  info  = psb_success_
  call psb_erractionsave(err_act)

  if (present(rowscale)) then
    rowscale_ = rowscale
  else
    rowscale_ = .true.
  end if

  ma  = a%get_nrows()
  na  = a%get_ncols()
  nza = a%get_nzeros()

  select type(a)
  type is (psb_ls_csr_sparse_mat)

    call psb_ensure_size(nr+1,a%irp,info)

    if (present(b)) then
      mb = b%get_nrows()
      nb = b%get_ncols()
      nzb = b%get_nzeros()

      select type (b)
      type is (psb_ls_csr_sparse_mat)
        call psb_ensure_size(nza+nzb,a%ja,info)
        call psb_ensure_size(nza+nzb,a%val,info)
        a%ja(nza+1:nza+nzb)  = b%ja(1:nzb)
        a%val(nza+1:nza+nzb) = b%val(1:nzb)
        do i=1, min(nr-ma,mb)
          a%irp(ma+i+1) =  a%irp(ma+i) + b%irp(i+1) - b%irp(i)
        end do
        do j=i,nr-ma
          a%irp(ma+j+1) = a%irp(ma+j)
        end do
      class default

        write(psb_err_unit,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
      end select
      call a%set_ncols(max(na,nb))

    else

      do j=ma+2,nr+1
        a%irp(j) = a%irp(j-1)
      end do

    end if

    call a%set_nrows(nr)

  type is (psb_ls_coo_sparse_mat)
    nza = a%get_nzeros()

    if (present(b)) then
      mb = b%get_nrows()
      nb = b%get_ncols()
      nzb = b%get_nzeros()
      call a%ensure_size(nza+nzb)

      select type(b)
      type is (psb_ls_coo_sparse_mat)

        if (rowscale_) then
          do j=1,nzb
            if ((ma + b%ia(j)) <= nr) then
              nza = nza + 1
              a%ia(nza)  = ma + b%ia(j)
              a%ja(nza)  = b%ja(j)
              a%val(nza) = b%val(j)
            end if
          enddo
        else
          do j=1,nzb
            if ((ma + b%ia(j)) <= nr) then
              nza = nza + 1
              a%ia(nza)  = b%ia(j)
              a%ja(nza)  = b%ja(j)
              a%val(nza) = b%val(j)
            end if
          enddo
        endif
        call a%set_nzeros(nza)

      type is (psb_ls_csr_sparse_mat)

        do i=1, min(nr-ma,mb)
          do jb = b%irp(i), b%irp(i+1)-1
            nza = nza + 1
            a%val(nza) = b%val(jb)
            a%ia(nza)  = ma + i
            a%ja(nza)  = b%ja(jb)
          end do
        end do
        call a%set_nzeros(nza)

      class default
        write(psb_err_unit,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'

      end select

      call a%set_ncols(max(na,nb))
    endif

    call a%set_nrows(nr)

  class default
    info = psb_err_unsupported_format_
    ch_err=a%get_fmt()
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999

  end select

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return

end subroutine psb_lsbase_rwextd


text:
!
! SignalFile_Utility
!
! Module containing routines to create and manipulate signal files
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 15-May-2007
!                       paul.vandelst@ssec.wisc.edu
!

MODULE SignalFile_Utility

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE File_Utility
  USE Message_Handler
  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: Create_SignalFile

CONTAINS

!################################################################################
!################################################################################
!##                                                                            ##
!##                         ## PUBLIC MODULE ROUTINES ##                       ##
!##                                                                            ##
!################################################################################
!################################################################################

!--------------------------------------------------------------------------------
!
! NAME:
!       Create_Signal_File
!
! PURPOSE:
!       Function to create a "signal" file that can be checked for in
!       processing scripts. The presence of a signal file generally
!       indicates "something" completed successfully.
!
! CALLING SEQUENCE:
!       Error_Status = Create_SignalFile( Rootname               , &  ! Input
!                                         Message_Log=Message_Log  )  ! Error messaging
!
! INPUT ARGUMENTS:
!       Rootname:     Root of the name of the signal file to create. The
!                     created signal file will be this name appended with
!                     the string ".signal"
!                     UNITS:      N/A
!                     TYPE:       CHARACTER(*)
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: INTENT(IN)
!
!
! OPTIONAL INPUT ARGUMENTS:
!       Message_Log:  Character string specifying a filename in which any
!                     messages will be logged. If not specified, or if an
!                     error occurs opening the log file, the default action
!                     is to output messages to standard output.
!                     UNITS:      N/A
!                     TYPE:       CHARACTER(*)
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: INTENT(IN), OPTIONAL
!
! FUNCTION RESULT:
!       Error_Status: The return value is an integer defining the error
!                     status. The error codes are defined in the
!                     Message_Handler module.
!                     If == SUCCESS the signal file creation was successful
!                        == FAILURE an unrecoverable error occured
!                     UNITS:      N/A
!                     TYPE:       INTEGER
!                     DIMENSION:  Scalar
!
! SIDE EFFECTS:
!       If the signal file already exists, it is overwritten.
!
!--------------------------------------------------------------------------------

  FUNCTION Create_SignalFile( Rootname   , &  ! Input
                              Message_Log) &  ! Error messaging
                            RESULT ( Error_Status )
    ! Arguments
    CHARACTER(*),           INTENT(IN) :: Rootname
    CHARACTER(*), OPTIONAL, INTENT(IN) :: Message_Log
    ! Function result
    INTEGER :: Error_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Create_SignalFile'
    ! Local variables
    INTEGER :: FileID
    INTEGER :: IO_Status

    ! Set up
    ! ------
    Error_Status = SUCCESS

    ! Create the signal file
    ! ----------------------
    ! Get a free unit number
    FileID = Get_Lun()
    IF ( FileID < 0 ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Error obtaining file unit for '//&
                            TRIM(Rootname)//' signal file.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Open the file
    OPEN( FileID, FILE  =TRIM(Rootname)//'.signal', &
                  STATUS='REPLACE', &
                  FORM  ='FORMATTED', &
                  ACCESS='SEQUENTIAL', &
                  IOSTAT=IO_Status )
    IF ( IO_Status /= 0 ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Error opening '//&
                            TRIM(Rootname)//' signal file.', &
                            Error_Status, &
                            Message_Log=Message_Log )
      RETURN
    END IF

    ! Write the Rootname to the signal file
    ! -------------------------------------
    WRITE( FileID, FMT   ='(a)',  &
                   IOSTAT=IO_Status ) TRIM(Rootname)
    IF ( IO_Status /= 0 ) THEN
      Error_Status = FAILURE
      CALL Display_Message( ROUTINE_NAME, &
                            'Error writing to '//&
                            TRIM(Rootname)//' signal file.', &
                            Error_Status, &
                            Message_Log = Message_Log )
      CLOSE( FileID, STATUS='DELETE' )
    END IF

    ! Close the signal file
    ! ---------------------
    CLOSE( FileID, STATUS='KEEP' )

  END FUNCTION Create_SignalFile

END MODULE SignalFile_Utility

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
      SUBROUTINE AB09HY( N, M, P, A, LDA, B, LDB, C, LDC, D, LDD,
     $                   SCALEC, SCALEO, S, LDS, R, LDR, IWORK,
     $                   DWORK, LDWORK, BWORK, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To compute the Cholesky factors Su and Ru of the controllability
C     Grammian P = Su*Su' and observability Grammian Q = Ru'*Ru,
C     respectively, satisfying
C
C            A*P  + P*A' +  scalec^2*B*B'   = 0,       (1)
C
C            A'*Q + Q*A  +  scaleo^2*Cw'*Cw = 0,       (2)
C
C     where
C            Cw = Hw - Bw'*X,
C            Hw = inv(Dw)*C,
C            Bw = (B*D' + P*C')*inv(Dw'),
C            D*D' = Dw*Dw' (Dw upper triangular),
C
C     and, with Aw = A - Bw*Hw, X is the stabilizing solution of the
C     Riccati equation
C
C            Aw'*X + X*Aw + Hw'*Hw + X*Bw*Bw'*X = 0.   (3)
C
C     The P-by-M matrix D must have full row rank. Matrix A must be
C     stable and in a real Schur form.
C
C     ARGUMENTS
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of state-space representation, i.e.,
C             the order of the matrix A.  N >= 0.
C
C     M       (input) INTEGER
C             The number of system inputs.  M >= 0.
C
C     P       (input) INTEGER
C             The number of system outputs.  M >= P >= 0.
C
C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
C             The leading N-by-N part of this array must contain the
C             stable state dynamics matrix A in a real Schur canonical
C             form.
C
C     LDA     INTEGER
C             The leading dimension of array A.  LDA >= MAX(1,N).
C
C     B       (input) DOUBLE PRECISION array, dimension (LDB,M)
C             The leading N-by-M part of this array must contain the
C             input/state matrix B, corresponding to the Schur matrix A.
C
C     LDB     INTEGER
C             The leading dimension of array B.  LDB >= MAX(1,N).
C
C     C       (input) DOUBLE PRECISION array, dimension (LDC,N)
C             The leading P-by-N part of this array must contain the
C             state/output matrix C, corresponding to the Schur
C             matrix A.
C
C     LDC     INTEGER
C             The leading dimension of array C.  LDC >= MAX(1,P).
C
C     D       (input) DOUBLE PRECISION array, dimension (LDD,M)
C             The leading P-by-M part of this array must
C             contain the full row rank input/output matrix D.
C
C     LDD     INTEGER
C             The leading dimension of array D.  LDD >= MAX(1,P).
C
C     SCALEC  (output) DOUBLE PRECISION
C             Scaling factor for the controllability Grammian in (1).
C
C     SCALEO  (output) DOUBLE PRECISION
C             Scaling factor for the observability Grammian in (2).
C
C     S       (output) DOUBLE PRECISION array, dimension (LDS,N)
C             The leading N-by-N upper triangular part of this array
C             contains the Cholesky factor Su of the cotrollability
C             Grammian P = Su*Su' satisfying (1).
C
C     LDS     INTEGER
C             The leading dimension of array S.  LDS >= MAX(1,N).
C
C     R       (output) DOUBLE PRECISION array, dimension (LDR,N)
C             The leading N-by-N upper triangular part of this array
C             contains the Cholesky factor Ru of the observability
C             Grammian Q = Ru'*Ru satisfying (2).
C
C     LDR     INTEGER
C             The leading dimension of array R.  LDR >= MAX(1,N).
C
C     Workspace
C
C     IWORK   INTEGER array, dimension 2*N
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK and DWORK(2) contains RCOND, the reciprocal
C             condition number of the U11 matrix from the expression
C             used to compute X = U21*inv(U11). A small value RCOND
C             indicates possible ill-conditioning of the Riccati
C             equation (3).
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= MAX( 2, N*(MAX(N,M,P)+5),
C                            2*N*P+MAX(P*(M+2),10*N*(N+1) ) ).
C             For optimum performance LDWORK should be larger.
C
C     BWORK   LOGICAL array, dimension 2*N
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             = 1:  the state matrix A is not stable or is not in a
C                   real Schur form;
C             = 2:  the reduction of Hamiltonian matrix to real Schur
C                   form failed;
C             = 3:  the reordering of the real Schur form of the
C                   Hamiltonian matrix failed;
C             = 4:  the Hamiltonian matrix has less than N stable
C                   eigenvalues;
C             = 5:  the coefficient matrix U11 in the linear system
C                   X*U11 = U21, used to determine X, is singular to
C                   working precision;
C             = 6:  the feedthrough matrix D has not a full row rank P.
C
C     CONTRIBUTORS
C
C     A. Varga, German Aerospace Center, Oberpfaffenhofen, May 2000.
C     D. Sima, University of Bucharest, May 2000.
C     V. Sima, Research Institute for Informatics, Bucharest, May 2000.
C     Based on the RASP routines SRGRO and SRGRO1, by A. Varga, 1992.
C
C     REVISIONS
C
C     V. Sima, Research Institute for Informatics, Bucharest, Oct. 2001.
C
C     KEYWORDS
C
C     Minimal realization, model reduction, multivariable system,
C     state-space model, state-space representation,
C     stochastic balancing.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION ZERO, ONE, TWO
      PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
C     .. Scalar Arguments ..
      INTEGER          INFO, LDA, LDB, LDC, LDD, LDR, LDS, LDWORK, M, N,
     $                 P
      DOUBLE PRECISION SCALEC, SCALEO
C     .. Array Arguments ..
      INTEGER          IWORK(*)
      DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
     $                 DWORK(*), R(LDR,*), S(LDS,*)
      LOGICAL          BWORK(*)
C     .. Local Scalars ..
      INTEGER          I, IERR, KBW, KCW, KD, KDW, KG, KQ, KS, KTAU, KU,
     $                 KW, KWI, KWR, LW, N2, WRKOPT
      DOUBLE PRECISION RCOND, RTOL
C     .. External Functions ..
      DOUBLE PRECISION DLANGE, DLAMCH
      EXTERNAL         DLANGE, DLAMCH
C     .. External Subroutines ..
      EXTERNAL         DGEMM, DGERQF, DLACPY, DORGRQ, DSYRK, DTRMM,
     $                 DTRSM, SB02MD, SB03OU, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC        ABS, DBLE, INT, MAX, MIN
C     .. Executable Statements ..
C
      INFO = 0
      LW   = MAX( 2, N*( MAX( N, M, P ) + 5 ),
     $            2*N*P + MAX( P*(M + 2), 10*N*(N + 1) ) )
C
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( M.LT.0 ) THEN
         INFO = -2
      ELSE IF( P.LT.0 .OR. P.GT.M ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
         INFO = -9
      ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
         INFO = -11
      ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
         INFO = -15
      ELSE IF( LDR.LT.MAX( 1, N ) ) THEN
         INFO = -17
      ELSE IF( LDWORK.LT.LW ) THEN
         INFO = -20
      END IF
C
      IF( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'AB09HY', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      SCALEC = ONE
      SCALEO = ONE
      IF( MIN( N, M, P ).EQ.0 ) THEN
         DWORK(1) = TWO
         DWORK(2) = ONE
         RETURN
      END IF
C
C     Solve for Su the Lyapunov equation
C                                      2
C     A*(Su*Su') + (Su*Su')*A' + scalec *B*B' = 0 .
C
C     Workspace:  need   N*(MAX(N,M) + 5);
C                 prefer larger.
C
      KU   = 1
      KTAU = KU + N*MAX( N, M )
      KW   = KTAU + N
C
      CALL DLACPY( 'Full', N, M, B, LDB, DWORK(KU), N )
      CALL SB03OU( .FALSE., .TRUE., N, M, A, LDA, DWORK(KU), N,
     $             DWORK(KTAU), S, LDS, SCALEC, DWORK(KW),
     $             LDWORK - KW + 1, IERR )
      IF( IERR.NE.0 ) THEN
         INFO = 1
         RETURN
      ENDIF
      WRKOPT = INT( DWORK(KW) ) + KW - 1
C
C     Allocate workspace for Bw' (P*N), Cw (P*N), Q2 (P*M),
C     where Q2 = inv(Dw)*D.
C     Workspace:  need   2*N*P + P*M.
C
      KBW  = 1
      KCW  = KBW  + P*N
      KD   = KCW  + P*N
      KDW  = KD   + P*(M - P)
      KTAU = KD   + P*M
      KW   = KTAU + P
C
C     Compute an upper-triangular Dw such that D*D' = Dw*Dw', using
C     the RQ-decomposition of D: D = [0 Dw]*( Q1 ).
C                                           ( Q2 )
C     Additional workspace:  need 2*P; prefer P + P*NB.
C
      CALL DLACPY( 'F', P, M, D, LDD, DWORK(KD), P )
      CALL DGERQF( P, M, DWORK(KD), P, DWORK(KTAU), DWORK(KW),
     $             LDWORK-KW+1, IERR )
      WRKOPT = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
C
C     Check the full row rank of D.
C
      RTOL = DBLE( M ) * DLAMCH( 'E' ) *
     $       DLANGE( '1', P, M, D, LDD, DWORK )
      DO 10 I = KDW, KDW+P*P-1, P+1
         IF( ABS( DWORK(I) ).LE.RTOL ) THEN
            INFO = 6
            RETURN
         END IF
   10 CONTINUE
C                    -1
C     Compute Hw = Dw  *C.
C
      CALL DLACPY( 'F', P, N, C, LDC, DWORK(KCW), P )
      CALL DTRSM( 'Left', 'Upper', 'No-transpose', 'Non-unit', P, N,
     $            ONE, DWORK(KDW), P, DWORK(KCW), P )
C
C     Compute Bw' = inv(Dw)*(D*B' + C*Su*Su').
C
C     Compute first Hw*Su*Su' in Bw'.
C
      CALL DLACPY( 'F', P, N, DWORK(KCW), P, DWORK(KBW), P )
      CALL DTRMM( 'Right', 'Upper', 'No-transpose', 'Non-unit', P, N,
     $            ONE, S, LDS, DWORK(KBW), P )
      CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit', P, N,
     $            ONE, S, LDS, DWORK(KBW), P )
C
C     Compute Q2 = inv(Dw)*D, as the last P lines of the orthogonal
C     matrix ( Q1 ) from the RQ decomposition of D.
C            ( Q2 )
C     Additional workspace:  need P; prefer P*NB.
C
      CALL DORGRQ( P, M, P, DWORK(KD), P, DWORK(KTAU), DWORK(KW),
     $             LDWORK-KW+1, IERR )
      WRKOPT = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
C
C     Compute Bw' <- Bw' + Q2*B'.
C
      CALL DGEMM( 'No-transpose', 'Transpose', P, N, M, ONE,
     $            DWORK(KD), P, B, LDB, ONE, DWORK(KBW), P )
C
C     Compute Aw = A - Bw*Hw in R.
C
      CALL DLACPY( 'F', N, N, A, LDA, R, LDR )
      CALL DGEMM( 'Transpose', 'No-transpose', N, N, P, -ONE,
     $            DWORK(KBW), P, DWORK(KCW), P, ONE, R, LDR )
C
C     Allocate storage to solve the Riccati equation (3) for
C     G(N*N), Q(N*N), WR(2N), WI(2N), S(2N*2N), U(2N*2N).
C
      N2  = N + N
      KG  = KD
      KQ  = KG  + N*N
      KWR = KQ  + N*N
      KWI = KWR + N2
      KS  = KWI + N2
      KU  = KS  + N2*N2
      KW  = KU  + N2*N2
C
C     Compute G = -Bw*Bw'.
C
      CALL DSYRK( 'Upper', 'Transpose', N, P, -ONE, DWORK(KBW), P, ZERO,
     $            DWORK(KG), N )
C
C     Compute Q = Hw'*Hw.
C
      CALL DSYRK( 'Upper', 'Transpose', N, P, ONE, DWORK(KCW), P, ZERO,
     $            DWORK(KQ), N )
C
C     Solve
C
C        Aw'*X + X*Aw + Q - X*G*X = 0,
C
C     with Q =  Hw'*Hw  and  G = -Bw*Bw'.
C     Additional workspace: need   6*N;
C                           prefer larger.
C
      CALL SB02MD( 'Continuous', 'None', 'Upper', 'General', 'Stable',
     $             N, R, LDR, DWORK(KG), N, DWORK(KQ), N, RCOND,
     $             DWORK(KWR), DWORK(KWI), DWORK(KS), N2,
     $             DWORK(KU), N2, IWORK, DWORK(KW), LDWORK-KW+1,
     $             BWORK, INFO )
      IF( INFO.NE.0 )
     $   RETURN
      WRKOPT = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
C
C     Compute Cw = Hw - Bw'*X.
C
      CALL DGEMM ( 'No-transpose', 'No-transpose', P, N, N, -ONE,
     $              DWORK(KBW), P, DWORK(KQ), N, ONE, DWORK(KCW), P )
C
C     Solve for Ru the Lyapunov equation
C                                      2
C     A'*(Ru'*Ru) + (Ru'*Ru)*A + scaleo  * Cw'*Cw = 0 .
C
C     Workspace:  need   N*(MAX(N,P) + 5);
C                 prefer larger.
C
      KTAU = KCW  + N*MAX( N, P )
      KW   = KTAU + N
C
      CALL SB03OU( .FALSE., .FALSE., N, P, A, LDA, DWORK(KCW), P,
     $             DWORK(KTAU), R, LDR, SCALEO, DWORK(KW),
     $             LDWORK - KW + 1, IERR )
      WRKOPT = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
C
C     Save optimal workspace and RCOND.
C
      DWORK(1) = WRKOPT
      DWORK(2) = RCOND
C
      RETURN
C *** Last line of AB09HY ***
      END

ext: f

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


text:
! tl_ad_tensor.f90
!
!> Tangent Linear (TL) and Adjoint (AD) model versions of MAOOAM.
!> Tensors definition module
!
!> @copyright
!> 2016 Lesley De Cruz & Jonathan Demaeyer.
!> See LICENSE.txt for license information.
!
!---------------------------------------------------------------------------!
!
!> @remark
!> The routines of this module should be called only after
!> params::init_params() and aotensor_def::init_aotensor() have been called !
!
!---------------------------------------------------------------------------!

MODULE tl_ad_tensor

  !-----------------------------------------------------!
  !                                                     !
  ! Preamble and variables declaration                  !
  !                                                     !
  !-----------------------------------------------------!

  USE params, only:ndim
  USE aotensor_def, only:aotensor
  USE tensor, only:coolist,jsparse_mul,jsparse_mul_mat,sparse_mul3,simplify
  IMPLICIT NONE

  PRIVATE

  !> Epsilon to test equality with 0
  REAL(KIND=8), PARAMETER :: real_eps = 2.2204460492503131e-16

  !> Vector used to count the tensor elements
  INTEGER, DIMENSION(:), ALLOCATABLE :: count_elems

  !> Tensor representation of the Tangent Linear tendencies.
  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: tltensor

  !> Tensor representation of the Adjoint tendencies.
  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: adtensor

  PUBLIC :: init_tltensor,init_adtensor,init_adtensor_ref,ad,tl,jacobian_mat

  !-----------------------------------------------------!
  !                                                     !
  ! End of preamble                                     !
  !                                                     !
  !-----------------------------------------------------!

CONTAINS

  !-----------------------------------------------------!
  !                                                     !
  ! Function declarations :                             !
  !                                                     !
  ! These functions should be called only after         !
  ! init_params and init_aotensor                       !
  !                                                     !
  !-----------------------------------------------------!

  !-----------------------------------------------------!
  !                                                     !
  ! Jacobian functions                                  !
  !                                                     !
  !-----------------------------------------------------!

  !> Compute the Jacobian of MAOOAM in point ystar.
  !> @param ystar array with variables in which the jacobian should be evaluated.
  !> @return Jacobian in coolist-form (table of tuples {i,j,0,value})
  FUNCTION jacobian(ystar)
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: ystar
    TYPE(coolist), DIMENSION(ndim) :: jacobian
    CALL jsparse_mul(aotensor,ystar,jacobian)
  END FUNCTION jacobian

  !> Compute the Jacobian of MAOOAM in point ystar.
  !> @param ystar array with variables in which the jacobian should be evaluated.
  !> @return Jacobian in matrix form
  FUNCTION jacobian_mat(ystar)
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: ystar
    REAL(KIND=8), DIMENSION(ndim,ndim) :: jacobian_mat
    CALL jsparse_mul_mat(aotensor,ystar,jacobian_mat)
  END FUNCTION jacobian_mat

  !-----------------------------------------------------!
  !                                                     !
  ! Tangent linear model functions                      !
  !                                                     !
  !-----------------------------------------------------!

  !> Routine to initialize the TL tensor
  SUBROUTINE init_tltensor
    INTEGER :: i
    INTEGER :: AllocStat
    ALLOCATE(tltensor(ndim),count_elems(ndim), STAT=AllocStat)
    IF (AllocStat /= 0) STOP "*** Not enough memory ! ***"
    count_elems=0
    CALL compute_tltensor(tl_add_count)

    DO i=1,ndim
       ALLOCATE(tltensor(i)%elems(count_elems(i)), STAT=AllocStat)
       IF (AllocStat /= 0) STOP "*** Not enough memory ! ***"
    END DO

    DEALLOCATE(count_elems, STAT=AllocStat)
    IF (AllocStat /= 0) STOP "*** Deallocation problem ! ***"

    CALL compute_tltensor(tl_coeff)

    CALL simplify(tltensor)

  END SUBROUTINE init_tltensor

  !> Routine to compute the TL tensor from the original MAOOAM one
  !> @param func subroutine used to do the computation
  SUBROUTINE compute_tltensor(func)
    EXTERNAL :: func
    INTERFACE
       SUBROUTINE func(i,j,k,v)
         INTEGER, INTENT(IN) :: i,j,k
         REAL(KIND=8), INTENT(IN) :: v
       END SUBROUTINE func
    END INTERFACE
    INTEGER :: i,j,k,n,ne
    REAL(KIND=8) :: v
    DO i=1,ndim
       ne=aotensor(i)%nelems
       DO n=1,ne
          j=aotensor(i)%elems(n)%j
          k=aotensor(i)%elems(n)%k
          v=aotensor(i)%elems(n)%v
          CALL func(i,j,k,v)
       ENDDO
    ENDDO
  END SUBROUTINE compute_tltensor

  !> Subroutine used to count the number of TL tensor entries
  !> @param i tensor \f$i\f$ index
  !> @param j tensor \f$j\f$ index
  !> @param k tensor \f$k\f$ index
  !> @param v value that will be added
  SUBROUTINE tl_add_count(i,j,k,v)
    INTEGER, INTENT(IN) :: i,j,k
    REAL(KIND=8), INTENT(IN)  :: v
    IF (ABS(v) .ge. real_eps) THEN
       IF (j /= 0) count_elems(i)=count_elems(i)+1
       IF (k /= 0) count_elems(i)=count_elems(i)+1
    ENDIF
  END SUBROUTINE tl_add_count

  !> Subroutine used to compute the TL tensor entries
  !> @param i tensor \f$i\f$ index
  !> @param j tensor \f$j\f$ index
  !> @param k tensor \f$k\f$ index
  !> @param v value to add
  SUBROUTINE tl_coeff(i,j,k,v)
    INTEGER, INTENT(IN) :: i,j,k
    REAL(KIND=8), INTENT(IN) :: v
    INTEGER :: n
    IF (.NOT. ALLOCATED(tltensor)) STOP "*** tl_coeff routine : tensor not yet allocated ***"
    IF (.NOT. ALLOCATED(tltensor(i)%elems)) STOP "*** tl_coeff routine : tensor not yet allocated ***"
    IF (ABS(v) .ge. real_eps) THEN
       IF (j /=0) THEN
          n=(tltensor(i)%nelems)+1
          tltensor(i)%elems(n)%j=j
          tltensor(i)%elems(n)%k=k
          tltensor(i)%elems(n)%v=v
          tltensor(i)%nelems=n
       END IF
       IF (k /=0) THEN
          n=(tltensor(i)%nelems)+1
          tltensor(i)%elems(n)%j=k
          tltensor(i)%elems(n)%k=j
          tltensor(i)%elems(n)%v=v
          tltensor(i)%nelems=n
       END IF
    END IF
  END SUBROUTINE tl_coeff

  !-----------------------------------------------------!
  !                                                     !
  ! Adjoint model functions                             !
  !                                                     !
  !-----------------------------------------------------!

  !> Routine to initialize the AD tensor
  SUBROUTINE init_adtensor
    INTEGER :: i
    INTEGER :: AllocStat
    ALLOCATE(adtensor(ndim),count_elems(ndim), STAT=AllocStat)
    IF (AllocStat /= 0) STOP "*** Not enough memory ! ***"
    count_elems=0
    CALL compute_adtensor(ad_add_count)

    DO i=1,ndim
       ALLOCATE(adtensor(i)%elems(count_elems(i)), STAT=AllocStat)
       IF (AllocStat /= 0) STOP "*** Not enough memory ! ***"
    END DO

    DEALLOCATE(count_elems, STAT=AllocStat)
    IF (AllocStat /= 0) STOP "*** Deallocation problem ! ***"

    CALL compute_adtensor(ad_coeff)

    CALL simplify(adtensor)

  END SUBROUTINE init_adtensor

  !> Subroutine to compute the AD tensor from the original MAOOAM one
  !> @param func subroutine used to do the computation
  SUBROUTINE compute_adtensor(func)
    EXTERNAL :: func
    INTERFACE
       SUBROUTINE func(i,j,k,v)
         INTEGER, INTENT(IN) :: i,j,k
         REAL(KIND=8), INTENT(IN) :: v
       END SUBROUTINE func
    END INTERFACE
    INTEGER :: i,j,k,n,ne
    REAL(KIND=8) :: v
    DO i=1,ndim
       ne=aotensor(i)%nelems
       DO n=1,ne
          j=aotensor(i)%elems(n)%j
          k=aotensor(i)%elems(n)%k
          v=aotensor(i)%elems(n)%v
          CALL func(i,j,k,v)
       ENDDO
    ENDDO
  END SUBROUTINE compute_adtensor

  !> Subroutine used to count the number of AD tensor entries
  !> @param i tensor \f$i\f$ index
  !> @param j tensor \f$j\f$ index
  !> @param k tensor \f$k\f$ index
  !> @param v value that will be added
  SUBROUTINE ad_add_count(i,j,k,v)
    INTEGER, INTENT(IN) :: i,j,k
    REAL(KIND=8), INTENT(IN)  :: v
    IF ((ABS(v) .ge. real_eps).AND.(i /= 0)) THEN
       IF (k /= 0) count_elems(k)=count_elems(k)+1
       IF (j /= 0) count_elems(j)=count_elems(j)+1
    ENDIF
  END SUBROUTINE ad_add_count

  ! Subroutine used to compute the AD tensor entries
  !> @param i tensor \f$i\f$ index
  !> @param j tensor \f$j\f$ index
  !> @param k tensor \f$k\f$ index
  !> @param v value to add
  SUBROUTINE ad_coeff(i,j,k,v)
    INTEGER, INTENT(IN) :: i,j,k
    REAL(KIND=8), INTENT(IN) :: v
    INTEGER :: n
    IF (.NOT. ALLOCATED(adtensor)) STOP "*** ad_coeff routine : tensor not yet allocated ***"
    IF ((ABS(v) .ge. real_eps).AND.(i /=0)) THEN
       IF (k /=0) THEN
          IF (.NOT. ALLOCATED(adtensor(k)%elems)) STOP "*** ad_coeff routine : tensor not yet allocated ***"
          n=(adtensor(k)%nelems)+1
          adtensor(k)%elems(n)%j=i
          adtensor(k)%elems(n)%k=j
          adtensor(k)%elems(n)%v=v
          adtensor(k)%nelems=n
       END IF
       IF (j /=0) THEN
          IF (.NOT. ALLOCATED(adtensor(j)%elems)) STOP "*** ad_coeff routine : tensor not yet allocated ***"
          n=(adtensor(j)%nelems)+1
          adtensor(j)%elems(n)%j=i
          adtensor(j)%elems(n)%k=k
          adtensor(j)%elems(n)%v=v
          adtensor(j)%nelems=n
       END IF
    END IF
  END SUBROUTINE ad_coeff

  !-----------------------------------------------------!
  !                                                     !
  ! Adjoint model functions : Alternate method          !
  !                           using the TL tensor       !
  !                           (must be initialized      !
  !                           before)                   !
  !                                                     !
  !-----------------------------------------------------!

  !> Alternate method to initialize the AD tensor from the TL
  !> tensor
  !> @remark The #tltensor must be initialised before using this method.
  SUBROUTINE init_adtensor_ref
    INTEGER :: i
    INTEGER :: AllocStat
    ALLOCATE(adtensor(ndim),count_elems(ndim), STAT=AllocStat)
    IF (AllocStat /= 0) STOP "*** Not enough memory ! ***"
    count_elems=0
    CALL compute_adtensor_ref(ad_add_count_ref)

    DO i=1,ndim
       ALLOCATE(adtensor(i)%elems(count_elems(i)), STAT=AllocStat)
       IF (AllocStat /= 0) STOP "*** Not enough memory ! ***"
    END DO

    DEALLOCATE(count_elems, STAT=AllocStat)
    IF (AllocStat /= 0) STOP "*** Deallocation problem ! ***"

    CALL compute_adtensor_ref(ad_coeff_ref)

    CALL simplify(adtensor)

  END SUBROUTINE init_adtensor_ref

  !> Alternate subroutine to compute the AD tensor from the TL one
  !> @param func subroutine used to do the computation
  SUBROUTINE compute_adtensor_ref(func)
    EXTERNAL :: func
    INTERFACE
       SUBROUTINE func(i,j,k,v)
         INTEGER, INTENT(IN) :: i,j,k
         REAL(KIND=8), INTENT(IN) :: v
       END SUBROUTINE func
    END INTERFACE
    INTEGER :: i,j,k,n,ne
    REAL(KIND=8) :: v
    IF (.NOT. ALLOCATED(tltensor)) STOP "*** compute_adtensor_ref routine : tensor TL should be allocated first***"
    DO i=1,ndim
       ne=tltensor(i)%nelems
       DO n=1,ne
          j=tltensor(i)%elems(n)%j
          k=tltensor(i)%elems(n)%k
          v=tltensor(i)%elems(n)%v
          CALL func(i,j,k,v)
       ENDDO
    ENDDO
  END SUBROUTINE compute_adtensor_ref

  !> Alternate subroutine used to count the number of AD tensor entries
  !> from the TL tensor
  !> @param i tensor \f$i\f$ index
  !> @param j tensor \f$j\f$ index
  !> @param k tensor \f$k\f$ index
  !> @param v value that will be added
  SUBROUTINE ad_add_count_ref(i,j,k,v)
    INTEGER, INTENT(IN) :: i,j,k
    REAL(KIND=8), INTENT(IN)  :: v
    IF ((ABS(v) .ge. real_eps).AND.(j /= 0)) count_elems(j)=count_elems(j)+1
  END SUBROUTINE ad_add_count_ref

  !> Alternate subroutine used to compute the AD tensor entries from
  !> the TL tensor
  !> @param i tensor \f$i\f$ index
  !> @param j tensor \f$j\f$ index
  !> @param k tensor \f$k\f$ index
  !> @param v value to add
  SUBROUTINE ad_coeff_ref(i,j,k,v)
    INTEGER, INTENT(IN) :: i,j,k
    REAL(KIND=8), INTENT(IN) :: v
    INTEGER :: n
    IF (.NOT. ALLOCATED(adtensor)) STOP "*** ad_coeff_ref routine : tensor not yet allocated ***"
    IF ((ABS(v) .ge. real_eps).AND.(j /=0)) THEN
       IF (.NOT. ALLOCATED(adtensor(j)%elems)) STOP "*** ad_coeff_ref routine : tensor not yet allocated ***"
       n=(adtensor(j)%nelems)+1
       adtensor(j)%elems(n)%j=i
       adtensor(j)%elems(n)%k=k
       adtensor(j)%elems(n)%v=v
       adtensor(j)%nelems=n
    END IF
  END SUBROUTINE ad_coeff_ref

  !-----------------------------------------------------!
  !                                                     !
  ! Tendencies                                          !
  !                                                     !
  !-----------------------------------------------------!

  !> Tendencies for the AD of MAOOAM in point ystar for perturbation deltay.
  !> @param t time
  !> @param ystar vector with the variables (current point in trajectory)
  !> @param deltay vector with the perturbation of the variables at time t
  !> @param buf vector (buffer) to store derivatives.
  SUBROUTINE ad(t,ystar,deltay,buf)
    REAL(KIND=8), INTENT(IN) :: t
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: ystar,deltay
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: buf
    CALL sparse_mul3(adtensor,deltay,ystar,buf)
  END SUBROUTINE ad

  !> Tendencies for the TL of MAOOAM in point ystar for perturbation deltay.
  !> @param t time
  !> @param ystar vector with the variables (current point in trajectory)
  !> @param deltay vector with the perturbation of the variables at time t
  !> @param buf vector (buffer) to store derivatives.
  SUBROUTINE tl(t,ystar,deltay,buf)
    REAL(KIND=8), INTENT(IN) :: t
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(IN) :: ystar,deltay
    REAL(KIND=8), DIMENSION(0:ndim), INTENT(OUT) :: buf
    CALL sparse_mul3(tltensor,deltay,ystar,buf)
  END SUBROUTINE tl

END MODULE tl_ad_tensor

max_stars_repo_name: solepomies/MAOOAM

max_issues_repo_name: solepomies/MAOOAM

max_forks_repo_name: solepomies/MAOOAM


! Bugs item #877956, was opened at 2004-01-15 17:59
! Message generated for change (Tracker Item Submitted) made by Item Submitter
! You can respond by visiting:
! https://sourceforge.net/tracker/?func=detail&atid=106556&aid=877956&group_id=6556
!
! Category: None
! Group: None
! Status: Open
! Resolution: None
! Priority: 5
! Submitted By: Randy Hood (randy762)
! Assigned to: Nobody/Anonymous (nobody)
! Summary: Broken Fortran variable listing after =-1
!
! Initial Comment:
! When I run ctags v5.5.2 on Redhat Linux 9 with the command
!
! ctags --kinds-Fortran=v  -x test.f90
!
! where test.f90 is
! ----------------------------------
PROGRAM test

 IMPLICIT NONE
 INTEGER :: cm1 =-1, c2 = 2

END PROGRAM test
! -------------------------------------
!
! I only get this one line of output
!
! cm1              variable      4 test.f90        INTEGER :: cm1 =-1, c2 = 2
!
! If I change one line of test.f90 so that it is now
! ----------------------------------------
PROGRAM test

 IMPLICIT NONE
 INTEGER :: cm1 = -1, c2 = 2

END PROGRAM test
! -----------------------------------------
! and run the command
!
! ctags --kinds-Fortran=v  -x test.f90
!
! I get this correct output
!
! c2               variable      4 test.f90        INTEGER :: cm1 = -1, c2 = 2
! cm1              variable      4 test.f90        INTEGER :: cm1 = -1, c2 = 2
!
! ----------------------------------------------------------------------
! You can respond by visiting:
! https://sourceforge.net/tracker/?func=detail&atid=106556&aid=877956&group_id=6556


! RUN: %S/test_errors.sh %s %flang %t
! Check that if constructs only accept scalar logical expressions.
! TODO: expand the test to check this restriction for more types.

INTEGER :: I
LOGICAL, DIMENSION (2) :: B

!ERROR: Must be a scalar value, but is a rank-1 array
if ( B ) then
  a = 1
end if

!ERROR: Must be a scalar value, but is a rank-1 array
if ( B ) then
  a = 2
else
  a = 3
endif

!ERROR: Must be a scalar value, but is a rank-1 array
if ( B ) then
  a = 4
!ERROR: Must be a scalar value, but is a rank-1 array
else if( B ) then
  a = 5
end if

!ERROR: Must be a scalar value, but is a rank-1 array
if ( B ) then
  a = 6
!ERROR: Must be a scalar value, but is a rank-1 array
else if( B ) then
  a = 7
!ERROR: Must be a scalar value, but is a rank-1 array
elseif( B ) then
  a = 8
end if

!ERROR: Must be a scalar value, but is a rank-1 array
if ( B ) then
  a = 9
!ERROR: Must be a scalar value, but is a rank-1 array
else if( B ) then
  a = 10
else
  a = 11
end if

!ERROR: Must be a scalar value, but is a rank-1 array
if ( B ) then
  a = 12
!ERROR: Must be a scalar value, but is a rank-1 array
else if( B ) then
  a = 13
!ERROR: Must be a scalar value, but is a rank-1 array
else if( B ) then
  a = 14
end if

!ERROR: Must have LOGICAL type, but is INTEGER(4)
if ( I ) then
  a = 1
end if

!ERROR: Must have LOGICAL type, but is INTEGER(4)
if ( I ) then
  a = 2
else
  a = 3
endif

!ERROR: Must have LOGICAL type, but is INTEGER(4)
if ( I ) then
  a = 4
!ERROR: Must have LOGICAL type, but is INTEGER(4)
else if( I ) then
  a = 5
end if

!ERROR: Must have LOGICAL type, but is INTEGER(4)
if ( I ) then
  a = 6
!ERROR: Must have LOGICAL type, but is INTEGER(4)
else if( I ) then
  a = 7
!ERROR: Must have LOGICAL type, but is INTEGER(4)
elseif( I ) then
  a = 8
end if

!ERROR: Must have LOGICAL type, but is INTEGER(4)
if ( I ) then
  a = 9
!ERROR: Must have LOGICAL type, but is INTEGER(4)
else if( I ) then
  a = 10
else
  a = 11
end if

!ERROR: Must have LOGICAL type, but is INTEGER(4)
if ( I ) then
  a = 12
!ERROR: Must have LOGICAL type, but is INTEGER(4)
else if( I ) then
  a = 13
!ERROR: Must have LOGICAL type, but is INTEGER(4)
else if( I ) then
  a = 14
end if

!ERROR: Must have LOGICAL type, but is REAL(4)
if (f()) then
  a = 15
end if

contains
  real function f()
    f = 1.0
  end
end


text:
C
C
*+ do_gate_map

       subroutine do_gate_map(map_array,status)
C      ----------------------------------------
C
C Gate a map with itself
C
C Given:
C   map data
       real*4        map_array(*)
C Returned:
C   error status word
       integer       status
C
C All map values below a specified gate are set to a BLANKing value.
*-
       include '/mrao/include/iolib_functions.inc'
C
       integer    minirt(8), n, iout
       real*4     gate_low, gate_high, value, blank_value
       integer    imap, imapo, ip_mapo

C check status on entry
       if (status.ne.0) return
C read map entry
       call map_getmap('Map : ','Default-Map','READ',imap,status)
C allocate input map to output map
       call map_alloc_toout(imap,'DIRECT',
     *                      map_array,imapo,ip_mapo,status)
       if (status.ne.0) goto 999

C read gate low and high
       call io_enqout(iout)
       call io_getr('Gate-low  : ','-1.0E+30',gate_low,status)
       call io_getr('Gate-high : ','1.0E+30',gate_high,status)
       call ennull(blank_value,status)
       value = blank_value
       call io_getr('Value to replace with [BLANK] : ',
     *           ' ',value,status)
       if (status.ne.0) goto 999

C blank the map region and update the computing redtape
       call enminirt(minirt,status)
       do n=0,minirt(5)*minirt(6)-1
         if ( map_array(ip_mapo+n).lt.gate_low .or.
     *        map_array(ip_mapo+n).gt.gate_high ) then
           map_array(ip_mapo+n) = value
         else
           map_array(ip_mapo+n) = map_array(ip_mapo+n)
         end if
       end do

C Update the TEXT of the redtape
       call adredt('created','GATE-MAP',status)

C Tidy up, make output map new current map etc.
999    call map_end_alloc(imapo,map_array,status)

C check STATUS value
       call cmd_err(status,'GATE-MAP','Failed ')
       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2019 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!     This subroutines enables to calculate the reynolds number correction after:
!     "Discharge coeffcients for incompressible non-cavitating flowthrough long orifices"
!     A. Lichtarowicz, R.K duggins and E. Markland
!     Journal  Mechanical Engineering Science , vol 7, No. 2, 1965
!
!     author: Yannick Muller
!
      subroutine cd_lichtarowicz(cd,cdu,reynolds,amod,bdh)
!
      implicit none
!
      real*8 cdu,reynolds,amod,bdh,eps,A1,cd_diff,cd0,cd
!
      cd0=cdu
      cd_diff=1.d0
!
      do
!
         if(cd_diff.lt.1.d-3) exit
!
         cd=cd0
         A1=20/(reynolds*dsqrt(1.d0-Amod**2))*(1.d0+2.25d0*bdh)
         eps=(0.005d0*bdh)/(1.d0+7.5d0*(log10(0.00015d0*reynolds*
     &        dsqrt(1.d0-Amod**2)/cd))**2)

         cd=((-1/cdu+eps)+dsqrt((1/cdu-eps)**2.d0+4.d0*A1))/(2*A1)
!
         cd_diff=dabs(cd-cd0)
!
         cd0=cd
!
       enddo
!
         return
         end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      SUBROUTINE MPY3P (Z,IZ,DZ)
C*****
C    PERFORMS MULTIPLICATION AND SUMMATION.
C*****
      DOUBLE PRECISION DZ(1),DFACT
C
C
C
      INTEGER CODE,PREC
      INTEGER ZPNTRS
C
C
C
      DIMENSION Z(1),IZ(1)
C
C
C SUBROUTINE CALL PARAMETERS
      COMMON / MPY3CP / ITRL,ICORE,N,NCB,M,DUM1(3),ZPNTRS(22),LAEND,
     1                  DUM2(6),KA,KB
C FILES
      COMMON / MPY3TL / FILEA(7),FILEB(7),FILEE(7),FILEC(7),SCR1,SCR2,
     1                  SCR,LKORE,CODE,PREC,LCORE,SCR3(7),BUF1,BUF2,
     2                  BUF3,BUF4,E
C
C
C
      EQUIVALENCE     (FACT,DFACT)
C OPEN CORE POINTERS
      EQUIVALENCE     (IPOINT,ZPNTRS(3)),      (IACOLS,ZPNTRS(5)),
     *                (ITRANS,ZPNTRS(7)),      (IC,ZPNTRS(9)),
     *                (IBCOLS,ZPNTRS(11)),     (IAKJ,ZPNTRS(21))
C*****
C    LOOP FOR ACCUMULATING SUMS.
C*****
      KJ = IAKJ + KA - 1
      KJ2 = (IAKJ - 1)/2 + KA
      KB = IBCOLS + PREC*((KB - 1)*N - 1)
      IF (CODE .EQ. 2 .OR. ICORE .EQ. 1) GO TO 100
C*****
C    A(T)BA CASE.
C*****
      LP = IPOINT - 1
      DO 90 L=1,N
C CALCULATE FACTOR = B(LK)*A(KJ) TO BE MULTIPLIED TO NON-ZERO TERMS IN
C LTH COLUMN OF A(T)
      KB = KB + PREC
      LP = LP + 1
      IF (IZ(LP) .EQ. 0) GO TO 90
      IF (PREC .EQ. 2) GO TO 10
      IF (Z(KB) .EQ. 0.0) GO TO 90
      FACT = Z(KB)*Z(KJ)
      GO TO 20
   10 KB2 = (KB + 1)/2
      IF (DZ(KB2) .EQ. 0.0D0) GO TO 90
      DFACT = DZ(KB2)*DZ(KJ2)
   20 I1 = IZ(LP)
      IF (L .EQ. N) GO TO 40
C ACCUMULATE SUMS FOR NON-ZERO TERMS IN COLUMN L OF A(T)
      L1 = L + 1
      LLP = LP
      DO 30 LL=L1,N
      LLP = LLP + 1
      IF (IZ(LLP) .NE. 0) GO TO 50
   30 CONTINUE
   40 I2 = LAEND
      GO TO 60
   50 I2 = IZ(LLP) - 1
   60 IAC = IACOLS + I1 - 2
      IF (PREC .EQ. 2) GO TO 80
C SINGLE PRECISION CASE
      IAT = ITRANS + I1 - 2
      DO 70 I=I1,I2
      IAC = IAC + 1
      IAT = IAT + 1
      II = IC + IZ(IAC) - 1
   70 Z(II) = Z(II) + Z(IAT)*FACT
      GO TO 90
C DOUBLE PRECISION CASE
   80 IAT = (ITRANS - 3)/2 + I1
      DO 85 I=I1,I2
      IAC = IAC + 1
      IAT = IAT + 1
      II = (IC - 1)/2 + IZ(IAC)
   85 DZ(II) = DZ(II) + DZ(IAT)*DFACT
      III = (IC - 1)/2 + 1
   90 CONTINUE
      GO TO 999
C*****
C    BA CASE.
C*****
  100 IF (PREC .EQ. 2) GO TO 140
C SINGLE PRECISION CASE
      II = IC - 1
      DO 130 I=1,N
      II = II + 1
      KB = KB + 1
      IF (Z(KB) .EQ. 0.0) GO TO 130
      Z(II) = Z(II) + Z(KB)*Z(KJ)
  130 CONTINUE
      GO TO 999
C DOUBLE PRECISION CASE
  140 II = (IC - 1)/2
      KB = (KB + 1)/2
      DO 150 I=1,N
      II = II + 1
      KB = KB + 1
      IF (DZ(KB) .EQ. 0.0D0) GO TO 150
      DZ(II) = DZ(II) + DZ(KB)*DZ(KJ2)
  150 CONTINUE
C
  999 RETURN
      END


text:
      subroutine runC_0000i(k,l,i1,DetGr,f,Gtwiddle,Gtt,
     . Shat4zz,Shat5zzzz,S0000i,Shat5zz,N0)
      implicit none
C     Fixes C0000i
C     known C00i,C0000  corrections of order Delta C00iii
      include 'TRconstants.f'
      include 'pvCnames.f'
      include 'pvCv.f'
      include 'Carraydef.f'
      include 'Carrays.f'
      integer ep,N0,k,l,n,m,i1,np
      parameter(np=2)
      double precision DetGr,Gtwiddle(np,np),Gtt(np,np,np,np),
     . f(np)
      double complex Shat4zz(np,z1max,-2:0),Shat5zzzz(np,-2:0),
     . S0000i(np,-2:0),Shat5zz(np,z2max,-2:0),bit,pole

      do ep=-2,0
      bit=czip
      do n=1,np
      do m=1,np
      bit=bit
     . +Gtt(k,n,l,m)*(f(n)*Shat4zz(m,i1,ep)
     . +2d0*(delta(n,i1)*Shat5zzzz(m,ep))
     . -f(n)*f(m)*Cv(czzi(i1)+N0,ep)
     . -2d0*(f(n)*delta(m,i1)+f(m)*delta(n,i1))*Cv(cc0000+N0,ep))
      enddo
      enddo
      pole=czip
      if (ep .gt. -2)
     . pole=-4d0*Gtwiddle(k,l)*Cv(czzzzi(i1)+N0,ep-1)
      Cv(czzzzi(i1)+N0,ep)=
     . -(pole
     . +DetGr*Cv(czziii(z3(k,l,i1))+N0,ep)
     . -Gtwiddle(k,l)*S0000i(i1,ep)
     . -Gtwiddle(1,l)*Shat5zz(1,z2(k,i1),ep)
     . -Gtwiddle(2,l)*Shat5zz(2,z2(k,i1),ep)
     . +Gtwiddle(k,l)
     . *(Shat5zz(1,z2(1,i1),ep)
     .  +Shat5zz(2,z2(2,i1),ep))
     . +bit)/(14d0*Gtwiddle(k,l))

c      if (ep .eq. -1) then
c      write(6,*) 'C0000i ',ep,k,l,i1,
c     . pole
c      write(6,*) 'C0000i ',ep,k,l,i1,
c     . +DetGr*Cv(czziii(z3(k,l,i1))+N0,ep)
c      write(6,*) 'C0000i ',ep,k,l,i1,
c     . -Gtwiddle(k,l)*S0000i(i1,ep)
c      write(6,*) 'C0000i ',ep,k,l,i1,
c     . -Gtwiddle(1,l)*Shat5zz(1,z2(k,i1),ep)
c      write(6,*) 'C0000i ',ep,k,l,i1,
c     . -Gtwiddle(2,l)*Shat5zz(2,z2(k,i1),ep)
c      write(6,*) 'C0000i ',ep,k,l,i1,
c     . +Gtwiddle(k,l)
c     . *(Shat5zz(1,z2(1,i1),ep)
c     .  +Shat5zz(2,z2(2,i1),ep))
c      write(6,*) 'C0000i ',ep,k,l,i1,
c     . +bit
c      endif

      enddo
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
	SUBROUTINE GD_FLTM  ( filnms, nfiles, maxt,
     +			      ngdftm, gdtlst, iret )
C************************************************************************
C* GD_FLTM								*
C*									*
C* This subroutine returns all the times present in a list of grid	*
C* files.								*
C*									*
C* GD_FLTM  ( FILNMS, NFILES, MAXT, NGDFTM, GDTLST, IRET )		*
C*									*
C* Input parameters:							*
C*	FILNMS(LLMXGT)	CHAR*256	Grid file names array		*
C*	NFILES		INTEGER		Number of grid filename		*
C*	MAXT  		INTEGER		Maximum number of times allowed	*
C*									*
C* Output parameters:							*
C*	NGDFTM		INTEGER		Number of times			*
C*	GDTLST(*)	CHAR*		List of GEMPAK times		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -4 = file not open		*
C*					 -6 = read/write error		*
C**									*
C* Log:									*
C* D.W.Plummer/NCEP	 8/98						*
C* S. Jacobs/NCEP	 9/98	Changed maxtim to LLMXGT in GD_GTIM call*
C* T. Piper/GSC		 3/99	Corrected prolog			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'GMBDTA.CMN'
	INCLUDE 	'grdcmn.cmn'
C*
	CHARACTER	filnms(LLMXGT)*256
	CHARACTER*(*)	gdtlst(*)
	CHARACTER*20 	timarr (LLMXGT)
C
	REAL		anlblk(128), rnvblk(256)
C-----------------------------------------------------------------------
	iret = 0
	ngdftm = 0
C
C*	    Get all times from the files returned.
C
	    ngdftm = 0
	    DO  nf = 1, nfiles
C
	        CALL ST_LSTR ( filnms(nf), lf, iret )
	        CALL GD_OFIL ( filnms(nf)(:lf), .false., .true.,
     +			       igdfln, navsz, rnvblk, ianlsz, anlblk,
     +			       ihdrsz, maxgrd, iret )
C
	        IF ( iret .eq. 0 )  THEN
C
	            CALL GD_GTIM ( igdfln, LLMXGT, timarr, nt, iret )
C
C*	            Loop through the times returned.
C
	            DO  i = 1, nt
C
C*	                Convert the times from integers to characters.
C
		        IF ( ngdftm .lt. maxt )  THEN
	                    ngdftm = ngdftm + 1
	                    CALL TG_ITOC  ( igdatm (1,i,igdfln),
     +				            gdtlst (ngdftm), ier )
		        END IF
	            END DO
C
	            CALL GD_CLOS ( igdfln, ier )
C
	        END IF
C
	    END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! { dg-do compile }
!
! PR 82932: [OOP] ICE in update_compcall_arglist, at fortran/resolve.c:5837
!
! Contributed by Janus Weil <janus@gcc.gnu.org>

module m

   implicit none

   type, abstract :: AT
   contains
      procedure(init_ifc),    deferred :: sinit
      procedure(missing_ifc), deferred :: missing
      generic :: init    => sinit
   end type

   abstract interface
      subroutine init_ifc(data)
         import AT
         class(AT) :: data
      end subroutine
      subroutine missing_ifc(data)
         import AT
         class(AT) :: data
      end subroutine
   end interface

end module

program p

   use m

   implicit none

   type, extends(AT) :: ET  ! { dg-error "must be ABSTRACT" }
   contains
      procedure :: sinit
   end type

   type(ET) :: c
   call c%init()

end


text:
	SUBROUTINE ST_UNPR  ( string, lenin, outstr, lenout, iret )
C************************************************************************
C* ST_UNPR								*
C*									*
C* This subroutine eliminates substrings of unprintable characters.	*
C* Substrings of control characters, i.e., characters less than a 	*
C* blank, are replaced by a single blank.  Characters greater than	*
C* '}' (CHAR (126)) are replaced by '~' (CHAR (127)).  This subroutine	*
C* can be used to replace control characters such as CR and LF with	*
C* a single blank.  Invalid characters in the ASCII character set	*
C* are replaced by '~' so that the lengths of fields in the record	*
C* will remain unchanged.  The input and output strings may be the	*
C* same variable.							*
C*									*
C* ST_UNPR  ( STRING, LENIN, OUTSTR, LENOUT, IRET )			*
C*									*
C* Input parameters:							*
C*	STRING		CHAR*		Input string			*
C*	LENIN		INTEGER		Length of input string		*
C*									*
C* Output parameters:							*
C*	OUTSTR		CHAR*		Output string			*
C*	LENOUT		INTEGER		Length of output string		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	12/87						*
C* M. desJardins/GSFC	 6/88	Documentation				*
C* J. Whistler/SSAI	 5/91	Put tilda definition in GEMPRM.PRM	*
C* D. Kidwell/NCEP	 6/04	Added lenomx for bug if ip > LEN(outstr)*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	CHARACTER*(*)	string, outstr
C*
	LOGICAL		bl
	CHARACTER	c*1
C------------------------------------------------------------------------
	iret   = 0
	bl     = .true.
	ip     = 0
	lenomx = LEN ( outstr )
C
C*	Check each character to see if it is a control character,
C*	a blank or out of range.
C
	DO  i = 1, lenin
	    IF ( ip .lt. lenomx ) THEN
	        c = string (i:i)
	        IF  ( ( c .ge. '!' ) .and. ( c .le. '}' ) )  THEN
C
C*		    Add non blanks to the string.
C
		    ip = ip + 1
		    outstr ( ip: ip ) = c
		    bl = .false.
	          ELSE IF  ( c .gt. '}' )  THEN
C
C*		    Replace out of range characters with a tilda.
C
		    ip = ip + 1
		    outstr ( ip : ip ) = CHTLDA
		    bl = .false.
C
C*		    Add a blank to the list only if the last character
C*		    is not already a blank.
C
	          ELSE IF ( .not. bl ) THEN
		    bl = .true.
		    ip = ip + 1
		    outstr ( ip : ip ) = ' '
	        END IF
	    END IF
	END DO
C
C*	Blank out remainder of string.
C
 	IF ( ip .lt. lenomx ) outstr ( ip+1 : lenomx ) = ' '
C
C*	Set output length.  If last character was a blank, eliminate it.
C
	IF  ( ( ip .gt. 0 ) .and. ( outstr (ip:ip) .eq. ' ' ) ) THEN
	    lenout = ip - 1
	  ELSE
	    lenout = ip
	END IF
C*
	RETURN
	END


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ############################################################
c     ##                                                        ##
c     ##  subroutine katom  --  atom type parameter assignment  ##
c     ##                                                        ##
c     ############################################################
c
c
c     "katom" assigns an atom type definitions to each atom in
c     the structure and processes any new or changed values
c
c
      subroutine katom
      implicit none
      include 'sizes.i'
      include 'atmtyp.i'
      include 'atoms.i'
      include 'couple.i'
      include 'inform.i'
      include 'iounit.i'
      include 'katoms.i'
      include 'keys.i'
      integer i,k,next
      integer cls,atn,lig
      real*8 wght
      logical header
      character*3 symb
      character*20 keyword
      character*24 notice
      character*120 record
      character*120 string
c
c
c     process keywords containing atom type parameters
c
      header = .true.
      do i = 1, nkey
         next = 1
         record = keyline(i)
         call gettext (record,keyword,next)
         call upcase (keyword)
         if (keyword(1:5) .eq. 'ATOM ') then
            k = 0
            cls = 0
            symb = ' '
            notice = ' '
            atn = 0
            wght = 0.0d0
            lig = 0
            call getnumb (record,k,next)
            call getnumb (record,cls,next)
            if (cls .eq. 0)  cls = k
            atmcls(k) = cls
            call gettext (record,symb,next)
            call getstring (record,notice,next)
            string = record(next:120)
            read (string,*,err=40,end=40)  atn,wght,lig
            if (k.ge.1 .and. k.le.maxtyp) then
               if (header .and. .not.silent) then
                  header = .false.
                  write (iout,10)
   10             format (/,' Additional Atom Type Parameters :',
     &                    //,5x,'Type  Class  Symbol  Description',
     &                       15x,'Atomic',4x,'Mass',3x,'Valence',/)
               end if
               symbol(k) = symb
               describe(k) = notice
               atmnum(k) = atn
               weight(k) = wght
               ligand(k) = lig
               if (.not. silent) then
                  write (iout,20)  k,cls,symb,notice,atn,wght,lig
   20             format (2x,i6,1x,i6,5x,a3,3x,a24,i6,f11.3,i6)
               end if
            else if (k .ge. maxtyp) then
               write (iout,30)
   30          format (/,' KATOM   --  Too many Atom Types;',
     &                    ' Increase MAXTYP')
               abort = .true.
            end if
   40       continue
         end if
      end do
c
c     transfer atom type values to individual atoms
c
      do i = 1, n
         k = type(i)
         if (k .eq. 0) then
            class(i) = 0
            atomic(i) = 0
            mass(i) = 0.0d0
            valence(i) = 0
            story(i) = 'Undefined Atom Type     '
         else
            if (symbol(k) .ne. '   ')  name(i) = symbol(k)
            class(i) = atmcls(k)
            atomic(i) = atmnum(k)
            mass(i) = weight(k)
            valence(i) = ligand(k)
            story(i) = describe(k)
         end if
      end do
c
c     process keywords containing atom types for specific atoms
c
      header = .true.
      do i = 1, nkey
         next = 1
         record = keyline(i)
         call gettext (record,keyword,next)
         call upcase (keyword)
         if (keyword(1:5) .eq. 'ATOM ') then
            k = 0
            symb = ' '
            notice = ' '
            atn = 0
            wght = 0.0d0
            lig = 0
            call getnumb (record,k,next)
            call getnumb (record,cls,next)
            call gettext (record,symb,next)
            call getstring (record,notice,next)
            string = record(next:120)
            read (string,*,err=70,end=70)  atn,wght,lig
            if (k.lt.0 .and. k.ge.-n) then
               if (header .and. .not.silent) then
                  header = .false.
                  write (iout,50)
   50             format (/,' Additional Atom Types for',
     &                       ' Specific Atoms :',
     &                    //,5x,'Atom  Class  Symbol  Description',
     &                       15x,'Atomic',4x,'Mass',3x,'Valence',/)
               end if
               k = -k
               if (cls .eq. 0)  cls = k
               class(k) = cls
               name(k) = symb
               story(k) = notice
               atomic(k) = atn
               mass(k) = wght
               valence(k) = lig
               if (.not. silent) then
                  write (iout,60)  k,cls,symb,notice,atn,wght,lig
   60             format (2x,i6,1x,i6,5x,a3,3x,a24,i6,f11.3,i6)
               end if
            end if
   70       continue
         end if
      end do
c
c     check for presence of undefined atom types or classes
c
      header = .true.
      do i = 1, n
         k = type(i)
         cls = class(i)
         if (k.lt.1 .or. k.gt.maxtyp
     &          .or. cls.lt.1 .or. cls.gt.maxclass) then
            abort = .true.
            if (header) then
               header = .false.
               write (iout,80)
   80          format (/,' Undefined Atom Types or Classes :',
     &                 //,' Type',10x,'Atom Number',5x,'Atom Type',
     &                    5x,'Atom Class',/)
            end if
            write (iout,90)  i,k,cls
   90       format (' Atom',12x,i5,10x,i5,10x,i5)
         end if
      end do
c
c     check the number of atoms attached to each atom
c
      header = .true.
      do i = 1, n
         if (n12(i) .ne. valence(i)) then
            if (header) then
               header = .false.
               write (iout,100)
  100          format (/,' Atoms with an Unusual Number of Attached',
     &                    ' Atoms :',
     &                 //,' Type',11x,'Atom Name',6x,'Atom Type',7x,
     &                    'Expected',4x,'Found',/)
            end if
            write (iout,110)  i,name(i),type(i),valence(i),n12(i)
  110       format (' Valence',7x,i5,'-',a3,8x,i5,10x,i5,5x,i5)
         end if
      end do
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
module chemistry_module

  use amrex_fort_module, only : amrex_real
  implicit none

  integer, save :: nelements   ! number of elements
  integer, save :: nspecies    ! number of species
  integer, save :: nreactions  ! number of reactions

  logical, save :: chemistry_initialized = .false.

  integer, parameter :: L_elem_name = 3 ! Each element name has at most 3 characters
  character*(L_elem_name), allocatable, save :: elem_names(:)

  integer, parameter :: L_spec_name = 8 ! Each species name has at most 8 characters
  character*(L_spec_name), allocatable, save :: spec_names(:)

  real(amrex_real), allocatable, save :: molecular_weight(:), inv_mwt(:)

  real(amrex_real), save :: Ru, Ruc, Patm

contains

end module chemistry_module

max_stars_repo_name: marient/PelePhysics

max_issues_repo_name: marient/PelePhysics

max_forks_repo_name: marient/PelePhysics