text:
c     ****************************************************************
c     *                                                              *
c     *                      subroutine ouhstr                       *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 3/22/21 rhd                *
c     *                                                              *
c     *     drive output of element strains/stresses to printed      *
c     *     hardcopy and/or to packets files                         *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine ouhstr( do_stress, wide, eform, prec,
     &                   noheader, out_packet_now, element_list,
     &                   num_list_entries )
c
      use global_data, only : ltmstp, iprops, lprops, elblks
      use main_data, only: packet_file_no, cohesive_ele_types,
     &                     elems_to_blocks, bar_types, link_types
c
      implicit none
c
      integer, intent(in) :: num_list_entries
      integer, intent(in) :: element_list(num_list_entries)
      logical, intent(in) :: wide, eform, prec, do_stress, noheader,
     &                       out_packet_now
c
c             local declarations
c
      integer :: local_count, index, elem, elem_type, num_enodes,
     &           num_int_points, output_loc, additional,
     &           cohes_pkt_offset, pkt_type, lnum, pgnum, lbltyp,
     &           span, blk, rel_elem, int_order, mat_type, totdof,
     &           num_short_stress, num_short_strain, cohesive_type,
     &           felem, num_enode_dof
      character(len=8) :: strlbl(30), hedtyp*30
      real :: dumr
      double precision :: dumd
      character :: dums
      logical :: bbar_flg, geo_non_flg, long_out_flg, solid_elem,
     &           nodpts_flg, center_output, cohesive_elem,
     &           at_intpts, at_enodes, at_center, bar_elem, link_elem
c
c             for packet output, loop over the user specified
c             element list to determine how many data records
c             will be in the packet. the packet contains
c             as many elements an in the list with data
c             at integration points, element nodes, or center.
c
      if( out_packet_now ) then
c
        local_count = 0
        do index = 1, num_list_entries
          elem           = element_list(index)
          elem_type      = iprops(1,elem)
          num_enodes     = iprops(2,elem)
          num_int_points = iprops(6,elem)
          output_loc     = iprops(12,elem)
          cohesive_elem  = cohesive_ele_types(elem_type)
          at_intpts      = output_loc .eq. 1
          at_enodes      = output_loc .eq. 2
          at_center      = output_loc .eq. 3
          additional     = 0
          cohes_pkt_offset = 0
          if( cohesive_elem ) then
             additional = num_int_points
             if( at_center )  additional = 1
             cohes_pkt_offset = 16
          else   ! solid element
             if( at_enodes  ) additional = num_enodes
             if( at_center ) additional = 1
             if( at_intpts ) additional = num_int_points
          end if
          local_count = local_count + additional
         end do
         pkt_type = 15 + cohes_pkt_offset
         if( .not. do_stress ) pkt_type = pkt_type + 1
         write(packet_file_no) pkt_type, local_count, ltmstp, 0
c
      end if
c
c             initialize parameters controlling hardcopy output.
c
      lnum               = 56  ! forces printing of headers
      pgnum              = 0
      lbltyp             = 0   ! tracks last printed labels
      strlbl(1:30)       = ' '
c
c             pass over the use list of elements for output.
c             call lower routines to compute secondary (derivable)
c             output values, then output values to hardcopy and/or
c             packet.
c
      do index = 1, num_list_entries
        elem              = element_list(index)
        span              = 1
        blk               = elems_to_blocks(elem,1)
        rel_elem          = elems_to_blocks(elem,2)
        felem             = elblks(1,blk)
        elem_type         = iprops(1,elem)
        int_order         = iprops(5,elem)
        mat_type          = iprops(25,elem)
        num_enodes        = iprops(2,elem)
        num_enode_dof     = iprops(4,elem)
        totdof            = num_enodes * num_enode_dof
        geo_non_flg       = lprops(18,elem)
        bbar_flg          = lprops(19,elem)
        num_int_points    = iprops(6,elem)
        long_out_flg      = lprops(16,elem)
        output_loc        = iprops(12,elem)
        nodpts_flg        = output_loc .eq. 2
        center_output     = output_loc .eq. 3
        if( out_packet_now ) long_out_flg = .true.
        num_short_stress  = 11
        num_short_strain  = 7
        cohesive_type     = iprops(27,elem)
        cohesive_elem     = cohesive_ele_types(elem_type)
        bar_elem          = bar_types(elem_type)
        link_elem         = link_types(elem_type)
c
c             duplicate necessary element block data.
c
        call oudups( span, elem, num_int_points, geo_non_flg,
     &               do_stress, cohesive_elem )
c
c             compute the element block stress/strain data to be
c             output. set to get new header labels if needed.
c             solids and interface-cohesive treated separately.
c
c             lbltyp = 1 for solids.
c                    = 2 for bars
c                    = 3 for links
c                    = 5 +  cohesive_type for cohesive elements
c
        if( cohesive_elem ) then
           if( lbltyp .ne. 5 + cohesive_type  ) then
            lbltyp = 5 + cohesive_type
            lnum = 56 ! forces new page & label printing
          end if
          call ouhrks_cohesive( elem, blk, felem, elem_type,
     &       int_order, num_int_points, num_enodes, center_output,
     &       do_stress, mat_type, cohesive_type, wide, eform,
     &       prec, lnum, pgnum, lbltyp, strlbl, hedtyp,  noheader,
     &       out_packet_now, geo_non_flg )
          cycle
        end if
c
        if( bar_elem ) then
          if( lbltyp .ne. 2  ) then
            lbltyp = 2
            lnum = 56 ! forces new page & label printing
          end if
          call ouhrks_bar( elem, blk, felem, elem_type, num_enodes,
     &        do_stress, mat_type, wide, eform, prec, lnum, pgnum,
     &        lbltyp, strlbl, hedtyp,  noheader, out_packet_now,
     &        geo_non_flg )
           cycle
        end if
c
        if( link_elem ) then
          if( lbltyp .ne. 3  ) then
            lbltyp = 3
            lnum = 56 ! forces new page & label printing
          end if
          call ouhrks_link( elem, blk, felem, elem_type, num_enodes,
     &        do_stress, mat_type, wide, eform, prec, lnum, pgnum,
     &        lbltyp, strlbl, hedtyp,  noheader, out_packet_now,
     &        geo_non_flg )
           cycle
        end if
c
c             process solid element
c
        if( lbltyp .ne. 1 ) then
          lbltyp = 1
          lnum = 56 ! forces new page & label printing
        end if
c
        call ouhrks( span, blk, felem, elem_type, int_order,
     &                num_int_points, num_enodes,
     &                geo_non_flg, long_out_flg,
     &                nodpts_flg, do_stress, mat_type,
     &                center_output, num_short_stress,
     &                num_short_strain )
c
       call ouhprn(  elem, elem_type, num_int_points,
     &                num_enodes, long_out_flg,
     &                nodpts_flg, do_stress, wide, eform, prec,
     &                lnum, pgnum, lbltyp, strlbl, hedtyp,
     &                num_short_stress, num_short_strain,
     &                center_output, geo_non_flg, noheader,
     &                out_packet_now )
c
      end do
c
      return
      end

c     ****************************************************************
c     *                                                              *
c     *                      subroutine ou_gastr                     *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 03/22/21 rhd               *
c     *                                                              *
c     ****************************************************************
c
c
c               moved here to stop compilers from complaining about
c               history_global is history....(blk)%ptr(location).
c               compilers don't like this type of old-style Fortran
c               array passing. it thinks ptr might be an array section.
c
c               by putting routine here and with separate compiles, the
c               cont's see what we're doing (it is legal!)
c
c
      subroutine ou_gastr( history_local, history_global, ngp,
     &                     mxhist, mxngp, hist_size, span, mxvl )
      implicit none
c
c               parameter declarations
c
      integer, intent(in) :: ngp, mxhist, mxngp, hist_size, mxvl, span
      double precision, intent(in) :: history_global(hist_size,ngp,span)
      double precision, intent(out) :: history_local(mxvl,mxhist,mxngp)
c
c               local declarations
c
      integer i, j, k
c
      if( ngp .ne. 8 ) then
        do k = 1, ngp
         do  j = 1, hist_size
!DIR$ VECTOR ALIGNED
            do  i = 1, span
               history_local(i,j,k) = history_global(j,k,i)
            end do
         end do
        end do
        return
      end if
c
c                number of gauss points = 8, unroll.
c
      do  j = 1, hist_size
!DIR$ VECTOR ALIGNED
        do  i = 1, span
            history_local(i,j,1) = history_global(j,1,i)
            history_local(i,j,2) = history_global(j,2,i)
            history_local(i,j,3) = history_global(j,3,i)
            history_local(i,j,4) = history_global(j,4,i)
            history_local(i,j,5) = history_global(j,5,i)
            history_local(i,j,6) = history_global(j,6,i)
            history_local(i,j,7) = history_global(j,7,i)
            history_local(i,j,8) = history_global(j,8,i)
        end do
      end do
c
      return
      end


text: Sheep can be found at the Animal Science sheep barn west of main campus.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c
c     integer function length(strg)
c
c     logical function fileOK(name,l)
c
c     integer function next_key(list,nlist,iof)
c
c     subroutine strip(strg)
c     logical stripOK(strg)
c
c     subroutine replace(strg,dummy,txt,l,n)
c
c     subroutine presskey(key)
c
c     subroutine write_s(A)
c     subroutine write_u(A)
c     subroutine write_4(X)
c

c==============================================================================

      integer function length(strg)
c------------------------------------------------------
c     return length of strg
c
      implicit none
      integer j
      character strg*(*)
      j = len(strg)
c      do while ((j.gt.0).and.(.not.(
c     *  ((strg(j:j).ge.'a').and.(strg(j:j).le.'z')).or.
c     *  ((strg(j:j).ge.'A').and.(strg(j:j).le.'Z')).or.
c     *  ((strg(j:j).ge.'1').and.(strg(j:j).le.'9')).or.
c     *   (strg(j:j).eq.'0').or.
c     *   (strg(j:j).eq.'.').or.
c     *   (strg(j:j).eq.'_'))))
      do while (j.gt.0.and.(strg(j:j).eq.' '.or.strg(j:j).eq.''))
        j = j - 1
      enddo
      length = j
      return
      end
c==============================================================================

      logical function fileOK(name,mxl)
c---------------------------------------------------------
c     check whether name is an admissible file name
c
      implicit none
      integer mxl, length, i, ll
      character name*(*)

      ll = length(name)
      if (ll.gt.mxl.or.ll.lt.1) then
        fileOK = .false.
      else
        fileOK = .true.
        do i=1, ll
          if (.not.(
     *      ((name(i:i).ge.'a').and.(name(i:i).le.'z')).or.
     *      ((name(i:i).ge.'A').and.(name(i:i).le.'Z')).or.
     *      ((name(i:i).ge.'1').and.(name(i:i).le.'9')).or.
     *       (name(i:i).eq.'0').or.
     *       (name(i:i).eq.'.').or.
     *       (name(i:i).eq.'_')))
     *      fileOK = .false.
        enddo
      endif

      return
      end
c==============================================================================

      integer function next_key(list,nlist,iof)
c----------------------------------------------------
c     read next key from file and return id number
c
      implicit none
      integer   nlist, iof, ios, length, i, l
      character strg*4, list(nlist)*4
      logical   pcomp

      i = nlist + 1
      read(iof,'(A)',iostat=ios) strg
      do while ((ios.eq.0).and.(i.gt.nlist))
        i = 1
        l = length(list(i))
        do while (.not.(pcomp(strg,list(i),l)).and.(i.le.nlist))
          i = i + 1
          l = length(list(i))
        enddo
        if (i.gt.nlist) read(iof,'(A)',iostat=ios) strg
      enddo

      if (ios.eq.0) then
        next_key = i
      else
        next_key = 0
      endif

      return
      end
c==============================================================================

      subroutine strip(strg)
c----------------------------------------------------
c     remove comments and empty parts at beginning and end from strg
c
      implicit none
      integer   l, ll, length
      character strg*(*)

c...  remove comments
      ll = length(strg)
      l  = 1
      do while ((strg(l:l).ne.'!').and.(strg(l:l).ne.'#').and.(l.le.ll))
        l = l + 1
      enddo
      strg = strg(1:l-1)
      ll = l - 1

c...  remove empty space at the beginning
      l = 1
      do while (strg(l:l).eq.' '.and.l.le.ll)
        l = l + 1
      enddo
      if (l.gt.ll) then
        strg = ''
      else
        strg = strg(l:ll)
      endif

c...  remove empty space at the end
      strg = strg(1:length(strg))

      return

      end
c==============================================================================

      logical function stripOK(strg)
c----------------------------------------------------
c     remove comments from strg
c     return .false. if clean strg is empty
c
      implicit none
      integer   length
      character strg*(*)

      stripOK = .true.

      call strip(strg)

      if (length(strg).eq.0) stripOK = .false.

      return

      end
c==============================================================================

      logical function mypcomp(a,b)
c----------------------------------------------------
c     compare two character strings a and b
c
      implicit  none
      integer   la, lb, length
      logical   pcomp
      character a*(*), b*(*)

      mypcomp = .false.

      la = length(a)
      lb = length(b)

c      write(*,*) la, lb
c      write(*,*) a(1:la)
c      write(*,*) b(1:lb)

      if (la.ne.lb) return

      mypcomp = pcomp(a,b,la)

      end
c==============================================================================

      logical function pcomp(a,b,n)
