text:
      subroutine configs_and_props_inc_chooser()
c For a given nFKSprocess, it fills the c_configs_inc common block with
c the configs.inc information (i.e. IFOREST(), SPROP(), TPRID() and
c MAPCONFIG())
      implicit none
      include 'nexternal.inc'
      include 'coupl.inc'
      double precision ZERO
      parameter (ZERO=0d0)
      include 'maxparticles.inc'
      include 'maxconfigs.inc'
      integer i,j,k
      INTEGER NFKSPROCESS
      COMMON/C_NFKSPROCESS/NFKSPROCESS
      integer iforest(2,-max_branch:-1,lmaxconfigs)
      integer sprop(-max_branch:-1,lmaxconfigs)
      integer tprid(-max_branch:-1,lmaxconfigs)
      integer mapconfig(0:lmaxconfigs)
      common/c_configs_inc/iforest,sprop,tprid,mapconfig
      double precision prmass(-max_branch:nexternal,lmaxconfigs)
      double precision prwidth(-max_branch:-1,lmaxconfigs)
      integer prow(-max_branch:-1,lmaxconfigs)
      common/c_props_inc/prmass,prwidth,prow
      double precision pmass(nexternal)
      logical firsttime
      data firsttime /.true./
      include 'configs_and_props_decl.inc'
      save mapconfig_d,iforest_d,sprop_d,tprid_d,pmass_d,pwidth_d,pow_d
      include "pmass.inc"
c
      if (max_branch_used.gt.max_branch) then
         write (*,*) 'ERROR in configs_and_props_inc_chooser:'/
     $        /' increase max_branch',max_branch,max_branch_used
         stop
      endif
      if (lmaxconfigs_used.gt.lmaxconfigs) then
         write (*,*) 'ERROR in configs_and_propsinc_chooser:'/
     $        /' increase lmaxconfigs' ,lmaxconfigs,lmaxconfigs_used
         stop
      endif

C the configurations and propagators infos are read at the first
C evaluation
      if (firsttime) then
        call read_configs_and_props_info(mapconfig_d,iforest_d,sprop_d,
     1                                   tprid_d,pmass_d,pwidth_d,pow_d)
        firsttime = .false.
      endif
c
c Fill the arrays of the c_configs_inc and c_props_inc common
c blocks. Some of the information might not be available in the
c configs_and_props_info.inc include file, but there is no easy way of skipping
c it. This will simply fill the common block with some bogus
c information.
      do i=0,MAPCONFIG_D(nFKSprocess,0)
         mapconfig(i)=MAPCONFIG_D(nFKSprocess,i)
         if (i.ne.0) then
            do j=-max_branch_used,-1
               do k=1,2
                  iforest(k,j,i)=IFOREST_D(nFKSprocess,k,j,i)
               enddo
               sprop(j,i)=SPROP_D(nFKSprocess,j,i)
               tprid(j,i)=TPRID_D(nFKSprocess,j,i)
               prmass(j,i)=PMASS_D(nFKSprocess,j,i)
               prwidth(j,i)=PWIDTH_D(nFKSprocess,j,i)
               prow(j,i)=POW_D(nFKSprocess,j,i)
            enddo
c for the mass, also fill for the external masses
            prmass(0,i)=0d0
            do j=1,nexternal
               prmass(j,i)=pmass(j)
            enddo
         endif
      enddo
c
      return
      end

      subroutine fks_inc_chooser()
c For a given nFKSprocess, it fills the c_fks_inc common block with the
c fks.inc information
      implicit none
      include 'nexternal.inc'
      include 'fks_info.inc'
      integer i,j
      INTEGER NFKSPROCESS
      COMMON/C_NFKSPROCESS/NFKSPROCESS
      integer fks_j_from_i(nexternal,0:nexternal)
     &     ,particle_type(nexternal),pdg_type(nexternal)
      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
      double precision particle_charge(nexternal), particle_charge_born(nexternal-1)
      common /c_charges/particle_charge
      common /c_charges_born/particle_charge_born
      integer i_fks,j_fks
      common/fks_indices/i_fks,j_fks
      integer i_type,j_type,m_type
      double precision ch_i,ch_j,ch_m
      integer particle_type_born(nexternal-1)
      common /c_particle_type_born/particle_type_born
      logical particle_tag(nexternal)
      common /c_particle_tag/particle_tag
      logical particle_tag_born(nexternal-1)
      common /c_particle_tag/particle_tag_born
      logical need_color_links, need_charge_links
      common /c_need_links/need_color_links, need_charge_links
      integer extra_cnt, isplitorder_born, isplitorder_cnt
      common /c_extra_cnt/extra_cnt, isplitorder_born, isplitorder_cnt
      include 'orders.inc'
      logical split_type(nsplitorders)
      common /c_split_type/split_type
      logical is_aorg(nexternal)
      common /c_is_aorg/is_aorg
      logical is_charged(nexternal)
      common /c_is_charged/is_charged
