text: array A = [3] U32 format "{} units"

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:
use plantfem

implicit none

type(Soybean_) :: soy
integer(int32) :: i
real(real64)   :: length, coord(3)

! Initialize soybean
call soy%init(config="Tutorial/obj/realSoybeanConfig.json")

! change node-length by its height
do i=1, size(soy%stem)
    if( soy%stem(i)%femdomain%mesh%empty()) exit
    coord = soy%stem(i)%getCoordinate("A")
    length = coord(3)*0.20d0 + 0.010d0
    print *, length
    call soy%stem(i)%grow(length=length)
enddo

! update geoetry
call soy%update()

! move soy
call soy%move(x=0.0d0)

! export as vtk
call soy%vtk(name="soy3")

end

max_stars_repo_name: kazulagi/plantFEM


text: locate type Fpp.ToCpp.Phases at "ToCpp.fpp"


text:
C     LAST UPDATE 17/02/97
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE FLIEB(Q,RN,LL,F,CHI)
      IMPLICIT NONE
C
C Purpose: This subroutine calculates the LIEBLER DIBLOCK COPOLYMER
C          L.LEIBLER, Macromolecule 13(1980)1602-1617 eqns IV-2 to IV-8
C          R.K.Heenan 13/1/97
C          Modified for FIT 17/02/97 by RCD.
C          Q=X, DIBLK=Y, A(1)=
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Parameters:
C

c==== NOTE MAY NOT BE VERY ROBUST TO OVERFLOWS, NON-PHYSICAL PARAMS ETC. !!!
C
      real*4 LL
      LT=1
      UU=0.0
C====  compute Rg**2
      RR = RN*LL*LL/6.0
      BIGF = DEBGAUSS(LT,Q,RR,UU)
      FM = 1.0-F
      D  = (F**2)* DEBGAUSS(LT,Q,RR*F ,UU)
      DM = (FM**2)*DEBGAUSS(LT,Q,RR*FM,UU)
      DEN = D*DM -0.25*(BIGF -D -DM)**2
      DIBLK = RN/( BIGF/DEN  -2.0*CHI*RN)
      RETURN
      END
c
c
      function DEBGAUSS(LTYP,QQ,RR,UU)
C==== P(Q) for Debye Gaussian coil.  Polydisp if LTYP >10 (i.e. 11 or 21)
c
C==== NOTE input RR is RG SQUARED not RG, but QQ is just Q
c
c==== carefully modified Feb 1992 by R.K.Heenan to avoid over & underflows
C==== 12/5/94 star polymer equation added as LTYP=31, fo rwhich UU is then
c==== the number of arms F
c
c==== polynomial expansion is used at small QR
c
c==== Written by R.K.Heenan, RAL
c
      QR=RR*(QQ**2)
      U=0.0
      FSTAR=1.0
      IF(LTYP.EQ.11)U=ABS(UU)
      IF(LTYP.EQ.31)FSTAR=AMAX1(1.0,ABS(UU))
C====  note U and FSTAR cannot be used simultaneously !
      Y=QR/( (1.0+2.0*U)*(3.0 -2.0/FSTAR) )
C==== polydisp version
      IF(U.GT.1.0E-02.AND.Y.GT.0.01)DEBGAUSS=
     >         2.0*( (1.0+U*Y)**(-1.0/U) +Y -1.0)/( (1.0+U)*(Y**2) )
C==== monodisp version
C2345 789012345678901234567890123456789012345678901234567890123456789012
      IF(U.LE.1.0E-02.AND.Y.GT.0.01)THEN
        AA=(EXPSPEC(-Y)-1.0)/(FSTAR*Y)
        DEBGAUSS= 2.*(AA + 1./FSTAR)/Y
        IF(FSTAR.GT.1.0)DEBGAUSS=DEBGAUSS+ FSTAR*(FSTAR-1.0)*AA*AA
      END IF
C==== EXPSPEC avoids underflows as QR**2 can be large
C==== low Q expansion for ALL versions
      IF(Y.LE.0.01)DEBGAUSS= 1.0 - QR*(1.0+ (0.25 - 0.75*U)*Y )/3.0
      return
      end
c
      FUNCTION EXPSPEC(A)
C==== expspec avoids underflows invoking the error handler
c==== the -88.5 is a machine dependent number
c
c==== The usual EXP( ) function may be fine - depends whether
c==== a particular compiler gives underflow errors
cc
c==== Written by R.K.Heenan, RAL
c
      EXPSPEC=0.0
      IF(A.GT.-88.5.AND.A.LT.88.0)THEN
         EXPSPEC=EXP(A)
      ELSE IF(A.GE.85.2)THEN
         EXPSPEC=1.0E37
      END IF
C==== NOTE EXP(88.02)=1.684996E38
C==== NOTE EXP(-88.7)=3.006635E-39
      RETURN
      END
C

ext: f


text:
        subroutine stat_3d(s,nx_s,nx_e,ny_s,ny_e,nz_s,nz_e)!,field_name)

!  subroutine to stat a  3-d array with ferret

!       character field_name(20)

        integer nx_s,nx_e,ny_s,ny_e,nz_s,nz_e
        integer i,j,k

        real s(nx_s:nx_e,ny_s:ny_e,nz_s:nz_e), s_ave, sum , std, missing_value
        integer number

        missing_value = -1.0E+34

          sum = 0.0d0
          number= 0

       do k = nz_s,nz_e
       do j = ny_s,ny_e
       do i = nx_s,nx_e

          if ( s(i,j,k) .gt. 0.98 * missing_value) then

             sum = s(i,j,k) + sum
             number= number+1

          endif

       enddo;enddo;enddo

!         write(6,*) 'in stat', sum

          s_ave = 0.0d0
          if( number .gt. 0) s_ave = sum / real(number)
          sum = 0.0d0

       do k = nz_s, nz_e
       do j = ny_s, ny_e
       do i = nx_s, nx_e

          if ( s(i,j,k) .gt. 0.99* missing_value)  then
             sum = (s(i,j,k)-s_ave)**2 + sum
          endif

       enddo;enddo;enddo

          std = 0.0d0
          if( number .gt. 0.0) std = sqrt( sum /  real(number) )

!         write (6,*) 'number of bad values;', (nx_e-nx_s+1)*(ny_e-ny_s+1)*(nz_e-nz_s+1)-number

          write(6,'(2d18.10,i10)')   s_ave, std , (nx_e-nx_s+1)*(ny_e-ny_s+1)*(nz_e-nz_s+1)-number!, 's_ave, std, # bad:',

        return
        end

max_stars_repo_name: diioriod/PlumeModel-Hotcross

max_issues_repo_name: diioriod/PlumeModel-Hotcross

max_forks_repo_name: diioriod/PlumeModel-Hotcross


text:
! $Id$
!
! testing Xpomp display facilities in Fortran 77 / HPF
!
! (c) Ronan.Keryell@cri.ensmp.fr 1996
!
      program fractal

      implicit none

! XPOMP library interface:
      include 'xpomp_graphic_F.h'

! also tells that a file of stubs is needed:
!fcd$ stubs xpomp_stubs.f

! Algorithm parameters:
      integer n_iteration
      parameter (n_iteration = 256)
      real*8 big_value
      parameter (big_value = 100)
! distorsion = 2 to have normal Mandelbrot:
      real*8 distorsion
      parameter (distorsion = 2)

! Size of the iteration space:
      integer x_size, y_size
      parameter(x_size = 800)
      Parameter(y_size = 800)

! The zooming ratio to display this iteration space:
      integer x_display_zoom, y_display_zoom
      parameter(x_display_zoom = 1)
      parameter(y_display_zoom = 1)

! The zooming factor used when zooming with the button:
      real*8 zooming_factor
      parameter (zooming_factor = 4)

      integer x_display_size, y_display_size
      parameter(x_display_size = x_size*x_display_zoom)
      parameter(y_display_size = y_size*y_display_zoom)

      integer counter
      integer k, x, y
      integer display
      character image(0:x_size-1, 0:y_size-1)

      real*8 zoom, xcenter, ycenter
      real*8 zr, zrp, zi, cr, ci, d
      integer status, button, state
      integer X0, Y0, X1, Y1

! Some HPF distributions:
! cyclic would make more sense as far as load balancing is concerned
! However HPFC would not be very good at it...
! No communication -> block distribution on the second dimension
!   for best cache effects
!
!hpf$ processors pe(number_of_processors())
!hpf$ distribute image(*,block) onto pe

! Initialize XPOMP
      call xpompf open display(x_display_size, y_display_size, display)
      call xpompf set color map(display, 1, 1, 0, 0, status)

      call xpompf show usage

      print *, 'Mouse button 1 to zoom in, 2 to recenter, 3 to zoom out'
      print *, 'A mouse + Shift: restart.'

! Initial position
      xcenter = 0
      ycenter = 0
      zoom = 5

! Main loop:
      counter = 0
 10   continue

!fcd$ time
! Compute a fractal image:
!hpf$ independent, new(cr,ci,zr,zi,d,zrp,k)
      do y = 0, y_size - 1
         ci = ycenter + (y - y_size/2)*zoom/y_size
!hpf$    independent
         do x = 0, x_size - 1
            cr = xcenter + (x - x_size/2)*zoom/x_size
            zr = cr
            zi = ci
            do k = 0, n_iteration - 1
               d = zr*zr + zi*zi
               zrp = zr*zr - zi*zi + cr
               zi = distorsion*zr*zi + ci
               zr = zrp
               if (d .gt. big_value) goto 300
            enddo
 300        continue
            image(x, y) = CHAR(k)
         enddo
      enddo

! Display the image:
      call xpompf flash(display, image, x_size, y_size,
     &     0, 0, x_display_zoom, y_display_zoom, status)

!fcd$ end time ('Computation of one image')

! Wait for user interaction:
      call xpompf wait mouse(display, x, y, state, button)

! Recenter
      xcenter = xcenter + (x/x_display_zoom - x_size/2)*zoom/x_size
      ycenter = ycenter + (y/y_display_zoom - y_size/2)*zoom/y_size

      if (state .eq. 1) then

! A button with shift -> restart from the begining:
       xcenter = 0
       ycenter = 0
       zoom = 5
      else if (button .eq. 1) then

! Zoom in
         zoom = zoom/zooming_factor
! Display the frame that will be zoomed in:
         X0 = x - x_display_size/zooming_factor/2
         Y0 = y - y_display_size/zooming_factor/2
         X1 = x + x_display_size/zooming_factor/2
         Y1 = y + y_display_size/zooming_factor/2
         call xpompf draw frame(display, ' ', 255, -1,
     &        X0, Y0, X1, Y1, -181, status)
      else if (button .eq .3) then

! Zoom out
         zoom = zoom*zooming_factor
      endif

      print *, 'Position (', xcenter, ',', ycenter, '), zoom =',
     &     zoom, ' button =', button, ', state =', state

! hpfc does not like infinite loop without an exit point...
! some bug related to the management of unstructured to be investigated
! I added this counter and exit as a temporary fix - FC.
      counter = counter+1
      if (counter.gt.100000) goto 20
      goto 10
 20   continue
      end

! end of it


text:
      real function gcum(x)
c-----------------------------------------------------------------------
c
c Evaluate the standard normal cdf given a normal deviate x.  gcum is
c the area under a unit normal curve to the left of x.  The results are
c accurate only to about 5 decimal places.
c
c
c-----------------------------------------------------------------------

      implicit integer (i-n)
      implicit double precision (a-h,o-z)

      z = x
      if(z.lt.0.) z = -z
      t    = 1./(1.+ 0.2316419*z)
      gcum = t*(0.31938153   + t*(-0.356563782 + t*(1.781477937 +
     +       t*(-1.821255978 + t*1.330274429))))
      e2   = 0.
c
c  6 standard deviations out gets treated as infinity:
c
      if(z.le.6.) e2 = exp(-z*z/2.)*0.3989422803
      gcum = 1.0- e2 * gcum
      if(x.ge.0.) return
      gcum = 1.0 - gcum
      return
      end

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


text:
! This is free and unencumbered software released into the public domain.
! SPDX-Identifier: Unlicense

!> Provides interfaces and implementations for hashing algorithms.
module reference_waterhash
   use, intrinsic :: iso_c_binding, only : c_loc, c_long
   implicit none
   private

   public :: c_waterhash

   !> Hash function to obtain a hashed value of the provided variable
   interface c_waterhash
      !> Interface to the waterhash implementation in C
      pure function c_waterhash(ptr, len, seed) result(val) bind(C)
         use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long
         type(c_ptr), value, intent(in) :: ptr
         integer(c_int), value, intent(in) :: len
         integer(c_long), value, intent(in) :: seed
         integer(c_int) :: val
      end function c_waterhash
      module procedure :: c_hash_character
   end interface c_waterhash