c----------------------------------------------------
c     compare two character strings a and b
c
      implicit  none
      integer   n, inc, i, ia, ib, iau, izu

      character a*(*),b*(*)

      pcomp = .false.

      iau = ichar('A')
      izu = ichar('Z')

      inc = iau - ichar('a')

      do i=1, n

        ia = ichar(a(i:i))
        ib = ichar(b(i:i))

        if (ia.ge.iau.and.ia.le.izu) ia = ia - inc
        if (ib.ge.iau.and.ib.le.izu) ib = ib - inc

        if (ia.ne.ib) return
      end do

      pcomp = .true.

      end
c==============================================================================

      subroutine replace(strg,dummy,txt,l,n)
c------------------------------------------------------
c     replace dummy(1..n) in strg by txt(1..n)
c
      implicit none
      integer   l, n, length,
     *          ls, lt, ld, i, il, j, k
      character strg*(*), dummy*(*), txt*(*)

      ls = length(strg)

      do i=1, n

        il = (i - 1) * l

        ld = length(dummy(il+1:il+l))
        lt = length(txt(il+1:il+l))

        j = 1

        do while (j.le.ls-ld)

          do while ((strg(j:j-1+ld).ne.dummy(il+1:il+l))
     *      .and.(j.le.ls-ld+1))
            j = j + 1
          enddo

          if (j.le.ls-ld+1) then

            ls = ls - ld + lt

            if (lt.gt.ld) then
              do k=ls, j+lt, -1
                strg(k:k) = strg(k-lt+ld:k-lt+ld)
              enddo
            elseif (ld.gt.lt) then
              do k=j+lt, ls
                strg(k:k) = strg(k+ld-lt:k+ld-lt)
              enddo
            endif

            do k=j, j+lt-1
              strg(k:k) = txt(il+k-j+1:il+k-j+1)
            enddo

            if (lt-ld.lt.0) then
              do k=ls+1, ls+ld-lt
                strg(k:k) = ''
              enddo
            endif

          endif

        enddo

      enddo

      return
      end
c==============================================================================

      subroutine presskey(key)
c------------------------------------------------------
c     wait until 'key' has been pressed
c
c     key = 1  -> ENTER
c           2  -> Esc
c           3  -> x - ENTER
c           4  -> q - ENTER
c
      implicit none

      integer   key

      character ch

      save

      read(*,'(a)') ch

      return
      end
c==============================================================================

      subroutine write_s(A)
      implicit none
      double precision A(6)
      write(*,'(/3(1x,f6.3))') A(1), A(4), A(6)
      write(*,'( 3(1x,f6.3))') A(4), A(2), A(5)
      write(*,'( 3(1x,f6.3))') A(6), A(5), A(3)
      return
      end

      subroutine write_u(A)
      implicit none
      double precision A(3,3)
      write(*,'(/3(1x,f8.3))') A(1,1), A(1,2), A(1,3)
      write(*,'( 3(1x,f8.3))') A(2,1), A(2,2), A(2,3)
      write(*,'( 3(1x,f8.3))') A(3,1), A(3,2), A(3,3)
      return
      end

      subroutine write_4(X)
      implicit none
      integer          i, j
      double precision X(6,6)
      write(*,*)
      do i=1, 6
        write(*,'(6(1x,f8.3))') (X(i,j),j=1,6)
      enddo
      return
      end

ext: for

max_stars_repo_name: chennachaos/stabfem

max_issues_repo_name: chennachaos/stabfem

max_forks_repo_name: chennachaos/mpap


text:
Describe Users/RobynLee here.

See also Users/ManagerSilverstone


module ECRad_python

contains

subroutine pre_initialize_ECRad(ecrad_verbose, dstf, ray_tracing, ecrad_Bt_ripple, &
                                rhopol_max_spline_knot, ecrad_weak_rel, &
                                ecrad_ratio_for_third_harmonic, ecrad_tau_ignore, &
                                ecrad_modes, reflec_X_mode, reflec_O_mode, ece_1O_flag, &
                                ecrad_max_points_svec, N_pts_BPD, &
                                ecrad_O2X_mode_conversion, &
                                rhopol_scal_te, rhopol_scal_ne, &
                                ecrad_ds_large, ecrad_ds_small, ecrad_R_shift, &
                                ecrad_z_shift, &
                                ecrad_N_ray, ecrad_N_freq, log_flag, N_vessel, vessel_R, vessel_z, &
                                f, df, R, phi, z, tor, pol, dist_foc, width, pol_coeff)
! Everything that is absolutely static in time is done over here
use mod_ECRad,      only: pre_initialize_ECRad_f2py
implicit none
character(*), intent(in)        :: dstf
real(kind=8), intent(in)      :: rhopol_max_spline_knot, ecrad_ratio_for_third_harmonic, ecrad_tau_ignore, &
                                              reflec_X_mode, reflec_O_mode, ecrad_O2X_mode_conversion, &
                                              rhopol_scal_te, rhopol_scal_ne, &
                                              ecrad_ds_large, ecrad_ds_small, ecrad_R_shift, ecrad_z_shift
integer, intent(in)   :: ecrad_modes, ecrad_max_points_svec, N_pts_BPD, ecrad_N_ray, &
                         ecrad_N_freq,ece_1O_flag