c
      i_fks=fks_i_D(nFKSprocess)
      j_fks=fks_j_D(nFKSprocess)
      extra_cnt = extra_cnt_D(nFKSprocess)
      isplitorder_born = isplitorder_born_D(nFKSprocess)
      isplitorder_cnt = isplitorder_cnt_D(nFKSprocess)
      need_color_links=need_color_links_d(nFKSprocess)
      need_charge_links=need_charge_links_d(nFKSprocess)
      do i=1,nexternal
         if (fks_j_from_i_D(nFKSprocess,i,0).ge.0 .and.
     &        fks_j_from_i_D(nFKSprocess,i,0).le.nexternal) then
            do j=0,fks_j_from_i_D(nFKSprocess,i,0)
               fks_j_from_i(i,j)=fks_j_from_i_D(nFKSprocess,i,j)
            enddo
         else
            write (*,*) 'ERROR in fks_inc_chooser'
            stop
         endif
         particle_type(i)=particle_type_D(nFKSprocess,i)
         particle_tag(i)=particle_tag_D(nFKSprocess,i)
         particle_charge(i)=particle_charge_D(nFKSprocess,i)
         pdg_type(i)=pdg_type_D(nFKSprocess,i)
         ! is_aorg is true if the particle can induce soft singularities
         ! currencly it is based on the pdg id (photon or gluon)
         is_aorg(i) = abs(pdg_type(i)).eq.21.or.pdg_type(i).eq.22
         ! is_charged is true if the particle has any color or electric
         ! charge
         is_charged(i) = particle_type(i).ne.1.or.
     &                   particle_charge(i).ne.0d0
      enddo
      do i=1,nexternal
         if (i.lt.min(i_fks,j_fks)) then
            particle_type_born(i)=particle_type(i)
            particle_charge_born(i)=particle_charge(i)
            particle_tag_born(i)=particle_tag(i)
         elseif (i.gt.max(i_fks,j_fks)) then
            particle_type_born(i-1)=particle_type(i)
            particle_charge_born(i-1)=particle_charge(i)
            particle_tag_born(i-1)=particle_tag(i)
         elseif (i.eq.min(i_fks,j_fks)) then
            i_type=particle_type(i_fks)
            j_type=particle_type(j_fks)
            ch_i=particle_charge(i_fks)
            ch_j=particle_charge(j_fks)
            particle_tag_born(i) = particle_tag(j_fks)
            call get_mother_col_charge(i_type,ch_i,j_type,ch_j,m_type,ch_m)
            particle_type_born(i)=m_type
            particle_charge_born(i)=ch_m
         elseif (i.ne.max(i_fks,j_fks)) then
            particle_tag_born(i) = particle_tag(i)
            particle_type_born(i)=particle_type(i)
            particle_charge_born(i)=particle_charge(i)
         endif
      enddo

      do i = 1, nsplitorders
         split_type(i) = split_type_d(nFKSprocess,i)
      enddo
      return
      end

      subroutine leshouche_inc_chooser()
c For a given nFKSprocess, it fills the c_leshouche_inc common block with the
c leshouche.inc information
      implicit none
      include 'nexternal.inc'
      include 'genps.inc'
      integer i,j,k
      INTEGER NFKSPROCESS
      COMMON/C_NFKSPROCESS/NFKSPROCESS
      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
     &     icolup(2,nexternal,maxflow),niprocs
      common /c_leshouche_inc/idup,mothup,icolup,niprocs
      logical firsttime
      data firsttime /.true./
      include 'leshouche_decl.inc'
      common/c_leshouche_idup_d/ idup_d
      save mothup_d, icolup_d, niprocs_d

c
      if (maxproc_used.gt.maxproc) then
         write (*,*) 'ERROR in leshouche_inc_chooser: increase maxproc',
     &        maxproc,maxproc_used
         stop
      endif
      if (maxflow_used.gt.maxflow) then
         write (*,*) 'ERROR in leshouche_inc_chooser: increase maxflow',
     &        maxflow,maxflow_used
         stop
      endif

      if (firsttime) then
        call read_leshouche_info(idup_d,mothup_d,icolup_d,niprocs_d)
        firsttime = .false.
      endif

      niprocs=niprocs_d(nFKSprocess)
      do j=1,niprocs
         do i=1,nexternal
            IDUP(i,j)=IDUP_D(nFKSprocess,i,j)
            MOTHUP(1,i,j)=MOTHUP_D(nFKSprocess,1,i,j)
            MOTHUP(2,i,j)=MOTHUP_D(nFKSprocess,2,i,j)
         enddo
      enddo
c
      do j=1,maxflow_used
         do i=1,nexternal
            ICOLUP(1,i,j)=ICOLUP_D(nFKSprocess,1,i,j)
            ICOLUP(2,i,j)=ICOLUP_D(nFKSprocess,2,i,j)
         enddo
      enddo
c
      return
      end

      subroutine read_configs_and_props_info(mapconfig_d,iforest_d
     $     ,sprop_d,tprid_d,pmass_d,pwidth_d,pow_d)