contains

   !> Hash a character value
   pure function c_hash_character(var, seed) result(val)
      !> Variable to hash
      character(len=*), target, intent(in) :: var
      !> Seed for the hashing function
      integer(c_long), intent(in) :: seed
      !> Hash value
      integer :: val

      val = c_waterhash(c_loc(var), len(var), seed)

   end function c_hash_character

end module reference_waterhash

max_stars_repo_name: awvwgk/waterhash

max_issues_repo_name: awvwgk/waterhash

max_forks_repo_name: awvwgk/waterhash


text:
      SUBROUTINE RWHOUT
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CSTVAL/ TSTART,ZL0(10),ZL1(10),ZL2(10),ZLA(10),TTST
      COMMON/VARBLS/ DEPEND(150),DERIV(150)
      WRITE(50)TTST,(DEPEND(I),I=1,15),(DERIV(I),I=1,15)
      RETURN
      END

ext: for


text:
!==============================================================================!
  subroutine Grid_Mod_Create_Levels(grid)
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
!-----------------------------------[Locals]-----------------------------------!
  integer :: lev, n_cells, n_faces, c, s
!==============================================================================!

  !--------------------------------!
  !                                !
  !   Allocate memory for levels   !
  !                                !
  !--------------------------------!

  ! At each level, allocate memory for individual data
  do lev = grid % n_levels + 1, 1, -1

    ! Tentative number of cells and faces
    n_cells = grid % level(lev) % n_cells
    n_faces = grid % n_faces  ! could be cells times a factor

    if(lev .eq. 1) then
      n_cells = grid % n_cells
      n_faces = grid % n_faces
    end if

    ! Connectivity with finest grid
    allocate(grid % level(lev) % cell(grid % n_cells))
    allocate(grid % level(lev) % face(grid % n_faces))

    ! Intra-level connectivity and geometry
    allocate(grid % level(lev) % faces_c(2, n_faces))
    allocate(grid % level(lev) % n_finest_cells(n_cells))
    allocate(grid % level(lev) % xc     (   n_cells))
    allocate(grid % level(lev) % yc     (   n_cells))
    allocate(grid % level(lev) % zc     (   n_cells))
  end do

  !---------------------------------------------------------!
  !                                                         !
  !   Initialize arrays for level zero and coarser levels   !
  !                                                         !
  !---------------------------------------------------------!

  !------------------------------------!
  !     Cells and faces on level 1     !
  !   (Level 1 is the original grid)   !
  !------------------------------------!

  ! Cells
  grid % level(1) % n_cells = grid % n_cells
  do c = 1, grid % n_cells
    grid % level(1) % cell(c) = c
  end do

  ! Faces
  grid % level(1) % n_faces = grid % n_faces
  do s = 1, grid % n_faces
    grid % level(1) % face(s) = s
    grid % level(1) % faces_c(1, s) = grid % faces_c(1, s)
    grid % level(1) % faces_c(2, s) = grid % faces_c(2, s)
  end do

  !-------------------------------------!
  !   Cells and faces on other levels   !
  !-------------------------------------!
  do lev = 2, grid % n_levels + 1

    ! Cells
    do c = 1, grid % n_cells
      grid % level(lev) % cell(c) = 1
    end do

    ! Faces
    do s = 1, grid % n_faces
      grid % level(lev) % face(s) = 0
    end do
  end do

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
c     Subroutine to export month station data
c     AJ_Kettle, Dec2/2017

      SUBROUTINE export_tab_sn3a_stndata(s_directory_output,s_filename1,
     +  s_date,
     +  l_sd_stnrecord,s_sd_stke,s_sd_stid,
     +  l_zip_use,s_stnnum,
     +  s_arch_mon_commonyear_st,s_arch_mon_commonmon_st,
     +  s_arch_mon_commonyear_en,s_arch_mon_commonmon_en,
     +  s_arch_geog_hgt_m,s_arch_geog_lat_deg,s_arch_geog_lon_deg,
     +  i_arch_lmonrec_common,
     +  f_arch_mon_airt_ngood,f_arch_mon_wspd_ngood,
     +  f_arch_mon_slpres_ngood,f_arch_mon_ppt_ngood)

      IMPLICIT NONE
c************************************************************************
c     Declare variables

      CHARACTER(LEN=200)  :: s_directory_output
      CHARACTER(LEN=200)  :: s_filename1
      CHARACTER(LEN=8)    :: s_date

      INTEGER             :: l_sd_stnrecord
      CHARACTER(LEN=5)    :: s_sd_stke(100)
      CHARACTER(LEN=5)    :: s_sd_stid(100)

      INTEGER             :: l_zip_use

      CHARACTER(LEN=5)    :: s_stnnum(1100)

      CHARACTER(LEN=2)    :: s_arch_mon_commonmon_st(1100)
      CHARACTER(LEN=2)    :: s_arch_mon_commonmon_en(1100)
      CHARACTER(LEN=4)    :: s_arch_mon_commonyear_st(1100)
      CHARACTER(LEN=4)    :: s_arch_mon_commonyear_en(1100)

      CHARACTER(LEN=8)    :: s_arch_geog_hgt_m(1100)
      CHARACTER(LEN=8)    :: s_arch_geog_lat_deg(1100)
      CHARACTER(LEN=8)    :: s_arch_geog_lon_deg(1100)

      INTEGER             :: i_arch_lmonrec_common(1100)
      REAL                :: f_arch_mon_airt_ngood(1100)
      REAL                :: f_arch_mon_wspd_ngood(1100)
      REAL                :: f_arch_mon_slpres_ngood(1100)
      REAL                :: f_arch_mon_ppt_ngood(1100)

      CHARACTER(LEN=4)    :: s_year_st
      CHARACTER(LEN=2)    :: s_month_st
      CHARACTER(LEN=4)    :: s_year_en
      CHARACTER(LEN=2)    :: s_month_en

      CHARACTER(LEN=18)   :: s_datestamp_st
      CHARACTER(LEN=18)   :: s_datestamp_en

c      CHARACTER(LEN=4)    :: s_code_relh          !38
      CHARACTER(LEN=4)    :: s_code_ppt           !44
c      CHARACTER(LEN=4)    :: s_code_snowcover     !53
c      CHARACTER(LEN=4)    :: s_code_snow_daytot   !55
      CHARACTER(LEN=4)    :: s_code_slpres          !58
      CHARACTER(LEN=4)    :: s_code_airtk          !85
c      CHARACTER(LEN=4)    :: s_code_wdir          !106
      CHARACTER(LEN=4)    :: s_code_wspd          !107

      CHARACTER(LEN=21)   :: s_code_allvar

      CHARACTER(LEN=300)  :: s_pathandname1

      INTEGER             :: i,j,k,ii,jj,kk
c************************************************************************
      print*,'just entered export_tab_sn3a_stndata'

c     Open file for export
      s_pathandname1=TRIM(s_directory_output)//TRIM(s_filename1)

      OPEN(UNIT=1,
     +  FILE=TRIM(s_pathandname1),
     +  FORM='formatted',
     +  STATUS='NEW',ACTION='WRITE')

      WRITE(UNIT=1,FMT=3009)
     + 'DWD - daily station information; 81 stations in record      '

c     SPACE
      WRITE(UNIT=1,FMT=3009)
     + '                                                            '

      WRITE(unit=1,FMT=3029) 'AJ_Kettle ',s_date
3029  FORMAT(t1,a10,t12,a8)

c     SPACE
      WRITE(UNIT=1,FMT=3009)
     + '                                                            '

      WRITE(UNIT=1,FMT=3009)
     + 'Main program: dd1.f                                         '
      WRITE(UNIT=1,FMT=3009)
     + 'Subroutine:   export_tab_sn3a_stndata.f                     '
3009  FORMAT(a60)

c     SPACE
      WRITE(UNIT=1,FMT=3009)
     + '                                                            '

c     SPACE
      WRITE(UNIT=1,FMT=3009)
     + 'Notes:                                                      '
      WRITE(UNIT=1,FMT=3009)
     + 'A:   no formal information in daily files on aut/man status '
      WRITE(UNIT=1,FMT=3009)
     + 'B:   44=ppt; 58=sea level pressure; 85=air temperature;     '
      WRITE(UNIT=1,FMT=3009)
     + '     107=wind speed                                         '
      WRITE(UNIT=1,FMT=3009)
     + 'C:   Deutscher Wetterdienst (DWD) & precursors              '
      WRITE(UNIT=1,FMT=3009)
     + 'D:   Walter Koelschtzky                                     '
      WRITE(UNIT=1,FMT=3009)
     + '     Head of Central Sales Management                       '
      WRITE(UNIT=1,FMT=3009)
     + '     Customer Relations Management                          '
      WRITE(UNIT=1,FMT=3009)
     + '     Deutscher Wetterdienst                                 '
      WRITE(UNIT=1,FMT=3009)
     + '     Business Unit Climate and Environment                  '
      WRITE(UNIT=1,FMT=3009)
     + '     Postbox 10 04 65                                       '
      WRITE(UNIT=1,FMT=3009)
     + '     D-63004 Offenbach, Germany                             '
      WRITE(UNIT=1,FMT=3009)
     + '     Phone: +49 (0)69 8062-4417                             '
      WRITE(UNIT=1,FMT=3009)
     + '     Fax:   +49 (0)69 8062-4499                             '
      WRITE(UNIT=1,FMT=3009)
     + '     Email: walter.koelschtzky@dwd.de                       '

c     SPACE
      WRITE(UNIT=1,FMT=3009)
     + '                                                            '

c     table headers
      WRITE(UNIT=1,FMT=3011)
     + 'Index','Orig_','alt_asl  ','Latitude  ','Longitude ','WMO  ',
     + 'N_mon ',
     + 'Report time stamp start     ','Report time stamp end       ',
     + 'Man','Observed_variables        ','Ope','Con'
      WRITE(UNIT=1,FMT=3011)
     + '     ','ID   ','         ','          ','          ','ID   ',
     + '      ',
     + '(monthly:yyyy-mm)','(monthly:yyyy-mm)',
     + 'Aut','see Note B                ','Ins','Add'
      WRITE(UNIT=1,FMT=3011)
     + '+---+','+---+','+-------+','+--------+','+--------+','+---+',
     + '+----+',
     + '+---------------+','+---------------+',
     + '+-+','+------------------------+','+-+','+-+'
3011  FORMAT(t1,a5,t7,a5,t13,a9,t23,a10,t34,a10,t45,a5,
     +  t51,a6,
     +  t58,a17,t76,a17,
     +  t94,a3,t98,a14,t113,a3,t117,a3)

c     Cycle through 81 subdaily records
      DO i=1,l_sd_stnrecord
c      Cycle through long daily file list
       DO j=1,l_zip_use

c       Output data if there is match
        IF (s_sd_stid(i).EQ.s_stnnum(j)) THEN

         s_year_st =s_arch_mon_commonyear_st(j)
         s_month_st=s_arch_mon_commonmon_st(j)
         s_year_en =s_arch_mon_commonyear_en(j)
         s_month_en=s_arch_mon_commonmon_en(j)
         s_datestamp_st='(Monthly:'//s_year_st//'-'//s_month_st//')'
         s_datestamp_en='(Monthly:'//s_year_en//'-'//s_month_en//')'
c****
c        Initialize all codes to null values
c         s_code_relh       =''   !38
         s_code_ppt        =''   !44
c         s_code_snowcover  =''   !53
c         s_code_snow_daytot=''   !55
         s_code_slpres     =''   !58
         s_code_airtk      =''   !85
c         s_code_wdir       =''   !106
         s_code_wspd       =''   !107

         IF (f_arch_mon_ppt_ngood(j).GT.0.0) THEN
          s_code_ppt        ='44,'
         ENDIF
         IF (f_arch_mon_slpres_ngood(j).GT.0.0) THEN
          s_code_slpres     ='58,'
         ENDIF
         IF (f_arch_mon_airt_ngood(j).GT.0.0) THEN
          s_code_airtk      ='85,'
         ENDIF
         IF (f_arch_mon_wspd_ngood(j).GT.0.0) THEN
          s_code_wspd       ='107,'
         ENDIF

         s_code_allvar=
     +    TRIM(s_code_ppt)//TRIM(s_code_slpres)//
     +    TRIM(s_code_airtk)//TRIM(s_code_wspd)
c****
       WRITE(UNIT=1,FMT=3013)
     +  i,s_sd_stid(i),s_arch_geog_hgt_m(j),
     +  s_arch_geog_lat_deg(j),s_arch_geog_lon_deg(j),
     +  s_sd_stke(i),
     +  i_arch_lmonrec_common(j),
     +  s_datestamp_st,s_datestamp_en,
     +  'A  ',s_code_allvar,'C  ','D  '