logical, intent(in)           :: ecrad_verbose, ecrad_Bt_ripple, ray_tracing, ecrad_weak_rel, log_flag
integer, intent(in)    :: N_vessel
real(kind=8), dimension(:), intent(in) :: vessel_R, vessel_z ! vessel contour
real(kind=8), dimension(:), intent(in), optional :: f, df, R, phi, z, tor, pol, dist_foc, width, pol_coeff
call pre_initialize_ECRad_f2py(ecrad_verbose, dstf, ray_tracing, ecrad_Bt_ripple, &
                               rhopol_max_spline_knot, ecrad_weak_rel, &
                               ecrad_ratio_for_third_harmonic, ecrad_tau_ignore, &
                               ecrad_modes, reflec_X_mode, reflec_O_mode, ece_1O_flag, &
                               ecrad_max_points_svec, N_pts_BPD ,& ! (modes = 1 -> pure X-mode, 2 -> pure O-mode, 3 both modes and filter
                               ecrad_O2X_mode_conversion, & ! mode conversion ratio from O-X due to wall reflections
                               ! Scaling of rhop axis for shifting on ne or Te
                               ! Every rhop value obtained in ray tracing will be multiplied by the corresponding scaling value when evaluating Te/ne
                               rhopol_scal_te, rhopol_scal_ne, &
                               ecrad_ds_large, ecrad_ds_small, ecrad_R_shift, &    ! Allows shifting the equilbrium - moves entire flux matrix
                               ecrad_z_shift, &    ! Allows shifting the equilbrium - moves entire flux matrix
                               ecrad_N_ray, ecrad_N_freq, log_flag, N_vessel, vessel_R, vessel_z, &
                               f, df, R, phi, z, tor, pol, dist_foc, width, pol_coeff)
end subroutine pre_initialize_ECRad

subroutine reset_ECRad()
use mod_ECRad,      only: clean_up_ECRad
implicit none
  call clean_up_ECRad()
end subroutine reset_ECRad

subroutine initialize_ECRad(N_ch, N_Te_spline_knots, N_ne_spline_knots, &
                            R, z, rhop, Br, Bt, Bz, R_ax, z_ax, rhopol_out)
! Hence, to keep the structure similiar all initizalization is performed here
! Initializations that depend on time are done here
use mod_ECRad,        only: initialize_ECRad_f2py
implicit none
integer, intent(in)                        		     :: N_ch
integer, intent(in)                                :: N_Te_spline_knots, N_ne_spline_knots
real(kind=8), intent(in), optional                 :: R_ax, z_ax
real(kind=8), dimension(:), intent(in)             :: R, z
real(kind=8), dimension(:,:), intent(in)           :: rhop, Br, Bt, Bz
real(kind=8), dimension(N_ch), intent(out)         :: rhopol_out
call initialize_ECRad_f2py(N_Te_spline_knots, N_ne_spline_knots, &
                           R, z, rhop, Br, Bt, Bz, R_ax, z_ax)
end subroutine initialize_ECRad

subroutine initialize_ECRad_2D_profs(N_ch, N_Te_spline_knots, N_ne_spline_knots, &
                            R, z, T_e, n_e, rhop, Br, Bt, Bz, R_ax, z_ax)
! Hence, to keep the structure similiar all initizalization is performed here
! Initializations that depend on time are done here
use mod_ECRad,        only: initialize_ECRad_f2py
implicit none
integer, intent(in)                                :: N_ch
integer, intent(in)                                :: N_Te_spline_knots, N_ne_spline_knots
real(kind=8), intent(in), optional                 :: R_ax, z_ax
real(kind=8), dimension(:), intent(in)             :: R, z
real(kind=8), dimension(:,:), intent(in)           :: rhop, Br, Bt, Bz, T_e, n_e
call initialize_ECRad_f2py(N_Te_spline_knots, N_ne_spline_knots, &
                           R, z, rhop, Br, Bt, Bz, R_ax, z_ax, T_e, n_e)
end subroutine initialize_ECRad_2D_profs

subroutine initialize_ECRad_3D(N_ch, N_Te_spline_knots, N_ne_spline_knots, &
                               equilibrium_file, equilibrium_type, use_mesh, &
                               use_symmetry, B_ref, s_plus, s_max, &
                               interpolation_acc, fourier_coeff_trunc, &
                               h_mesh, delta_phi_mesh, vessel_filename, &
                               rhopol_out)
! Hence, to keep the structure similiar all initizalization is performed here
! Initializations that depend on time are done here
use mod_ECRad,        only: initialize_ECRad_3D_f2py
implicit none
integer, intent(in)      :: N_ch
integer(4), intent(in) :: N_Te_spline_knots, N_ne_spline_knots
CHARACTER(*), intent(in) :: equilibrium_file
CHARACTER(*), intent(in) :: equilibrium_type
logical, intent(in)      :: use_mesh, use_symmetry
real(8), intent(in)  :: B_ref, s_plus, s_max, h_mesh, delta_phi_mesh, &
                            interpolation_acc, fourier_coeff_trunc
CHARACTER(*), intent(in) :: vessel_filename
real(8), dimension(N_ch), intent(out)  :: rhopol_out
call initialize_ECRad_3D_f2py(N_Te_spline_knots, N_ne_spline_knots, &
                              equilibrium_file, equilibrium_type, use_mesh, &
                              use_symmetry, B_ref, s_plus, s_max, &
                              interpolation_acc, fourier_coeff_trunc, &
                              h_mesh, delta_phi_mesh, vessel_filename, &
                              rhopol_out)
end subroutine initialize_ECRad_3D

subroutine set_ECRad_FFP_dist(rho, u, pitch, f)
  use mod_ECRad,        only: set_ECRad_FFP_dist_f2py
  implicit none
  real(kind=8), dimension(:), intent(in)             :: rho, u, pitch
  real(kind=8), dimension(:,:,:), intent(in)         :: f
  call set_ECRad_FFP_dist_f2py(rho, u, pitch, f)
end subroutine set_ECRad_FFP_dist

subroutine set_ECRad_GENE_dist(rho, vpar, mu, f_0, g)
  use mod_ECRad,        only: set_ECRad_GENE_dist_f2py
  implicit none
  real(kind=8), dimension(:), intent(in)             :: rho, vpar, mu
  real(kind=8), dimension(:,:), intent(in)           :: f_0
  real(kind=8), dimension(:,:,:), intent(in)         :: g
  call set_ECRad_GENE_dist_f2py(rho, vpar, mu, f_0, g)
end subroutine set_ECRad_GENE_dist

subroutine make_rays_ECRad(N_ch, rhop_knots_ne, n_e, rhop_knots_Te, T_e, rhop_res)
use mod_ECRad,        only: make_rays_ECRad_f2py
implicit none
integer, intent(in)                        :: N_ch
real(kind=8), dimension(:), intent(in) :: rhop_knots_ne, n_e, rhop_knots_Te, T_e
real(kind=8), dimension(N_ch),  intent(out) :: rhop_res
integer                          :: idiag, ich
call make_rays_ECRad_f2py(rhop_knots_ne=rhop_knots_ne, n_e=n_e, rhop_knots_Te=rhop_knots_Te, &
					                T_e=T_e,  rhop_res=rhop_res)
end subroutine make_rays_ECRad

subroutine make_rays_ECRad_2D(N_ch, rhop_res)
use mod_ECRad,        only: make_rays_ECRad_F2py
implicit none
integer, intent(in)                        :: N_ch
real(kind=8), dimension(N_ch),  intent(out) :: rhop_res
integer                          :: idiag, ich
! For 2D profiles Te and ne are already in place
call make_rays_ECRad_f2py(rhop_res=rhop_res)
end subroutine make_rays_ECRad_2D

subroutine make_rays_ECRad_spline(N_ch, rhop_knots_ne, n_e, n_e_dx2, rhop_knots_Te, T_e, T_e_dx2, &
                                  rhop_res)
use mod_ECRad,        only: make_rays_ECRad_f2py
implicit none
integer, intent(in)                        :: N_ch
real(kind=8), dimension(:), intent(in) :: rhop_knots_ne, n_e, n_e_dx2, rhop_knots_Te, T_e, T_e_dx2
real(kind=8), dimension(N_ch),  intent(out) :: rhop_res
integer                          :: idiag, ich
call make_rays_ECRad_f2py(rhop_knots_ne, n_e, n_e_dx2, rhop_knots_Te, T_e, T_e_dx2, rhop_res)
end subroutine make_rays_ECRad_spline

subroutine make_dat_model_ECRad(N_ch, rhop_knots_ne, n_e, rhop_knots_Te, T_e, &
                                        ne_rhop_scal, reflec_X_new, &
                                        reflec_O_new, ece_fm_flag_ch, rp_min, &
                                        dat_model_ece, tau, set_grid_dynamic, verbose)
use mod_ECRad,        only: make_dat_model_ece_ECRad_f2py
implicit none
integer, intent(in)             :: N_ch
real(kind=8), dimension(:), intent(in)  :: rhop_knots_ne, n_e, rhop_knots_Te, T_e
real(kind=8),               intent(in)  :: ne_rhop_scal, reflec_X_new, reflec_O_new, rp_min
logical,     dimension(:), intent(in)  :: ece_fm_flag_ch
real(kind=8), intent(out) :: dat_model_ece(N_ch)
real(kind=8), intent(out)  :: tau(N_ch)
logical,      intent(in)     :: verbose
logical, intent(in)          :: set_grid_dynamic
integer                         :: ich
call make_dat_model_ece_ECRad_f2py(rhop_knots_ne=rhop_knots_ne, n_e=n_e, rhop_knots_Te=rhop_knots_Te, &
								                   T_e=T_e, ne_rhop_scal=ne_rhop_scal, reflec_X_new=reflec_X_new, & ! in
                                   reflec_O_new=reflec_O_new, ece_fm_flag_ch=ece_fm_flag_ch, rp_min=rp_min, &
                                   dat_model_ece=dat_model_ece, tau=tau, set_grid_dynamic=set_grid_dynamic, &
                                   verbose=verbose)
end subroutine make_dat_model_ECRad

subroutine make_TRad_direct()
use mod_ECRad,        only: make_ece_rad_temp, set_for_single_eval
implicit none
integer :: idiag
idiag = 1
! For single call. For usage with fixed Te/ne profiles.
call set_for_single_eval()
call make_ece_rad_temp()
end subroutine make_TRad_direct

subroutine make_dat_model_ECRad_spline(N_ch, rhop_knots_ne, n_e, n_e_dx2, rhop_knots_Te, T_e, T_e_dx2, &
                                        ne_rhop_scal, reflec_X_new, &
                                        reflec_O_new, ece_fm_flag_ch, rp_min, &
                                        dat_model_ece, tau, set_grid_dynamic, verbose)
use mod_ECRad,        only: make_dat_model_ece_ECRad_f2py
implicit none
integer, intent(in)             :: N_ch
real(kind=8), dimension(:), intent(in)  :: rhop_knots_ne, n_e, rhop_knots_Te, T_e
real(kind=8),               intent(in)  :: ne_rhop_scal, reflec_X_new, reflec_O_new, rp_min
real(kind=8), dimension(:), intent(in)  :: n_e_dx2, T_e_dx2
logical,     dimension(:), intent(in)   :: ece_fm_flag_ch
real(kind=8), intent(out) :: dat_model_ece(N_ch)
real(kind=8), intent(out)  :: tau(N_ch)
logical,      intent(in)     :: verbose
logical, intent(in)          :: set_grid_dynamic
integer                         :: ich
  call make_dat_model_ece_ECRad_f2py(rhop_knots_ne, n_e, n_e_dx2, rhop_knots_Te, T_e, T_e_dx2, &
                                ne_rhop_scal, reflec_X_new, & ! in
                                reflec_O_new, ece_fm_flag_ch, rp_min, &
                                dat_model_ece, tau, set_grid_dynamic, verbose)
end subroutine make_dat_model_ECRad_spline

subroutine update_Te_ne_ECRad(rhop_knots_ne, n_e, n_e_dx2, rhop_knots_Te, T_e, T_e_dx2)
use mod_ECRad,               only: update_Te_ne
  implicit none
  real(kind=8), dimension(:), intent(in)   :: rhop_knots_ne, n_e, rhop_knots_Te, T_e
  real(kind=8), dimension(:), intent(in), optional   :: n_e_dx2, T_e_dx2
  call update_Te_ne(rhop_knots_ne, n_e, n_e_dx2, rhop_knots_Te, T_e, T_e_dx2)
end subroutine update_Te_ne_ECRad

subroutine get_Trad_resonances_basic(imode, N_ch, Trad, tau, s_res, R_res, z_res, rho_res)
use mod_ECRad,        only: get_Trad_resonances_basic_f2py
implicit None
integer, intent(in)   :: imode, N_ch
! <= 0 -> average over modes
! > 0 corresponding to imode in ECRad
real(kind=8), dimension(N_ch), intent(out) :: Trad, tau, s_res, R_res, z_res, rho_res
  call get_Trad_resonances_basic_f2py(imode, N_ch, Trad, tau, s_res, R_res, z_res, rho_res)
end subroutine get_Trad_resonances_basic

subroutine get_Trad_resonances_extra_output(imode, N_ch, Trad_secondary, tau_secondary, &
                                            rel_s_res, rel_rho_res, &
                                            rel_R_res, rel_z_res, &
                                            rel_s_res_secondary, rel_rho_res_secondary, &
                                            rel_R_res_secondary, rel_z_res_secondary)
use mod_ECRad,        only: get_Trad_resonances_extra_output_f2py
implicit None
integer, intent(in)   :: imode, N_ch
! <= 0 -> average over modes
! > 0 corresponding to imode in ECRad
real(kind=8), dimension(N_ch), intent(out) :: Trad_secondary, tau_secondary, &
                                              rel_s_res, rel_rho_res, &
                                              rel_R_res, rel_z_res, &
                                              rel_s_res_secondary, rel_rho_res_secondary, &
                                              rel_R_res_secondary, rel_z_res_secondary

  call get_Trad_resonances_extra_output_f2py(imode, Trad_secondary, tau_secondary, &
                                             rel_s_res, rel_rho_res, &
                                             rel_R_res, rel_z_res, &
                                             rel_s_res_secondary, rel_rho_res_secondary, &
                                             rel_R_res_secondary, rel_z_res_secondary)
end subroutine get_Trad_resonances_extra_output

subroutine get_BPD(ich, imode, N_BPD, rho, BPD, BPD_second)
use mod_ECRad,        only: get_BPD_f2py
implicit None
integer, intent(in)   :: ich, imode, N_BPD
real(kind=8), dimension(N_BPD), intent(out) :: rho, BPD, BPD_second
  call get_BPD_f2py(ich, imode, rho, BPD, BPD_second)
end subroutine get_BPD

subroutine get_ray_length(ich, imode, ir, N_LOS)
use mod_ECRad,        only: get_ray_length_f2py
implicit None
integer, intent(in)   :: ich, imode, ir
integer, intent(out)  :: N_LOS
integer               :: idiag
  call get_ray_length_f2py(ich, imode, ir, N_LOS)
end subroutine get_ray_length

subroutine get_ray_data(ich, imode, ir, N_LOS, s, x, y, z, Nx, Ny, Nz, &
                        Bx, By, Bz, rho, T_e, n_e, theta, N_cold, H, v_g_perp, &
                        Trad, Trad_secondary, em, em_secondary, ab, ab_secondary, T, &
                        T_secondary, BPD, BPD_secondary)
use mod_ECRad,        only: get_ray_data_f2py
implicit None
integer, intent(in)   :: ich, imode, ir, N_LOS
real(kind=8), dimension(N_LOS), intent(out) :: s, x, y, z, Nx, Ny, Nz, Bx, By, Bz, rho, &
                                               T_e, n_e, theta, N_cold, H, v_g_perp, &
                                               Trad, Trad_secondary, em, em_secondary, &
                                               ab, ab_secondary, T, T_secondary, BPD, BPD_secondary
  call get_ray_data_f2py(ich, imode, ir, s, x, y, z, Nx, Ny, Nz, &
                        Bx, By, Bz, rho, T_e, n_e, theta, N_cold, H, v_g_perp, &
                        Trad, Trad_secondary, em, em_secondary, ab, ab_secondary, T, &
                        T_secondary, BPD, BPD_secondary)
end subroutine get_ray_data

subroutine get_mode_weights(N_ch, imode, pol_coeff, pol_coeff_secondary)
  use mod_ECRad,        only:  get_mode_weights_f2py
  implicit none
    integer, intent(in) :: N_ch, imode
    real(kind=8), dimension(N_ch), intent(out) :: pol_coeff, pol_coeff_secondary
    call get_mode_weights_f2py(N_ch, imode, pol_coeff, pol_coeff_secondary)
end subroutine get_mode_weights

subroutine get_weights(N_ray, N_freq, ich,  ray_weights, freq_weights)
use mod_ECRad,        only: get_weights_f2py
implicit none
  integer, intent(in)   :: N_ray, N_freq, ich
  real(kind=8), dimension(N_ray), intent(out) :: ray_weights
  real(kind=8), dimension(N_freq), intent(out) :: freq_weights
  call get_weights_f2py(ich, ray_weights, freq_weights)
end subroutine get_weights

!subroutine make_BPD_w_res_ch_ECRad(pnts_BPD, idiag, ich, rhop_knots_ne, n_e, n_e_dx2, rhop_knots_Te, T_e, T_e_dx2, &
!                                   ne_rhop_scal, reflec_X_new, &
!                                   reflec_O_new, rp_min, &
!                                   rhop, BPD, rhop_res_warm)
!use mod_ECRad,        only: make_BPD_w_res_ch
!implicit none
!integer, intent(in)             :: pnts_BPD
!integer, intent(in)             :: idiag, ich
!real(kind=8), dimension(:), intent(in)  :: rhop_knots_ne, n_e, rhop_knots_Te, T_e, n_e_dx2, T_e_dx2
!real(kind=8),               intent(in)  :: ne_rhop_scal, reflec_X_new, reflec_O_new, rp_min
!real(kind=8),  intent(out) :: rhop(pnts_BPD), BPD(pnts_BPD)
!real(kind=8), intent(out)               :: rhop_res_warm
!real(kind=8), dimension(:), allocatable :: rhop_ECRad_out, BPD_ECRad_out
!call make_BPD_w_res_ch(idiag, ich, rhop_knots_ne, n_e, n_e_dx2, rhop_knots_Te, T_e, T_e_dx2, &
!                       ne_rhop_scal, reflec_X_new, & ! in
!                       reflec_O_new, rp_min, &
!                       rhop_ECRad_out, BPD_ECRad_out, rhop_res_warm)
!rhop(1:pnts_BPD) = rhop_ECRad_out(1:pnts_BPD)
!BPD(1:pnts_BPD) = BPD_ECRad_out(1:pnts_BPD)
!deallocate(rhop_ECRad_out, BPD_ECRad_out)
!end subroutine make_BPD_w_res_ch_ECRad

end module ECRad_python


text:
http://cce.ucdavis.edu//content/view/15/33/ Wildlife Society meets in Haring Hall Haring 2016. They are an affiliate group of the Campus Center for the Environment.

Most recent activities seems to be handled through their http://www.facebook.com/reqs.php#/group.php?gid2223054541 facebook group.


text:
	SUBROUTINE NGD_NLST ( alias, namstr, lenstr, iret )
C************************************************************************
C* NGD_NLST								*
C*									*
C* This routine returns a list of names given GRID information.  The 	*
C* alias is used to locate the grid data files on disk.  The times are	*
C* returned in a single string separated by a semi-colon (;).		*
C*									*
C* NGD_NLST ( ALIAS, NAMSTR, LENSTR, IRET )				*
C*									*
C* Input parameters:							*
C*	ALIAS		CHAR*		Alias for GRID data		*
C*									*
C* Output parameters:							*
C*	NAMSTR		CHAR*		String containing all names	*
C*	LENSTR		INTEGER		Length of string 		*
C*	IRET		INTEGER		Return code			*
C*					  -3 = no files found		*
C*					  -6 = invalid template		*
C**									*
C* Log:									*
C* S. Jacobs/NCEP	 3/01	Created					*
C* A. Hardy/SAIC         2/02   Changed call FL_SCND			*
C* T. Lee/SAIC		 9/04	Replaced FL_TMPL with CTB_DTGET		*
C* A. Hardy/NCEP	11/04	Added calls to ST_RNUL			*
C* m.gamazaychikov/SAIC 12/04   Added ion flag to CTB_DTGET CS          *
C* m.gamazaychikov/SAIC 04/06   Added idtmch flag to CTB_DTGET CS       *
C* F. J. Yen/NCEP        4/08   Added bin mins & mstrct to CTB_DTGET CSC*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	alias, namstr
C*
	CHARACTER	path*(MXFLSZ), tmplt*(MXFLSZ),
     +			names(MXNMFL)*(MXFLSZ),
     +			namarr(MXNMFL)*(MXFLSZ),
     +			fnull*(MXFLSZ)
C------------------------------------------------------------------------
	iret = 0
        nexp = MXNMFL
C
C*	Get the info for the given alias.
C
	CALL ST_NULL ( alias, fnull, na, ier )
	path  = ' '
	tmplt = ' '
	CALL CTB_DTGET ( fnull, path, tmplt, ic, is, if, ir, ii, ion,
     +          ihb, mnb, iha, mna, mstrct, idtmch, ier )
	CALL ST_RNUL ( path, path, lens, ier )
	CALL ST_RNUL ( tmplt, tmplt, len, ier )
C
C*	Find the "*" in the template. This is where the name starts.
C
	CALL ST_NOCC ( tmplt, '*', 1, ipos, ier )
	IF  ( ( ier .ne. 0 ) .or. ( ipos .eq. 0 ) )  THEN
	    iret = -6
	    RETURN
	END IF
C
C*	Get all of the file names. Return if there are none.
C
	CALL FL_SCND ( path, tmplt, 1, nexp, names, nn, ier )
	IF  ( nn .eq. 0 )  THEN
	    iret = -3
	    RETURN
	END IF
C
C*	Keep only the "name" part of the file name.
C
	DO  i = 1, nn
	    namarr(i) = names(i)(ipos:)
	END DO
C
C*	Sort and get only the unique names.
C
	CALL ST_SORT ( 2, nn, namarr, nout, namarr, ier )
C
C*	Construct a single string from the array of names.
C
	CALL ST_LSTC ( namarr, nout, ';', namstr, ier )
	CALL ST_LSTR ( namstr, lenstr, ier )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Subroutine aran9(qran, ndim)
  Dimension qran(10)
  Common /sedvax/num1
  Save
  Do i = 1, ndim
    qran(i) = ranart(num1)
  End Do
  Return
End Subroutine aran9

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 03.03/06                              *
C  *    TEST TITLE : Retrieving open structure             *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

C conflict resolution
C                maintain    abandon     update
      INTEGER    PCRMNT,     PCRABA,     PCRUPD
      PARAMETER (PCRMNT = 0, PCRABA = 1, PCRUPD = 2)

      INTEGER    ARID, ARNM, SPCSTR(20), SPCSIZ

      CALL INITGL ('03.03/06')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)