C read the various information from the configs_and_props_info.dat file
      implicit none
      integer i,j,k
      integer ndau, idau, dau, id
      character *200 buff
      double precision get_mass_from_id, get_width_from_id
      include 'configs_and_props_decl.inc'
      mapconfig_d=0
      iforest_d=0
      sprop_d=0
      tprid_d=0
      pmass_d=0d0
      pwidth_d=0d0
      pow_d=0
      open(unit=78, file='configs_and_props_info.dat', status='old')
      do while (.true.)
        read(78,'(a)',end=999) buff
        if (buff(:1).eq.'#') cycle
        if (buff(:1).eq.'C') then
        ! mapconfig
        ! C  i   j   k -> MAPCONFIG_D(i,j)=k
          read(buff(2:),*) i,j,k
          mapconfig_d(i,j) = k
        else if (buff(:1).eq.'F') then
        ! iforest
        ! after the first line there are as many lines
        ! as the daughters
        ! F  i   j   k  ndau
        ! D dau_1
        ! D ...
        ! D dau_ndau        -> IFORREST_D(i,idau,i,k)=dau_idau
          read(buff(2:),*) i,j,k,ndau
          do idau=1,ndau
            read(78,'(a)') buff
            if (buff(:1).ne.'D') then
              write(*,*) 'ERROR #1 in read_configs_and_props_info',
     1                    i,j,k,ndau,buff
              stop
            endif
            read(buff(2:),*) dau
            iforest_d(i,idau,j,k) = dau
          enddo
        else if (buff(:1).eq.'S') then
        ! sprop
        ! S  i   j   k  id -> SPROP_D(i,j,k)=id
          read(buff(2:),*) i,j,k,id
          sprop_d(i,j,k) = id
        else if (buff(:1).eq.'T') then
        ! tprid
        ! T  i   j   k  id -> TPRID_D(i,j,k)=id
          read(buff(2:),*) i,j,k,id
          tprid_d(i,j,k) = id
        else if (buff(:1).eq.'M') then
        ! pmass and pwidth
          read(buff(2:),*) i,j,k,id
        ! M  i   j   k  id -> gives id of particle of which
        ! the mass/width is stored in PMASS/WIDTH_D(i,j,k)
          pmass_d(i,j,k) = get_mass_from_id(id)
          pwidth_d(i,j,k) = get_width_from_id(id)
        else if (buff(:1).eq.'P') then
        ! pow
        ! P  i   j   k  id -> POW_D(i,j,k)=id
          read(buff(2:),*) i,j,k,id
          pow_d(i,j,k) = id
        endif
      enddo
 999  continue
      close(78)

      return
      end

      subroutine read_leshouche_info(idup_d,mothup_d,icolup_d,niprocs_d)
C read the various information from the configs_and_props_info.dat file
      implicit none
      include "nexternal.inc"
      integer itmp_array(nexternal)
      integer i,j,k,l
      character *200 buff
      include 'leshouche_decl.inc'
      include 'nFKSconfigs.inc'
      include 'fks_info.inc'
      include 'born_maxamps.inc'
      integer idup(nexternal,maxproc)
      integer mothup(2,nexternal,maxproc)
      integer icolup(2,nexternal,maxflow)
      include 'born_leshouche.inc'
      if (fks_configs.eq.1) then
         if (pdg_type_d(1,fks_i_d(1)).eq.-21) then
c SPECIAL for [LOonly=QCD] process. Simply use the information from the
c born_leshouche.inc file.
            do j=1,maxproc
               do i=1,nexternal-1
                  idup_d(1,i,j)=idup(i,j)
                  mothup_d(1,1,i,j)=mothup(1,i,j)
                  mothup_d(1,2,i,j)=mothup(2,i,j)
               enddo
               idup_d(1,nexternal,j)=-21
               mothup_d(1,1,nexternal,j)=mothup(1,fks_j_d(1),j)
               mothup_d(1,2,nexternal,j)=mothup(2,fks_j_d(1),j)
            enddo
            do j=1,maxflow
               do i=1,nexternal-1
                  icolup_d(1,1,i,j)=icolup(1,i,j)
                  icolup_d(1,2,i,j)=icolup(2,i,j)
               enddo
               icolup_d(1,1,nexternal,j)=-99999 ! should not be used
               icolup_d(1,2,nexternal,j)=-99999
            enddo
            niprocs_d(1)=maxproc_used
            return
         endif
      endif

      open(unit=78, file='leshouche_info.dat', status='old')
      do while (.true.)
        read(78,'(a)',end=999) buff
        if (buff(:1).eq.'#') cycle
        if (buff(:1).eq.'I') then
        ! idup
        ! I  i   j   id1 ..idn -> IDUP_D(i,k,j)=idk
          read(buff(2:),*) i,j,(itmp_array(k),k=1,nexternal)
          do k=1,nexternal
            idup_d(i,k,j)=itmp_array(k)
          enddo
          niprocs_d(i)=j
        else if (buff(:1).eq.'M') then
        ! idup
        ! I  i   j   l   id1 ..idn -> MOTHUP_D(i,j,k,l)=idk
          read(buff(2:),*) i,j,l,(itmp_array(k),k=1,nexternal)
          do k=1,nexternal
            mothup_d(i,j,k,l)=itmp_array(k)
          enddo
        else if (buff(:1).eq.'C') then
        ! idup
        ! I  i   j   l   id1 ..idn -> ICOLUP_D(i,j,k,l)=idk
          read(buff(2:),*) i,j,l,(itmp_array(k),k=1,nexternal)
          do k=1,nexternal
            icolup_d(i,j,k,l)=itmp_array(k)
          enddo
        endif
      enddo
 999  continue
      close(78)

      return
      end

      subroutine get_mother_col_charge(i_type, ch_i, j_type, ch_j,
     $     m_type, ch_m)