3013   FORMAT(t1,i2,t7,a5,t13,a9,t23,a10,t34,a10,t45,a5,
     +  t51,i6,
     +  t58,a18,t76,a18,
     +  t96,a3,t98,a14,t113,a3,t117,a3)
c****

        ENDIF
       ENDDO
      ENDDO

      CLOSE(UNIT=1)
c************************************************************************
      print*,'just leaving export_tab_sn3a_stndata'

      RETURN
      END

max_stars_repo_name: ajkettle/glamod-nuim

max_issues_repo_name: ajkettle/glamod-nuim

max_forks_repo_name: ajkettle/glamod-nuim


!
! CRTM_AntennaCorrection
!
! Module containgin routines to apply the antenna correction to the
! RTSolution brightness temperatures.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 20-Jul-2007
!                       paul.vandelst@noaa.gov
!

MODULE CRTM_AntennaCorrection

  ! ------------------
  ! Environment set up
  ! ------------------
  ! Module use statements
  USE Type_Kinds              , ONLY: fp
  USE CRTM_Parameters         , ONLY: TSPACE
  USE CRTM_GeometryInfo_Define, ONLY: CRTM_GeometryInfo_type, &
                                      CRTM_GeometryInfo_GetValue
  USE CRTM_RTSolution_Define  , ONLY: CRTM_RTSolution_type
  USE CRTM_SpcCoeff           , ONLY: SC
  ! Disable all implicit typing
  IMPLICIT NONE

  ! --------------------
  ! Default visibilities
  ! --------------------
  ! Everything private by default
  PRIVATE
  ! Module procedures
  PUBLIC :: CRTM_Compute_AntCorr
  PUBLIC :: CRTM_Compute_AntCorr_TL
  PUBLIC :: CRTM_Compute_AntCorr_AD

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*),  PARAMETER :: MODULE_VERSION_ID = &
  '$Id$'

CONTAINS

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

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Compute_AntCorr
!
! PURPOSE:
!       Subroutine to compute the antenna correction and apply it the RTSolution
!       brightness temperatures
!
! CALLING SEQUENCE:
!       CALL CRTM_Compute_AntCorr( GeometryInfo, &
!                                  SensorIndex , &
!                                  ChannelIndex, &
!                                  RTSolution    )
!
! INPUTS:
!       GeometryInfo:    Structure containing the view geometry information.
!                        UNITS:      N/A
!                        TYPE:       CRTM_GeometryInfo_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       SensorIndex:     Sensor index id. This is a unique index associated
!                        with a (supported) sensor used to access the
!                        shared coefficient data for a particular sensor.
!                        See the ChannelIndex argument.
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       ChannelIndex:    Channel index id. This is a unique index associated
!                        with a (supported) sensor channel used to access the
!                        shared coefficient data for a particular sensor's
!                        channel.
!                        See the SensorIndex argument.
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       RTSolution:      Structure containing the radiative transfer solution.
!                        UNITS:      N/A
!                        TYPE:       CRTM_RTSolution_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN OUT)
!:sdoc-:
!------------------------------------------------------------------------------

  SUBROUTINE CRTM_Compute_AntCorr( gI, &  ! Input
                                   n , &  ! Input
                                   l , &  ! Input
                                   RT  )  ! In/Output
    ! Arguments
    TYPE(CRTM_GeometryInfo_type), INTENT(IN)     :: gI
    INTEGER                     , INTENT(IN)     :: n  ! SensorIndex
    INTEGER                     , INTENT(IN)     :: l  ! ChannelIndex
    TYPE(CRTM_RTSolution_type)  , INTENT(IN OUT) :: RT
    ! Local variables
    INTEGER :: iFOV

    ! Get the FOV index value
    CALL CRTM_GeometryInfo_GetValue( gI, iFOV = iFOV )

    ! Compute the antenna temperature
    ! Note the earth temperature is used as a proxy
    ! for the platform temperature.
    RT%Brightness_Temperature = SC(n)%AC%A_earth(   iFOV,l)*RT%Brightness_Temperature + &
                                SC(n)%AC%A_platform(iFOV,l)*RT%Brightness_Temperature + &
                                SC(n)%AC%A_space(   iFOV,l)*TSPACE

  END SUBROUTINE CRTM_Compute_AntCorr

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Compute_AntCorr_TL
!
! PURPOSE:
!       Subroutine to compute the tangent-linear antenna correction and
!       apply it to the tangent-linear RTSolution brightness temperatures.
!
! CALLING SEQUENCE:
!       CALL CRTM_Compute_AntCorr_TL( GeometryInfo, &
!                                     SensorIndex , &
!                                     ChannelIndex, &
!                                     RTSolution_TL )
!
! INPUTS:
!       GeometryInfo:    Structure containing the view geometry information.
!                        UNITS:      N/A
!                        TYPE:       CRTM_GeometryInfo_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       SensorIndex:     Sensor index id. This is a unique index associated
!                        with a (supported) sensor used to access the
!                        shared coefficient data for a particular sensor.
!                        See the ChannelIndex argument.
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       ChannelIndex:    Channel index id. This is a unique index associated
!                        with a (supported) sensor channel used to access the
!                        shared coefficient data for a particular sensor's
!                        channel.
!                        See the SensorIndex argument.
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       RTSolution_TL:   Structure containing the tangent-linear radiative
!                        transfer solution.
!                        UNITS:      N/A
!                        TYPE:       CRTM_RTSolution_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN OUT)
!:sdoc-:
!------------------------------------------------------------------------------

  SUBROUTINE CRTM_Compute_AntCorr_TL( gI   , &  ! Input
                                      n    , &  ! Input
                                      l    , &  ! Input
                                      RT_TL  )  ! In/Output
    ! Arguments
    TYPE(CRTM_GeometryInfo_type), INTENT(IN)     :: gI
    INTEGER                     , INTENT(IN)     :: n  ! SensorIndex
    INTEGER                     , INTENT(IN)     :: l  ! ChannelIndex
    TYPE(CRTM_RTSolution_type)  , INTENT(IN OUT) :: RT_TL
    ! Local variables
    INTEGER :: iFOV
    REAL(fp) :: c

    ! Get the FOV index value
    CALL CRTM_GeometryInfo_GetValue( gI, iFOV = iFOV )

    ! Compute the tangent linear antenna temperature
    ! Note the A_platform term has to be included even
    ! though the earth temperature is used as a proxy
    ! for the platform temperature.
    c = SC(n)%AC%A_earth(iFOV,l) + SC(n)%AC%A_platform(iFOV,l)
    RT_TL%Brightness_Temperature = c * RT_TL%Brightness_Temperature

  END SUBROUTINE CRTM_Compute_AntCorr_TL

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Compute_AntCorr_AD
!
! PURPOSE:
!       Subroutine to compute the adjoint antenna correction and
!       apply it to the adjoint RTSolution brightness temperatures.
!
! CALLING SEQUENCE:
!       CALL CRTM_Compute_AntCorr_AD( GeometryInfo, &
!                                     SensorIndex , &
!                                     ChannelIndex, &
!                                     RTSolution_AD )
!
! INPUTS:
!       GeometryInfo:    Structure containing the view geometry information.
!                        UNITS:      N/A
!                        TYPE:       CRTM_GeometryInfo_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       SensorIndex:     Sensor index id. This is a unique index associated
!                        with a (supported) sensor used to access the
!                        shared coefficient data for a particular sensor.
!                        See the ChannelIndex argument.
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       ChannelIndex:    Channel index id. This is a unique index associated
!                        with a (supported) sensor channel used to access the
!                        shared coefficient data for a particular sensor's
!                        channel.
!                        See the SensorIndex argument.
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       RTSolution_AD:   Structure containing the adjoint radiative
!                        transfer solution.
!                        UNITS:      N/A
!                        TYPE:       CRTM_RTSolution_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN OUT)
!:sdoc-:
!------------------------------------------------------------------------------

  SUBROUTINE CRTM_Compute_AntCorr_AD( gI   , &  ! Input
                                      n    , &  ! Input
                                      l    , &  ! Input
                                      RT_AD  )  ! In/Output
    ! Arguments
    TYPE(CRTM_GeometryInfo_type), INTENT(IN)     :: gI
    INTEGER                     , INTENT(IN)     :: n  ! SensorIndex
    INTEGER                     , INTENT(IN)     :: l  ! ChannelIndex
    TYPE(CRTM_RTSolution_type)  , INTENT(IN OUT) :: RT_AD
    ! Local variables
    INTEGER :: iFOV
    REAL(fp) :: c

    ! Get the FOV index value
    CALL CRTM_GeometryInfo_GetValue( gI, iFOV = iFOV )

    ! Compute the adjoint of the antenna temperature
    c = SC(n)%AC%A_earth(iFOV,l) + SC(n)%AC%A_platform(iFOV,l)
    RT_AD%Brightness_Temperature = c * RT_AD%Brightness_Temperature

  END SUBROUTINE CRTM_Compute_AntCorr_AD

END MODULE CRTM_AntennaCorrection


text:
!>
!! @copyright Copyright (c) 2016-, Issam SAID <said.issam@gmail.com>
!! All rights reserved.
!!
!! Redistribution and use in source and binary forms, with or without
!! modification, are permetted 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 inthe
!!    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 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
!! INCLUDING, BUT NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY AND FITNESS
!! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
!! HOLDER OR ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
!! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
!! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
!! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
!! LIABILITY, WETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
!! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
!! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!!
!! @file m_ezcu_mem.f90
!! @author Issam SAID
!! @brief This file implements the Fortran interface of the ezcu memory
!! manipulation routines.
!!
!! @details All the ezcu memory management routines are implemented in Fortran
!! for arrays preallocated in stack, for allocatables and for pointers. In each
!! case the arrays can be 1 dimensional, 2 dimensional or 3 dimensional.
!! @see ezcu/mem.h
!<
module m_ezcu_mem
    use, intrinsic :: iso_c_binding
    use m_ezcu_flags
    use m_ezcu_types
    use m_ezcu_mem_interface
    use m_ezcu_mem_wrap
    use m_ezcu_mem_update
end module m_ezcu_mem

max_stars_repo_name: issamsaid/ezCU

max_issues_repo_name: issamsaid/ezCU

max_forks_repo_name: issamsaid/ezCU


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

      SUBROUTINE GETDUM (I, DUMMY, LEN)
C***********************************************************************

C  SUBROUTINE GETDUM = GETS AN INTEGER INTO A DUMMY CHARACTER STRING

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

      CHARACTER*72 DUMMY

      DUMMY = ' '
      IF (I .LT. -9999) THEN
         WRITE(DUMMY(1:6),10050)I
         LEN = 6
      ELSEIF (I .LT. -999) THEN
         WRITE (DUMMY(1:5), 10040) I
         LEN = 5
      ELSEIF (I .LT. -99) THEN
         WRITE (DUMMY(1:4), 10030) I
         LEN = 4
      ELSEIF (I .LT. -9) THEN
         WRITE (DUMMY(1:3), 10020) I
         LEN = 3
      ELSEIF (I .LT. 0) THEN
         WRITE (DUMMY(1:2), 10010) I
         LEN = 6
      ELSEIF(I .LT. 10) THEN
         WRITE (DUMMY(1:1), 10000) I
         LEN = 1
      ELSEIF(I.LT.100)THEN
         WRITE (DUMMY(1:2), 10010) I
         LEN = 2
      ELSEIF (I .LT. 1000) THEN
         WRITE (DUMMY(1:3), 10020) I
         LEN = 3
      ELSEIF (I .LT. 10000) THEN
         WRITE (DUMMY(1:4), 10030) I
         LEN = 4
      ELSE
         WRITE (DUMMY(1:5), 10040) I
         LEN = 5
      ENDIF
      RETURN

10000 FORMAT (I1)
10010 FORMAT (I2)
10020 FORMAT (I3)
10030 FORMAT (I4)
10040 FORMAT (I5)
10050 FORMAT (I6)

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
c=======================================================================
c FFTLog
c    computes the discrete Fast Fourier Transform
c    or Fast Hankel Transform (of arbitrary real index)
c    of a periodic logarithmic sequence.
c
c---Start of revision history-------------------------------------------
c
c Original version 27 March 1999.
c
c 5 July 1999:
c Bug fix in singular transform.
c Store amp, arg of last element of wsave separately, for even n,
c enlarging wsave by 1, for even n and q != 0.