C  open new archive file, with arid = archive identifier
      CALL AVARNM (ARNM)
      ARID = 6
      CALL POPARF (ARID, ARNM)

C  set retrieval conflict resolution flag = UPDATE
      CALL PSCNRS (PCRABA, PCRUPD)

C  set up structures:
C     archive : 102-A, 104-A (3 elements), with 102-A invoking 104-A
C     CSS     : 101-C, 104-C (4 elements)
      CALL PDAS
      CALL SETSTR (104, 'L911,L912,L913')
      CALL SETSTR (102, 'E104')
      CALL PDASAR (ARID)
      CALL PARAST (ARID)
      CALL PDAS
      CALL SETSTR (101, 'L931')
      CALL SETSTR (104, 'L931,L932,L933,L934')

C  <open structure> 104
      CALL POPST (104)
C  <retrieve all structures>
      CALL PRAST (ARID)

      CALL SETMSG ('11 16 18', 'When <retrieve all structures> '   //
     1             'causes a structure in the archive to replace ' //
     2             'the open CSS structure, the structure should ' //
     3             'remain open, with the element pointer '        //
     4             'indicating the last element.')

C  pass/fail depending on (104-A is open structure and
C                          element pointer = 3)
      CALL RETOPN (104, 3, '67,911,67,912,67,913')

C  set up structures:
C     CSS     : 101-C, 104-C (4 elements)
      CALL PDAS
      CALL SETSTR (101, 'L931')
      CALL SETSTR (104, 'L931,L932,L933,L934')

C  <open structure> 101
      CALL POPST (101)
C  <retrieve structures> 102, 101
      CALL SETVS  ('102,101', SPCSTR, SPCSIZ)
      CALL ERRCTL (.TRUE.)
      CALL PREST  (ARID, SPCSIZ, SPCSTR)
      CALL ERRCTL (.FALSE.)

      CALL SETMSG ('9 17 18', 'When <retrieve structures> specifies ' //
     1             'a non-existent archive structure to replace the ' //
     2             'open CSS structure, the structure should become ' //
     3             'empty but remain open, with the element pointer ' //
     4             'indicating zero.')

C  pass/fail depending on (101-E is open structure and
C                          element pointer = 0)
      CALL RETOPN (101, 0, ' ')

C  set up structures:
C     CSS     : 101-C, 104-C (2 elements)
      CALL PDAS
      CALL SETSTR (101, 'L931')
      CALL SETSTR (104, 'L931,L932')

C  <open structure> 104
      CALL POPST (104)
C  <retrieve structure networks> 102
      CALL SETVS  ('102', SPCSTR, SPCSIZ)
      CALL PRESN  (ARID, SPCSIZ, SPCSTR)

      CALL SETMSG ('10 16 18', 'When <retrieve structure networks> '  //
     1             'causes an implicitly specified structure in the ' //
     2             'archive to replace the open CSS structure, the '  //
     3             'structure should remain open, with the element '  //
     4             'pointer indicating the last element.')

C  pass/fail depending on (104-A is open structure and
C                          element pointer = 3)
      CALL RETOPN (104, 3, '67,911,67,912,67,913')

666   CONTINUE
C close archive and PHIGS
      CALL PCLARF (ARID)
      CALL ENDIT
      END


!
!     ------------------------------------------------------------
!        Z Z
!     ------------------------------------------------------------
!
      REAL(KIND(0.0D0)) FUNCTION ZZ (I1, I2, I3, I4, K)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE1=>DOUBLE
      use nel_C
!
!...Translated by Pacific-Sierra Research 77to90  4.3E  00:56:55  11/16/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE tk_I
      USE uk_I
      USE sn_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER  :: I1
      INTEGER  :: I2
      INTEGER  :: I3
      INTEGER  :: I4
      INTEGER  :: K
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER :: NWD = 128
      INTEGER, PARAMETER :: NOD = 220
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      REAL(DOUBLE1) :: C, C1, C2
!-----------------------------------------------

!      POINTER(IQP,P(NOD,1)),(IQN,N(1)),(IQL,L(1)),(IQAZ,AZ(1)),
!     :       (IQMAX,MAX_(1))
!      COMMON/NEL/IQP,IQN,IQL,IQAZ,IQMAX,IQ(7)
      C = 2*K*(K + 1)
      ZZ = -C*(TK(I1,I2,I3,I4,K + 1) - TK(I1,I2,I3,I4,K - 1))
!
      C1 = L(I1)*(L(I1)+1) - L(I3)*(L(I3)+1) - K*(K + 1)
      ZZ = ZZ - C1*(UK(I1,I2,I3,I4,K + 1) - UK(I1,I2,I3,I4,K - 1))
!
      C2 = L(I2)*(L(I2)+1) - L(I4)*(L(I4)+1) - K*(K + 1)
      ZZ = ZZ - C2*(UK(I2,I1,I4,I3,K + 1) - UK(I2,I1,I4,I3,K - 1))
!
      C = C1*C2/2
      C1 = C*(K - 2)/K/(K + K - 1)
      C2 = C*(K + 3)/(K + 1)/(K + K + 3)
      ZZ = ZZ - C1*(SN(I1,I2,I3,I4,K - 2) + SN(I2,I1,I4,I3,K - 2)) + C2*(SN(I1,&
         I2,I3,I4,K) + SN(I2,I1,I4,I3,K))
      RETURN
      END FUNCTION ZZ


text:
      PROGRAM xpolcof
C     driver for routine polcof
      INTEGER NP
      REAL PI
      PARAMETER(NP=5,PI=3.141593)
      INTEGER i,j,nfunc
      REAL f,sum,x,xa(NP),ya(NP),coeff(NP)
      do 15 nfunc=1,2
        if (nfunc.eq.1) then
          write(*,*) 'Sine function from 0 to PI'
          do 11 i=1,NP
            xa(i)=i*PI/NP
            ya(i)=sin(xa(i))
11        continue
        else if (nfunc.eq.2) then
          write(*,*) 'Exponential function from 0 to 1'
          do 12 i=1,NP
            xa(i)=1.0*i/NP
            ya(i)=exp(xa(i))
12        continue
        else
          stop
        endif
        call polcof(xa,ya,NP,coeff)
        write(*,*) '    coefficients'
        write(*,'(1x,6f12.6)') (coeff(i),i=1,NP)
        write(*,'(1x,t10,a1,t20,a4,t29,a10)')
     *       'x','f(x)','polynomial'
        do 14 i=1,10
          if (nfunc.eq.1) then
            x=(-0.05+i/10.0)*PI
            f=sin(x)
          else if (nfunc.eq.2) then
            x=-0.05+i/10.0
            f=exp(x)
          endif
          sum=coeff(NP)
          do 13 j=NP-1,1,-1
            sum=coeff(j)+sum*x
13        continue
          write(*,'(1x,3f12.6)') x,f,sum
14      continue
        write(*,*) '***********************************'
        write(*,*) 'Press RETURN'
        read(*,*)
15    continue
      END

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xpolcof.for


text:
module modules_05_mod
implicit none
private
integer, parameter, public :: a = 5
integer, parameter :: b = 6
integer, parameter :: c = 7
public :: c
end module

program modules_05
use modules_05_mod, only: a, c
print *, a, c
end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
! **********************************************************************
! *** LAMELLAR DI-BLOCK COPOLYMER + NANODIMERS *************************
! *** Dynamical Data Analysis                  *************************
! *** Author: J. Javier Burgos Mármol -- 2021  *************************
! **********************************************************************
!BSD 3-Clause License
!
!Copyright (c) 2021, J. Javier Burgos-Mármol (Github username: jjavier-bm)
!All rights reserved.
!
!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. Neither the name of the copyright holder nor the names of its
!   contributors may be used to endorse or promote products derived from
!   this software without specific prior 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 COPYRIGHT HOLDER OR 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.
PROGRAM LamellarDynamics
IMPLICIT NONE
!---DATA------------------------------------
!CONSTANTS
REAL*8, PARAMETER :: PI=3.141592654d0, TWOPI=6.283185307d0, HALFPI=1.570796327d0
!PARAMETERS
INTEGER, PARAMETER :: Nchains = 12696, LchainsA = 15, LchainsB = 15	!Number of Chains and length of each of the blocks
REAL*8, PARAMETER :: Mbeadchains=1.d0 ! Mass of each bead in the chains
INTEGER, PARAMETER :: Npart = 1080 !Number of nano-dimers
REAL*8, PARAMETER :: NPdiam=3.d0, RbondEqNP=2.5d0 ! NP bead diameter, NP bond eq. distance (mass derived from here)
REAL*8, PARAMETER :: delta_t=0.001d0
INTEGER, PARAMETER :: normaldir=3
INTEGER, PARAMETER :: NtimestepsperFile=1000000                 !Number of timesteps per file (power of 10)
INTEGER, PARAMETER :: Nfiles=100                                !Number of statistical repetitions (files, power of 10)
!OTHER PARAMETERS
INTEGER, PARAMETER :: Lchains=LchainsA+LchainsB
REAL*8, PARAMETER :: Vlens=(PI*(NPdiam-RbondEqNP)**2*(RbondEqNP**2+2.d0*RbondEqNP*NPdiam)/(12.d0*RbondEqNP)) ! Overlapping volume (lens shaped)
REAL*8, PARAMETER :: MbeadNpart=0.74d0*(NPdiam**3/6.d0-Vlens/2.d0) ! Mass of Nano-dimer's beads (scalates as FCC packed)
REAL*8, PARAMETER :: Mass=dfloat(NChains*Lchains)*Mbeadchains+dfloat(2*Npart)*MbeadNpart
INTEGER, PARAMETER :: finaltimestep=NtimestepsperFile*Nfiles !Final timestep printed (in a log scale)
INTEGER, PARAMETER :: Nconfigs=int(dlog10(dfloat(finaltimestep)))*9&
                               +finaltimestep/(10**(int(dlog10(dfloat(finaltimestep))))) !Total number of configurations
INTEGER, PARAMETER :: NconfigsperFile=int(dlog10(dfloat(NtimestepsperFile)))*9&
                                      +NtimestepsperFile/(10**(int(dlog10(dfloat(NtimestepsperFile))))) !Number of configurations per file
!OTHER VARIABLES
REAL*8, DIMENSION(:) :: RChain(Nchains,0:Lchains,3) !Positions of chains ! 0 is for the CoM
REAL*8, DIMENSION(:) :: Rpart(Npart,0:2,3) !Positions of Nano-dimers ! 0 is for the CoM
REAL*8, DIMENSION(:) :: R0Chain(Nchains,0:Lchains,3,0:Nfiles) !Positions of chains at t=0! 0 is for the CoM
REAL*8, DIMENSION(:) :: R0part(Npart,0:2,3,0:Nfiles) !Positions of Nano-dimers at t=0! 0 is for the CoM
REAL*8, DIMENSION(:) :: RCoM(3) !Position of Centre of Mass
REAL*8, DIMENSION(:) :: R0CoM(3,0:Nfiles) !Position of Centre of Mass at t=0
REAL*8, DIMENSION(:) :: pos(3), vel(3), jvol(6) ! jvol(stresstensor*volume): xx,yy,zz,xy,xz,yz
INTEGER, DIMENSION(:) :: t(0:Nconfigs), avconf(Nconfigs)
INTEGER ::  i, j, k, l, conf, conf0, d, d2, f, f2, id, idmol, idbead, tipo
REAL*8, DIMENSION(:) :: msd(Nconfigs,6,0:4) ! 1=ChainMonomers, 2=Chain Centre beads, 3=Chain End beads, 4=Chain CoMs, 5=Particle Beads, 6=Particle CoMs; 0=3D, 1=x,2=y,3=z,4=2D plane
REAL*8, DIMENSION(:) :: avdir(0:Nconfigs), integral(Nconfigs)
!***********************************************************************************************************************************************
R0chain=0.d0
R0part=0.d0
R0CoM=0.d0
msd=0.d0
avdir=0.d0
avdir(0)=1.d0
!print*, Nconfigs, NConfigsPerFile
DO conf=0,Nconfigs
    IF (conf<10) THEN
        t(conf)=conf
    ELSE
        t(conf)=t(conf-9)*10
    END IF
    IF (conf>=0 .and. conf<=NconfigsPerFile) avconf(conf)=Nfiles
    IF (conf>NconfigsperFile) avconf(conf)=0