C Given the type (color representation) and charges of i and j, return
C the type and charges of the mother particle
      implicit none
      integer i_type, j_type, m_type
      double precision ch_i, ch_j, ch_m
      include 'nexternal.inc'
      integer i_fks,j_fks

      common/fks_indices/i_fks,j_fks

      if (abs(i_type).eq.abs(j_type) .and.
     &    abs(ch_i).eq.abs(ch_j) .and.
     &    abs(i_type).gt.1) then
        ! neutral color octet splitting
         m_type=8
         ch_m = 0d0
         if ( (j_fks.le.nincoming .and.
     &         abs(i_type).eq.3 .and. j_type.ne.i_type) .or.
     &        (j_fks.gt.nincoming .and.
     &         abs(i_type).eq.3 .and. j_type.ne.-i_type)) then
            write(*,*)'Flavour mismatch #1col in get_mother_col_charge',
     &                 i_fks,j_fks,i_type,j_type
            stop
         endif
      elseif (abs(i_type).eq.abs(j_type) .and.
     &        dabs(ch_i).eq.dabs(ch_j).and.abs(i_type).eq.1) then
        ! neutral color singlet splitting
         m_type=1
         ch_m = 0d0
         if ( (j_fks.le.nincoming .and.
     &         dabs(ch_i).gt.0d0 .and. ch_j.ne.ch_i) .or.
     &        (j_fks.gt.nincoming .and.
     &         dabs(ch_i).gt.0d0 .and. ch_j.ne.-ch_i)) then
            write(*,*)'Flavour mismatch #1chg in get_mother_col_charge',
     &                 i_fks,j_fks,ch_i,ch_j
            stop
         endif
      elseif ((abs(i_type).eq.3 .and. j_type.eq.8) .or.
     &        (dabs(ch_i).gt.0d0 .and. ch_j.eq.0d0) ) then
         if(j_fks.le.nincoming)then
            m_type=-i_type
            if (m_type.eq.-1) m_type=1
            ch_m = -ch_i
         else
            write(*,*) 'Error in get_mother_col_charge: (i,j)=(q,g)'
            stop
         endif
      elseif ((i_type.eq.8 .and. abs(j_type).eq.3) .or.
     &         (ch_i.eq.0d0 .and. dabs(ch_j).gt.0d0) ) then
         m_type=j_type
         ch_m= ch_j
!     special for processes without a proper NLO contribution (i.e., PDG(i_fks)=-21)
      elseif (i_type.eq.8 .and. j_type.eq.1 .and. ch_i.eq.0d0 .and.
     $        ch_j.eq.0d0) then
         m_type=0
         ch_m=0d0
         continue
      else
         write(*,*)'Flavour mismatch #2 in get_mother_col_charge',
     &      i_type,j_type,m_type
         stop
      endif
      return
      end

      subroutine set_pdg(ict,iFKS)
c fills the pdg and pdg_uborn variables. It uses only the 1st IPROC. For
c the pdg_uborn (the PDG codes for the underlying Born process) the PDG
c codes of i_fks and j_fks are combined to give the PDG code of the
c mother and the extra (n+1) parton is given the PDG code of the gluon.
      use weight_lines
      implicit none
      include 'nexternal.inc'
      include 'fks_info.inc'
      include 'genps.inc'
      integer k,ict,iFKS
      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
     $     icolup(2,nexternal,maxflow),niprocs
      common /c_leshouche_inc/idup,mothup,icolup,niprocs
      include 'orders.inc'
      do k=1,nexternal
         pdg(k,ict)=idup(k,1)
      enddo
      do k=1,nexternal
         if (k.lt.fks_j_d(iFKS)) then
            pdg_uborn(k,ict)=pdg(k,ict)
         elseif(k.eq.fks_j_d(iFKS)) then
            if ( abs(pdg(fks_i_d(iFKS),ict)) .eq.
     &           abs(pdg(fks_j_d(iFKS),ict)) .and.
     &           abs(pdg(fks_i_d(iFKS),ict)).ne.21.and.
     &           abs(pdg(fks_i_d(iFKS),ict)).ne.22) then