c 17 July 1999:
c Reorganized wsave more logically.
c drffti now gets first, not last, 2*n+15 elements of wsave.
c fhti now treats last element of wsave with others, not by itself;
c wsave is enlarged from original, for even n,
c by 1 for q = 0, and by 2 for q != 0.
c
c 18 July 1999:
c A backward transform (dir = -1) is the same as
c a forward transform with q -> -q [and rk -> 1/rk]
c for any kr if n is odd,
c but only for low-ringing kr if n is even.
c Original falsely stated that transforms were same irrespective of kr.
c
c 19 Dec 1999:
c Slightly spiffed up interactive option to change kr to low-ringing kr.
c
c 13 Mar 2000:
c Made g77-safe:
c 1. Checked for any automatic variables that needed explicitly
c    `save'ing - actually there weren't any.
c    g77 is unusual in that it does NOT save variables not in
c    common or data statements.
c 2. All double precision constants explicitly suffixed d0.
c    g77 is unusual in that it does NOT automatically promote
c    real constants such as 1.1 to double precision 1.1d0.
c
c---End of revision history---------------------------------------------
c
c For more information about FFTLog, see
c
c http://casa.colorado.edu/~ajsh/FFTLog/
c
c Andrew J S Hamilton March 1999
c email: Andrew.Hamilton@Colorado.EDU
c
c Refs:	Talman J. D., 1978, J. Comp. Phys., 29, 35
c	Hamilton A. J. S., 2000, MNRAS, 312, 257
c	( http://xxx.lanl.gov/abs/astro-ph/9905191 )
c
c FFTLog uses the NCAR suite of FFT routines,
c and a modified version of the complex Gamma function
c from the gamerf package at
c http://momonga.t.u-tokyo.ac.jp/~ooura/gamerf.html .
c The original gamerf copyright statement states:
c   Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
c   You may use, copy, modify this code for any purpose and
c   without fee. You may distribute this ORIGINAL package.
c
c Permission to distribute the modified gamma function code
c with the FFTLog package has been granted
c (email from Takuya Ooura to Andrew Hamilton dated 16 March 1999).
c
c-----------------------------------------------------------------------
c Note: email programs may zap the 8-bit characters in the program
c       comments by converting them into 7-bit characters.
c       If this happens, get an unzapped copy from
c       http://casa.colorado.edu/~ajsh/FFTLog/ .
c       The fortran code itself contains only 7-bit characters;
c       the 8-bit characters occur only in comments.
c
c       If you are on a UNIX machine and 8-bit characters are not
c       showing up properly (they may appear as \nnn where n are
c       integers), then try putting
c          setenv LC_CTYPE iso_8859_1
c       in your .login file.
c
c-----------------------------------------------------------------------
c FFTLog computes a discrete version of
c the Hankel Transform (= Fourier-Bessel Transform)
c with a power law bias (k r)^q
c
c          infinity
c           /           q
c    ã(k) = | a(r) (k r)  J  (k r) k dr
c           /              mu
c          0
c
c          infinity
c           /           -q
c    a(r) = | ã(k) (k r)   J  (k r) r dk
c           /               mu
c          0
c
c where J_mu is the Bessel function of order mu.
c The index mu may be any real number, positive or negative.
c
c The input array a_j is a periodic sequence of length n,
c uniformly logarithmically spaced with spacing dlnr
c    a_j = a(r_j)   at   r_j = r_c exp[(j-j_c) dlnr]
c centred about the point r_c.  The central index j_c = (n+1)/2
c is 1/2 integral if n is even.
c Similarly, the output array ã_j is a periodic sequence of length n,
c also uniformly logarithmically spaced with spacing dlnr
c    ã_j = ã(k_j)   at   k_j = k_c exp[(j-j_c) dlnr]
c centred about the point k_c.
c
c The centre points r_c and k_c of the periodic intervals may be
c chosen arbitrarily; but it would be normal to choose the product
c kr = k_c r_c = k_j r_(n+1-j) = k_(n+1-j) r_j
c to be about 1 (or 2, or pi, to taste).
c
c The FFTLog algorithm is (see Hamilton 2000):
c 1. FFT the input array a_j to obtain the Fourier coefficients c_m ;
c 2. Multiply c_m by
c       u_m = (kr)^[- i 2 m pi/(n dlnr)] U_mu[q + i 2 m pi/(n dlnr)]
c    where
c       U_mu(x) = 2^x Gamma[(mu+1+x)/2] / Gamma[(mu+1-x)/2]
c    to obtain c_m u_m ;
c 3. FFT c_m u_m back to obtain the discrete Hankel transform ã_j .
c
c-----------------------------------------------------------------------
c The Fourier sine and cosine transforms
c
c                     infinity
c                      /
c    Ã(k) = sqrt(2/pi) | A(r) sin(k r) dr
c                      /
c                     0
c
c                     infinity
c                      /
c    Ã(k) = sqrt(2/pi) | A(r) cos(k r) dr
c                      /
c                     0
c
c may be regarded as special cases of the Hankel transform
c with mu = 1/2 and -1/2 since
c
c    sqrt(2/pi) sin(x) = sqrt(x) J   (x)
c                                 1/2
c
c    sqrt(2/pi) cos(x) = sqrt(x) J    (x)
c                                 -1/2
c
c The Fourier transforms may be done by making the substitutions
c                 q-(1/2)                      -q-(1/2)
c    A(r) = a(r) r          and   Ã(k) = ã(k) k
c
c and Hankel transforming a(r) with a power law bias (k r)^q
c
c          infinity
c           /           q
c    ã(k) = | a(r) (k r)  J    (k r) k dr
c           /              ±1/2
c          0
c
c Different choices of power law bias q lead to different discrete
c Fourier transforms of A(r), because the assumption of periodicity of
c a(r) = A(r) r^[-q+(1/2)] is different for different q.
c
c If A(r) is a power law, A(r) proportional to r^[q-(1/2)],
c then applying a bias q yields a discrete Fourier transform Ã(k)
c that is exactly equal to the continuous Fourier transform,
c because then a(r) is a constant, which is a periodic function.
c
c-----------------------------------------------------------------------
c The Hankel transform
c
c          infinity
c           /
c    Ã(k) = | A(r) J  (k r) k dr
c           /       mu
c          0
c
c may be done by making the substitutions
c                 q                      -q
c    A(r) = a(r) r    and   Ã(k) = ã(k) k
c
c and Hankel transforming a(r) with a power law bias (k r)^q
c
c          infinity
c           /           q
c    ã(k) = | a(r) (k r)  J  (k r) k dr
c           /              mu
c          0
c
c Different choices of power law bias q lead to different discrete
c Hankel transforms of A(r), because the assumption of periodicity of
c a(r) = A(r) r^-q is different for different q.
c
c If A(r) is a power law, A(r) proportional to r^q,
c then applying a bias q yields a discrete Hankel transform Ã(k)
c that is exactly equal to the continuous Hankel transform,
c because then a(r) is a constant, which is a periodic function.
c
c-----------------------------------------------------------------------
c ------------------------
c There are five routines:
c ------------------------
c Comments in the subroutines contain further details.
c
c (1) subroutine fhti(n,mu,q,dlnr,kr,kropt,wsave,ok)
c     is an initialization routine.
c
c (2) subroutine fftl(n,a,rk,dir,wsave)
c     computes the discrete Fourier sine or cosine transform
c     of a logarithmically spaced periodic sequence.
c     This is a driver routine that calls fhtq.
c
c (3) subroutine fht(n,a,dir,wsave)
c     computes the discrete Hankel transform
c     of a logarithmically spaced periodic sequence.
c     This is a driver routine that calls fhtq.
c
c (4) subroutine fhtq(n,a,dir,wsave)
c     computes the biased discrete Hankel transform
c     of a logarithmically spaced periodic sequence.
c     This is the basic FFTLog routine.
c
c (5) real*8 function krgood(mu,q,dlnr,kr)
c     takes an input kr and returns the nearest low-ringing kr.
c     This is an optional routine called by fhti.
c
c=======================================================================
c THE FFTLog CODE
c=======================================================================
      subroutine fhti(n,mu,q,dlnr,kr,kropt,wsave,ok)
      integer n,kropt
      logical ok
      real*8 mu,q,dlnr,kr,wsave(*)
c
c fhti initializes the working array wsave
c used by fftl, fht, and fhtq, and by NCAR routines drfftf and drfftb.
c fhti need be called once, whereafter fftl, fht, or fhtq may be
c called many times, as long as n, mu, q, dlnr, and kr remain unchanged.
c fhti should be called each time n, mu, q, dlnr, or kr is changed.
c The work array wsave should not be changed between calls to fftl,
c fht, or fhtq.
c
c If you are using the g77 fortran compiler, which by default does not
c save locally declared variables (unlike most other fortran compilers),
c then you should explicitly
c       save wsave
c in the calling program.
c
c  Input: n = number of points in the array to be transformed;
c             n may be any positive integer, but the NCAR FFT routines
c             run fastest if n is a product of small primes 2, 3, 5.
c	  mu = index of J_mu in Hankel transform;
c	       mu may be any real number, positive or negative.
c	  q = exponent of power law bias;
c	      q may be any real number, positive or negative.
c             If in doubt, use q = 0, for which case the Hankel
c             transform is orthogonal, i.e. self-inverse,
c             provided also that, for n even, kr is low-ringing.
c             Non-zero q may yield better approximations to the
c             continuous Hankel transform for some functions.
c         dlnr = separation between natural log of points;
c		 dlnr may be positive or negative.
c	  kr = k_c r_c where c is central point of array
c              = k_j r_(n+1-j) = k_(n+1-j) r_j .
c	       Normally one would choose kr to be about 1
c	       (or 2, or pi, to taste).
c         kropt = 0 to use input kr as is;
c                 1 to change kr to nearest low-ringing kr, quietly;
c                 2 to change kr to nearest low-ringing kr, verbosely;
c                 3 for option to change kr interactively.
c Output: wsave = working array used by fftl, fht, and fhtq,
c                 and by NCAR FFT routines drfftf and drfftb.
c                 wsave should be dimensioned at least:
c                 for q = 0 (unbiased transform):
c                    2*n+2*(n/2)+18
c                       [= 3*n+18 if n even, 3*n+17 if n odd]
c                 for q != 0 (biased transform):
c                    2*n+3*(n/2)+19
c                       [= (7*n)/2+19 if n even, (7*n)/2+18 if n odd]
c                 The first 2*n+15 elements of wsave are used by
c                 the NCAR FFT routines.
c         ok = .true. if all went ok;
c              .false. if not; currently only occurs if interactive
c                      response causes exit.
c
c        parameters
      real*8 PI
      parameter (PI=3.141592653589793238462643383279502884197d0)
      real*8 ZERO,ONE,TWO
      parameter (ZERO=0.d0,ONE=1.d0,TWO=2.d0)
      real*8 ROUND
      parameter (ROUND=1.d-15)
c        externals
      integer lnblnk,stblnk
      real*8 krgood
      complex*16 cdgamma
c        local (automatic) variables
      character*1 go
      character*64 temp
      integer l,m
      real*8 amp,arg,d,ln2,ln2kr,xm,xp,y
      complex*16 zm,zp
c
c--------adjust kr
c        keep kr as is
      if (kropt.eq.0) then
        continue
c        change kr to low-ringing kr quietly
      elseif (kropt.eq.1) then
        kr=krgood(mu,q,dlnr,kr)
c        change kr to low-ringing kr verbosely
      elseif (kropt.eq.2) then
        d=krgood(mu,q,dlnr,kr)
        if (abs(kr/d-ONE).gt.ROUND) then
          kr=d
          write (*,'(" kr changed to",g24.16)') kr
        endif
c        option to change kr to low-ringing kr interactively
      else
        d=krgood(mu,q,dlnr,kr)
        if (abs(kr/d-ONE).gt.ROUND) then
c        fortran demonstrates its inferiority to C
          write (*,'(" change kr = ",$)')
          write (temp,*) kr
          write (*,'(a,$)') temp(stblnk(temp):lnblnk(temp))
          write (*,'(" to low-ringing kr = ",$)')
          write (temp,*) d
          write (*,'(a,$)') temp(stblnk(temp):lnblnk(temp))
          write (*,'("? [CR,y=yes, n=no, x=exit]: ",$)')
          read (*,'(a1)') go
          if (go.eq.' '.or.go.eq.'y'.or.go.eq.'Y') then
            kr=d
            write (*,'(" kr changed to",g24.16)') kr
          elseif (go.eq.'n'.or.go.eq.'N') then
            print *,'kr left unchanged at',kr
          else
            print *,'exit'
            goto 300
          endif
        endif
      endif
c--------return if n is <= 0
      if (n.le.0) goto 200
c--------initialize normal FFT
      call drffti(n,wsave)
c        drfft uses first 2*n+15 elements of wsave
      l=2*n+15
c--------put q, dlnr, kr in next 3 elements of wsave
      wsave(l+1)=q
      wsave(l+2)=dlnr
      wsave(l+3)=kr
c        so far used 2*n+18 elements of wsave
      l=l+3
c--------rest of wsave used by fhtq - unbiased case (q = 0)
      if (q.eq.ZERO) then
        ln2kr=log(TWO/kr)
        xp=(mu+ONE)/TWO
        d=PI/(n*dlnr)
        do m=1,n/2
c        y = m pi/(n dlnr)
          y=m*d
          zp=dcmplx(xp,y)
          zp=cdgamma(zp,1)
c        Argument of kr^(-2 i y) U_mu(2 i y)
          arg=TWO*(ln2kr*y+dimag(zp))
          wsave(l+2*m-1)=cos(arg)
          wsave(l+2*m)=sin(arg)
        enddo
c Altogether 2*n+18 + 2*(n/2) = 2*((3*n)/2)+18 elements used for q = 0,
c which is 3*n+18 for even n, 3*n+17 for odd n.
c--------rest of wsave used by fhtq - biased case (q != 0)
      else
        ln2=log(TWO)
        ln2kr=log(TWO/kr)
        xp=(mu+ONE+q)/TWO
        xm=(mu+ONE-q)/TWO
c........first element of rest of wsave
        y=ZERO
c        case where xp or xm is a negative integer
        if ((anint(xp).eq.xp.and.anint(xp).le.ZERO)
     *    .or.(anint(xm).eq.xm.and.anint(xm).le.ZERO)) then
c        case where xp and xm are both negative integers
c U_mu(q) = 2^q Gamma[xp]/Gamma[xm] is finite in this case
          if ((anint(xp).eq.xp.and.anint(xp).le.ZERO)
     *      .and.(anint(xm).eq.xm.and.anint(xm).le.ZERO)) then
c        Amplitude and Argument of U_mu(q)
            amp=exp(ln2*q)
            if (xp.gt.xm) then
              do m=1,nint(xp-xm)
                amp=amp*(xm+m-ONE)
              enddo
            elseif (xp.lt.xm) then
              do m=1,nint(xm-xp)
                amp=amp/(xp+m-ONE)
              enddo
            endif
            arg=anint(xp+xm)*PI
c        one of xp or xm is a negative integer
          else
c Transformation is singular if xp is -ve integer,
c and inverse transformation is singular if xm is -ve integer,
c but transformation may be well-defined if sum_j a_j = 0,
c as may well occur in physical cases.
c Policy is to drop the potentially infinite constant in the transform.
            if (anint(xp).eq.xp.and.anint(xp).le.ZERO) then
              print *,'fhti: (mu+1+q)/2 =',nint(xp),
     &          ' is -ve integer, yields singular transform:'
              print *,'   transform will omit additive constant that is
     &generically infinite,'
              print *,'   but that may be finite or zero if'
              print *,'   the sum of the elements of the input array a_j
     & is zero.'
            else
              print *,'fhti: (mu+1-q)/2 =',nint(xm),
     &          ' is -ve integer, yields singular inverse transform:'
              print *,'   inverse transform will omit additive constant
     &that is generically'
              print *,'   infinite, but that may be finite or zero if'
              print *,'   the sum of the elements of the input array a_j
     & is zero.'
            endif
            amp=ZERO
            arg=ZERO
          endif
c        neither xp nor xm is a negative integer
        else
          zp=dcmplx(xp,y)
          zm=dcmplx(xm,y)
          zp=cdgamma(zp,1)
          zm=cdgamma(zm,1)
c        Amplitude and Argument of U_mu(q)
          amp=exp(ln2*q+dble(zp)-dble(zm))
c        note +Im(zm) to get conjugate value below real axis
          arg=dimag(zp)+dimag(zm)
        endif
c        cos(arg) = ±1, sin(arg) = 0
        wsave(l+1)=amp*cos(arg)
c........remaining elements of wsave
        d=PI/(n*dlnr)
        do m=1,n/2
c        y = m pi/(n dlnr)
          y=m*d
          zp=dcmplx(xp,y)
          zm=dcmplx(xm,y)
          zp=cdgamma(zp,1)
          zm=cdgamma(zm,1)
c        Amplitude and Argument of kr^(-2 i y) U_mu(q + 2 i y)
          amp=exp(ln2*q+dble(zp)-dble(zm))
          arg=2*ln2kr*y+dimag(zp)+dimag(zm)
          wsave(l+3*m-1)=amp
          wsave(l+3*m)=cos(arg)
          wsave(l+3*m+1)=sin(arg)
        enddo
c Altogether 2*n+18 + 3*(n/2)+1 elements used for q != 0,
c which is (7*n)/2+19 for even n, (7*n)/2+18 for odd n.
c For even n, the very last element of wsave
c [i.e. wsave(l+3*m+1)=sin(arg) for m=n/2] is not used within FFTLog;
c if a low-ringing kr is used, this element should be zero.
c The last element is computed in case somebody wants it.
      endif
  200 ok=.true.
      return
c
c--------error returns
  300 ok=.false.
      return
      end
c
c=======================================================================
      subroutine fftl(n,a,rk,dir,wsave)
      integer n,dir
      real*8 a(n),rk,wsave(*)
c
c This is a driver routine that calls fhtq.
c
c fftl computes a discrete version of the Fourier
c sine (if mu = 1/2) or cosine (if mu = -1/2) transform
c
c                     infinity
c                      /
c    Ã(k) = sqrt(2/pi) | A(r) sin(k r) dr
c                      /
c                     0
c
c                     infinity
c                      /
c    Ã(k) = sqrt(2/pi) | A(r) cos(k r) dr
c                      /
c                     0
c
c by making the substitutions
c                 q-(1/2)                      -q-(1/2)
c    A(r) = a(r) r          and   Ã(k) = ã(k) k
c
c and applying a biased Hankel transform to a(r).
c
c The steps are:
c 1. a(r) = A(r) r^[-dir*(q-.5)]
c 2. call fhtq to transform a(r) -> ã(k)
c 3. Ã(k) = ã(k) k^[-dir*(q+.5)]
c
c fhti must be called before the first call to fftl,
c with mu = 1/2 for a sine transform,
c or mu = -1/2 for a cosine transform.
c
c A call to fftl with dir=1 followed by
c a call to fftl with dir=-1 (and rk unchanged), or vice versa,
c leaves the array a unchanged.
c
c  Input: n = length of a array.
c         rk = r_c/k_c
c              = r_j/k_j (a constant, the same constant for any j);
c              rk is not (necessarily) the same quantity as kr.
c              rk is used only to multiply the output array by
c              sqrt(rk)^dir, so if you want to do the normalization
c              later, or you don't care about the normalization,
c              you can set rk = 1.
c	  dir = 1 for forward transform,
c		-1 for backward transform.
c               A backward transform (dir = -1) is the same as
c               a forward transform with q -> -q and rk -> 1/rk,
c               for any kr if n is odd,
c               for low-ringing kr if n is even.
c	  wsave = working array set up by fhti.
c Input/Output:
c         a on  input is the array A(r) to transform:
c             a(j) is A(r_j) at r_j = r_c exp[(j-jc) dlnr]
c                where jc = (n+1)/2 = central index of array.
c         a on output is the transformed array Ã(k):
c             a(j) is Ã(k_j) at k_j = k_c exp[(j-jc) dlnr].
c
c        parameters
      real*8 ONE,TWO,HALF
      parameter (ONE=1.d0,TWO=2.d0,HALF=ONE/TWO)
c        local (automatic) variables
      integer j,l
      real*8 dlnr,jc,kr,lnkr,lnrk,q
c
      l=2*n+15
      q=wsave(l+1)
      dlnr=wsave(l+2)
      kr=wsave(l+3)
c........a(r) = A(r) (r/rc)^[-dir*(q-.5)]
c        centre point of array
      jc=dble(n+1)/TWO
      do j=1,n
        a(j)=a(j)*exp(-dir*(q-HALF)*(j-jc)*dlnr)
      enddo
c........transform a(r) -> ã(k)
      call fhtq(n,a,dir,wsave)
c........Ã(k) = ã(k) k^[-dir*(q+.5)] rc^[-dir*(q-.5)]
c        = ã(k) (k/kc)^[-dir*(q+.5)] (kc rc)^(-dir*q) (rc/kc)^(dir*.5)
      lnkr=log(kr)
      lnrk=log(rk)
      do j=1,n
        a(j)=a(j)*exp(-dir*((q+HALF)*(j-jc)*dlnr+q*lnkr-HALF*lnrk))
      enddo
      return
      end
c
c=======================================================================
      subroutine fht(n,a,dir,wsave)
      integer n,dir
      real*8 a(n),wsave(*)
c
c This is a driver routine that calls fhtq.
c
c fht computes a discrete version of the Hankel transform
c
c          infinity
c           /
c    Ã(k) = | A(r) J  (k r) k dr
c           /       mu
c          0
c
c by making the substitutions
c                 q                      -q
c    A(r) = a(r) r    and   Ã(k) = ã(k) k
c
c and applying a biased Hankel transform to a(r).
c
c The steps are:
c 1. a(r) = A(r) r^(-dir*q)
c 2. call fhtq to transform a(r) -> ã(k)
c 3. Ã(k) = ã(k) k^(-dir*q)
c
c fhti must be called before the first call to fht.
c
c A call to fht with dir=1 followed by
c a call to fht with dir=-1, or vice versa,
c leaves the array a unchanged.
c
c  Input: n = length of a array.
c	  dir = 1 for forward transform,
c		-1 for backward transform.
c               A backward transform (dir = -1) is the same as
c               a forward transform with q -> -q,
c               for any kr if n is odd,
c               for low-ringing kr if n is even.
c	  wsave = working array set up by fhti.
c Input/Output:
c         a on  input is the array A(r) to transform:
c             a(j) is A(r_j) at r_j = r_c exp[(j-jc) dlnr]
c                where jc = (n+1)/2 = central index of array.
c         a on output is the transformed array Ã(k):
c             a(j) is Ã(k_j) at k_j = k_c exp[(j-jc) dlnr].
c
c        parameters
      real*8 ZERO,TWO
      parameter (ZERO=0.d0,TWO=2.d0)
c        local (automatic) variables
      integer j,l
      real*8 dlnr,jc,kr,lnkr,q
c
      l=2*n+15
      q=wsave(l+1)
      dlnr=wsave(l+2)
      kr=wsave(l+3)
c........a(r) = A(r) (r/rc)^(-dir*q)
      if (q.ne.ZERO) then
c        centre point of array
        jc=dble(n+1)/TWO
        do j=1,n
          a(j)=a(j)*exp(-dir*q*(j-jc)*dlnr)
        enddo
      endif
c........transform a(r) -> ã(k)
      call fhtq(n,a,dir,wsave)
c........Ã(k) = ã(k) (k rc)^(-dir*q)
c             = ã(k) (k/kc)^(-dir*q) (kc rc)^(-dir*q)
      if (q.ne.ZERO) then
        lnkr=log(kr)
        do j=1,n
          a(j)=a(j)*exp(-dir*q*((j-jc)*dlnr+lnkr))
        enddo
      endif
      return
      end
c
c=======================================================================
      subroutine fhtq(n,a,dir,wsave)
      integer n,dir
      real*8 a(n),wsave(*)
c
c This is the basic FFTLog routine.
c
c fhtq computes a discrete version of the biased Hankel transform
c
c          infinity
c           /           q
c    ã(k) = | a(r) (k r)  J  (k r) k dr
c           /              mu
c          0
c
c fhti must be called before the first call to fhtq.
c
c A call to fhtq with dir=1 followed by
c a call to fhtq with dir=-1, or vice versa,
c leaves the array a unchanged.
c
c  Input: n = length of a array.
c	  dir = 1 for forward transform,
c		-1 for backward transform.
c               A backward transform (dir = -1) is the same as
c               a forward transform with q -> -q,
c               for any kr if n is odd,
c               for low-ringing kr if n is even.
c Input/Output:
c         a on  input is the periodic array a(r) to transform:
c             a(j) is a(r_j) at r_j = r_c exp[(j-jc) dlnr]
c                 where jc = (n+1)/2 = central index of array.
c         a on output is the transformed periodic array ã(k):
c             a(j) is ã(k_j) at k_j = k_c exp[(j-jc) dlnr].
c
c        parameters
      real*8 ZERO
      parameter (ZERO=0.d0)
c        local (automatic) variables
      integer l,m
      real*8 ai,ar,q
c
      l=2*n+15
      q=wsave(l+1)
      l=l+3
c--------normal FFT
      call drfftf(n,a,wsave)
c--------unbiased (q = 0) transform
      if (q.eq.ZERO) then
c........multiply by
c        (kr)^[- i 2 m pi/(n dlnr)] U_mu[i 2 m pi/(n dlnr)]
        do m=1,(n-1)/2
          ar=a(2*m)
          ai=a(2*m+1)
          a(2*m)=ar*wsave(l+2*m-1)-ai*wsave(l+2*m)
          a(2*m+1)=ar*wsave(l+2*m)+ai*wsave(l+2*m-1)
        enddo
c        problematical last element, for even n
        if (mod(n,2).eq.0) then
          ar=wsave(l+n-1)
c        forward transform: multiply by real part
c Why?  See http://casa.colorado.edu/~ajsh/FFTLog/index.html#ure
          if (dir.eq.1) then
            a(n)=a(n)*ar
c        backward transform: divide by real part
          elseif (dir.eq.-1) then
c Real part ar can be zero for maximally bad choice of kr.
c This is unlikely to happen by chance, but if it does,
c policy is to let it happen.
c For low-ringing kr, imaginary part ai is zero by construction,
c and real part ar is guaranteed nonzero.
            a(n)=a(n)/ar
          endif
        endif
c--------biased (q != 0) transform
      else
c........multiply by
c        (kr)^[- i 2 m pi/(n dlnr)] U_mu[q + i 2 m pi/(n dlnr)]
c        phase
        do m=1,(n-1)/2
          ar=a(2*m)
          ai=a(2*m+1)
          a(2*m)=ar*wsave(l+3*m)-ai*wsave(l+3*m+1)
          a(2*m+1)=ar*wsave(l+3*m+1)+ai*wsave(l+3*m)
        enddo
c        forward transform: multiply by amplitude
        if (dir.eq.1) then
          a(1)=a(1)*wsave(l+1)
          do m=1,(n-1)/2
            a(2*m)=a(2*m)*wsave(l+3*m-1)
            a(2*m+1)=a(2*m+1)*wsave(l+3*m-1)
          enddo
c        backward transform: divide by amplitude
        elseif (dir.eq.-1) then
c        amplitude of m=0 element
          ar=wsave(l+1)
          if (ar.eq.ZERO) then
c Amplitude of m=0 element can be zero for some mu, q combinations
c (singular inverse); policy is to drop potentially infinite constant.
            a(1)=ZERO
          else
            a(1)=a(1)/ar
          endif
c        remaining amplitudes should never be zero
          do m=1,(n-1)/2
            a(2*m)=a(2*m)/wsave(l+3*m-1)
            a(2*m+1)=a(2*m+1)/wsave(l+3*m-1)
          enddo
        endif
c        problematical last element, for even n
        if (mod(n,2).eq.0) then
          m=n/2
          ar=wsave(l+3*m)*wsave(l+3*m-1)
c        forward transform: multiply by real part
          if (dir.eq.1) then
            a(n)=a(n)*ar
c        backward transform: divide by real part
          elseif (dir.eq.-1) then
c Real part ar can be zero for maximally bad choice of kr.
c This is unlikely to happen by chance, but if it does,
c policy is to let it happen.
c For low-ringing kr, imaginary part ai is zero by construction,
c and real part ar is guaranteed nonzero.
            a(n)=a(n)/ar
          endif
        endif
      endif
c--------normal FFT back
      call drfftb(n,a,wsave)
c--------reverse the array
c        and at the same time undo the FFTs' multiplication by n
      do m=1,n/2
        ar=a(m)
        a(m)=a(n+1-m)/n
        a(n+1-m)=ar/n
      enddo
      if (mod(n,2).eq.1) then
        m=(n+1)/2
        a(m)=a(m)/n
      endif
      return
      end
c
c=======================================================================
      real*8 function krgood(mu,q,dlnr,kr)
      real*8 mu,q,dlnr,kr
c
c Use of this routine is optional.
c
c Choosing kr so that
c     (kr)^(- i pi/dlnr) U_mu(q + i pi/dlnr)
c is real may reduce ringing of the discrete Hankel transform,
c because it makes the transition of this function across the period
c boundary smoother.
c
c  Input: mu = index of J_mu in Hankel transform.
c	  q = exponent of power law bias.
c         dlnr = separation between natural log of points.
c	  kr = suggested value of kr.
c Output: krgood = low-ringing value of kr nearest to input kr.
c                  ln(krgood) is always within dlnr/2 of ln(kr).
c
c        parameters
      real*8 PI
      parameter (PI=3.141592653589793238462643383279502884197d0)
      real*8 ZERO,ONE,TWO
      parameter (ZERO=0.d0,ONE=1.d0,TWO=2.d0)
c        externals
      complex*16 cdgamma
c        local (automatic) variables
      real*8 arg,iarg,xm,xp,y
      complex*16 zm,zp
c
      krgood=kr
      if (dlnr.eq.ZERO) return
      xp=(mu+ONE+q)/TWO
      xm=(mu+ONE-q)/TWO
      y=PI/(TWO*dlnr)
      zp=dcmplx(xp,y)
      zm=dcmplx(xm,y)
      zp=cdgamma(zp,1)
      zm=cdgamma(zm,1)
c        low-ringing condition is that following should be integral
      arg=log(TWO/kr)/dlnr+(dimag(zp)+dimag(zm))/PI
      iarg=anint(arg)
c        if should ensure arg = +-Infinity dealt with correctly
      if (arg.ne.iarg) then
c        low-ringing kr
        krgood=kr*exp((arg-iarg)*dlnr)
      endif
      return
      end
c
c-----------------------------------------------------------------------
      integer function stblnk(string)
      character*(*) string
c
c        externals
      integer lnblnk
c        local (automatic) variables
      integer i
c *
c * Return index of first non-blank character in string.
c * If string is all blank, returned index is 1+lnblnk(string).
c *
      do i=1,lnblnk(string)
        if (string(i:i).ne.' ') goto 120
      enddo
  120 stblnk=i
      return
      end
c


text:
      subroutine BMG2_SymStd_SETUP_comm(COMM, NOL, MyProc, &
           Min_GSz, mp)&
           BIND(C, NAME='MPI_BMG2_SymStd_SETUP_comm')

      USE ModInterface
      use message_passing
      IMPLICIT NONE

! ==============================================================
! Include Statements
! ==============================================================

      INCLUDE 'mpif.h'

      INCLUDE 'BMG_parameters_f90.h'

! ==============================================================
! Input Variables
! ==============================================================

      integer(c_int), VALUE :: NOL, MyProc, Min_GSz

      integer :: COMM(2,2:NOL)
      type(c_ptr) :: mp

! ==============================================================
! Local Variables
! ==============================================================

      INTEGER TempProc

      INTEGER INSIZE, GSIZE, HSIZE

      INTEGER ACTIVE, N

      INTEGER IERR

! ==============================================================
! Initialize local variables
! ==============================================================

      call MPI_COMM_SIZE(COMM(1,NOL), INSIZE, IERR)

      N = NOL
      TempProc = MyProc
      ACTIVE = 1

! ==============================================================
! Call comm_coarsen to build communicators
! ==============================================================

      IF (INSIZE .GE. 2*Min_GSz) THEN

10       CONTINUE

         !
         ! If there are any process left in the current
         ! level communicator try to coasren.
         ! This is just a while-loop
         !
         IF (ACTIVE .EQ. 1) THEN

            call comm_coarsen(COMM(1,N), INSIZE, &
     &                        COMM(2,N), COMM(1,N-1),&
     &                        HSIZE, Min_GSz,&
     &                        TempProc, ACTIVE, mp)

            INSIZE = HSIZE
            N = N - 1

            goto 10

         ELSE

            goto 20

         END IF

20    CONTINUE

      END IF

      !
      ! Set COMM(2,1) equal to COMM(1,1).  This facilitates
      ! the loops used in the up and down portions of the
      ! linesolve.
      !

      COMM(2,2) = COMM(1,2)

      END

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

      SUBROUTINE comm_coarsen(INCOMM, INSIZE, &
     &                        GCOMM, HCOMM,&
     &                        NGPS, Min_GSz, &
     &                        MyProc, ACTIVE, mp)

      use message_passing
      IMPLICIT NONE

      INCLUDE 'mpif.h'

      INTEGER INCOMM, INSIZE
      INTEGER GCOMM, HCOMM
      INTEGER NGPS, Min_GSz
      INTEGER ACTIVE

      INTEGER GSIZE, HSIZE

      INTEGER MPI_MY_ID, MyProc
      INTEGER HFLAG, HRANK
      INTEGER MyGroup, MyGroupRank
      INTEGER LBRK, RBRK, REM, CLUMP

      INTEGER EXITFLAG, K, IERR

      type(c_ptr) :: mp

      !
      ! Compute the number of new groups to form
      !
      NGPS = INSIZE / Min_GSz
      REM = mod(INSIZE, Min_GSz)

      !
      ! This computers CLUMP = ceiling(REM/NGPS)
      !
      IF (real(REM)/NGPS .GT. REM/NGPS) THEN
         CLUMP = REM/NGPS + 1
      ELSe
         CLUMP = REM/NGPS
      END IF

      !
      ! Initialize counters and flags
      !
      LBRK = 0
      RBRK = 0
      EXITFLAG = 0
      K = 1

      !
      ! Determien the size of each group taking
      ! load balancing issues into account
      !
30    CONTINUE

      !
      ! This is just a while-loop
      !
      IF (EXITFLAG .NE. 1) THEN

         IF (REM .GE. CLUMP) THEN
            GSIZE = Min_GSz + CLUMP
            REM = REM - CLUMP
         ELSE
            GSIZE = Min_GSz + REM
            REM = 0
         END IF

         RBRK = RBRK + GSIZE

         !
         ! if local process fits in this
         ! group set group and local
         ! rank flags
         !
         IF (MyProc .LE. RBRK) THEN
            MyGroup = K
            MyGroupRank = MyProc - LBRK - 1
            EXITFLAG = 1
         ELSE
            K = K + 1
         END IF

         LBRK = RBRK

         goto 30

      ELSE

         goto 40

      END IF

40    CONTINUE

      !
      ! Form group communicators
      !
      call cedar_comm_split(mp, INCOMM, MyGroup, MyGroupRank,&
           GCOMM, IERR)

      !
      ! If I'm the first process in my group designate
      ! me the HEAD of my group.  Also check to see if
      ! I will be part of any other communicators at
      ! lower levels.  If I won't be then set ACTIVE
      ! to 0
      !
      IF (MyGroupRank .EQ. 0) THEN

         HFLAG = 1

         IF (NGPS .GE. 2*Min_GSz) THEN
            ACTIVE = 1
         ELSE
            ACTIVE = 0
         END IF

      ELSE

         HFLAG = MPI_UNDEFINED
         ACTIVE = 0

      END IF

      !
      ! Identify the rank of this process relative
      ! to its new subgroup
      !
      HRANK = MyGroup - 1

      !
      ! Form communicator with all of the HEAD nodes
      ! in it.  Non-head nodes get assigned to a
      ! communicator with type MPI_PROC_NULL
      !
      call cedar_comm_split(mp, INCOMM, HFLAG, HRANK, &
           HCOMM, IERR)

      !
      ! Update the rank that this process will have
      ! in the new head communicator
      !
      MyProc = MyGroup

      END

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar


text:
      program synthesis01

C     Goal: debug the code synthesis in initializer()
C     Check that a call graph of this routine can be obtained

      call foo(x)

      end


text:
C***********************************************************************
C***********************************************************************
C
C   Version:        0.3
C   Last modified:  December 27, 1994
C   Authors:        Joseph W.H. Liu
C
C   Mathematical Sciences Section, Oak Ridge National Laboratory
C
C***********************************************************************
C***********************************************************************
C--- SPARSPAK-A (ANSI FORTRAN) RELEASE III --- NAME = MMDINT
C  (C)  UNIVERSITY OF WATERLOO   JANUARY 1984
C***********************************************************************
C***********************************************************************
C***     MMDINT ..... MULT MINIMUM DEGREE INITIALIZATION     ***********
C***********************************************************************
C***********************************************************************
C
C     PURPOSE - THIS ROUTINE PERFORMS INITIALIZATION FOR THE
C        MULTIPLE ELIMINATION VERSION OF THE MINIMUM DEGREE
C        ALGORITHM.
C
C     INPUT PARAMETERS -
C        NEQNS  - NUMBER OF EQUATIONS.
C        (XADJ,ADJNCY) - ADJACENCY STRUCTURE.
C
C     OUTPUT PARAMETERS -
C        (DHEAD,DFORW,DBAKW) - DEGREE DOUBLY LINKED STRUCTURE.
C        QSIZE  - SIZE OF SUPERNODE (INITIALIZED TO ONE).
C        LLIST  - LINKED LIST.
C        MARKER - MARKER VECTOR.
C
C***********************************************************************
C
      SUBROUTINE  MMDINT ( NEQNS, XADJ, ADJNCY, DHEAD, DFORW,
     1                     DBAKW, QSIZE, LLIST, MARKER )
C
C***********************************************************************
C
         INTEGER    ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) ,
     1              LLIST(1) , MARKER(1), QSIZE(1)
         INTEGER    XADJ(1)
         INTEGER    FNODE , NDEG  , NEQNS , NODE
C
C***********************************************************************
C
         DO  100  NODE = 1, NEQNS
             DHEAD(NODE) = 0
             QSIZE(NODE) = 1
             MARKER(NODE) = 0
             LLIST(NODE) = 0
  100    CONTINUE
C        ------------------------------------------
C        INITIALIZE THE DEGREE DOUBLY LINKED LISTS.
C        ------------------------------------------
         DO  200  NODE = 1, NEQNS
             NDEG = XADJ(NODE+1) - XADJ(NODE) + 1
             FNODE = DHEAD(NDEG)
             DFORW(NODE) = FNODE
             DHEAD(NDEG) = NODE
             IF  ( FNODE .GT. 0 )  DBAKW(FNODE) = NODE
             DBAKW(NODE) = - NDEG
  200    CONTINUE
         RETURN
C
      END

max_stars_repo_name: Pangqiyuangh/SeIInv

max_issues_repo_name: Pangqiyuangh/SeIInv

max_forks_repo_name: Pangqiyuangh/SeIInv


text:
program ch0503
    ! Integer Division and Real Assignment
    implicit none

    integer :: i, j, k
    real :: answer

    i = 5
    j = 2
    k = 4
    answer = i/j*k
    print *, i, j, k, answer
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
module q5_fit
  use solver
  implicit none

  interface
    function RV2R(x)
      !! [Real] -> Real
      real, intent(in) :: x(:)
      real :: RV2R
    end function RV2R

    function RV2RV(x)
      !! [Real] -> [Real]
      real, intent(in) :: x(:)
      real :: RV2RV(size(x))
    end function RV2RV

    pure function pRVR2R(c, x)
      !! [Real], Real -> Real
      real, intent(in) :: c(:), x
      real :: pRVR2R
    end function pRVR2R

    pure function pRVR2RV(c, x)
      !! [Real], Real -> [Real]
      real, intent(in) :: c(:), x
      real :: pRVR2RV(size(c))
    end function pRVR2RV
  end interface

contains

  subroutine q5()
    !! Fit coeffs by finding minimum of chi-square via Levenberg-Marquardt
    integer, parameter :: n = 3
    real, dimension(:), allocatable :: x, y, y_fit
    real, dimension(n + 2) :: coeffs
    integer :: i

    call read_data("data.dat", x, y)
    allocate(y_fit(size(x)))

    print "(a)", "5. Data fit:"
    print *

    print "(a)", "   Levenberg-Marquardt:"
    coeffs = 0.0
    coeffs = levenberg_marquardt(x, y, model, dmodel, &
        c0=coeffs, lambda0=0.7, tol=0.1, maxiter=100)
    print "('   Fit params (0..n): ', 4f9.3)", coeffs(1:n+1)
    print "('   Fit params (const):', f9.3)",  coeffs(n+2)
    print "('   Loss:              ', f9.0)", loss(coeffs)
    print *

    forall (i=1:size(x)) y_fit(i) = model(coeffs, x(i))
    call write_csv("results_levenberg_marquardt.csv", x, y, y_fit)

    print "(a)", "   Gradient Descent:"
    coeffs = 0.0
    coeffs = gradient_descent(loss, dloss, c0=coeffs, tol=0.1, maxiter=2000)
    print "('   Fit params (0..n): ', 4f9.3)", coeffs(1:n+1)
    print "('   Fit params (const):', f9.3)",  coeffs(n+2)
    print "('   Loss:              ', f9.0)", loss(coeffs)
    print *

    forall (i=1:size(x)) y_fit(i) = model(coeffs, x(i))
    call write_csv("results_gradient_descent.csv", x, y, y_fit)

    print "(a)", "Evidently, Levenberg-Marquardt requires far fewer iterations"
    print "(a)", "than Gradient Descent. It also results in a better fit,"
    print "(a)", "with less tweaking needed -- GD does not converge for poorly"
    print "(a)", "chosen step sizes and step decay rates."
    print *

    deallocate(x, y, y_fit)

  contains

    function loss(c)
      real, intent(in) :: c(:)
      real :: loss

      forall (i=1:size(x)) y_fit(i) = model(c, x(i))
      loss = 0.5 * sum((y - y_fit)**2)
    end function

    function dloss(c)
      real, intent(in) :: c(:)
      real :: dloss(size(c))
      real :: df_dc(size(c), size(x))

      forall (i=1:size(x)) y_fit(i) = model(c, x(i))
      forall (i=1:size(x)) df_dc(:, i) = dmodel(c, x(i))
      forall (i=1:size(c)) dloss(i) = sum((y_fit - y) * df_dc(i, :))
    end function

  end subroutine

  subroutine read_data(filename, x, y)
    !! Read raw (x,y) data from .dat file
    character(len=*) :: filename
    real, dimension(:), allocatable :: x, y
    integer, parameter :: ifh = 1
    integer, parameter :: max_lines = 1048576
    integer :: i, stat
    real, dimension(max_lines) :: x_, y_

    ! Read data
    open(unit=ifh, file=filename, action="read", status="old", &
      access="sequential", form="formatted")
    do i = 1, max_lines
      read(ifh, *, iostat=stat) x_(i), y_(i)
      if (stat < 0) exit
    end do
    close(ifh)

    ! Allocate properly-sized arrays
    i = i - 1
    allocate(x(i), y(i))
    x = x_(1:i)
    y = y_(1:i)
  end subroutine

  subroutine write_csv(filename, x, y, y_fit)
    !! Output a csv file with plottable data
    real, dimension(:) :: x, y, y_fit
    character(len=*) :: filename
    integer, parameter :: ofh = 2
    integer :: i

    open(unit=ofh, file=filename, action="write", status="replace")
    write(ofh, *) "x, y, f(x)"
    do i = 1, size(x)
      write(ofh, *) x(i), ", ", y(i), ", ", y_fit(i)
    end do
    close(ofh)
  end subroutine

  function gradient_descent(f, df, c0, tol, maxiter) result(c)
    !! Minimize f using df starting at c0
    procedure(RV2R) :: f
    procedure(RV2RV) :: df
    real, intent(in) :: c0(:), tol
    integer, intent(in) :: maxiter
    real :: c(size(c0)), step, tol_new, loss, loss_new
    integer :: iter

    step = 1e-5
    c = c0
    loss = f(c)

    do iter = 1, maxiter
      c = c - step * df(c)
      step = step * 0.995

      loss_new = f(c)
      tol_new = abs(loss_new - loss)
      loss = loss_new

      if (tol_new <= tol) exit
    end do

    print "('   Iterations:        ', i9)", iter - 1
  end function

  function levenberg_marquardt(x, y, f, df, c0, lambda0, tol, maxiter) result(c)
    !! Minimize f w.r.t. chi-square for given data x, y points starting at c0
    procedure(pRVR2R) :: f
    procedure(pRVR2RV) :: df
    real, intent(in) :: x(:), y(:), c0(:), lambda0, tol
    integer, intent(in) :: maxiter
    real, dimension(size(c0)) :: c, c_new, dloss
    real, dimension(size(x)) :: r, r_new
    real, dimension(size(c0), size(c0)) :: g, JTJ
    real :: J(size(x), size(c0))
    real :: JT(size(c0), size(x))
    real :: loss, loss_new, lambda, tol_new
    logical :: invalidated
    integer :: i, iter
    real, parameter :: lambda_up = 1.1, lambda_down = 1 / 1.1

    lambda = lambda0
    c = c0
    forall (i=1:size(x)) r(i) = y(i) - f(c, x(i))
    loss = 0.5 * sum(r**2)
    invalidated = .true.

    do iter = 1, maxiter
      ! Compute J, JT, JTJ, and loss-gradient
      ! but only if coefficients have changed
      if (invalidated) then
        forall (i=1:size(x)) JT(:, i) = df(c, x(i))
        J = transpose(JT)
        JTJ = matmul(JT, J)
        dloss = matmul(JT, r)
        ! invalidated = .false.  ! NOTE this seems to cause weird issues
      end if

      ! Compute g
      g = JTJ
      forall (i=1:size(c0)) g(i, i) = (1.0 + lambda) * g(i, i)

      ! Compute new coefficients, new residuals, and new loss
      c_new = c + solve(g, dloss)
      forall (i=1:size(x)) r_new(i) = y(i) - f(c_new, x(i))
      loss_new = 0.5 * sum(r_new**2)
      tol_new = abs(loss_new - loss)

      ! If loss is smaller, accept new coefficients
      if (loss_new < loss) then
        c = c_new
        r = r_new
        loss = loss_new
        lambda = lambda * lambda_down
        invalidated = .true.
      else
        lambda = lambda * lambda_up
      end if

      if (tol_new <= tol) exit
    end do

    print "('   Iterations:        ', i9)", iter - 1
  end function

  elemental function basis(a, x)
    !! An even basis... for a more civilized age
    integer, intent(in) :: a
    real, intent(in) :: x
    real :: basis
    real, parameter :: pi = 3.14159265358979323846264338327950288419716939937510

    basis = cos(2 * pi * a * x)
  end function

  pure function model(c, x)
    !! Non-linear model that we will fit our data to
    !! f(x) = exp(sum(c_a * b_a(x))) + const
    real, intent(in) :: c(:), x
    real :: model, b(size(c) - 1)
    integer :: a, n
    n = size(c) - 1

    forall (a=1:n) b(a) = basis(a - 1, x)
    model = exp(sum(c(1:n) * b)) + c(n+1)
  end function

  pure function dmodel(c, x)
    !! Derivative of non-linear model with respect to coefficients
    !! df(x)/dc_k = exp(sum(c_a * b_a(x))) * b_k(x)
    real, intent(in) :: c(:), x
    real :: dmodel(size(c)), b(size(c) - 1), factor
    integer :: a, n
    n = size(c) - 1

    forall (a=1:n) b(a) = basis(a - 1, x)
    factor = exp(sum(c(1:n) * b))
    dmodel(1:n) = factor * b
    dmodel(n+1) = 1.0
  end function

end module q5_fit

max_stars_repo_name: YodaEmbedding/experiments

max_issues_repo_name: YodaEmbedding/experiments

max_forks_repo_name: YodaEmbedding/experiments


text:
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!  BIGSTICK configuration-interaction shell-model code
!
!  file BDIAGH.F
      module diagh

      implicit none
! diagonal interaction from single-particle energies

      type blockhspe
         real, pointer :: hspe(:)
      end type blockhspe

      real, pointer :: pspeME(:),nspeME(:)

      real, allocatable, target :: pspe_obs(:),nspe_obs(:)
      real, allocatable, target :: pspe_h(:),nspe_h(:)

contains

!
!  Diagonal part of Hamiltonian:
!  Hamiltonian routines for BIGSTICK
!
!  CALLED BY:
!      MAIN
!      lanczos_output
!      setup4obsmaster
!      pocc_compute_spocc
!      particle_occupation_p_orig
!
   subroutine makespeme(it,obschar)

!   use diagh
   use basis
   use coupledmatrixelements
   use interaction_deformed
   implicit none
!   include 'binterfaces.inc'

   integer it
   character(1) obschar
   real, pointer :: speME(:)
   real, pointer :: xspe(:)
   type (blockhspe), pointer :: lbspe(:), rbspe(:)
   integer :: aerr

   if(.not.call_spe)return

   select case (obschar)

   case ('H','J')

   if(it == 1)then
     if(.not.allocated(pspe_h)) then
        allocate(pspe_h( nXsd(it) ), stat=aerr)
        if(aerr /= 0) call memerror("makespeme 1")
     end if
     speME => pspe_h
     xspe  => pspe
   else
     if(.not.allocated(nspe_h)) then
        allocate(nspe_h( nXsd(it) ), stat=aerr)
        if(aerr /= 0) call memerror("makespeme 2")
     end if
     speME => nspe_h
     xspe  => nspe
   endif

   case ('T')
   if(it == 1)then
     if(.not.allocated(pspe_obs)) then
        allocate(pspe_obs( nXsd(it) ), stat=aerr)
        if(aerr /= 0) call memerror("makespeme 3")
     end if
     speME => pspe_obs
     xspe  => pspe
   else
     if(.not.allocated(nspe_obs)) then
        allocate(nspe_obs( nXsd(it) ), stat=aerr)
        if(aerr /= 0) call memerror("makespeme 4")
     end if
     speME => nspe_obs
     xspe  => nspe
   endif

   end select
   speME = 0.0
   if(deformed .and. obschar=='H')then
      call makehaikuspes_deformed(-it,lbspe)
      call makehaikuspes_deformed(it,rbspe)
   else
      call makehaikuspes(-it,lbspe,xspe)
      call makehaikuspes(it,rbspe,xspe)
   end if
   call combinehaikuspes(it,lbspe,rbspe,speME)
   deallocate(lbspe,rbspe)
   return
   end subroutine makespeme
!========================================================
!
!
   subroutine makehaikuspes(ith,bspe,xspe)

   use bitstuff
!   use diagh
   use blocks
   use haiku_info
   use spstate

   implicit none
   integer ith
   type (blockhspe), pointer :: bspe(:)
   real, pointer :: xspe(:)

!------------------------------------
   integer iblock
   integer i
   integer isps
   integer j,iword,imask
   integer iorb
   integer, pointer :: hsd(:,:)
   integer :: aerr

   allocate(bspe( hblock(ith)%nhblocks ), stat=aerr)
   if(aerr /= 0) call memerror("makehaikuspes 1")
   do iblock = 1,hblock(ith)%nhblocks
     allocate( bspe( iblock)%hspe ( hblock(ith)%list(iblock)%nhsd ), stat=aerr)
     if(aerr /= 0) call memerror("makehaikuspes 10")
     bspe(iblock)%hspe = 0.0
     do i = 1,hblock(ith)%list(iblock)%nhsd
        hsd => hblock(ith)%list(iblock)%hsd
        do isps = 1,nhsps(ith)
          iword = (isps-1)/max_bit_word+1
          j = mod(isps-1,max_bit_word)+1
          imask = ishft(1,j-1)
          if ( iand(imask,hsd(iword,i)) /=0 ) then
              iorb = hspsqn(ith,isps)%orb
              bspe(iblock)%hspe(i) = bspe(iblock)%hspe(i) + xspe(iorb)
          end if

        enddo ! iop
     enddo ! i
   enddo  ! iblock
   return
   end subroutine makehaikuspes
!
!========================================================

   subroutine combinehaikuspes(it,lbspe,rbspe,speME)

   use sectors
!   use diagh
   use blocks

   implicit none
   integer it
   type (blockhspe), pointer :: lbspe(:),rbspe(:)
   real, pointer :: speME(:)

!------------------------------------
   integer is
   integer iblock
   integer rblock,lblock
   integer(8) :: xsdstart
   integer nradd,nladd
   integer ladd,radd
   integer(8) :: ix
   do is = 1,nsectors(it)
     do iblock = 1,xsd(it)%sector(is)%nhblocks
        rblock = xsd(it)%sector(is)%rhblock(iblock)
        lblock = xsd(it)%sector(is)%lhblock(iblock)
        xsdstart = xsd(it)%sector(is)%blockstart(iblock)-1

!---------------- always make left address the outermost loop (a convention)
        nradd = hblock(it)%list(rblock)%nhsd
        nladd = hblock(-it)%list(lblock)%nhsd

        do ladd = 1,nladd
            do radd = 1,nradd
                ix = xsdstart + radd+(ladd-1)*nradd
                speME(ix) = lbspe(lblock)%hspe(ladd)+rbspe(rblock)%hspe(radd)
            enddo  ! radd
        enddo  ! ladd

     enddo ! iblock

  enddo ! is

  return
  end subroutine combinehaikuspes
!=======================================================
end module diagh

max_stars_repo_name: noctildon/pyBigstick

max_issues_repo_name: noctildon/pyBigstick

max_forks_repo_name: noctildon/pyBigstick


text:
*----------------------------------------------------------------------*
      subroutine set_Zcontrib(flist,ansatz,approx,
     &     irdag,irbreve,
     &     idx_intm,idx_op,nop,op_info)
*----------------------------------------------------------------------*
*     set Z type contributions to B matrix
*       R12 P12 f R12 = R12 [O1+O2-2O1O2] R12tilde
*     skip for HY1
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'def_contraction.h'
      include 'mdef_operator_info.h'
      include 'def_formula_item.h'
      include 'def_formula.h'

      integer, parameter ::
     &     ntest = 00

      type(formula_item), intent(inout), target ::
     &     flist
      character*(*) ::
     &     approx
      integer, intent(in) ::
     &     ansatz,nop,irdag,irbreve,
     &     idx_intm,idx_op(nop)
      type(operator_info), intent(inout) ::
     &     op_info

      type(formula_item), pointer ::
     &     flist_pnt, flist_pnt0

      integer ::
     &     idx_1, idx_2, idx_prj, njoined_intm

      if (approx(1:1).eq.'A') return
      if (approx(4:6).eq.'GBC') return

      if (irdag.gt.nop.or.
     &    irbreve.gt.nop) then
        write(lulog,*) 'idx: ',irdag,irbreve
        write(lulog,*) 'nop: ',nop
        call quit(1,'set_Zcontrib',
     &         'not enough operators on input list')
      end if

      njoined_intm = op_info%op_arr(idx_intm)%op%njoined

      if (idx_op(irdag).le.0.or.
     &    idx_op(irbreve).le.0) then
        write(lulog,*) 'idx: ',idx_op(irdag),idx_op(irbreve)
        call quit(1,'set_Zcontrib',
     &         'operator(s) not on input list')
      end if

      !----------------------------------!
      ! - R12+ Rbreve                    !
      !----------------------------------!
      idx_1 = idx_op(irdag)
      idx_2 = idx_op(irbreve)

      ! go to end of list
      flist_pnt => flist
      do while(associated(flist_pnt%next))
        flist_pnt => flist_pnt%next
      end do
      flist_pnt0 => flist_pnt
      ! set projector
      idx_prj = 5
      if (ansatz.gt.1) idx_prj = 6
      if (njoined_intm.eq.1) then
c        call expand_op_product2(flist_pnt,idx_intm,
c     &       -1d0,4,3,
c     &       (/idx_intm,-idx_1,idx_2,idx_intm/),
c     &       (/1       ,2     ,3    ,1       /),
c     &       -1, -1,
c     &       (/2,3/),1,  ! force additional contraction
c     &       0,0,
c     &       (/2,3,2,idx_prj/),1, ! def. of projector
c     &       .false.,op_info)
c        flist_pnt => flist
c        do while(associated(flist_pnt%next))
c          flist_pnt => flist_pnt%next
c        end do
c        ! this gives the terms with ONLY the projector as contraction
        call expand_op_product2(flist_pnt,idx_intm,
     &       -1d0,4,3,
     &       (/idx_intm,-idx_1,idx_2,idx_intm/),
     &       (/1       ,2     ,3    ,1       /),
     &       -1, -1,
     &       0,0,
     &       0,0,
     &       (/2,3,2,idx_prj/),1, ! def. of projector
     &       .false.,op_info)
      else if (njoined_intm.eq.2) then
        call expand_op_product2(flist_pnt,idx_intm,
     &       -1d0,6,3,
     &       (/idx_intm,-idx_1,idx_intm,idx_intm,idx_2,idx_intm/),
     &       (/1       ,2     ,1       ,1       ,3    ,1       /),
     &       -1, -1,
     &       0,0,
     &       (/2,6, 1,5/),2,    ! avoid cross contrib. to external lines
     &       (/2,5,2,idx_prj/),1, ! def. of projector
     &       .false.,op_info)
      else
        call quit(1,'set_Zcontrib','unexpected: njoined_intm>2')
      end if

      if (ntest.ge.100) then
        write(lulog,*) 'Z contribution'
        call print_form_list(lulog,flist_pnt,op_info)
      end if

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
subroutine OpenFileDom(DataFolder,namda,lun_pdata_dom)
  use   typre
  use   Mod_iofile
  implicit none

  character(150), intent(in) :: DataFolder,namda
  integer(ip), intent(inout)    :: lun_pdata_dom
  character(150) :: fil_pdata_dom, fil_outpu

  fil_pdata_dom = trim(DataFolder)//'/'//adjustl(trim(namda))//'.dom.dat'
  call iofile(zero,lun_pdata_dom,fil_pdata_dom,fil_pdata_dom,'old')

end subroutine OpenFileDom

subroutine CloseFileDom(lun_pdata_dom)
  use typre
  use Mod_iofile
  implicit none

  integer(ip), intent(inout)    :: lun_pdata_dom
  call iofile(two,lun_pdata_dom,'','DOMAIN','old')

end subroutine CloseFileDom

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
! Check whether nested WHEREs work
program where_16
   integer :: a(9)
   integer :: b(9)
   integer :: c(9)

   a = (/0, 0, 0, 1, 1, 1, 2, 2, 2/)
   b = (/0, 1, 2, 0, 1, 2, 0, 1, 2/)
   c = (/0, 0, 0, 0, 0, 0, 0, 0, 0/)

   where (a .eq. 0)
     where (b .eq. 0)
       c = 1
     else where (b .eq. 1)
       c = 2
     else where
       c = 3
     endwhere
   elsewhere (a .eq. 1)
     where (b .eq. 0)
       c = 4
     else where (b .eq. 1)
       c = 5
     else where
       c = 6
     endwhere
   elsewhere
     where (b .eq. 0)
       c = 7
     else where (b .eq. 1)
       c = 8
     else where
       c = 9
     endwhere
   endwhere
   if (any (c .ne. (/1, 2, 3, 4, 5, 6, 7, 8, 9/))) &
     call abort
end program

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
\ ***************************************************************************
\ File: ATMega328p.f -- Definitions for the Atmel328P
\ MCU: ATMega328P
\ ***************************************************************************

-atmega328p
marker -atmega328p

$23 constant PINB  \ PortB Input Register
$24 constant DDRB  \ PortB Dara Direction Register
$25 constant PORTB \ PortB Output Register
$26 constant PINC  \ PortC Input Register
$27 constant DDRC  \ PortC Data Direction Register
$28 constant PORTC \ PortC Output Register
$29 constant PIND  \ PortD Input Register
$30 constant DDRD  \ PortD Data Direction Register
$31 constant PORTD \ PortD Output Register

%10000000 constant PIN7
%01000000 constant PIN6
%00100000 constant PIN5
%00010000 constant PIN4
%00001000 constant PIN3
%00000100 constant PIN2
%00000010 constant PIN1
%00000001 constant PIN0

max_stars_repo_name: Philbywhizz/ff-experiments

max_issues_repo_name: Philbywhizz/ff-experiments

max_forks_repo_name: Philbywhizz/ff-experiments


text:
! { dg-do compile }
!
! PR fortran/30940
program test
implicit none
interface
  subroutine foobar(a)
     character(len=1),dimension(4) :: a
  end subroutine foobar
  subroutine arr(a)
     character(len=1),dimension(1,2,1,2) :: a
  end subroutine arr
end interface

  call foobar( [ "bar" ]) ! { dg-warning "contains too few elements" }
  call foobar( ["ba ","r33"])
  call arr( [ "bar" ]) ! { dg-warning "contains too few elements" }
  call arr( reshape(["b","a","r","3"], [2,2]))
  call arr( reshape(["b","a"], [1,2])) ! { dg-warning "contains too few elements" }
  call arr( reshape(["b","a"], [2,1])) ! { dg-warning "contains too few elements" }
end program test

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


! Compiled with pr83149_a.f90
! { dg-do run }
! { dg-options "-fno-whole-file" }
! { dg-compile-aux-modules "pr83149_a.f90" }
! { dg-additional-sources pr83149_a.f90 }
!
! Contributed by Neil Carlson  <neil.n.carlson@gmail.com>
!
  use mod
  string = 'fubar'
  select case (get_string())
    case ('fubar')
    case default
      stop 1
  end select
end


text:
      SUBROUTINE WCNPC (WCX,WCY, NPCX,NPCY, NPCPWC)

C  WCNPC converts a 2D point in WC to the equivalent point in NPC
C  within view#1 for the picture area, which maps 0:1,0:1 (WC) to
C  PICL:PICR, PICB:PICT (NPC).  It also returns the NPC/WC ratio.

      COMMON /SCRFMT/ PICL,PICR,PICB,PICT, PRTL,PRTR,PRTB,PRTT,
     1                ECHL,ECHR,ECHB,ECHT
      REAL            PICL,PICR,PICB,PICT, PRTL,PRTR,PRTB,PRTT,
     1                ECHL,ECHR,ECHB,ECHT

      REAL    WCX,WCY, NPCX,NPCY, NPCPWC

      NPCX = WCX * (PICR-PICL) + PICL
      NPCY = WCY * (PICT-PICB) + PICB
      NPCPWC = PICT-PICB

      END