END DO !i N configs
DO f=1,Nfiles
    open(100+f)
    !print*, 'Reading file fort.',100+f
    IF (f==1) conf0=0
    IF (f>1) conf0=1
    DO conf=conf0, NconfigsperFile
        !print*, conf
        Rchain=0.d0
        Rpart=0.d0
        RCoM=0.d0
        DO i=1,9
            READ(100+f,*)
        END DO !i
        ! READ THE COORDINATES
        DO i=1,Nchains*Lchains+Npart*2
            READ(100+f,*) id, pos(1), pos(2), pos(3), vel(1), vel(2), vel(3),&
                         jvol(1), jvol(2), jvol(3), jvol(4), jvol(5), jvol(6)
            IF (id<=Nchains*Lchains) THEN
                idmol=int((id-1)/Lchains)+1
                idbead=id-(idmol-1)*Lchains
                DO d=1,3
                    Rchain(idmol,idbead,d)=pos(d)
                    Rchain(idmol,0,d)=Rchain(idmol,0,d)+pos(d)*MbeadChains
                    RCoM(d)=RCoM(d)+pos(d)*MbeadChains
                END DO !dims
            ELSE
                idmol=int((id-Nchains*Lchains-1)/2)+1
                idbead=id-Nchains*Lchains-(idmol-1)*2
                DO d=1,3
                    Rpart(idmol,idbead,d)=pos(d)
                    Rpart(idmol,0,d)=Rpart(idmol,0,d)+pos(d)*MbeadNpart
                    RCoM(d)=RCoM(d)+pos(d)*MbeadNpart
                END DO !dims
            END IF !id
        END DO !i
        ! COMPUTE CoM
        DO d=1,3
            DO i=1,Nchains
                Rchain(i,0,d)=Rchain(i,0,d)/(dfloat(Lchains)*MbeadChains)
                IF (conf==0 .or. conf==NconfigsperFile) THEN
                    DO j=0,Lchains
                        IF (conf==0) R0chain(i,j,d,0)=Rchain(i,j,d)
                        IF (conf==NconfigsperFile) R0chain(i,j,d,f)=Rchain(i,j,d)
                    END DO !j Lchains
                END IF !terminal confs
            END DO !i Nchains
            DO i=1,Npart
                Rpart(i,0,d)=Rpart(i,0,d)/(2.d0*MbeadNpart)
                IF (conf==0 .or. conf==NconfigsperFile) THEN
                    DO j=0,2
                        IF (conf==0) R0part(i,j,d,0)=Rpart(i,j,d)
                        IF (conf==NconfigsperFile) R0part(i,j,d,f)=Rpart(i,j,d)
                    END DO !j dimer beads and CoM
                END IF !terminal confs
            END DO !i Npart
            RCoM(d)=RCoM(d)/Mass
            IF (conf==0) R0CoM(d,0)=RCoM(d)
            IF (conf==NconfigsperFile) R0CoM(d,f)=RCoM(d)
        END DO !dims
        IF (conf>0) THEN
            ! Compute MSD and director
            DO i=1,Nchains
                DO d=1,3
                    DO j=1,Lchains
                        msd(conf,1,d)=msd(conf,1,d)+(Rchain(i,j,d)-R0chain(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !All chain beads, 1D
                        msd(conf,1,0)=msd(conf,1,0)+(Rchain(i,j,d)-R0chain(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 ! 3D
                        IF (d/=normaldir) msd(conf,1,4)=msd(conf,1,4)+(Rchain(i,j,d)-R0chain(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !2D
                        IF (j==LchainsA .or. j==LchainsA+1) THEN !Central beads
                            msd(conf,2,d)=msd(conf,2,d)+(Rchain(i,j,d)-R0chain(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !1D
                            msd(conf,2,0)=msd(conf,2,0)+(Rchain(i,j,d)-R0chain(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !3D
                            IF (d/=normaldir) msd(conf,2,4)=msd(conf,2,4)+(Rchain(i,j,d)-R0chain(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !2D
                        ELSE IF (j==1 .or. j==Lchains) THEN !End beads
                            msd(conf,3,d)=msd(conf,3,d)+(Rchain(i,j,d)-R0chain(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !1D
                            msd(conf,3,0)=msd(conf,3,0)+(Rchain(i,j,d)-R0chain(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !3D
                            IF (d/=normaldir) msd(conf,3,4)=msd(conf,3,4)+(Rchain(i,j,d)-R0chain(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !2D
                        END IF!bead type
                    END DO !j Lchains
                    msd(conf,4,d)=msd(conf,4,d)+(Rchain(i,0,d)-R0chain(i,0,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !Chain CoMs, 1D
                    msd(conf,4,0)=msd(conf,4,0)+(Rchain(i,0,d)-R0chain(i,0,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !3D
                    IF (d/=normaldir) msd(conf,4,4)=msd(conf,4,4)+(Rchain(i,0,d)-R0chain(i,0,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !2D
                END DO !dims
            END DO !i Nchains
            DO i=1,Npart
                DO d=1,3
                    DO j=1,2
                        msd(conf,5,d)=msd(conf,5,d)+(Rpart(i,j,d)-R0part(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !All ND beads, 1D
                        msd(conf,5,0)=msd(conf,5,0)+(Rpart(i,j,d)-R0part(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 ! 3D
                        IF (d/=normaldir) msd(conf,5,4)=msd(conf,5,4)+(Rpart(i,j,d)-R0part(i,j,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !2D
                    END DO !dimer bead
                    msd(conf,6,d)=msd(conf,6,d)+(Rpart(i,0,d)-R0part(i,0,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !NP CoMs, 1D
                    msd(conf,6,0)=msd(conf,6,0)+(Rpart(i,0,d)-R0part(i,0,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !3D
                    IF (d/=normaldir) msd(conf,6,4)=msd(conf,6,4)+(Rpart(i,0,d)-R0part(i,0,d,f-1)-RCoM(d)+R0CoM(d,f-1))**2 !2D
                    avdir(conf)=avdir(conf)+(Rpart(i,2,d)-Rpart(i,1,d))*(R0part(i,2,d,f-1)-R0part(i,1,d,f-1))&
                                           /(dsqrt((Rpart(i,2,1)-Rpart(i,1,1))**2&
                                                  +(Rpart(i,2,2)-Rpart(i,1,2))**2&
                                                  +(Rpart(i,2,3)-Rpart(i,1,3))**2)&
                                            *dsqrt((R0part(i,2,1,f-1)-R0part(i,1,1,f-1))**2&
                                                  +(R0part(i,2,2,f-1)-R0part(i,1,2,f-1))**2&
                                                  +(R0part(i,2,3,f-1)-R0part(i,1,3,f-1))**2))
                END DO !dims
            END DO !Npart
        END IF !conf
    END DO !conf
    close(100+f)
    !print*, '    --> file read'
    !print*, 'Accumulators'
    ! Multi-file statistics (only for conf>NconfigsperFile)
    DO f2=1,f-1
        conf=0
        DO l=NconfigsperFile+1,Nconfigs
            IF ((f-(f2-1))*NtimestepsperFile==t(l)) conf=l
        END DO !l loop
        IF (conf>0) THEN
            avconf(conf)=avconf(conf)+1
            DO d=1,3
                DO i=1,Nchains
                    DO j=1,Lchains
                        msd(conf,1,d)=msd(conf,1,d)+&
                                    +(Rchain(i,j,d)-R0chain(i,j,d,f2-1)&
                                    -RCoM(d)+R0CoM(d,f2-1))**2 !All chain beads, 1D
                        msd(conf,1,0)=msd(conf,1,0)+&
                                    +(Rchain(i,j,d)-R0chain(i,j,d,f2-1)&
                                    -RCoM(d)+R0CoM(d,f2-1))**2 !All chain beads, 3D
                        IF (d/=normaldir) msd(conf,1,4)=msd(conf,1,4)+&
                                    +(Rchain(i,j,d)-R0chain(i,j,d,f2-1)&
                                    -RCoM(d)+R0CoM(d,f2-1))**2 !All chain beads, 2D
                        IF (j==LchainsA .or. j==LchainsA+1) THEN !Central beads
                            msd(conf,2,d)=msd(conf,2,d)+&
                                        +(Rchain(i,j,d)-R0chain(i,j,d,f2-1)&
                                        -RCoM(d)+R0CoM(d,f2-1))**2 !1D
                            msd(conf,2,0)=msd(conf,2,0)+&
                                        +(Rchain(i,j,d)-R0chain(i,j,d,f2-1)&
                                        -RCoM(d)+R0CoM(d,f2-1))**2 !3D
                            IF (d/=normaldir) msd(conf,2,4)=msd(conf,2,4)+&
                                        +(Rchain(i,j,d)-R0chain(i,j,d,f2-1)&
                                        -RCoM(d)+R0CoM(d,f2-1))**2 !2D
                        ELSE IF (j==1 .or. j==Lchains) THEN !End beads
                            msd(conf,3,d)=msd(conf,3,d)+&
                                        +(Rchain(i,j,d)-R0chain(i,j,d,f2-1)&
                                        -RCoM(d)+R0CoM(d,f2-1))**2 !1D
                            msd(conf,3,0)=msd(conf,3,0)+&
                                        +(Rchain(i,j,d)-R0chain(i,j,d,f2-1)&
                                        -RCoM(d)+R0CoM(d,f2-1))**2 !3D
                            IF (d/=normaldir) msd(conf,3,4)=msd(conf,3,4)+&
                                        +(Rchain(i,j,d)-R0chain(i,j,d,f2-1)&
                                        -RCoM(d)+R0CoM(d,f2-1))**2 !2D
                        END IF !bead type
                    END DO !j Lchains
                    msd(conf,4,d)=msd(conf,4,d)+&
                                +(Rchain(i,0,d)-R0chain(i,0,d,f2-1)&
                                -RCoM(d)+R0CoM(d,f2-1))**2 !Chain CoMs, 1D
                    msd(conf,4,0)=msd(conf,4,0)+&
                                +(Rchain(i,0,d)-R0chain(i,0,d,f2-1)&
                                -RCoM(d)+R0CoM(d,f2-1))**2 !Chain CoMs, 3D
                    IF (d/=normaldir) msd(conf,4,4)=msd(conf,4,4)+&
                                +(Rchain(i,0,d)-R0chain(i,0,d,f2-1)&
                                -RCoM(d)+R0CoM(d,f2-1))**2 !Chain CoMs, 2D
                END DO !i Nchains
                DO i=1, Npart
                    DO j=1, 2
                        msd(conf,5,d)=msd(conf,5,d)&
                                +(Rpart(i,j,d)-R0part(i,j,d,f2-1)-RCoM(d)+R0CoM(d,f2-1))**2 !All ND beads, 1D
                        msd(conf,5,0)=msd(conf,5,0)&
                                +(Rpart(i,j,d)-R0part(i,j,d,f2-1)-RCoM(d)+R0CoM(d,f2-1))**2 ! 3D
                        IF (d/=normaldir) msd(conf,5,4)=msd(conf,5,4)&
                                +(Rpart(i,j,d)-R0part(i,j,d,f2-1)-RCoM(d)+R0CoM(d,f2-1))**2 !2D
                    END DO !j dimer bead
                    msd(conf,6,d)=msd(conf,6,d)&
                            +(Rpart(i,0,d)-R0part(i,0,d,f2-1)-RCoM(d)+R0CoM(d,f2-1))**2 !NP CoMs, 1D
                    msd(conf,6,0)=msd(conf,6,0)&
                            +(Rpart(i,0,d)-R0part(i,0,d,f2-1)-RCoM(d)+R0CoM(d,f2-1))**2 !3D
                    IF (d/=normaldir) msd(conf,6,4)=msd(conf,6,4)&
                            +(Rpart(i,0,d)-R0part(i,0,d,f2-1)-RCoM(d)+R0CoM(d,f2-1))**2 !2D
                    avdir(conf)=avdir(conf)&
                            +(Rpart(i,2,d)-Rpart(i,1,d))*(R0part(i,2,d,f2-1)-R0part(i,1,d,f2-1))&
                            /(dsqrt((Rpart(i,2,1)-Rpart(i,1,1))**2&
                                   +(Rpart(i,2,2)-Rpart(i,1,2))**2&
                                   +(Rpart(i,2,3)-Rpart(i,1,3))**2)&
                             *dsqrt((R0part(i,2,1,f2-1)-R0part(i,1,1,f2-1))**2&
                                   +(R0part(i,2,2,f2-1)-R0part(i,1,2,f2-1))**2&
                                   +(R0part(i,2,3,f2-1)-R0part(i,1,3,f2-1))**2))
                END DO !i Npart
            END DO !dims
        END IF !conf
    END DO !f2 files
    !print*, '    --> accumulation finished'
END DO !f

!print*, 'Averaging'
! DIVIDE BY NUMBER OF TIMES COMPUTED
DO i=1,6
    IF (i==1) d2=Nchains*Lchains
    IF (i==2 .or. i==3) d2=Nchains*2
    IF (i==4) d2=Nchains
    IF (i==5) d2=Npart*2
    IF (i==6) d2=Npart
    DO conf=1,Nconfigs
        IF (Npart>0 .and. i==6) avdir(conf)=avdir(conf)/dfloat(d2*avconf(conf))
        DO d=0,4
            msd(conf,i,d)=msd(conf,i,d)/dfloat(d2*avconf(conf))
        END DO !dims
    END DO !conf Nconfigs
END DO !i types
!print*, '    --> average done'
!print*, 'Integrating autocorrelation'
! Trapezoid rule
IF (Npart>0) THEN
    integral=0.d0
    DO i=1,Nconfigs
        DO j=1,i
            integral(i)=integral(i)+0.5d0*(avdir(j)+avdir(j-1))*delta_t*dfloat(t(j)-t(j-1))
        END DO !j 1,i
    END DO !i Nconfigs
!print*, '    --> integration done'
!print*, 'Printing outputs'
!OUTPUT MSD, ROT
open(501,file='msd_allmonomers.dat')
open(502,file='msd_centralmonomers.dat')
open(503,file='msd_endingmonomers.dat')
open(504,file='msd_chainsCM.dat')
IF (Npart>0) THEN
    open(505,file='msd_NPbeads.dat')
    open(506,file='msd_NPCM.dat')
END IF
tipo=4
IF (Npart>0) tipo=6
DO i=1,tipo
    WRITE(500+i,*) '# time(tau) msd_x msd_y msd_z msd_plane msd_3d'
    DO conf=1, Nconfigs
        WRITE(500+i,500) dfloat(t(conf))*delta_t, msd(conf,i,1), msd(conf,i,2),&
            msd(conf,i,3), msd(conf,i,4), msd(conf,i,0)
    END DO !Nconfigs
    close(500+i)
END DO !Ntypes
500 FORMAT (F14.6,2X,F20.12,2X,F20.12,2X,F20.12,2X,F20.12,2X,F20.12)
open(601,file='avdir_autocorr_Drot.dat')
WRITE(601,*) '# time(tau) avdir ln_avdir integral_avdir(tau) D_rot(tau^-1)'
DO i=1,Nconfigs
    WRITE(601,600) dfloat(t(i))*delta_t, avdir(i), dlog(avdir(i)), integral(i), 0.5d0/integral(i)
END DO
600 FORMAT (F14.3,2X,F20.12,2X,F20.12,2X,F20.12,2X,F20.12)
close(601)

END PROGRAM LamellarDynamics

max_stars_repo_name: jjavier-bm/LamAnalysis

max_issues_repo_name: jjavier-bm/LamAnalysis

max_forks_repo_name: jjavier-bm/LamAnalysis


text:
passive component C {
}
instance c: C base id 0x100
topology T {
  connections C {
    c.out -> c.in
  }
}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
      SUBROUTINE BMG3_SymStd_SETUP_cg_boxmg(&
     &                II, JJ, KK, NGx, NGy, NGz,&
     &                iGs, jGs, kGs,&
     &                BMG_iPARMS, BMG_rPARMS, BMG_IOFLAG,&
     &                SO, NStncl, NOG,&
     &                BMG_iWORK_CS, NBMG_iWORK_CS,&
     &                BMG_rWORK_CS, NBMG_rWORK_CS,&
     &                WS, NMSGr, &
     &                NProcI, NProcJ, NProcK, NProc, MyProc,&
     &                ProcGrid, ProcCoord, LocArrSize, MPICOMM&
     &                ) BIND(C, NAME='BMG3_SymStd_SETUP_cg_boxmg')

! ======================================================================
!  --------------------
!   DESCRIPTION:
!  --------------------
!
!     BMG3_SymStd_SETUP_cg_boxmg sets up the stencil, and workspace,
!     to use the serial version of BoxMG for the coarse-grid solve.
!
! ======================================================================
! $license_flag$
! ======================================================================
!  --------------------
!   INPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   INPUT/OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   OUTPUT:
!  --------------------
!
!
!
! ======================================================================
!  --------------------
!   LOCAL:
!  --------------------
!
!
! ======================================================================

      IMPLICIT NONE

! ----------------------------
!     Includes
!
      INCLUDE 'BMG_stencils_f90.h'
      INCLUDE 'BMG_constants_f90.h'
      INCLUDE 'BMG_workspace_f90.h'
      INCLUDE 'BMG_parameters_f90.h'

      INCLUDE 'mpif.h'

! ----------------------------
!     Argument Declarations
! ----------------------------

      !
      !  Global/Local indexing
      !
      INTEGER II, JJ, KK,&
     &        NGx, NGy, NGz,&
     &        iGs, jGs, kGs

      !
      !  Stencil neighbours
      !
      INTEGER  NStncl, NOG
      REAL*8   SO(II+1,JJ+1,KK+1,NStncl)

      !
      !  Processor grid
      !
      INTEGER NProcI, NProcJ, NProcK, NProc, MyProc
      INTEGER ProcCoord(3, NProc),&
     &        ProcGrid(NProcI,NProcJ,NProcK)

      !
      !  MPI REAL buffer space
      !
      INTEGER NMSGr
      REAL*8  WS(NMSGr)

      !
      !  MPI communicator
      !
      INTEGER  MPICOMM

      !
      !
      !
      INTEGER LocArrSize(3,*)

      !
      !
      !
      INTEGER BMG_iPARMS(NBMG_iPARMS)
      REAL*8  BMG_rPARMS(NBMG_rPARMS)
      LOGICAL BMG_IOFLAG(NBMG_IOFLAG)

      !
      !  Coarse-grid solve workspace
      !
      INTEGER  NBMG_iWORK_CS, NBMG_rWORK_CS
      INTEGER  BMG_iWORK_CS(NBMG_iWORK_CS)
      REAL*8   BMG_rWORK_CS(NBMG_rWORK_CS)

! ----------------------------
!     Local Declarations
! ----------------------------

      INTEGER i, i1, i2, ibw, j, j1, k, k1, kl, KK_t, n
      INTEGER IIG, JJG, KKG, proc, KKMAX, P1, P2, P3
      INTEGER P1SUM, P2SUM, P3SUM, PXP1, PYP2, PZP3, p_WS
      INTEGER INT_TEMP, INT_TEMP1, INT_TEMP2, INT_TEMP3
      INTEGER LARGESTNODES, IERR, MPI_IERROR

      INTEGER i_WS, iGs_ws, jGs_ws, kGs_ws, NLx_ws, NLy_ws, NLz_ws,&
     &        Proc_ws, ProcI_ws, ProcJ_ws, ProcK_ws

      INTEGER p_pWORK, p_iPARMS, p_iWORK, p_iWORK_PL, &
     &        p_rPARMS, p_rWORK, p_rWORK_PL
      INTEGER NBMG_SER_iWORK, NBMG_SER_iWORK_PL, &
     &        NBMG_SER_rWORK, NBMG_SER_rWORK_PL

      INTEGER NOGm_SER, NOG_SER, NFm_SER, NSOm_SER

      INTEGER NF, NC, NCI, NSO, NSOR, NCBW, NCU,&
     &        p_CI, p_CSO, p_CU, p_iGRD, p_Q,&
     &        p_RES, p_SO, p_SOR, p_U

      INTEGER pSI, pSR
      LOGICAL BMG_SER_IOFLAG(NBMG_SER_IOFLAG),&
     &        BMG_SER_InWORK(NBMG_SER_InWORK)

      INTEGER ifd, kg

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

      i1=ii-1
      j1=jj-1
      k1=kk-1

      i2=i1-1
!      n=i2*(j1-1)*(kk-2)

! ------------------------------------------------
!     Calculate the global number of points
! ------------------------------------------------

      !
      ! Global number in x
      !
      IIG=2
      DO i=1,NProcI
         IIG = IIG + LocArrSize(1,ProcGrid(i,1,1)) - 2
      ENDDO

      !
      ! Global number in y
      !
      JJG=2
      DO j=1,NProcJ
         JJG = JJG + LocArrSize(2,ProcGrid(1,j,1)) - 2
      ENDDO

      !
      ! Global number in z
      !
      KKG=2
      DO k=1,NProcK
         KKG = KKG + LocArrSize(3,ProcGrid(1,1,k)) - 2
      ENDDO

! ------------------------------------------------
!     Find the largest local array
! ------------------------------------------------

      !
      ! Loop over local x dimensions
      !
      INT_TEMP = LocArrSize(1,1) - 2
      INT_TEMP1 = INT_TEMP
      DO proc=2, NProcI*NProcJ*NProcK
         INT_TEMP = LocArrSize(1,proc) - 2
         IF (INT_TEMP.gt.INT_TEMP1) THEN
            INT_TEMP1 = INT_TEMP
         END IF
      END DO

      !
      ! Loop over local y dimensions
      !
      INT_TEMP = LocArrSize(2,1) - 2
      INT_TEMP2 = INT_TEMP
      DO proc=2, NProcI*NProcJ*NProcK
         INT_TEMP = LocArrSize(2,proc) - 2
         IF (INT_TEMP.gt.INT_TEMP2) THEN
            INT_TEMP2 = INT_TEMP
         END IF
      END DO

      !
      ! Loop over local z dimensions
      !
      INT_TEMP = LocArrSize(3,1) - 2
      INT_TEMP3 = INT_TEMP
      DO proc=2, NProcI*NProcJ*NProcK
         INT_TEMP = LocArrSize(3,proc) - 2
         IF (INT_TEMP.gt.INT_TEMP3) THEN
            INT_TEMP3 = INT_TEMP
         END IF
      END DO

      !
      ! Conservative: take largest from each
      !
      LARGESTNODES = INT_TEMP1 * INT_TEMP2 * INT_TEMP3 * 14 + 7

! ------------------------------------------------
!     Copy all information into the buffer
! ------------------------------------------------

      WS(1) = MyProc
      WS(2) = II
      WS(3) = JJ
      WS(4) = KK
      WS(5) = iGs
      WS(6) = jGs
      WS(7) = kGs
      INT_TEMP = 8

      IF ( NStncl.EQ.14 ) THEN
         !
         DO K = 2, KK - 1
            DO J = 2, JJ - 1
               DO I = 2, II - 1
                  WS(INT_TEMP)    = SO(I  ,J  ,k  ,kp)
                  WS(INT_TEMP+ 1) = SO(I  ,J  ,k  ,kpw)
                  WS(INT_TEMP+ 2) = SO(I+1,J  ,k  ,kpnw)
                  WS(INT_TEMP+ 3) = SO(I  ,J  ,k  ,kps)
                  WS(INT_TEMP+ 4) = SO(I  ,J  ,k  ,kpsw)
                  WS(INT_TEMP+ 5) = SO(I+1,J+1,k  ,kbne)
                  WS(INT_TEMP+ 6) = SO(I  ,J+1,k  ,kbn)
                  WS(INT_TEMP+ 7) = SO(I  ,J+1,k  ,kbnw)
                  WS(INT_TEMP+ 8) = SO(I+1,J  ,k  ,kbe)
                  WS(INT_TEMP+ 9) = SO(I  ,J  ,k  ,kb)
                  WS(INT_TEMP+10) = SO(I  ,J  ,k  ,kbw)
                  WS(INT_TEMP+11) = SO(I+1,J  ,k  ,kbse)
                  WS(INT_TEMP+12) = SO(I  ,J  ,k  ,kbs)
                  WS(INT_TEMP+13) = SO(I  ,J  ,k  ,kbsw)
                  INT_TEMP = INT_TEMP + 14
               END DO
            END DO
         END DO
         !
      ELSE IF ( NStncl.EQ.4 ) THEN
         !
         DO K = 2, KK - 1
            DO J = 2, JJ - 1
               DO I = 2, II - 1
                  WS(INT_TEMP)    = SO(I  ,J  ,k  ,kp)
                  WS(INT_TEMP+ 1) = SO(I  ,J  ,k  ,kpw)
                  WS(INT_TEMP+ 2) = rZERO
                  WS(INT_TEMP+ 3) = SO(I  ,J  ,k  ,kps)
                  WS(INT_TEMP+ 4) = rZERO
                  WS(INT_TEMP+ 5) = rZERO
                  WS(INT_TEMP+ 6) = rZERO
                  WS(INT_TEMP+ 7) = rZERO
                  WS(INT_TEMP+ 8) = rZERO
                  WS(INT_TEMP+ 9) = SO(I  ,J  ,k  ,kb)
                  WS(INT_TEMP+10) = rZERO
                  WS(INT_TEMP+11) = rZERO
                  WS(INT_TEMP+12) = rZERO
                  WS(INT_TEMP+13) = rZERO
                  INT_TEMP = INT_TEMP + 14
               END DO
            END DO
         END DO
         !
      ENDIF

! ------------------------------------------------
!     Send/Receive information to/from everybody
! ------------------------------------------------

      IF ( BMG_iPARMS(id_BMG3_CG_COMM).EQ.BMG_CG_ALLGATHER ) THEN

         CALL MPI_ALLGATHER ( &
     &        WS(1), LARGESTNODES, MPI_DOUBLE_PRECISION,&
     &        WS(LARGESTNODES+1), LARGESTNODES, MPI_DOUBLE_PRECISION,&
     &        MPICOMM, IERR&
     &        )

      ELSE

         CALL MPI_GATHER ( &
     &        WS(1), LARGESTNODES, MPI_DOUBLE_PRECISION,&
     &        WS(LARGESTNODES+1), LARGESTNODES, MPI_DOUBLE_PRECISION,&
     &        iZERO, MPICOMM, IERR&
     &        )

      ENDIF

! ------------------------------------------------
!     Setup BoxMG serial:
! ------------------------------------------------

      !
      ! Setup local pointers
      !
      p_pWORK  = BMG_iWORK_CS(ip_BMG_pWORK_CS)

      p_iPARMS = BMG_iWORK_CS(ip_BMG_iPARMS_CS)
      p_rPARMS = BMG_iWORK_CS(ip_BMG_rPARMS_CS)

      p_iWORK  = BMG_iWORK_CS(ip_BMG_iWORK_CS)
      p_rWORK  = BMG_iWORK_CS(ip_BMG_rWORK_CS)

      p_iWORK_PL = BMG_iWORK_CS(ip_BMG_iWORK_PL_CS)
      p_rWORK_PL = BMG_iWORK_CS(ip_BMG_rWORK_PL_CS)

      !
      ! Synchronize the default parameter values
      !
      CALL BMG3_SymStd_SETUP_cg_parms(&
     &          BMG_iPARMS, BMG_rPARMS, BMG_IOFLAG,&
     &          BMG_iWORK_CS(p_iPARMS), &
     &          BMG_rWORK_CS(p_rPARMS), &
     &          BMG_SER_IOFLAG,&
     &          NGx, NGy, NGz, NOG&
     &          )

      !
      ! Setup workspace dimensions
      !
      NBMG_SER_iWORK = BMG_iWORK_CS(id_BMG_iWORK_CS)
      NBMG_SER_rWORK = BMG_iWORK_CS(id_BMG_rWORK_CS)

      NBMG_SER_iWORK_PL = BMG_iWORK_CS(id_BMG_iWORK_PL_CS)
      NBMG_SER_rWORK_PL = BMG_iWORK_CS(id_BMG_rWORK_PL_CS)

      !
      ! Set the maximum number of grids
      !
      NOGm_SER = BMG_iWORK_CS(id_BMG_NOGm_CS)

      !
      ! Allocate multigrid related arrays in the workspace
      !
      DO i=1, NBMG_SER_InWORK
         BMG_SER_InWORK(i) = .FALSE.
      ENDDO

      BMG_SER_InWork(i_InWORK_SO)  = .TRUE.
      BMG_SER_InWork(i_InWORK_U)   = .TRUE.
      BMG_SER_InWork(i_InWORK_Q)   = .TRUE.
      BMG_SER_InWork(i_InWORK_RES) = .TRUE.

      !
      ! Initialize irrelevant data
      !
      NFm_SER = 1
      NSOm_SER = 1

      !
      ! Restore the free space counters for BMG_SER_[ir]WORK
      !
      pSI = p_iWORK
      pSR = p_rWORK

      !
      ! Get the workspace estimate ...
      !
      CALL BMG3_SER_SymStd_SETUP_PtrWork( &
     &          NGx-2, NGy-2, NGz-2,&
     &          BMG_iWORK_CS(p_iPARMS),&
     &          NOGm_SER, NFm_SER, NSOm_SER, &
     &          NBMG_SER_iWORK, NBMG_SER_rWORK,&
     &          NBMG_SER_iWORK_PL, NBMG_SER_rWORK_PL,&
     &          BMG_iWORK_CS(p_pWORK), BMG_SER_InWork,&
     &          pSR, pSI&
     &          )

      !
      !  Local pointers for the SERIAL call
      !
      p_SO  = BMG_iWORK_CS(p_pWORK+ip_SO-1)
      p_U   = BMG_iWORK_CS(p_pWORK+ip_U-1)
      p_Q   = BMG_iWORK_CS(p_pWORK+ip_Q-1)
      p_CI  = BMG_iWORK_CS(p_pWORK+ip_CI-1)
      p_RES = BMG_iWORK_CS(p_pWORK+ip_RES-1)
      p_SOR = BMG_iWORK_CS(p_pWORK+ip_SOR-1)
      p_CSO = BMG_iWORK_CS(p_pWORK+ip_CSO-1)
      p_CU  = BMG_iWORK_CS(p_pWORK+ip_CU-1)

      p_iGRD  = BMG_iWORK_CS(p_pWORK+ip_iG-1)

      !
      !  Local dimensional parameters for the serial call
      !
      NOG_SER  = BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_DIM_NOG-1)

      NF       = BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_DIM_NF-1)
      NC       = BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_DIM_NC-1)
      NSO      = BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_DIM_NSO-1)
      NCI      = BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_DIM_NCI-1)
      NSOR     = BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_DIM_NSOR-1)
      NCBW     = BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_DIM_NCBW-1)
      NCU      = BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_DIM_NCU-1)

      IF ( BMG_iPARMS(id_BMG3_CG_COMM).EQ.BMG_CG_ALLGATHER&
     &    .OR. MyProc.EQ.iONE ) THEN

         !
         !  Copy WS blocks into the SERIAL stencil
         !
         DO proc=1, NProcI*NProcJ*NProcK

            i_WS = proc * LARGESTNODES + 1

            Proc_ws = WS(i_WS)
            i_WS = i_WS + 1

            NLx_ws = WS(i_WS)
            NLy_ws = WS(i_WS+1)
            NLz_ws = WS(i_WS+2)
            i_WS = i_WS + 3

            iGs_ws = WS(i_WS)
            jGs_ws = WS(i_WS+1)
            kGs_ws = WS(i_WS+2)
            i_WS = i_WS + 3

            IF ( NLx_ws .NE. LocArrSize(1,Proc_ws)&
     &         .OR.  NLy_ws .NE. LocArrSize(2,Proc_ws) &
     &         .OR.  NLz_ws .NE. LocArrSize(3,Proc_ws)  ) THEN
               !
               IERR = IERR + 1
               !
               WRITE(*,*) 'Error: LocArrSize is inconsitent ... '
               WRITE(*,*)
               WRITE(*,*) ' Proc_ws = ', Proc_ws
               WRITE(*,*)
               WRITE(*,*) ' NLx_ws = ', NLx_ws
               WRITE(*,*) ' NLy_ws = ', NLy_ws
               WRITE(*,*) ' NLz_ws = ', NLz_ws
               WRITE(*,*)
               WRITE(*,*) ' NLx_ws = ', NLx_ws
               WRITE(*,*) ' NLy_ws = ', NLy_ws
               WRITE(*,*) ' NLz_ws = ', NLz_ws
            ENDIF

            ProcI_ws = ProcCoord(1,Proc_ws)
            ProcJ_ws = ProcCoord(2,Proc_ws)
            ProcK_ws = ProcCoord(3,Proc_ws)

            !
            ! Copy the current block of WS into the
            ! correct block of the SERIAL stencil.
            !
            CALL BMG3_SymStd_COPY_cg_WS_SO(&
     &                WS(i_WS), NLx_ws-2, NLy_ws-2, NLz_ws-2, &
     &                BMG_rWORK_CS(p_SO), NGx, NGy, NGz,&
     &                i_WS, iGs_ws, jGs_ws, kGs_ws&
     &                )

            !
         ENDDO

         !
         !  Override defaults (.FALSE.) for debugging:
         !
         BMG_SER_IOFLAG(iBMG3_SER_BUG_STENCIL_FG)  = .FALSE.
         BMG_SER_IOFLAG(iBMG3_SER_BUG_STENCIL_CG)  = .FALSE.
         BMG_SER_IOFLAG(iBMG3_SER_BUG_STENCIL_CG1) = .FALSE.

         !
         ! Call SERIAL BoxMG to perform the setup
         !
         BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_SETUP-1) = BMG_SER_SETUP_only

         CALL BMG3_SER_SymStd_SOLVE_boxmg( &
     &             NGx-2, NGy-2, NGz-2,&
     &             BMG_iWORK_CS(p_iPARMS), BMG_rWORK_CS(p_rPARMS),&
     &             BMG_SER_IOFLAG,&
     &             BMG_rWORK_CS(p_U), BMG_rWORK_CS(p_Q),&
     &             BMG_rWORK_CS(p_RES), NF, NC, &
     &             BMG_rWORK_CS(p_SO), NSO,&
     &             BMG_rWORK_CS(p_SOR), NSOR,&
     &             BMG_rWORK_CS(p_CI), NCI,&
     &             BMG_rWORK_CS(p_CSO),&
     &             BMG_rWORK_CS(p_CU), NCBW, NCU,&
     &             BMG_iWORK_CS(p_iGRD), NOGm_SER, NOG_SER,&
     &             BMG_iWORK_CS(p_iWORK_PL), NBMG_SER_iWORK_PL,&
     &             BMG_rWORK_CS(p_rWORK_PL), NBMG_SER_rWORK_PL&
     &             )

         BMG_iWORK_CS(p_iPARMS+id_BMG3_SER_SETUP-1) = BMG_SER_SETUP_none

      ENDIF

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

 500  FORMAT (/,'FATAL ERROR: BMG3_SymStd_SOLVE_cg_boxmg.f',/,5X,A)
 510  FORMAT (5X,A,1X,I3)

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

      RETURN
      END

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar


text:
Feel like youve been ripped off, Food Poisoning poisoned, or offended by a business owner? Youre not alone and chances are there are other people whove had similar fates. You have plenty of recourse to prevent these acts from repeating, either to you or to anybody.

See also Tenant Rights, and the Restaurant Guide for an explanation of policies from a business owners point of view.

 First Course of Action
An unspoken problem is an unknown problem. The best method of making your unpleasant experience known is by speaking up! Sometimes businesses are not aware theyre doing something illegal or unethical. If an employee can not solve your problem ask for the manager. In restaurants, depending on the situation, you should not be shy about asking for discounts or compensation when  something major goes wrong. Remember, if you dont ask, only the most benevolent of businesses will automatically compensate you. If the kitchen or establishment screwed up, make sure to still tipping tip the waitstaff appropriately for the work they performed.

 Health Issues
Bring your problem to the attention of the management and use your judgment to determine whether to take further steps. Finding a violation could be an isolated incident. Alternatively, it could be a symptom of a more systematic problem with maintaining health standards. This is not just a matter of cleanliness, its a matter of personal safety.  Call the http://www.dcn.davis.ca.us/YoloLINK/programs/pEnviroHealthServic4626.html Yolo County Health Department at (530)6668646. Take some responsibility if you dont lodge a complaint about the cockroach in your Tofu surprise, chances are no one else has wanted to either.

The following constitute health code violations that should be dealt with immediately. These are assessed during a health inspection: vermin and/or droppings, noncommercial grade pesticides, toxic chemicals on top of or near food, food stored between 41 and 135 degrees, employees not washing hands, employees washing soiled hands in dishwashing sink, and employees with hand wounds not wearing gloves or finger condoms. Concerns that may potentially be isolated incidents include meat not cooked to appropriate temperatures, foreign objects in food (e.g., bandaids, hair, etc.), and spoiled food.

See also Food Poisoning

 Receipt Issues
Always check your bill to ensure that you are only paying for services rendered. Ask for clarification on illegible receipts and either mark or put a cross through tip and total fields.  Businesses such as Gas Stations frequently charge ATMs ATM fees, but the legality of such transactions for credit cards is unclear and warrants further documentation. Some http://www.thecreditreportsite.com/credit_cards.asp#surcharges sources say its illegal to charge fees for credit cards, but http://www.ucop.edu/finmgt/banking/ccpayguide.pdf others say that credit card companies allow convenience fees at automated locations, but most do not allow surcharges. Needless to say, the vendor must disclose their convenience charge up front. Contact your credit card directly if you think a violation has occurred, and please share your experience with the Wiki.

Temporary overcharging (actually an authorization hold) is done at many restaurants, especially the corporate chain restaurants (Applebees, Chilis, TGIFridays...) to ensure that once your tip is added onto a bill it wont get declined. So if you swipe your card and go online that evening it will show the amount the server initially swiped for. When the server cashes out at the end of their shift, the correct tip is entered into the computer, and that is the charge that get run rather than the amount held initially. This may take a few days to post depending on the bank.

 You were ripped off
File a complaint with the Sacramento Valley Better Business Bureau Better Business Bureau or with the California Department of http://www.dca.ca.gov/complainthelp/citizencomplaint.html consumer affairs. You might also consider taking the business to Small Claims Court.

 How to avoid getting ripped off

Do not be embarrassed to ask how much something costs. For example,  Upscale Restaurants upscale restaurants bank on you ordering specials that are listed without a price. They know youre more willing to pay an exorbitant price after youve eaten the meal.

If a server offers you something you did not directly ask for, ask if it will cost you extra. With few exceptions Do you want rice with that means the same thing as Do you want fries with that, i.e. Give us another $1.50.


subroutine boundary_condition_c
  use variables
  implicit none
  real(8) :: tend

  ! Uc(1,0:Ny,1)=0.25d0
  !x方向境界条件
  !上下流端水位一定
  ! Uc(1:1,2:Ny-1,1)=0.1d0
  ! Uc(1:1,2:Ny-1,4)=0.3d0
  ! Uc(1:3,2:Ny-1,2)=0.430887d0
  select case(boundtype)
  case(1)
    !x方向上下流端の壁でuh=0
    Uc(0,0:Ny+1,2)=0.d0
    Ec(0,0:Ny+1,1)=0.d0
    Ec(0,0:Ny+1,3)=0.d0
    Ec(0,0:Ny+1,4)=0.d0
    Fc(0,0:Ny+1,2)=0.d0
    Uc(Nx+1,0:Ny+1,2)=0.d0
    Ec(Nx+1,0:Ny+1,1)=0.d0
    Ec(Nx+1,0:Ny+1,3)=0.d0
    Ec(Nx+1,0:Ny+1,4)=0.d0
    Fc(Nx+1,0:Ny+1,2)=0.d0
    !x方向側壁でuvhの差が0
    Fc(0:Nx+1,0,2)=Fc(0:Nx+1,1,2)
    Fc(0:Nx+1,Ny+1,2)=Fc(0:Nx+1,Ny,2)
    Ec(0:Nx+1,0,3)=Ec(0:Nx+1,1,3)
    Ec(0:Nx+1,Ny+1,3)=Ec(0:Nx+1,Ny,3)
    !y方向上下流端の壁でvh=0
    Uc(0:Nx+1,0,3)=0.0d0
    Ec(0:Nx+1,0,3)=0.0d0
    Fc(0:Nx+1,0,1)=0.0d0
    Fc(0:Nx+1,0,2)=0.0d0
    Fc(0:Nx+1,0,4)=0.0d0
    Uc(0:Nx+1,Ny+1,3)=0.0d0
    Ec(0:Nx+1,Ny+1,3)=0.0d0
    Fc(0:Nx+1,Ny+1,1)=0.0d0
    Fc(0:Nx+1,Ny+1,2)=0.0d0
    Fc(0:Nx+1,Ny+1,4)=0.0d0
    !y方向側壁でuvhの差が0
    Fc(0,0:Ny,2)=Fc(1,0:Ny,2)
    Fc(Nx+1,0:Ny,2)=Fc(Nx,0:Ny,2)
    Ec(0,0:Ny,3)=Ec(1,0:Ny,3)
    Ec(Nx+1,0:Ny,3)=Ec(Nx,0:Ny,3)
  case(2)
    !x方向上流端の壁でh,c=const
    Uc(0,1:Ny,1)=0.1d0
    Uc(0,1:Ny,4)=0.0d0
    !x方向上流端の壁でuh=0
    Uc(0,0:Ny+1,2)=0.d0
    Ec(0,0:Ny+1,1)=0.d0
    Ec(0,0:Ny+1,3)=0.d0
    Ec(0,0:Ny+1,4)=0.d0
    Fc(0,0:Ny+1,2)=0.d0
    !x方向側壁でuvhの差が0
    Fc(0:Nx+1,0,2)=Fc(0:Nx+1,1,2)
    Fc(0:Nx+1,Ny+1,2)=Fc(0:Nx+1,Ny,2)
    Ec(0:Nx+1,0,3)=Ec(0:Nx+1,1,3)
    Ec(0:Nx+1,Ny+1,3)=Ec(0:Nx+1,Ny,3)
    !y方向側壁でuvhの差が0
    Fc(0,0:Ny+1,2)=Fc(1,0:Ny+1,2)
    Fc(Nx+1,0:Ny+1,2)=Fc(Nx,0:Ny+1,2)
    Ec(0,0:Ny+1,3)=Ec(1,0:Ny+1,3)
    Ec(Nx+1,0:Ny+1,3)=Ec(Nx,0:Ny+1,3)
    !y方向上下流端の壁でvh=0
    Uc(0:Nx+1,0,3)=0.0d0
    Ec(0:Nx+1,0,3)=0.0d0
    Fc(0:Nx+1,0,1)=0.0d0
    Fc(0:Nx+1,0,2)=0.0d0
    Fc(0:Nx+1,0,4)=0.0d0
    Uc(0:Nx+1,Ny+1,3)=0.0d0
    Ec(0:Nx+1,Ny+1,3)=0.0d0
    Fc(0:Nx+1,Ny+1,1)=0.0d0
    Fc(0:Nx+1,Ny+1,2)=0.0d0
    Fc(0:Nx+1,Ny+1,4)=0.0d0
  case(3)
    !x方向上流端の壁でh,c=const
    Uc(0:1,1:Ny,1)=hup
    Uc(0:1,1:Ny,2)=qup/Bini
    ! Uc(1,1:Ny,2)=Uc(2,1:Ny,2)
    Uc(0:1,1:Ny,3)=0.d0
    Uc(0,1:Ny,4)=hup*cup !土砂ありの場合
    Uc(0,1:Ny,5)=Zbini(0,1:Ny) !土砂ありの場合
    !下流端
    Uc(Nx:Nx+1,1:Ny,1)=hdown
    Uc(Nx:Nx+1,1:Ny,2)=qup/Bini
    Uc(Nx:Nx+1,1:Ny,3)=0.d0
    ! Uc(Nx+1,1:Ny,3)=hup*cup
    ! Uc(Nx+1,1:Ny,3)=Zbini(Nx+1,1:Ny)
    ! Uc(Nx+1,1:Ny,4)=hup*cup

    ! Uc(Nx,1:Ny,1)=Uc(Nx-1,1:Ny,1)
    ! Uc(Nx,1:Ny,2)=Uc(Nx-1,1:Ny,2)
    ! Uc(Nx+1,1:Ny,1)=Uc(Nx,1:Ny,1)
    ! Uc(Nx+1,1:Ny,2)=Uc(Nx,1:Ny,2)
    ! Uc(Nx+1,1:Ny,3)=0.d0
    ! Uc(Nx+1,1:Ny,4)=hup*cup

    ! Uc(Nx+1,1:Ny,1)=hup
    ! Uc(Nx+1,1:Ny,2)=qup/Bini
    ! Uc(Nx+1,1:Ny,3)=0.d0

    !両岸の境界条件は，差分時に使う？
    ! Uc(0,0:Ny+1,4)=hup*cup !土砂ありの場合
    ! Ec(0,0:Ny+1,2)=(Qup/Bini/hup)**2*hup+0.5d0*g*hup
    ! Fc(0,0:Ny+1,3)=0.5d0*g*hup
    !x方向上流端の壁でuh=0
    ! Ec(0,0:Ny+1,1)=qup/Bini
    ! Ec(0,0:Ny+1,3)=0.d0
    ! Ec(0,0:Ny+1,4)=0.d0
    ! Fc(0,0:Ny+1,2)=0.d0
    !x方向側壁でuvhの差が0
    ! Fc(0:Nx+1,0,2)=Fc(0:Nx+1,1,2)
    ! Fc(0:Nx+1,Ny+1,2)=Fc(0:Nx+1,Ny,2)
    ! Ec(0:Nx+1,0,3)=Ec(0:Nx+1,1,3)
    ! Ec(0:Nx+1,Ny+1,3)=Ec(0:Nx+1,Ny,3)
    !y方向側壁でuvhの差が0
    ! Fc(0,0:Ny+1,2)=Fc(1,0:Ny+1,2)
    ! Fc(Nx+1,0:Ny+1,2)=Fc(Nx,0:Ny+1,2)
    ! Ec(0,0:Ny+1,3)=Ec(1,0:Ny+1,3)
    ! Ec(Nx+1,0:Ny+1,3)=Ec(Nx,0:Ny+1,3)
    !y方向上下流端の壁でvh=0
    ! Uc(0:Nx+1,0,3)=0.0d0
    ! Ec(0:Nx+1,0,3)=0.0d0
    ! Fc(0:Nx+1,0,1)=0.0d0
    ! Fc(0:Nx+1,0,2)=0.0d0
    ! Fc(0:Nx+1,0,4)=0.0d0
    ! Uc(0:Nx+1,Ny+1,3)=0.0d0
    ! Ec(0:Nx+1,Ny+1,3)=0.0d0
    ! Fc(0:Nx+1,Ny+1,1)=0.0d0
    ! Fc(0:Nx+1,Ny+1,2)=0.0d0
    ! Fc(0:Nx+1,Ny+1,4)=0.0d0
    !x方向上流端の壁でh,c=const,uh=0
    ! Uc(0,0:Ny+1,1)=hup
    ! Uc(0,0:Ny+1,2)=qup/Bini
    ! Uc(0,0:Ny+1,3)=0.0d0
    ! Uc(0,0:Ny+1,4)=hup*cup !土砂ありの場合
    ! Uc(0,0:Ny+1,5)=Uc(0,0:Ny+1,5)
    ! Ec(0,0:Ny+1,1)=Qup/Bini
    ! Ec(0,0:Ny+1,2)=(Qup/Bini/hup)**2*hup+0.5d0*g*hup
    ! Ec(0,0:Ny+1,3)=0.d0
    ! Ec(0,0:Ny+1,4)=Qup/Bini*Cup
    ! Ec(0,0:Ny+1,5)=0.d0
    ! Fc(0,0:Ny+1,1)=0.d0
    ! Fc(0,0:Ny+1,2)=0.d0
    ! Fc(0,0:Ny+1,3)=0.5d0*g*hup
    ! Fc(0,0:Ny+1,4)=0.d0
    ! Fc(0,0:Ny+1,5)=0.d0
    ! !x方向側壁(左岸)でuvhの差が0
    ! Fc(0:Nx+1,0,2)=Fc(0:Nx+1,1,2)
    ! Fc(0:Nx+1,Ny+1,2)=Fc(0:Nx+1,Ny,2)
    ! Ec(0:Nx+1,0,3)=Ec(0:Nx+1,1,3)
    ! Ec(0:Nx+1,Ny+1,3)=Ec(0:Nx+1,Ny,3)
    ! !y方向側壁でuvhの差が0
    ! Fc(0,0:Ny+1,2)=Fc(1,0:Ny+1,2)
    ! Fc(Nx+1,0:Ny+1,2)=Fc(Nx,0:Ny+1,2)
    ! Ec(0,0:Ny+1,3)=Ec(1,0:Ny+1,3)
    ! Ec(Nx+1,0:Ny+1,3)=Ec(Nx,0:Ny+1,3)
    ! !y方向上下流端の壁でvh=0
    ! Uc(0:Nx+1,0,3)=0.0d0
    ! Ec(0:Nx+1,0,3)=0.0d0
    ! Fc(0:Nx+1,0,1)=0.0d0
    ! Fc(0:Nx+1,0,2)=0.0d0
    ! Fc(0:Nx+1,0,4)=0.0d0
    ! Uc(0:Nx+1,Ny+1,3)=0.0d0
    ! Ec(0:Nx+1,Ny+1,3)=0.0d0
    ! Fc(0:Nx+1,Ny+1,1)=0.0d0
    ! Fc(0:Nx+1,Ny+1,2)=0.0d0
    ! Fc(0:Nx+1,Ny+1,4)=0.0d0
  case(4)
    !x方向上流端の壁でh,c=const
    Uc(0,1:Ny,1)=Uc(1,1:Ny,1) !hup
    Uc(0,1:Ny,2)=Uc(1,1:Ny,2) !qup/Bini
    Uc(0:1,1:Ny,3)=0.d0 !0.d0
    Uc(0,1:Ny,4)=Uc(1,1:Ny,4) !hup*cup !土砂ありの場合
    Uc(0:1,1:Ny,5)=Zbini(0:1,1:Ny) !土砂ありの場合
    !下流端
    Uc(Nx+1,1:Ny,1:5)=Uc(Nx,1:Ny,1:5)
  case(5)
    !x方向上流端の壁でh,c=const
    Tend=20.d0
    Uc(0,1:Ny,1:5)=Uc(1,1:Ny,1:5)
    ! Uc(0,1:Ny,1)=0.d0 !Uc(1,1:Ny,1) !hup
    ! Uc(0,1:Ny,2)=Uc(0,1:Ny,2) !qup/Bini
    ! Uc(0,1:Ny,3)=0.d0 !0.d0
    !下流端
    Uc(Nx+1,1:Ny,1:5)=Uc(Nx,1:Ny,1:5)
    if(T<Tend)then
      Uc(2:3,46:50,1)=hup !0.04d0
      Uc(2:3,46:50,2)=qup !1.d0/1000.d0
      Uc(2:3,46:50,3)=0.d0
      Uc(2:3,46:50,4)=hup*cup !0.04d0*0.3d0
      Uc(2:3,46:50,5)=Zbini(2:3,46:50)
    else
      Uc(2:3,46:50,1)=0.d0 !0.04d0
      Uc(2:3,46:50,2)=0.d0 !1.d0/1000.d0
      Uc(2:3,46:50,3)=0.d0
      Uc(2:3,46:50,4)=0.d0 !0.04d0*0.3d0
      Uc(2:3,46:50,5)=Zbini(2:3,46:50)
    endif
  case(6)
    Uc(0,1:Ny,1:5)=Uc(1,1:Ny,1:5)
    Uc(Nx+1,1:Ny,1:5)=Uc(Nx,1:Ny,1:5)
    Uc(1:Nx,0,1:5)=Uc(1:Nx,1,1:5)
    Uc(1:Nx,Ny+1,1:5)=Uc(1:Nx,Ny,1:5)
  case(7)
    if(rank_l(my_rank)<0)Uc(0,1:Ny,1:5)=Uc(1,1:Ny,1:5)
    if(rank_r(my_rank)<0)Uc(Nx+1,1:Ny,1:5)=Uc(Nx,1:Ny,1:5)
    if(rank_u(my_rank)<0)Uc(1:Nx,0,1:5)=Uc(1:Nx,1,1:5)
    if(rank_d(my_rank)<0)Uc(1:Nx,Ny+1,1:5)=Uc(1:Nx,Ny,1:5)
  end select

end subroutine boundary_condition_c