c gluon splitting:  g/a -> ff
               !!!pdg_uborn(k,ict)=21
               ! check if any extra cnt is needed
               if (extra_cnt_d(iFKS).eq.0) then
                  ! if not, assign photon/gluon depending on split_type
                  if (split_type_d(iFKS,qcd_pos)) then
                    pdg_uborn(k,ict)=21
                  else if (split_type_d(iFKS,qed_pos)) then
                    pdg_uborn(k,ict)=22
                  else
                    write (*,*) 'set_pdg ',
     &                'ERROR#1 in PDG assigment for underlying Born'
                    stop 1
                  endif
               else
                  ! if there are extra cnt's, assign the pdg of the
                  ! mother in the born (according to isplitorder_born_d)
                  if (isplitorder_born_d(iFKS).eq.qcd_pos) then
                    pdg_uborn(k,ict)=21
                  else if (isplitorder_born_d(iFKS).eq.qcd_pos) then
                    pdg_uborn(k,ict)=22
                  else
                    write (*,*) 'set_pdg ',
     &                'ERROR#2 in PDG assigment for underlying Born'
                    stop 1
                  endif
               endif
            elseif (abs(pdg(fks_i_d(iFKS),ict)).eq.21.or.
     &              abs(pdg(fks_i_d(iFKS),ict)).eq.22) then
c final state gluon radiation:  X -> Xg
               pdg_uborn(k,ict)=pdg(fks_j_d(iFKS),ict)
            elseif (pdg(fks_j_d(iFKS),ict).eq.21.or.
     &              pdg(fks_j_d(iFKS),ict).eq.22) then
c initial state gluon splitting (gluon is j_fks):  g -> XX
               pdg_uborn(k,ict)=-pdg(fks_i_d(iFKS),ict)
            else
               write (*,*)
     &          'set_pdg ERROR#3 in PDG assigment for underlying Born'
               stop 1
            endif
         elseif(k.lt.fks_i_d(iFKS)) then
            pdg_uborn(k,ict)=pdg(k,ict)
         elseif(k.eq.nexternal) then
            if (split_type_d(iFKS,qcd_pos)) then
              pdg_uborn(k,ict)=21  ! give the extra particle a gluon PDG code
            elseif (split_type_d(iFKS,qed_pos)) then
              pdg_uborn(k,ict)=22  ! give the extra particle a photon PDG code
            endif
         elseif(k.ge.fks_i_d(iFKS)) then
            pdg_uborn(k,ict)=pdg(k+1,ict)
         endif
      enddo
      return
      end


text:
!
! This file contains a code-snippet that is #included in the subroutine "update_U" in domain_mod. It allows us to make the
! variable "friction_power_depth" a parameter. For efficiency it's good for the code to know the power-law power in the friction
! term.
!

    !subroutine update_U(domain, dt)
    !    !!
    !    !! Update the values of domain%U (i.e. the main flow variables), based on the fluxes and sources in domain
    !    !!
    !    class(domain_type), intent(inout):: domain
    !    real(dp), intent(in):: dt !! Timestep to advance

        real(dp):: inv_cell_area_dt, depth, implicit_factor, dt_gravity, fs
        integer(ip):: j, i, kk

        domain%dt_last_update = dt

        !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(domain, dt)
        dt_gravity = dt * gravity
        !$OMP DO SCHEDULE(STATIC)
        do j = 1, domain%nx(2)
            ! For spherical coordiantes, cell area changes with y.
            ! For cartesian coordinates this could be moved out of the loop
            inv_cell_area_dt = dt / domain%area_cell_y(j)

            !$OMP SIMD PRIVATE(kk)
            do i = 1, domain%nx(1)
                ! Fluxes
                do kk = 1, 3
                    domain%U(i,j,kk) = domain%U(i,j,kk) - inv_cell_area_dt * ( &
                        (domain%flux_NS(i, j+1, kk) - domain%flux_NS(i, j, kk)) + &
                        (domain%flux_EW(i+1, j, kk) - domain%flux_EW(i, j, kk) ))
                end do
            end do

            !$OMP SIMD PRIVATE(depth, fs, implicit_factor)
            do i = 1, domain%nx(1)

                ! Velocity clipping
                depth = domain%depth(i,j)

                if (domain%U(i,j,STG) <= (domain%U(i,j,ELV) + minimum_allowed_depth)) then
                    domain%U(i,j,UH) = ZERO_dp
                    domain%U(i,j,VH) = ZERO_dp
                else
#ifndef NOFRICTION
                    ! Implicit friction slope update
                    ! U_new = U_last + U_explicit_update - dt_gravity*depth*friction_slope_multiplier*U_new -
                    !         dt * linear_friction_coeff * U_new

                    ! If we multiply this by UH or VH, we get the associated friction slope term
                    fs = domain%manning_squared(i,j) * &
                        sqrt(domain%velocity(i,j,UH) * domain%velocity(i,j,UH) + &
                             domain%velocity(i,j,VH) * domain%velocity(i,j,VH)) * &
                        max(depth, minimum_allowed_depth)**friction_power_depth
                        !exp(log(max(depth, minimum_allowed_depth))*NEG_SEVEN_ON_THREE_dp)

                    implicit_factor = ONE_dp/(ONE_dp + dt_gravity*max(depth, ZERO_dp)*fs + dt*domain%linear_friction_coeff)
                    !if(domain%U(i,j,STG) <= (domain%U(i,j,ELV) + minimum_allowed_depth)) implicit_factor = ZERO_dp
#else
                    implicit_factor = ONE_dp
#endif
                    ! Pressure gradients
                    domain%U(i,j,UH) = domain%U(i,j,UH) + inv_cell_area_dt * domain%explicit_source(i, j, UH)
                    ! Friction
                    domain%U(i,j,UH) = domain%U(i,j,UH) * implicit_factor

                    ! Here we add an extra source to take care of the j-1 pressure gradient addition, which
                    ! could not be updated in parallel (since j-1 might be affected by another OMP thread)
                    domain%U(i,j,VH) = domain%U(i,j,VH) + inv_cell_area_dt * &
                        (domain%explicit_source(i, j, VH) + domain%explicit_source_VH_j_minus_1(i, j+1))
                    ! Friction
                    domain%U(i,j,VH) = domain%U(i,j,VH) * implicit_factor
                end if
            end do
        end do
        !$OMP END DO
        !$OMP END PARALLEL
        domain%time = domain%time + dt

!    end subroutine

max_stars_repo_name: KayaGeoAus/ptha

max_issues_repo_name: KayaGeoAus/ptha

max_forks_repo_name: KayaGeoAus/ptha


text:
c     sub-program exchange
      program ffmod6
c     subroutine ffmod6 (iabs)

c     final calculations for various spectroscopies
c     (EXAFS, XANES, FPRIME, DANES, XES)
c     written by a.ankudinov 2000
c     modified by a.ankudinov 2001 for new I/O structure

c     INPUT: mod6.inp global.dat xsect.bin fms.bin list.dat and feff.pad
c     OUTPUT: xmu.dat (chi.dat for EXAFS)

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

      include '../HEADERS/dim.h'
      include '../HEADERS/parallel.h'
      integer elnes,ipmin,ipmax,ipstep  !KJ my variables 1-06
      integer absolu !KJ 3-06

      call par_begin
      if (worker) go to 400

c     sub-program exchange
      iabs = 1
c

c     open the log file, unit 11.  See subroutine wlog.
      open (unit=11, file='log6.dat', status='unknown', iostat=ios)
      call chopen (ios, 'log6.dat', 'feff')

c     read  input files
      call reff2x(mchi, ispec, ipr6, idwopt, critcw, s02, sig2g,
     1                   tk, thetad, mbconv, absolu,  !KJ added absolu 3-06
     1                   vrcorr, vicorr, alphat, thetae, nabs,
     4            elnes,ipmin,ipmax,ipstep)   !KJ added this line 1-06

      if (mchi .eq. 1)  then
         call wlog(' Calculating chi...')
         if (ispec.gt.0 .and. ispec.lt.3) then
c           using FMS+Paths method
            call ff2xmu (ispec, ipr6, idwopt, critcw, s02, sig2g,
     1                   tk, thetad, mbconv, absolu,  !KJ added absolu 3-06
     1                   vrcorr, vicorr, alphat, thetae, iabs, nabs,
     4            ipmin,ipmax,ipstep)   !KJ added this line  1-06
         elseif (ispec.eq.3 .or. ispec.eq.4) then
c           using FMS+Paths method
            call ff2afs ( ipr6, idwopt, critcw, s02, sig2g,
     1                   tk, thetad, mbconv, absolu,  !KJ added absolu 4-06
     1                   vrcorr, vicorr, alphat, thetae, iabs, nabs,
     4            ipmin,ipmax,ipstep)   !KJ added this line 4-06
         else
c           using MS Paths expansion
            call ff2chi (ispec, ipr6, idwopt, critcw, s02, sig2g,
     1           tk, thetad, mbconv, absolu, !KJ added absolu 3-06
     1           vrcorr, vicorr, alphat, thetae, iabs, nabs,
     4           ipmin,ipmax,ipstep) !KJ added this line 1-06
         endif
         call wlog(' Done with module 6: DW + final sum over paths.')
      endif
      close (unit=11)

  400 call par_barrier
      call par_end

c     sub-program exchange
      stop
c     return

      end

ext: f


text:
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c common macros for all routines
c define(good,0.0)
c define(goodpt,0.0)
c define(badpt,2.0)
c define(bad,2.0)
c define(yes,1)
c define(no,0)
c define(macheps,.0001)
c define(rinfinity, 10.e32)
c define(iinfinity, 1000000000)
c define(maxlv,12) - max. permitted levels of nesting
c define(maxcls,192) - max. number of clusters (corner dimensions)
c define(maxgr,192) - max. number of grids (node dimensions)
c define(rsize,12)
c define(nsize,17)
c define(maxstorage,allocsize)    allocsize set in calloc.i

c macros for the nodes of the grid tree
c for node (integer part of record)
c define(parent,1)
c define(bndryptr,2) - for perimeter pieces - phasing out, use as
c                      temp. storage for enlarged grids in adv. and ee.
c define(owner,3) - for use in multiprocessing version
c define(nestlevel,4)
c define(maxnumrow,5)
c define(maxnumcol,6)
c define(store1,7)
c define(store2,8)
c define(-------,9) - not used now
c define(levelptr,10)
c define(errptr,11)
c define(prevlevel,12)
c define(cbndryptr,13) - temporary storage for for coarsened perimeter
c                       pieces and errest temp.  storage
c define(permstore,14) - pointer to irreg array for each grid
c define(cfluxptr,15)
c define(ffluxptr,16)
c define (lstptr,17) - pointer to start of irregular list for that grid.

c for rnode (real part of record)
c define(corn1x,1)
c define(corn1y,2)
c define(corn2x,3)
c define(corn2y,4)
c define(corn3x,5)
c define(corn3y,6)
c define(corn4x,7)
c define(corn4y,8)
c define(hrow,9)
c define(hcol,10)
c define(ktime,11)
c define(timemult,12)

c tree - related macros:
c define(nextfree,2)
c define(null, 0)
c define(noinsect, 1)
c define(nil,  0)
c
c
c --------------------------------------------------------------
c
      subroutine advanc (level,nvar,dtlev,vtime,steady,lastout)
c
      implicit double precision (a-h,o-z)
      logical    rflag,vtime,steady
      logical    graf
      parameter  (maxgr = 192, maxlv=12)
      common /nodal/ rnode(12,maxgr),node(17,maxgr),lstart(maxlv),
     *  newstl(maxlv),
     *  listsp(maxlv),intrat(maxlv),tol,bzonex,bzoney,mstart,ndfree,
     *  lfine,iorder, mxnest,kcheck,lwidth,
     *  maxcl, graf, lhead
      common   /stats/   evol, rvol, rvoll(maxlv), lentot, lenmax
      include "calloc.i"
      data    rflag/.true./
      logical first/.false./, lastout
c
c  fill all boundary pieces of grids at this level. done
c  in this order to (a) help parallelize and (b) don't need
c  2 copies of finest level storage, since only needed
c  to fill boundary pieces for grids not yet advanced.
c  integrate all grids at 'level' by one step of its delta(t)
c
      mptr = lstart(level)
 3    continue
          maxi   = node(5,mptr)
          maxj   = node(6,mptr)
          maxip1 = maxi + 1
          maxjp1 = maxj + 1
          mitot  = maxi-1 + 2*lwidth
          mjtot  = maxj-1 + 2*lwidth
          locbig = igetsp(mitot*mjtot*nvar)
          node(2,mptr) = locbig
          locnew = node(7,mptr)
          xl   = rnode(1, mptr)
          yb   = rnode(2, mptr)
          xr   = rnode(7, mptr)
          yt   = rnode(4, mptr)
          hx   = rnode(9, mptr)
          hy   = rnode(10,mptr)
          time = rnode(12,mptr)
          intime = intcnt(level)
          call prem(alloc(locbig),alloc(locnew),nvar,maxip1,maxjp1,
     1              mitot,mjtot,lwidth)

c  #### read cart3d solution for restart
         if (first) then
            call readcart(alloc(locbig),mitot,mjtot)
         endif
          call bound(xl,xr,yb,yt,hx,hy,time,intime,level,nvar,lwidth,
     1               alloc(locbig),mitot,mjtot,alloc(node(14,mptr)),
     2               node(17,mptr))
        mptr = node(10, mptr)
        if (mptr .ne. 0) go to 3
c
      dtlev = 1.e32
      mptr  = lstart(level)
 5    continue
          locold = node(8, mptr)
          locnew = node(7, mptr)
          maxi   = node(5, mptr)
          maxj   = node(6, mptr)
          maxip1 = maxi + 1
          maxjp1 = maxj + 1
          delt   = possk(level)
          rnode(11,mptr) = delt
          hx     = rnode(9,mptr)
          hy     = rnode(10,mptr)
          time   = rnode(12,mptr)
c
          mitot  = maxi-1 + 2*lwidth
          mjtot  = maxj-1 + 2*lwidth
          locbig = node(2, mptr)
          locsc1 = igetsp(mitot*mjtot*nvar)
          locsc5 = igetsp(mitot*mjtot*nvar)
          locqx  = igetsp(mitot*mjtot*nvar)
          locqy  = igetsp(mitot*mjtot*nvar)
c
c  copy old soln. values into  next time steps soln. values
c  since integrator will overwrite it. only for grids not at
c  the finest level. finest level grids do not maintain copies
c  of old and new time solution values.
c
          if (level .lt. mxnest) then
             ntot   = maxip1 * maxjp1 * nvar
cdir$ ivdep
             do 10 i = 1, ntot
 10          alloc(locold + i - 1) = alloc(locnew + i - 1)
          endif
c
      xlow = rnode(1,mptr) - lwidth*hx
      ylow = rnode(2,mptr) - lwidth*hy
      rvol = rvol + (maxip1 - 2)*(maxjp1 - 2)
      rvoll(level) = rvoll(level) + (maxip1 - 2)*(maxjp1 - 2)
      locirr = node(14,mptr)
      locncount = locirr + mitot*mjtot
      locnumHoods = locncount + mitot*mjtot
c
      call method(alloc(locbig),alloc(locsc1),alloc(locsc5),
     1            alloc(locirr),mitot,mjtot,lwidth,
     2            delt,dtnew,node(17,mptr),hx,hy,rflag,
     3            iorder,xlow,ylow,mptr,vtime,steady,
     4            alloc(locqx),alloc(locqy),level,difmax,
     5            lastout,nvar,time)
       call postm(alloc(locbig),alloc(locnew),nvar,maxip1,maxjp1,
     1            mitot,mjtot,lwidth)
      if (node(15,mptr) .ne. 0)
     1   call fluxsv(mptr,alloc(locsc1),alloc(locsc5),
     2               alloc(node(15,mptr)),mitot,mjtot,
     3               nvar,listsp(level))
      if (node(16,mptr) .ne. 0) then
         lenbc = 2*((maxi-1)/intrat(level-1)+(maxj-1)/intrat(level-1))
         call fluxad(alloc(locsc1),alloc(locsc5),
     1               alloc(node(16,mptr)),mptr,mitot,mjtot,maxi,maxj,
     2               nvar,lenbc,intrat(level-1))
      endif
c
          call reclam(locsc1,mitot*mjtot*nvar)
          call reclam(locsc5,mitot*mjtot*nvar)
          call reclam(locqx, mitot*mjtot*nvar)
          call reclam(locqy, mitot*mjtot*nvar)
          call reclam(locbig,mitot*mjtot*nvar)
c
          if (vtime) then
             dtlev = dmin1(dtlev,dtnew)
          endif

c
          rnode(12,mptr)  = rnode(12,mptr)+delt
          mptr            = node(10, mptr)
          if (mptr .ne. 0) go to 5
c
      return
      end
c
c
      subroutine readcart(val,mitot,mjtot)
      implicit double precision(a-h,o-z)
      dimension val(mitot,mjtot,4)
      dimension q(4)

10    continue
      read(8,100,end=99) i,j,(q(k),k=1,4)
 100  format(23x,i3,3x,i3,4x,4e16.7)
c     write(6,*)i,j,(q(k),k=1,4)
c     q is in primitive vars
      val(i,j,1) = q(1)
      val(i,j,2) = q(1)*q(2)
      val(i,j,3) = q(1)*q(3)
      val(i,j,4) = q(4)/.4 + .5d0*q(1)*(q(2)*q(2)+q(3)*q(3))

      go to 10

99    return
      end


text:
!! The IMPLICIT NONE statement causes the compiler to regard S in the function
!! declaration as undeclared even though its type is declared in the body of
!! the function.  This is eerily similar to SR98771 which I reported (and NAG
!! fixed for 6.1) back in Aug '17.
!!
!! $ nagfor -c nag-20180203e.f90
!! NAG Fortran Compiler Release 6.2(Chiyoda) Build 6205
!! Error: nag-20180203e.f90, line 16: Implicit type for S
!!        detected at (@S
!! Error: nag-20180203e.f90, line 17: Symbol S has already been implicitly typed
!!        detected at S@<end-of-statement>

module mod
  implicit none
contains
  character(len(s)) function foo(s)
    character(*), intent(in) :: s
    foo = s
  end function
end module

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
Then, put whatever you want down here!  Note that http://local.google.com is sometimes a good and easy place to find times/phone. Better yet, look them up in the phone book, and ask them.

Please leave a link back to the media page in your writeup.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module Mod_lmn_ComputeAdvectionVelocity
   use typre
   use Mod_lmn_BaseElmope
   implicit none
   private
   public SetPointersAdvectionVelocity

   !External Procedures
   procedure() :: NULLSUB

   integer(ip), allocatable :: kfl_IsSet

contains

   subroutine SetPointersAdvectionVelocity(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
            kfl_IsSet = 1

            ProcPointer%ComputeAdvectionVelocity => ComputeAdvectionVelocity

            !Advection Velocity
            if (a%kfl_advec == 0) then
               ProcPointer%ComputeAdvectionVelocity => NULLSUB
            endif

            !NonLinear Subscales
            if (a%kfl_nolsg == 1) then
               call ConcatenateProcedures(ProcPointer%ComputeAdvectionVelocity,NonLinearSGSAdvectionVelocity)
            endif
         endif

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

      end select
   end subroutine

   !-----------------------------------------------------------------------
   !Actual Computations
   !AdvectionVelocity
   subroutine ComputeAdvectionVelocity
      implicit none

      gpadv = gpvel(:,1)
   end subroutine

   !Non-linear subscales
   subroutine NonLinearSGSAdvectionVelocity
      implicit none

      gpadv = gpadv + a%vesgs(ielem)%a(:,1,e%igaus)
   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


text:
program ej2

  !Programa que resuelve ecuaciones de segundo grado. Soporta raices complejas
  complex x_1, x_2, aux
  real a, b, c, disc

  print *, "Sea la ecuación con la forma ax² + bx + c = 0"
  print *, "Porfavor ingresa los valores de los coeficientes a, b y c"
  read *, a, b, c

  disc = b**2 - (4*a*c)

  if (disc < 0) then
    disc = - sqrt(abs(disc))
    aux = cmplx(0, disc)

  end if

  x_1 = (-b + aux)/(2*a)
  x_2 = (-b - aux)/(2*a)

  print *, "Las raices son: x = {", x_1, ", ", x_2, "}"

end program ej2

max_stars_repo_name: kotoromo/Fortran-Basico

max_issues_repo_name: kotoromo/Fortran-Basico

max_forks_repo_name: kotoromo/Fortran-Basico