text:
great_rsd(kk1,hh1).
great_rsd(k1,ii1).
great_rsd(ee1,ii1).
great_rsd(t1,aa1).
great_rsd(ff1,jj1).
great_rsd(v1,kk1).
great_rsd(f1,a1).
great_rsd(f1,cc1).
great_rsd(q1,ee1).
great_rsd(v1,x1).
great_rsd(w1,ii1).
great_rsd(t1,g1).
great_rsd(v1,h1).
great_rsd(v1,q1).
great_rsd(a1,u1).
great_rsd(p1,ff1).
great_rsd(p1,dd1).
great_rsd(ff1,l1).
great_rsd(m1,h1).
great_rsd(p1,v1).
great_rsd(k1,kk1).
great_rsd(v1,jj1).
great_rsd(a1,z1).
great_rsd(ii1,dd1).
great_rsd(aa1,x1).
great_rsd(u1,hh1).
great_rsd(ff1,u1).
great_rsd(g1,p1).
great_rsd(m1,u1).
great_rsd(z1,aa1).
great_rsd(h1,c1).
great_rsd(k1,t1).

ext: f

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
BMoney is allegedly Billy Durham, according to Greg Webb who witnessed him make a comment on a wikipage and apparently Steven Ostrowski who talked with him Durham on the phone

BMoney used a computer from the ASUCD Presidents office and is of interest to the USA frat. It is believed by Upsilon Sigma Alpha that he is an ASUCD official.

lol wat is this? This is great )

20071030 18:26:52 nbsp USA member Greg Webb knows that BMoney is Billy Durham.   Users/GregWebb

20071030 18:29:51 nbsp No, it IS Billy Durham.  Fact. Users/GregWebb

20071030 18:31:40 nbsp Is of interest? Hes an official? OH SNAP! What does that have to do with anything?  What the heck is wrong with you people?  Run Billy Run! Users/EdWins

20071030 18:33:21 nbsp And the computer was not even in the MU at all. Users/GregWebb

20071030 18:35:48 nbsp So you saw him do it?  Users/RonPayne

20071030 18:39:48 nbsp Yep. Users/GregWebb

20071030 18:41:00 nbsp And does Joe know about this?  Users/RonPayne

20071030 18:41:51 nbsp Cant you guys take this to AIM or save it for the next USA slumber party or something? O.o Users/EdWins
  lol slumber party

20071030 18:42:54 nbsp Probably not. Users/GregWebb

20071030 18:54:11 nbsp He said he did it on the phone and one of our members was right with him.  It seems obvious to us now.  Users/RonPayne

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      logical doipsgen
      integer ipsgen,maxipsgen
      common/ipsgen/doipsgen,ipsgen,maxipsgen
!$omp threadprivate(/ipsgen/)

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


!============================================================================
! CALCULO DE LA ENERGIA DE UNA DE LAS PARTICULAS DE LA CELDA
! POZO CUADRADO (SW)
!
! Autor: Martin Alejandro Paredes Sosa
!============================================================================

Subroutine EnergyPart(Rx1, Ry1, i, V)
  Use cte
  Implicit None
  Real :: V, VNew, Dist, Rx1, Rxd, Ry1, Rz1, Rzd, Ryd
  Integer :: i, j
  !INICIAR ENERGIA EN 0
  V = 0

  BuscarPart: Do j=1, N

     NoLaMisma: If(i .NE. j) Then

        Rxd = Rx1 - X(j)
        Ryd = Ry1 - Y(j)
	Rzd = Rz1 - Z(j)

        !CONDICION DE IMAGEN MINIMA (LOCALIZAR PARTICULAS EN CELDAS CERCANAS)
        Rxd = Rxd - BoxL*Anint(Rxd/BoxL)
        Ryd = Ryd - BoxL*Anint(Ryd/BoxL)
	Rzd = Rzd - BoxL*Anint(Rzd/BoxL)

        !INGRESANDO MODELO DE INTERACCON (DISCOS DUROS)
        Dist = sqrt( Rxd*Rxd + Ryd*Ryd + Rzd*Rzd )
        !If(Dist .LE. 1.0) Write(*,*) Dist, i,j

        ChecarInter: If(Dist .LT. RCut)  Then

           ChecarCercania: If (Dist .LE. 1.0) Then                    ! DESPUES DEL POZO
              VNew = 1.0E+10
           Else If( (Dist .GT. 1.0) .AND. (Dist .LT. Lambda)  ) Then  ! EL POZO
              VNew = -1.0 / TP
           Else                                                       ! ANTES DEL POZO
              VNew = 0
           End If ChecarCercania
           V = V + VNew
        End If ChecarInter

     End If NoLaMisma

  End Do BuscarPart

End Subroutine EnergyPart


text:
! Routine for determining step length bounds. It is used prior to the
! the golden section method for unidirectional search in n-dimensional
! problem.
! Igor Lopes, February 2015
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! ARGUMENTS
!! XO,FO        >   Initial point and corresponding function
!! ALPHL,ALPHU  <   Lower and upper limits
!! XMAX,XMIN    >   Side constraints: limits for domain
!! SDIR         >   Search direction vector
!! NDIM         >   Dimension of the problem (nr design variables)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE FINDBOUNDS(XO,FO,ALPHL,ALPHU,XMAX,XMIN,SDIR,NDIM)
    IMPLICIT NONE
    ! Parameters
    REAL(8) R0 /0.0D0/
    REAL(8) R1 /1.0D0/
    REAL(8) R2 /2.0D0/
    REAL(8) R5 /5.0D0/
    REAL(8) DELTA /1.0D-2/
    REAL(8) INF /1.0D-8/
    INTEGER MITER /50/
    ! Arguments
    INTEGER NDIM
    REAL(8),DIMENSION(NDIM) :: XO, XMAX, XMIN, SDIR
    REAL(8)     :: ALPHL, ALPHU, FO
    ! Locals
    REAL(8),DIMENSION(NDIM) :: XN
    REAL(8)     :: GOLD, ALPH1, AMAX, AMIN, F1, FU, FL, EVALFUNC, AUX
    INTEGER     :: I, NITER
    ! Initialize
    GOLD=(R1+DSQRT(R5))/R2
    AMAX=1.0E9
    AMIN=-1.0E9
    DO I=1,NDIM
        IF(DABS(SDIR(I)).LE.INF)GOTO 5
        IF(SDIR(I).GT.INF)THEN
            AUX=(XMAX(I)-XO(I))/SDIR(I)
            IF(AUX.LT.AMAX) AMAX=AUX
            AUX=(XMIN(I)-XO(I))/SDIR(I)
            IF(AUX.GT.AMIN) AMIN=AUX
        ELSE
            AUX=(XMIN(I)-XO(I))/SDIR(I)
            IF(AUX.LT.AMAX) AMAX=AUX
            AUX=(XMAX(I)-XO(I))/SDIR(I)
            IF(AUX.GT.AMIN) AMIN=AUX
        ENDIF
5       CONTINUE
    ENDDO
    IF(AMAX.LT.AMIN)THEN
        PAUSE
    ENDIF
    ! Begin algorithm
    ALPHU=DELTA*AMAX
    XN=XO+ALPHU*SDIR
    FU=EVALFUNC(XN,NDIM)
    IF(FU.GE.FO)THEN
        ! Upper limit is found... but lower limit is not
        GOTO 10
    ENDIF
    ALPHL=R0
    FL=FO
    DO
        ALPH1=ALPHU
        F1=FU
        ALPHU=ALPH1+(ALPH1-ALPHL)*GOLD
        IF(ALPHU.GT.AMAX) GOTO 15
        DO I=1,NDIM
            XN(I)=XO(I)+ALPHU*SDIR(I)
        ENDDO
        FU=EVALFUNC(XN,NDIM)
        IF(FU.GE.F1)GOTO 99
        ALPHL=ALPH1
        FL=F1
    ENDDO
    ! Find lower limit
10  ALPHL=DELTA*AMIN
    XN=XO+ALPHL*SDIR
    FL=EVALFUNC(XN,NDIM)
    IF(FL.GE.FO)GOTO 99 !XL is found
    DO
        ALPH1=ALPHL
        F1=FL
        ALPHL=ALPHL-(ALPHU-ALPH1)*GOLD
        IF(ALPHL.LT.AMIN)GOTO 20
        DO I=1,NDIM
            XN(I)=XO(I)+ALPHL*SDIR(I)
        ENDDO
        FL=EVALFUNC(XN,NDIM)
        IF(FL.GE.F1)GOTO 99
        ALPHU=ALPH1
        FU=F1
    ENDDO

15  ALPHU=AMAX
    DO I=1,NDIM
        XN(I)=XO(I)+ALPHU*SDIR(I)
    ENDDO
    FU=EVALFUNC(XN,NDIM)
    GOTO 99

20  ALPHL=AMIN
    DO I=1,NDIM
        XN(I)=XO(I)+ALPHL*SDIR(I)
    ENDDO
    FL=EVALFUNC(XN,NDIM)

99  CONTINUE
END SUBROUTINE

max_stars_repo_name: iarlopes/GPROPT

max_issues_repo_name: iarlopes/GPROPT

max_forks_repo_name: iarlopes/GPROPT


module io_module
implicit none

contains

 subroutine read_io_condition(id, file_name, var_name, units, fillval_name, fillval, var_long_name)

 ! This routine read a line of the file IO_CONDITIONS:
 !
 ! The entire line must not be more than 1000 characters long, and must have the format:
 !
 ! 'skipped    file_path    [var_name    [var_units]    [var_fillval_name    [var_fillval    [var_long_name]]]]'
 !
 ! where 'skipped', 'file_path', 'var_name', 'var_units', 'var_fillval_name' and 'var_long_name' are STRINGS,
 ! readable with 'read(..., fmt=*)' (ie, they must not contain string separators, unless the whole string
 ! is between quotes); and 'var_fillval' is FLOAT (double precision, but single is accepted too).
 !
 ! An instance between brackets [] means that it is optional: it will only be read if the subroutine
 ! is asked to. However, if it is the case and the instance is not present, an error will be raised,
 ! with the only expection of 'var_long_name' (a blank string will be returned if the code is unable
 ! to read it, and a warning message will be displayed).
 !
 ! To get a "latter" optional instance (like 'var_fillval'), the previous ones MUST be present, with
 ! the excpetion of 'var_units'.

 use utils, only: read_comment

 integer, intent(in):: id
 character(len=*), intent(out):: file_name
 character(len=*), intent(out), optional:: var_name, units, fillval_name, var_long_name
 double precision, intent(out), optional:: fillval
 character(len=1000) :: line, dummy
 double precision:: dummyfloat
 integer:: ierr

 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
 call read_comment(id)
 read(unit=id,fmt='(A)') line
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

 file_name = ''

 read(line,fmt=*,iostat=ierr) dummy, file_name
 !
 if (ierr/=0) then
   print *
   print *, 'error while reading file name in line:'
   print *, trim(line)
   read(line,fmt=*) dummy, file_name
 end if

 if (present(var_name)) then

   var_name = ''

   read(line,fmt=*,iostat=ierr) dummy, dummy, var_name
   !
   if (ierr/=0) then
     print *
     print *, 'error while reading variable name in line:'
     print *, trim(line)
     read(line,fmt=*) dummy, dummy, var_name
   end if

 end if

 if (present(units)) then

   units = ''
   !
   read(line,fmt=*,iostat=ierr) dummy, dummy, dummy, units
   !
   if (ierr/=0) then
     print *
     print *, 'error while reading variable units in line:'
     print *, trim(line)
     read(line,fmt=*) dummy, dummy, dummy, units
   end if

   if (present(fillval_name)) then

     fillval_name = ''
     !
     read(line,fmt=*,iostat=ierr) dummy, dummy, dummy, dummy, fillval_name
     !
     if (ierr/=0) then
       print *
       print *, 'error while reading variable fillval name in line:'
       print *, trim(line)
       read(line,fmt=*) dummy, dummy, dummy, dummy, fillval_name
     end if

   end if

   if (present(fillval)) then

     read(line,fmt=*,iostat=ierr) dummy, dummy, dummy, dummy, dummy, fillval
     !
     if (ierr/=0) then
       print *
       print *, 'error while reading variable fillvalue in line:'
       print *, trim(line)
       read(line,fmt=*) dummy, dummy, dummy, dummy, dummy, fillval
     end if

   end if

   if (present(var_long_name)) then

     read(line,fmt=*,iostat=ierr) dummy, dummy, dummy, dummy, dummy, dummyfloat, var_long_name
     !
     if (ierr/=0) then
       print *
       print *, 'Warning: unable to read variable long name in line'
       print *, trim(line)
       print *, 'Ignore "long_name" attribute'
       var_long_name = ''
     end if

   end if

 else ! if units not asked, skip it (ie, read one less instance):

   if (present(fillval_name)) then

     fillval_name = ''
     !
     read(line,fmt=*,iostat=ierr) dummy, dummy, dummy, fillval_name
     !
     if (ierr/=0) then
       print *
       print *, 'error while reading variable fillval name in line:'
       print *, trim(line)
       read(line,fmt=*) dummy, dummy, dummy, fillval_name
     end if

   end if

   if (present(fillval)) then

     read(line,fmt=*,iostat=ierr) dummy, dummy, dummy, dummy, fillval
     !
     if (ierr/=0) then
       print *
       print *, 'error while reading variable fillvalue in line:'
       print *, trim(line)
       read(line,fmt=*) dummy, dummy, dummy, dummy, fillval
     end if

   end if

   if (present(var_long_name)) then

     read(line,fmt=*,iostat=ierr) dummy, dummy, dummy, dummy, dummyfloat, var_long_name
     !
     if (ierr/=0) then
       print *
       print *, 'Warning: unable to read variable long name in line'
       print *, trim(line)
       print *, 'Ignore "long_name" attribute'
       var_long_name = ''
     end if

   end if

 end if

 end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 subroutine load_netcdf(filename,varname,VAR, units, missval_name,missval, x_axis_ref, y_axis_ref, expected_units)
 ! load variables only on the first two dimensions, the other dimensions will be loaded
 ! at the end of the dimension. ie: start = (/1,1,n3,n4,.../) ; stop = (/n1,n2,n3,n4,.../)
 ! The 2D-variable is then unravelled in a 1D variable, the first dimension
 ! being the inner-loop unravelled one.
 use netcdf

  character(len=*), intent(in):: filename, varname
  double precision, intent(out), dimension(:):: VAR
  character(len=*), intent(out), optional:: units
  character(len=*), intent(in),  optional:: missval_name
  double precision, intent(out), optional:: missval
  double precision, dimension(:), intent(in), optional:: x_axis_ref, y_axis_ref ! reference axis that current file should match
  character(len=*), intent(in), optional:: expected_units

  character(len=50):: missval_n
  integer:: fileid, varid, ierr, i, l, n
  integer, dimension(:), allocatable:: dimids, dimlen, stt, cnt
  character(len=50), dimension(:), allocatable:: dimname
  double precision, dimension(:), allocatable:: x, y

  integer, dimension(5):: ERROR_HANDLING_OPTION
  common /error/ ERROR_HANDLING_OPTION

!----------------!
! File openning: !
!----------------!

  ierr = nf90_open(filename, NF90_NOWRITE, fileid)
  call nf90check(ierr,'Error while openning input file '//filename)

!-------------------!
! Loading variable: !
!-------------------!

  ! get variable ID and number of dimension:
  ierr = nf90_inq_varid( fileid, varname, varid )
  call nf90check(ierr,'Error while getting identifiers of variable '//varname)
  ierr = nf90_inquire_variable( fileid, varid, ndims=n )
  call nf90check(ierr,'Error while getting number of dimensions of variable '//varname)

  if (n<2) then
    print *
    print *, 'Error while loading variable'//varname
    print *, 'Two dimension at least are expected'
    stop
  end if

  allocate(dimids(n))
  allocate(dimlen(n))
  allocate(stt(n))
  allocate(cnt(n))
  allocate(dimname(n))

  ! get variable shape:
  ierr = nf90_inquire_variable( fileid, varid, dimids=dimids )
  call nf90check(ierr,'Error while getting dimensions identifiers of variable '//varname)
  do i = 1,n
    ierr = nf90_inquire_dimension( fileid, dimids(i), dimname(i), dimlen(i) )
    call nf90check(ierr,'Error while getting dimensions lengths and names of variable '//varname)
  end do

  ! print report:
  print *
  print *, 'netCDF variable: '//varname
  print *, '  - Loaded dimension:'
  print *, '      x1: ',dimname(1)
  print *, '      x2: ',dimname(2)
  if (n>3) print *, '  - unloaded dimension (get end of dim):'
  do i = 3,n
    print *, '        ', dimname(i)
  end do

  ! select data range (except 2nd dimension):
  stt(1) = 1
  stt(3:n) = dimlen(3:n)
  cnt(1) = dimlen(1)
  cnt(2:n)=1

  ! variable loading:
  if (dimlen(1)*dimlen(2) /= size(VAR)) then
    print *
    write(*,'(A)') &
      ' Error while loading variable '//varname
    write(*,'(A,I0,A,I0,A,I0)') &
      ' Inconsistent shape, get: ',dimlen(1),'x',dimlen(2),' = ',dimlen(1)*dimlen(2)
    write(*,'(A,I0)') &
      ' while expected:          ',size(VAR)
    stop
  else
    do i=1,dimlen(2) ! loop on second dimension
      stt(2) = i
      ierr = nf90_get_var( fileid, varid, VAR( 1+dimlen(1)*(i-1) : dimlen(1)*i ), start=stt, count=cnt )
      !         dimension unravelling ---------^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    call nf90check(ierr,'Error while getting variable '//varname)
    end do
  end if

!----------------------------!
! Loading attribute "units": !
!----------------------------!

  if (present(units)) then
    ierr = nf90_get_att( fileid, varid, 'units', units)
    call nf90check(ierr,'Error while getting variable '//varname//' attribute "units"')
  end if

!-------------------------------------!
! Loading attributes "missing-value": !
!-------------------------------------!

  if (present(missval)) then
    if (present(missval_name)) then
      l=len(missval_name)
      missval_n(1:l) = missval_name
    else
      l=10
      missval_n(1:l) = '_FillValue'
    end if
    ierr = nf90_get_att( fileid, varid, missval_n(1:l), missval)
    call nf90check(ierr,'Error while getting variable '//varname//' attribute '//missval_n(1:l))
  end if

!--------------------------------------------!
! Check axis consistency with reference axis !
!--------------------------------------------!

  if (present(x_axis_ref) .and. present(y_axis_ref)) then
      ! Assume: 1st dim = x axis, 2nd dim = y axis
      allocate(x(dimlen(1)))
      allocate(y(dimlen(2)))
      ierr = nf90_inq_varid(fileid, dimname(1), varid)
      call nf90check(ierr,'Warning: unable to get identifier of variable '//trim(dimname(1)), kill=.false.)
      if (ierr==NF90_NOERR) then
          ierr = nf90_get_var(fileid, varid, x)
          call nf90check(ierr,'Error while getting variable '//trim(dimname(1)))
          ierr = nf90_inq_varid(fileid, dimname(2), varid)
          call nf90check(ierr,'Warning: unable to get identifier of variable '//trim(dimname(2)), kill=.false.)
          if (ierr==NF90_NOERR) then
              ierr = nf90_get_var(fileid, varid, y)
              call nf90check(ierr,'Error while getting variable '//trim(dimname(2)))
              call check_axis('', x, y, x_axis_ref, y_axis_ref, ERROR_HANDLING_OPTION(1))
          else
              print *, 'Cannot check axis consistency'
          end if
      else
          print *, 'Cannot check axis consistency'
      end if
      deallocate(x)
      deallocate(y)
  end if

!---------------!
! File closing: !
!---------------!

  ierr = nf90_close(fileid)
  call nf90check(ierr,'Error while closing input file '//filename)

  ! print report:
  if (present(expected_units)) then
      print *, '  - units: '//trim(units)//' (expected: '//trim(expected_units)//')'
  else
      print *, '  - units: '//trim(units)
  end if
  print *

 end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 subroutine load_netcdf_3D(filename,varname,VAR, units, missval_name,missval, expected_units)
 ! load variables only on the first three dimensions, the other dimensions will be loaded
 ! at the end of the dimension. ie: start = (/1,1,1,n4,.../) ; stop = (/n1,n2,n3,n4,.../)
 ! The first two dimensions of the 3D-variable are then unravelled, leading to a 2D variable.
 ! The first dimension being the inner-loop unravelled one.
 ! Then, the order of the 2 remaining dimensions is reversed.
 use netcdf

  character(len=*), intent(in):: filename, varname
  double precision, intent(out), dimension(:,:):: VAR
  character(len=*), intent(out), optional:: units
  character(len=*), intent(in),  optional:: missval_name
  double precision, intent(out), optional:: missval
  character(len=*), intent(in), optional:: expected_units

  character(len=50):: missval_n
  integer:: fileid, varid, ierr, i, j, l, n
  integer, dimension(:), allocatable:: dimids, dimlen, stt, cnt
  character(len=50), dimension(:), allocatable:: dimname

!----------------!
! File openning: !
!----------------!

  ierr = nf90_open(filename, NF90_NOWRITE, fileid)
  call nf90check(ierr,'Error while openning input file '//filename)

!-------------------!
! Loading variable: !
!-------------------!

  ! get variable ID and number of dimension:
  ierr = nf90_inq_varid( fileid, varname, varid )
  call nf90check(ierr,'Error while getting identifiers of variable '//varname)
  ierr = nf90_inquire_variable( fileid, varid, ndims=n )
  call nf90check(ierr,'Error while getting number of dimensions of variable '//varname)

  if (n<3) then
    print *
    print *, 'Error while loading variable '//varname
    print *, 'Three dimension at least are expected'
    stop
  end if

  allocate(dimids(n))
  allocate(dimlen(n))
  allocate(dimname(n))
  allocate(stt(n))
  allocate(cnt(n))

  ! get variable shape:
  ierr = nf90_inquire_variable( fileid, varid, dimids=dimids )
  call nf90check(ierr,'Error while getting dimensions identifiers of variable '//varname)
  do i = 1,n
    ierr = nf90_inquire_dimension( fileid, dimids(i), dimname(i), dimlen(i) )
    call nf90check(ierr,'Error while getting dimensions lengths and names of variable '//varname)
  end do

  ! print report:
  print *, 'netCDF variable: '//varname
  print *, '  - Loaded dimension:'
  print *, '      x1: ',dimname(1)
  print *, '      x2: ',dimname(2)
  print *, '      x3: ',dimname(3)
  if (n>4) print *, '  - unloaded dimension (get end of dim):'
  do i = 4,n
    print *, '        ', dimname(i)
  end do

  ! select data range (except "start" of 2nd and 3rd dimension):
  stt(1) = 1
  cnt(1) = dimlen(1)
  cnt(2:n) = 1
  stt(4:n) = dimlen(4:n)

  ! variable loading:
  if (dimlen(1)*dimlen(2) /= size(VAR,2) .or. dimlen(3) /= size(VAR,1) ) then
    print *
    write(*,'(A)') &
      ' Error while loading variable '//varname
    write(*,'(A,I0,A,I0,A,I0,A,I0,A,I0)') &
      ' Inconsistent shape, get: ',dimlen(1),'x',dimlen(2),'x',dimlen(3),' => ',dimlen(1)*dimlen(2),'x',dimlen(3)
    write(*,'(A,I0,A,I0)') &
      ' while expected:          ',size(VAR,2),'x',size(VAR,1)
    stop
  else
    do i=1,dimlen(2)
      do j=1,dimlen(3)
        stt(2) = i
        stt(3) = j
        !   3rd dimension becomes 1st dimension: v
        ierr = nf90_get_var( fileid, varid, VAR( j , 1+dimlen(1)*(i-1) : dimlen(1)*i ), start=stt, count=cnt )
        !     1st-2nd dimension unravelling ---------^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        call nf90check(ierr,'Error while getting variable '//varname)
      end do
    end do
  end if

!----------------------------!
! Loading attribute "units": !
!----------------------------!

  if (present(units)) then
    ierr = nf90_get_att( fileid, varid, 'units', units)
    call nf90check(ierr,'Error while getting variable '//varname//' attribute "units"')
  end if

!-------------------------------------!
! Loading attributes "missing-value": !
!-------------------------------------!

  if (present(missval)) then
    if (present(missval_name)) then
      l=len(missval_name)
      missval_n(1:l) = missval_name
    else
      l=10
      missval_n(1:l) = '_FillValue'
    end if
    ierr = nf90_get_att( fileid, varid, missval_n(1:l), missval)
    call nf90check(ierr,'Error while getting variable '//varname//' attribute '//missval_n(1:l))
  end if

!---------------!
! File closing: !
!---------------!

  ierr = nf90_close(fileid)
  call nf90check(ierr,'Error while closing input file '//filename)

  ! print report:
  if (present(expected_units)) then
      print *, '  - units: '//trim(units)//' (expected: '//trim(expected_units)//')'
  else
      print *, '  - units: '//trim(units)
  end if
  print *

 end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 subroutine load_netcdf_4D(filename,varname,VAR, units, missval_name,missval, expected_units)
 ! load variables only on the first four dimensions, the other dimensions will be loaded
 ! at the end of the dimension. ie: start = (/1,1,1,1,n5,.../) ; stop = (/n1,n2,n3,n4,n5,.../)
 ! The first two dimensions of the 4D-variable are then unravelled, leading to a 3D variable.
 ! The first dimension being the inner-loop unravelled one.
 ! Then, the order of the 3 remaining dimensions is reversed.
 use netcdf

  character(len=*), intent(in):: filename, varname
  double precision, intent(out), dimension(:,:,:):: VAR
  character(len=*), intent(out), optional:: units
  character(len=*), intent(in),  optional:: missval_name
  double precision, intent(out), optional:: missval
  character(len=*), intent(in), optional:: expected_units

  character(len=50):: missval_n
  integer:: fileid, varid, ierr, i, j, k, l, n
  integer, dimension(:), allocatable:: dimids, dimlen, stt, cnt
  character(len=50), dimension(:), allocatable:: dimname

!----------------!
! File openning: !
!----------------!

  ierr = nf90_open(filename, NF90_NOWRITE, fileid)
  call nf90check(ierr,'Error while openning input file '//filename)

!-------------------!
! Loading variable: !
!-------------------!

  ! get variable ID and number of dimension:
  ierr = nf90_inq_varid( fileid, varname, varid )
  call nf90check(ierr,'Error while getting identifiers of variable '//varname)
  ierr = nf90_inquire_variable( fileid, varid, ndims=n )
  call nf90check(ierr,'Error while getting number of dimensions of variable '//varname)

  if (n<4) then
    print *
    print *, 'Error while loading variable '//varname
    print *, 'Four dimension at least are expected'
    stop
  end if

  allocate(dimids(n))
  allocate(dimlen(n))
  allocate(dimname(n))
  allocate(stt(n))
  allocate(cnt(n))

  ! get variable shape:
  ierr = nf90_inquire_variable( fileid, varid, dimids=dimids )
  call nf90check(ierr,'Error while getting dimensions identifiers of variable '//varname)
  do i = 1,n
    ierr = nf90_inquire_dimension( fileid, dimids(i), dimname(i), dimlen(i) )
    call nf90check(ierr,'Error while getting dimensions lengths and names of variable '//varname)
  end do

  ! print report:
  print *, 'netCDF variable: '//varname
  print *, '  - Loaded dimension:'
  print *, '      x1: ',dimname(1)
  print *, '      x2: ',dimname(2)
  print *, '      x3: ',dimname(3)
  print *, '      x4: ',dimname(4)
  if (n>5) print *, '  - unloaded dimension (get end of dim):'
  do i = 5,n
    print *, '        ', dimname(i)
  end do

  ! select data range (except "start" of 2nd, 3rd and 4th dimension):
  stt(1) = 1
  cnt(1) = dimlen(1)
  cnt(2:n) = 1
  stt(5:n) = dimlen(5:n)

  ! variable loading:
  if (dimlen(1)*dimlen(2) /= size(VAR,3) .or. dimlen(3) /= size(VAR,2) .or. dimlen(4) /= size(VAR,1) ) then
    print *
    write(*,'(A)') &
      ' Error while loading variable '//varname
    write(*,'(A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0)') &
      ' Inconsistent shape, get: ',dimlen(1),'x',dimlen(2),'x',dimlen(3),'x',dimlen(4),&
      ' => ',dimlen(1)*dimlen(2),'x',dimlen(3),'x',dimlen(4)
    write(*,'(A,I0,A,I0,A,I0)') &
      ' while expected:          ',size(VAR,3),'x',size(VAR,2),'x',size(VAR,1)
    stop
  else
    do i=1,dimlen(2)
      do j=1,dimlen(3)
        do k=1,dimlen(4)
          stt(2) = i
          stt(3) = j
          stt(4) = k
          !   3rd dimension becomes 2nd dimension: ----|
          !   4th dimension becomes 1st dimension: v   v
          ierr = nf90_get_var( fileid, varid, VAR( k , j , 1+dimlen(1)*(i-1) : dimlen(1)*i ), start=stt, count=cnt )
          !         1st-2nd dimension unravelling ---------^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          call nf90check(ierr,'Error while getting variable '//varname)
        end do
      end do
    end do
  end if

!----------------------------!
! Loading attribute "units": !
!----------------------------!

  if (present(units)) then
    ierr = nf90_get_att( fileid, varid, 'units', units)
    call nf90check(ierr,'Error while getting variable '//varname//' attribute "units"')
  end if

!-------------------------------------!
! Loading attributes "missing-value": !
!-------------------------------------!

  if (present(missval)) then
    if (present(missval_name)) then
      l=len(missval_name)
      missval_n(1:l) = missval_name
    else
      l=10
      missval_n(1:l) = '_FillValue'
    end if
    ierr = nf90_get_att( fileid, varid, missval_n(1:l), missval)
    call nf90check(ierr,'Error while getting variable '//varname//' attribute '//missval_n(1:l))
  end if

!---------------!
! File closing: !
!---------------!

  ierr = nf90_close(fileid)
  call nf90check(ierr,'Error while closing input file '//filename)

  ! print report:
  if (present(expected_units)) then
      print *, '  - units: '//trim(units)//' (expected: '//trim(expected_units)//')'
  else
      print *, '  - units: '//trim(units)
  end if
  print *

 end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 subroutine load_netcdf_dimvar(filename,varnameX,varnameY,varnameZ,varX,varY,varZ,unitsX,unitsY,unitsZ)

 use netcdf

  character(len=*), intent(in):: filename
  character(len=*), intent(in), optional:: varnameX, varnameY, varnameZ
  double precision, intent(out), dimension(:), optional:: varX, varY, varZ
  character(len=*), intent(out), optional:: unitsX, unitsY, unitsZ
  integer:: fileid, varid, ierr

! File openning:

  ierr = nf90_open(filename, NF90_NOWRITE, fileid)
  call nf90check(ierr,'Error while openning input file '//filename)

!-----------------------------------!
! Loading variables and attributes: !
!-----------------------------------!

  if (present(varnameX)) then
    ierr = nf90_inq_varid( fileid, varnameX, varid )
    call nf90check(ierr,'Error while getting variable '//varnameX//' identifiant')
    if (present(varX)) then
      ierr = nf90_get_var( fileid, varid, varX )
      call nf90check(ierr,'Error while getting variable '//varnameX)
    end if
    if (present(unitsX)) then
      ierr = nf90_get_att( fileid, varid, 'units', unitsX)
      call nf90check(ierr,'Error while getting variable '//varnameX//' attribute "units"')
    end if
  end if
  if (present(varnameY)) then
    ierr = nf90_inq_varid( fileid, varnameY, varid )
    call nf90check(ierr,'Error while getting  variable '//varnameY//' identifiant')
    if (present(varY)) then
      ierr = nf90_get_var( fileid, varid, varY )
      call nf90check(ierr,'Error while getting variable '//varnameY)
    end if
    if (present(unitsY)) then
      ierr = nf90_get_att( fileid, varid, 'units', unitsY)
      call nf90check(ierr,'Error while getting variable '//varnameY//' attribute "units"')
    end if
  end if
  if (present(varnameZ)) then
    ierr = nf90_inq_varid( fileid, varnameZ, varid )
    call nf90check(ierr,'Error while getting  variable '//varnameZ//' identifiant')
    if (present(varZ)) then
      ierr = nf90_get_var( fileid, varid, varZ )
      call nf90check(ierr,'Error while getting variable '//varnameZ)
    end if
    if (present(unitsZ)) then
      ierr = nf90_get_att( fileid, varid, 'units', unitsZ)
      call nf90check(ierr,'Error while getting variable '//varnameZ//' attribute "units"')
    end if
  end if

!---------------!
! File closing: !
!---------------!

  ierr = nf90_close(fileid)
  call nf90check(ierr,'Error while closing input file '//filename)

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine check_landcells_generic(message, erratic_points, landarea, error_handling)

    character(len=*), intent(in):: message
    logical, dimension(:), intent(in):: erratic_points
    double precision, dimension(:), intent(inout):: landarea
    integer, intent(in):: error_handling
    logical, dimension(:), allocatable:: errormask
    double precision:: area_err, tot_landarea
    integer:: nerr, answer
    logical:: loop

    tot_landarea = sum(landarea)

    allocate(errormask(size(landarea)))

    errormask = (landarea>0 .and. erratic_points)

    nerr = count(errormask)
    area_err = sum(landarea, mask=errormask)

    ! If error found:
    if (nerr > 0) then

        print *
        write(*,'(A)')         ' '//trim(message)
        write(*,'(A,I0,A,I0)') '     Number of problematic continent cells:  ', nerr, ' / ', count(landarea>0)
        write(*,'(A,E14.7)')   '     Total area of those cells (km2):        ', 1d6*area_err ! Note: expect area from geoclim in 1e6km2
        write(*,'(A,E14.7)')   '     Which is a fraction of total land area: ', area_err/tot_landarea

        ! Error handling:
        select case (error_handling)
            case (-1) ! ask user interactively
                print *
                loop = .true.
                do while (loop)
                    loop = .false.
                    print *, 'Enter one of the following options:'
                    print *, '    0: abort the program'
                    print *, '    1: remove erratic points (set area=0)'
                    print *, '    2: ignore the issue and continue the execution'
                    read *, answer
                    select case (answer)
                        case (0); stop
                        case (1); where (errormask) landarea=0
                        case (2) ! do nothing
                        case default; loop=.true.
                    end select
                end do

            case (0) ! abort the program
                stop

            case (1) ! automatic correction
                print *
                print *, 'Automatic correction => remove erratic points (set area=0)'
                where (errormask) landarea=0

            case (2) ! do nothing

            case default
                print *
                print *, 'INTERNAL ERROR: illegal error handling option:', error_handling
                stop

        end select

    end if

    deallocate(errormask)

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine check_landcells(varname, fillval, landarea, error_handling, var1D, var2D, axis)

    character(len=*), intent(in):: varname
    double precision, intent(in):: fillval
    double precision, dimension(:), intent(inout):: landarea
    integer, intent(in):: error_handling
    double precision, dimension(:), intent(in), optional:: var1D
    double precision, dimension(:,:), intent(in), optional:: var2D
    integer, intent(in), optional:: axis
    integer:: loc_axis

    ! Default value for non-pixel axis (CO2 levels, lithology, ...)
    if (present(axis)) then
        loc_axis = axis
    else
        loc_axis = 2
        ! Default value => Climate variable convention: 1st axis is pixels, 2nd is CO2 level
    end if

    if (present(var1D)) then
        call check_landcells_generic('WARNING: found missing values on continental cells of variable "'//trim(varname)//'"', &
                                    (var1D==fillval), landarea, error_handling)
    end if

    if (present(var2D)) then
        call check_landcells_generic('WARNING: found missing values on continental cells of variable "'//trim(varname)//'"', &
                                    any(var2D==fillval, dim=loc_axis), landarea, error_handling)
    end if

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine check_invalid(which_var, landarea, error_handling, var1D, var2D, axis)

    character(len=*), intent(in):: which_var
    double precision, dimension(:), intent(inout):: landarea
    integer, intent(in):: error_handling
    double precision, dimension(:), intent(inout), optional:: var1D
    double precision, dimension(:,:), intent(inout), optional:: var2D
    integer, intent(in), optional:: axis
    logical, dimension(:), allocatable:: errormask
    double precision:: area_err, tot_landarea, example_invalid, replacement_value
    integer:: k, nerr, answer, loc_axis
    logical:: loop

    if (.not. (present(var1D) .or. (present(var2D)))) then
        print *
        print *, 'INTERNAL ERROR in subroutine "check_invalid" of module "io_module".'
        print *, 'No variable to check was given to the subroutine.'
        stop
    end if

    ! Default value for non-pixel axis (CO2 levels, lithology, ...)
    if (present(axis)) then
        loc_axis = axis
    else
        loc_axis = 2
        ! Default value => Climate variable convention: 1st axis is pixels, 2nd is CO2 level
    end if

    allocate(errormask(size(landarea)))

    select case (which_var)

        case ('runoff')
            ! Check if negative runoff
            if (present(var1D)) then
                errormask = (landarea>0 .and. var1D<0)
                example_invalid = minval(var1D, mask=errormask)
            else ! var2D
                errormask = (landarea>0 .and. any(var2D<0, dim=loc_axis))
                example_invalid = minval(minval(var2D, dim=loc_axis), mask=errormask)
            end if
            replacement_value = 0d0

        case ('slope')
            ! Check if negative or null slope
            if (present(var1D)) then
                errormask = (landarea>0 .and. var1D<=0)
                example_invalid = minval(var1D, mask=errormask)
                replacement_value = minval(var1D, mask=(.not. errormask))
            else ! var2D
                errormask = (landarea>0 .and. any(var2D<=0, dim=loc_axis))
                example_invalid = minval(minval(var2D, dim=loc_axis), mask=errormask)
                replacement_value = minval(minval(var2D, dim=loc_axis), mask=(.not. errormask))
            end if

        case ('lithology fraction')
            ! Check if sum of lithology is 1 (with 1d-6 accuracy)
            if (present(var1D)) then
                print *
                print *, 'INTERNAL ERROR in subroutine "check_invalid" of module "io_module".'
                print *, 'Unexpected 1-dimension lithology variable (should be nlitho x npixel)'
                stop
            else ! var2D
                errormask = (landarea>0 .and. abs(sum(var2D, dim=loc_axis) - 1d0)>1d-6)
                example_invalid = maxval(abs(sum(var2D, dim=loc_axis) - 1d0), mask=errormask)
            end if

        case default
            print *
            print *, 'INTERNAL ERROR in subroutine "check_invalid" of module "io_module".'
            print *, 'Unkown variable case "'//trim(which_var)//'"'
            stop

        end select

    tot_landarea = sum(landarea)

    nerr = count(errormask)
    area_err = sum(landarea, mask=errormask)

    ! If error found:
    if (nerr > 0) then

        print *
        select case (which_var)
            case ('runoff')
                write(*,'(A)')         ' WARNING: found negative runoff on continental cells.'
            case ('slope')
                write(*,'(A)')         ' WARNING: found null or negative slope on continental cells.'
            case ('lithology fraction')
                write(*,'(A)')         ' WARNING: found continental cells where sum of all lithology fractions is not 1.'
        end select
        write(*,'(A,I0,A,I0)') '     Number of problematic continent cells:  ', nerr, ' / ', count(landarea>0)
        write(*,'(A,E14.7)')   '     Total area of those cells (km2):        ', 1d6*area_err ! Note: expect area from geoclim in 1e6km2
        write(*,'(A,E14.7)')   '     Which is a fraction of total land area: ', area_err/tot_landarea
        if (which_var=='lithology fraction') then
            write(*,'(A,E14.7)')   '     Maximum deviation from 1 found: ', example_invalid
        else
            write(*,'(A,E14.7)')   '     Minimum '//trim(which_var)//' value found: ', example_invalid
        end if

        ! Error handling:
        select case (error_handling)
            case (-1) ! ask user interactively
                print *
                loop = .true.
                do while (loop)
                    loop = .false.
                    print *, 'Enter one of the following options:'
                    print *, '    0: abort the program'
                    print *, '    1: remove erratic points (set area=0)'
                    print *, '    2: ignore the issue and continue the execution'
                    select case (which_var)
                        case ('runoff')
                            print *, '    3: replace invalid values by 0'
                        case ('slope')
                            print *, '    3: replace invalid values by minimum positive slope:', replacement_value
                    end select
                    read *, answer
                    select case (answer)
                        case (0)
                            stop
                        case (1)
                            where (errormask) landarea=0
                        case (2)
                            ! do nothing
                        case (3)
                            if (which_var=='lithology fraction') then
                                loop=.true.
                            else
                                if (present(var1D)) then
                                    where (errormask) var1D=replacement_value
                                else ! var2D
                                    do k = 1,size(var2D,2) ! Runclimber convention: 1st dimension is pixels, 2nd is CO2 levels
                                        where (errormask) var2D(:,k)=replacement_value
                                    end do
                                end if
                            end if
                        case default
                            loop=.true.
                    end select
                end do

            case (0) ! abort the program
                stop

            case (1) ! automatic correction
                print *
                print *, 'Automatic correction => remove erratic points (set area=0)'
                where (errormask) landarea=0

            case (2) ! do nothing

            case (3)
                print *
                select case (which_var)
                    case ('runoff')
                        print *, 'Automatic correction => replace invalid values by 0'
                    case ('slope')
                        print *, 'Automatic correction => replace invalid values minimum positive slope'
                    case ('lithology fraction')
                        print *, 'Bad error handling option "3": No replacement value for invalid lithology fractions sum'
                        stop
                end select
                if (present(var1D)) then
                    where (errormask) var1D=replacement_value
                end if
                if (present(var2D)) then
                    do k = 1,size(var2D,2) ! Runclimber convention: 1st dimension is pixels, 2nd is CO2 levels
                        where (errormask) var2D(:,k)=replacement_value
                    end do
                end if

            case default
                print *
                print *, 'INTERNAL ERROR: illegal error handling option:', error_handling
                stop

        end select

    end if

    deallocate(errormask)

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine raise_axis_error(which_axis, nerr, axis_len, max_mismatch, error_handling)
    character(len=*), intent(in):: which_axis
    integer, intent(in):: nerr, axis_len
    double precision, intent(in):: max_mismatch
    integer, intent(in):: error_handling
    integer:: answer
    logical:: loop

    print *
    write(*,'(A)')         ' WARNING: found mismatch of '//which_axis//' axis with reference axis.'
    write(*,'(A,I0,A,I0)') '     Number of points concerned / axis length:  ', nerr, ' / ', axis_len
    write(*,'(A,E14.7)')   '     Maximum mismatch found (in axis units):    ', max_mismatch

    ! Error handling:
    select case (error_handling)
        case (-1) ! ask user interactively
            print *
            loop = .true.
            do while (loop)
                loop = .false.
                print *, 'Enter one of the following options:'
                print *, '    0: abort the program'
                print *, '    2: ignore the issue and continue the execution'
                read *, answer
                select case (answer)
                    case (0); stop
                    case (2) ! do nothing
                    case default; loop=.true.
                end select
            end do

        case (0) ! abort the program
            stop

        case (1) ! asked for automatic correction
            print *
            print *, 'Bad error handling option "1": cannot remove points with axis mismatch.'
            stop

        case (2) ! do nothing

        case default
            print *
            print *, 'INTERNAL ERROR: illegal error handling option:', error_handling
            stop

    end select

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine check_axis(filename, x, y, xref, yref, error_handling)
    character(len=*), intent(in):: filename
    double precision, dimension(:), intent(in):: x, y, xref, yref
    integer, intent(in):: error_handling
    double precision:: daxis, max_mismatch
    integer:: nerr
    !
    double precision, parameter:: axis_relat_accuracy = 1d-5

    ! x axis
    daxis = abs(x(2) - x(1))
    nerr = count( abs(x - xref)  >  daxis*axis_relat_accuracy )
    max_mismatch = maxval(abs(x - xref))
    if (nerr > 0) then
        call raise_axis_error(trim(filename)//' x', nerr, size(x), max_mismatch, error_handling)
    end if

    ! y axis
    daxis = abs(y(2) - y(1))
    nerr = count( abs(y - yref)  >  daxis*axis_relat_accuracy )
    max_mismatch = maxval(abs(y - yref))
    if (nerr > 0) then
        call raise_axis_error(trim(filename)//' y', nerr, size(y), max_mismatch, error_handling)
    end if

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine nf90check(ierr, message, kill)

    use netcdf
    integer, intent(in):: ierr
    character(len=*), intent(in):: message
    logical, optional, intent(in):: kill
    logical:: loc_kill

    if (present(kill)) then
      loc_kill = kill
    else
      loc_kill = .true.
    end if

    if (ierr/=NF90_NOERR) then
      print *
      print *, message
      print *, nf90_strerror(ierr)
      if (loc_kill) stop
    end if

  end subroutine

end module


text:
! { dg-do compile }
!
! PR 39996: Double typing of function results not detected
!
! Contributed by Janus Weil <janus@gcc.gnu.org>

  interface
    real function A ()
    end function
  end interface
  real :: A  ! { dg-error "already has basic type of" }

  real :: B
  interface
    real function B ()  ! { dg-error "already has basic type of" }
    end function  ! { dg-error "Expecting END INTERFACE statement" }
  end interface

  interface
    function C ()
      real :: C
    end function
  end interface
  real :: C  ! { dg-error "already has basic type of" }

  real :: D
  interface
    function D ()
      real :: D  ! { dg-error "already has basic type of" }
    end function
  end interface

  interface
    function E () result (s)
      real ::s
    end function
  end interface
  real :: E  ! { dg-error "already has basic type of" }

  real :: F
  interface
    function F () result (s)
      real ::s  ! { dg-error "already has basic type of" }
    end function F
  end interface

end

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
SUBROUTINE bee8(bee,coord,xi,eta,det)
!
! Analytical version of the bee matrix for an 8-node quad
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::coord(:,:),xi,eta
 REAL(iwp),INTENT(OUT)::bee(:,:),det
 REAL::x1,x2,x3,x4,x5,x6,x7,x8,y1,y2,y3,y4,y5,y6,y7,y8,xi2,xi3,et2,et3,   &
   xi2et2,xi1et1,xi2et1,xi1et2,xy12,xy13,xy14,xy15,xy16,xy17,xy18,xy23,   &
   xy24,xy25,xy26,xy27,xy28,xy34,xy35,xy36,xy37,xy38,xy45,xy46,xy47,xy48, &
   xy56,xy57,xy58,xy67,xy68,xy78,xcu,xsq,xon,ecu,esq,eon,x2e2,x1e1,x2e1,  &
   x1e2,cons,bot,top,dn1dx,dn2dx,dn3dx,dn4dx,dn5dx,dn6dx,dn7dx,dn8dx,     &
   dn1dy,dn2dy,dn3dy,dn4dy,dn5dy,dn6dy,dn7dy,dn8dy,zero=0.0_iwp,          &
   one=1.0_iwp,two=2.0_iwp,d3=3.0_iwp,d4=4.0_iwp,d5=5.0_iwp,d6=6.0_iwp,   &
   d8=8.0_iwp,pt125=0.125_iwp
 x1=coord(1,1)
 x2=coord(2,1)
 x3=coord(3,1)
 x4=coord(4,1)
 x5=coord(5,1)
 x6=coord(6,1)
 x7=coord(7,1)
 x8=coord(8,1)
 y1=coord(1,2)
 y2=coord(2,2)
 y3=coord(3,2)
 y4=coord(4,2)
 y5=coord(5,2)
 y6=coord(6,2)
 y7=coord(7,2)
 y8=coord(8,2)
 xi2=xi*xi
 xi3=xi2*xi
 et2=eta*eta
 et3=et2*eta
 xi2et2=xi2*et2
 xi1et1=xi*eta
 xi2et1=xi2*eta
 xi1et2=xi*et2
 xy12=x1*y2-y1*x2
 xy13=x1*y3-y1*x3
 xy14=x1*y4-y1*x4
 xy15=x1*y5-y1*x5
 xy16=x1*y6-y1*x6
 xy17=x1*y7-y1*x7
 xy18=x1*y8-y1*x8
 xy23=x2*y3-y2*x3
 xy24=x2*y4-y2*x4
 xy25=x2*y5-y2*x5
 xy26=x2*y6-y2*x6
 xy27=x2*y7-y2*x7
 xy28=x2*y8-y2*x8
 xy34=x3*y4-y3*x4
 xy35=x3*y5-y3*x5
 xy36=x3*y6-y3*x6
 xy37=x3*y7-y3*x7
 xy38=x3*y8-y3*x8
 xy45=x4*y5-y4*x5
 xy46=x4*y6-y4*x6
 xy47=x4*y7-y4*x7
 xy48=x4*y8-y4*x8
 xy56=x5*y6-y5*x6
 xy57=x5*y7-y5*x7
 xy58=x5*y8-y5*x8
 xy67=x6*y7-y6*x7
 xy68=x6*y8-y6*x8
 xy78=x7*y8-y7*x8
 xcu=-d8*xy48+d4*(-xy14+xy38+xy47+xy58)+two*(xy13+xy15-xy37-xy57)
 xsq=two*(-xy13+xy14-xy17+xy18+xy24-xy28-xy34+xy35-xy38-xy45+xy46+xy47-   &
   xy57+xy58+xy68-xy78)+one*(-xy12+xy16-xy23-xy25+xy27-xy36-xy56-xy67)
 xon=d8*xy48+two*(xy14-xy18+xy34-xy38-xy45-xy47-xy58-xy78)+xy12-xy16+xy23-&
   xy25+xy27+xy36-xy56-xy67
 ecu=-d8*xy26+d4*(xy16+xy25+xy36+xy27)+two*(-xy15-xy17-xy35-xy37)
 esq=two*(-xy12+xy13-xy16+xy17-xy23+xy24+xy25-xy27-xy28-xy35+xy36+xy46-   &
   xy56+xy57-xy67+xy68)-xy14+xy18-xy34+xy38-xy45-xy47-xy58-xy78
 eon=d8*xy26+two*( xy12-xy16-xy23-xy25-xy27-xy36-xy56+xy67)+xy14-xy18-    &
   xy34+xy38-xy45+xy47-xy58+xy78
 x2e2=d6*(xy24-xy28+xy46+xy68)+d3*(-xy12+xy13-xy14+xy16-xy17+xy18-xy23-   &
   xy25+xy27-xy34+xy35-xy36+xy38-xy45-xy47-xy56+xy57-xy58-xy67-xy78)
 x1e1=d8*(-xy24-xy28+xy46-xy68)+d6*(-xy12+xy18+xy23+xy34-xy45-xy56+xy67+  &
   xy78)+two*(xy14-xy16+xy25+xy27-xy36+xy38-xy47+xy58)
 x2e1=d8*(xy24+xy28+xy46-xy68)+d5*( xy17-xy18-xy34+xy35-xy45+xy78)+d4*    &
   (xy12-xy16-xy23-xy25-xy27-xy36-xy56+xy67)+d3*(-xy14+xy15+xy37-xy38-    &
   xy47+xy58)
 x1e2=d8*(-xy24+xy28+xy46+xy68)+d5*(xy12-xy13+xy23+xy57-xy56-xy67)+d4*    &
   (xy14-xy18+xy34-xy38-xy45-xy47-xy58-xy78)+d3*(-xy15+xy16+xy25-xy27-    &
   xy36+xy37)
 cons= two*(-xy24+xy28-xy46-xy68)
!
 bot=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 bot=bot+xi2et2*x2e2+xi1et1*x1e1
 bot=bot+xi2et1*x2e1+xi1et2*x1e2+cons
!
 det=pt125*bot
!
 xcu=two*y3-d4*y4+two*y5
 xsq=-y2-two*y3+two*y4+y6-two*y7+two*y8
 xon=y2+two*y4 -y6-two*y8
 ecu=-two*y5+d4*y6-two*y7
 esq=-two*y2+two*y3-y4-two*y6+two*y7+y8
 eon=two*y2+y4-two*y6-y8
 x2e2=-d3*y2+d3*y3-d3*y4+d3*y6-d3*y7+d3*y8
 x1e1=-d6*y2+two*y4-two*y6+d6*y8
 x2e1=d4*y2-d3*y4+d3*y5-d4*y6+d5*y7-d5*y8
 x1e2=d5*y2-d5*y3+d4*y4-d3*y5+d3*y6-d4*y8
 cons=zero
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn1dx=top/bot
!
 xcu=zero
 xsq=y1-y3+two*y4-y5+y7-two*y8
 xon=-y1+y3-y5+y7
 ecu=d4*y5-d8*y6+d4*y7
 esq=two*y1-two*y3+two*y4+two*y5-two*y7-two*y8
 eon=-two*y1-two*y3-two*y5+d8*y6-two*y7
 x2e2= d3*y1-d3*y3+d6*y4-d3*y5+d3*y7-d6*y8
 x1e1= d6*y1+d6*y3-d8*y4+two*y5+two*y7-d8*y8
 x2e1=-d4*y1-d4*y3+d8*y4-d4*y5-d4*y7+d8*y8
 x1e2=-d5*y1+d5*y3-d8*y4+d3*y5-d3*y7+d8*y8
 cons=-two*y4+two*y8
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn2dx=top/bot
!
 xcu=-two*y1-two*y7+d4*y8
 xsq=two*y1+y2-two*y4+two*y5-y6-two*y8
 xon=-y2+two*y4+y6-two*y8
 ecu=-two*y5+d4*y6-two*y7
 esq=-two*y1+two*y2-y4-two*y5+two*y6+y8
 eon=two*y2-y4-two*y6+y8
 x2e2=-d3*y1+d3*y2-d3*y4+d3*y5-d3*y6+d3*y8
 x1e1=-d6*y2+d6*y4-two*y6+two*y8
 x2e1=+d4*y2-d5*y4+d5*y5-d4*y6+d3*y7-d3*y8
 x1e2=d5*y1-d5*y2+d4*y4-d3*y6+d3*y7-d4*y8
 cons=zero
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn3dx=top/bot
!
 xcu=d4*y1+d4*y7-d8*y8
 xsq=-two*y1-two*y2+two*y3-two*y5+two*y6+two*y7
 xon=-two*y1-two*y3-two*y5-two*y7+d8*y8
 ecu=zero
 esq=y1-two*y2+y3-y5+two*y6-y7
 eon=-y1+y3-y5+y7
 x2e2=d3*y1-d6*y2+d3*y3-d3*y5+d6*y6-d3*y7
 x1e1=-two*y1+d8*y2-d6*y3-d6*y5+d8*y6-two*y7
 x2e1=d3*y1-d8*y2+d5*y3-d5*y5+d8*y6-d3*y7
 x1e2=-d4*y1+d8*y2-d4*y3-d4*y5+d8*y6-d4*y7
 cons=two*y2-two*y6
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn4dx=top/bot
!
 xcu=-two*y1-two*y7+d4*y8
 xsq=y2-two*y3+two*y4-y6-two*y7+two*y8
 xon=y2 +two*y4-y6-two*y8
 ecu=two*y1-d4*y2+two*y3
 esq=-two*y2+two*y3+y4-two*y6+two*y7-y8
 eon=two*y2+y4-two*y6-y8
 x2e2=d3*y2-d3*y3+d3*y4-d3*y6+d3*y7-d3*y8
 x1e1=-two*y2+d6*y4-d6*y6+two*y8
 x2e1=-d3*y1+d4*y2-d5*y3+d5*y4-d4*y6+d3*y8
 x1e2=d3*y1-d3*y2+d4*y4-d5*y6+d5*y7-d4*y8
 cons=zero
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn5dx=top/bot
!
 xcu=zero
 xsq=-y1+y3-two*y4+y5-y7+two*y8
 xon=y1-y3+y5-y7
 ecu=-d4*y1+d8*y2-d4*y3
 esq=two*y1-two*y3-two*y4+two*y5-two*y7+two*y8
 eon=two*y1-d8*y2+two*y3+two*y5+two*y7
 x2e2=-d3*y1+d3*y3-d6*y4+d3*y5-d3*y7+d6*y8
 x1e1=two*y1+two*y3-d8*y4+d6*y5+d6*y7-d8*y8
 x2e1=d4*y1+d4*y3-d8*y4+d4*y5+d4*y7-d8*y8
 x1e2=-d3*y1+d3*y3-d8*y4+d5*y5-d5*y7+d8*y8
 cons=two*y4-two*y8
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn6dx=top/bot
!
 xcu=two*y3-d4*y4+two*y5
 xsq=two*y1-y2-two*y4+two*y5+y6-two*y8
 xon=-y2+two*y4+y6-two*y8
 ecu=two*y1-d4*y2+two*y3
 esq=-two*y1+two*y2+y4-two*y5+two*y6-y8
 eon=two*y2-y4-two*y6+y8
 x2e2=d3*y1-d3*y2+d3*y4-d3*y5+d3*y6-d3*y8
 x1e1=-two*y2+two*y4-d6*y6+d6*y8
 x2e1=-d5*y1+d4*y2-d3*y3+d3*y4-d4*y6+d5*y8
 x1e2=d3*y2-d3*y3+d4*y4-d5*y5+d5*y6-d4*y8
 cons=zero
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn7dx=top/bot
!
 xcu=-d4*y3+d8*y4-d4*y5
 xsq=-two*y1+two*y2+two*y3-two*y5-two*y6+two*y7
 xon=two*y1+two*y3-d8*y4+two*y5+two*y7
 ecu=zero
 esq=-y1+two*y2-y3+y5-two*y6+y7
 eon=y1-y3+y5-y7
 x2e2=-d3*y1+d6*y2-d3*y3+d3*y5-d6*y6+d3*y7
 x1e1=-d6*y1+d8*y2-two*y3-two*y5+d8*y6-d6*y7
 x2e1=d5*y1-d8*y2+d3*y3-d3*y5+d8*y6-d5*y7
 x1e2=d4*y1-d8*y2+d4*y3+d4*y5-d8*y6+d4*y7
 cons=-two*y2+two*y6
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn8dx=top/bot
!
 y1=-x1
 y2=-x2
 y3=-x3
 y4=-x4
 y5=-x5
 y6=-x6
 y7=-x7
 y8=-x8
!
 xcu=two*y3-d4*y4+two*y5
 xsq=-y2-two*y3+two*y4+y6-two*y7+two*y8
 xon=y2+two*y4-y6-two*y8
 ecu=-two*y5+d4*y6-two*y7
 esq=-two*y2+two*y3-y4-two*y6+two*y7+y8
 eon=two*y2+y4-two*y6-y8
 x2e2=-d3*y2+d3*y3-d3*y4 +d3*y6-d3*y7+d3*y8
 x1e1=-d6*y2+two*y4-two*y6+d6*y8
 x2e1=d4*y2-d3*y4+d3*y5-d4*y6+d5*y7-d5*y8
 x1e2=d5*y2-d5*y3+d4*y4-d3*y5+d3*y6-d4*y8
 cons=zero
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn1dy=top/bot
!
 xcu=zero
 xsq=y1-y3+two*y4-y5+y7-two*y8
 xon=-y1+y3-y5+y7
 ecu=d4*y5-d8*y6+d4*y7
 esq=two*y1-two*y3+two*y4+two*y5-two*y7-two*y8
 eon=-two*y1-two*y3-two*y5+d8*y6-two*y7
 x2e2=d3*y1-d3*y3+d6*y4-d3*y5+d3*y7-d6*y8
 x1e1=d6*y1+d6*y3-d8*y4+two*y5+two*y7-d8*y8
 x2e1=-d4*y1-d4*y3+d8*y4-d4*y5-d4*y7+d8*y8
 x1e2=-d5*y1+d5*y3-d8*y4+d3*y5-d3*y7+d8*y8
 cons=-two*y4+two*y8
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn2dy=top/bot
!
 xcu=-two*y1-two*y7+d4*y8
 xsq=two*y1+y2-two*y4+two*y5-y6-two*y8
 xon=-y2+two*y4+y6-two*y8
 ecu=-two*y5+d4*y6-two*y7
 esq=-two*y1+two*y2-y4-two*y5+two*y6+y8
 eon=two*y2-y4-two*y6+y8
 x2e2=-d3*y1+d3*y2-d3*y4+d3*y5-d3*y6+d3*y8
 x1e1=-d6*y2+d6*y4-two*y6+two*y8
 x2e1=d4*y2-d5*y4+d5*y5-d4*y6+d3*y7-d3*y8
 x1e2=d5*y1-d5*y2+d4*y4-d3*y6+d3*y7-d4*y8
 cons=zero
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn3dy=top/bot
!
 xcu=d4*y1+d4*y7-d8*y8
 xsq=-two*y1-two*y2+two*y3-two*y5+two*y6+two*y7
 xon=-two*y1-two*y3-two*y5-two*y7+d8*y8
 ecu=zero
 esq=y1-two*y2+y3-y5+two*y6-y7
 eon=-y1+y3-y5+y7
 x2e2=d3*y1-d6*y2+d3*y3-d3*y5+d6*y6-d3*y7
 x1e1=-two*y1+d8*y2-d6*y3-d6*y5+d8*y6-two*y7
 x2e1=d3*y1-d8*y2+d5*y3-d5*y5+d8*y6-d3*y7
 x1e2=-d4*y1+d8*y2-d4*y3-d4*y5+d8*y6-d4*y7
 cons=two*y2-two*y6
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn4dy=top/bot
!
 xcu=-two*y1-two*y7+d4*y8
 xsq=y2-two*y3+two*y4-y6-two*y7+two*y8
 xon=y2+two*y4-y6-two*y8
 ecu=two*y1-d4*y2+two*y3
 esq=-two*y2+two*y3+y4-two*y6+two*y7-y8
 eon=two*y2+y4-two*y6-y8
 x2e2=d3*y2-d3*y3+d3*y4-d3*y6+d3*y7-d3*y8
 x1e1=-two*y2+d6*y4-d6*y6+two*y8
 x2e1=-d3*y1+d4*y2-d5*y3+d5*y4-d4*y6 +d3*y8
 x1e2=d3*y1-d3*y2+d4*y4-d5*y6+d5*y7-d4*y8
 cons=zero
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn5dy=top/bot
!
 xcu=zero
 xsq=-y1+y3-two*y4+y5-y7+two*y8
 xon=y1-y3+y5-y7
 ecu=-d4*y1+d8*y2-d4*y3
 esq=two*y1-two*y3-two*y4+two*y5-two*y7+two*y8
 eon=two*y1-d8*y2+two*y3+two*y5+two*y7
 x2e2=-d3*y1+d3*y3-d6*y4+d3*y5-d3*y7+d6*y8
 x1e1=two*y1+two*y3-d8*y4+d6*y5+d6*y7-d8*y8
 x2e1=d4*y1+d4*y3-d8*y4+d4*y5+d4*y7-d8*y8
 x1e2=-d3*y1+d3*y3-d8*y4+d5*y5-d5*y7+d8*y8
 cons=two*y4-two*y8
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn6dy=top/bot
!
 xcu=two*y3-d4*y4+two*y5
 xsq=two*y1-y2-two*y4+two*y5+y6-two*y8
 xon=-y2+two*y4+y6-two*y8
 ecu=two*y1-d4*y2+two*y3
 esq=-two*y1+two*y2+y4-two*y5+two*y6-y8
 eon=two*y2-y4-two*y6+y8
 x2e2=d3*y1-d3*y2+d3*y4-d3*y5+d3*y6-d3*y8
 x1e1=-two*y2+two*y4-d6*y6+d6*y8
 x2e1=-d5*y1+d4*y2-d3*y3+d3*y4-d4*y6+d5*y8
 x1e2=d3*y2-d3*y3+d4*y4-d5*y5+d5*y6-d4*y8
 cons=zero
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn7dy=top/bot
!
 xcu=-d4*y3+d8*y4-d4*y5
 xsq=-two*y1+two*y2+two*y3-two*y5-two*y6+two*y7
 xon=two*y1+two*y3-d8*y4+two*y5+two*y7
 ecu=zero
 esq=-y1+two*y2-y3+y5-two*y6+y7
 eon=y1-y3+y5-y7
 x2e2=-d3*y1+d6*y2-d3*y3+d3*y5-d6*y6+d3*y7
 x1e1=-d6*y1+d8*y2-two*y3-two*y5+d8*y6-d6*y7
 x2e1=d5*y1-d8*y2+d3*y3-d3*y5+d8*y6-d5*y7
 x1e2=d4*y1-d8*y2+d4*y3+d4*y5-d8*y6+d4*y7
 cons=-two*y2+two*y6
!
 top=xi3*xcu+xi2*xsq+xi*xon+et3*ecu+et2*esq+eta*eon
 top=top+xi2et2*x2e2+xi1et1*x1e1
 top=top+xi2et1*x2e1+xi1et2*x1e2+cons
!
 dn8dy=top/bot
!
 bee=zero
 bee(1,1)=dn1dx
 bee(1,3)=dn2dx
 bee(1,5)=dn3dx
 bee(1,7)=dn4dx
 bee(1,9)=dn5dx
 bee(1,11)=dn6dx
 bee(1,13)=dn7dx
 bee(1,15)=dn8dx
 bee(2,2)=dn1dy
 bee(2,4)=dn2dy
 bee(2,6)=dn3dy
 bee(2,8)=dn4dy
 bee(2,10)=dn5dy
 bee(2,12)=dn6dy
 bee(2,14)=dn7dy
 bee(2,16)=dn8dy
 bee(3,1)=dn1dy
 bee(3,3)=dn2dy
 bee(3,5)=dn3dy
 bee(3,7)=dn4dy
 bee(3,9)=dn5dy
 bee(3,11)=dn6dy
 bee(3,13)=dn7dy
 bee(3,15)=dn8dy
 bee(3,2)=dn1dx
 bee(3,4)=dn2dx
 bee(3,6)=dn3dx
 bee(3,8)=dn4dx
 bee(3,10)=dn5dx
 bee(3,12)=dn6dx
 bee(3,14)=dn7dx
 bee(3,16)=dn8dx
 RETURN
END SUBROUTINE bee8

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


text:
! Test copy-in / copy-out of non-contiguous variable passed as F77 array arguments.
! RUN: bbc -emit-fir %s -o - | FileCheck %s

! Nominal test
! CHECK-LABEL: func @_QPtest_assumed_shape_to_array(
! CHECK-SAME: %[[x:.*]]: !fir.box<!fir.array<?xf32>>{{.*}}) {
subroutine test_assumed_shape_to_array(x)
  real :: x(:)
! Creating temp
! CHECK:  %[[dim:.*]]:3 = fir.box_dims %[[x:.*]], %c0{{.*}} : (!fir.box<!fir.array<?xf32>>, index) -> (index, index, index)
! CHECK:  %[[temp:.*]] = fir.allocmem !fir.array<?xf32>, %[[dim]]#1 {uniq_name = ".copyinout"}

! Copy-in
! CHECK-DAG:  %[[shape:.*]] = fir.shape %[[dim]]#1 : (index) -> !fir.shape<1>
! CHECK-DAG:  %[[temp_load:.*]] = fir.array_load %[[temp]](%[[shape]]) : (!fir.heap<!fir.array<?xf32>>, !fir.shape<1>) -> !fir.array<?xf32>
! CHECK-DAG:  %[[x_load:.*]] = fir.array_load %[[x]] : (!fir.box<!fir.array<?xf32>>) -> !fir.array<?xf32>
! CHECK:  %[[copyin:.*]] = fir.do_loop %[[i:.*]] = %{{.*}} to %{{.*}} step %{{.*}} iter_args(%[[res:.*]] = %[[temp_load]]) -> (!fir.array<?xf32>) {
! CHECK:    %[[fetch:.*]] = fir.array_fetch %[[x_load]], %[[i]] : (!fir.array<?xf32>, index) -> f32
! CHECK:    %[[update:.*]] = fir.array_update %[[res]], %[[fetch]], %[[i]] : (!fir.array<?xf32>, f32, index) -> !fir.array<?xf32>
! CHECK:    fir.result %[[update]] : !fir.array<?xf32>
! CHECK:  }
! CHECK:  fir.array_merge_store %[[temp_load]], %[[copyin:.*]] to %[[temp]] : !fir.array<?xf32>, !fir.array<?xf32>, !fir.heap<!fir.array<?xf32>>

! CHECK:  %[[cast:.*]] = fir.convert %[[temp]] : (!fir.heap<!fir.array<?xf32>>) -> !fir.ref<!fir.array<?xf32>>
! CHECK:  fir.call @_QPbar(%[[cast]]) : (!fir.ref<!fir.array<?xf32>>) -> ()

! Copy-out

! CHECK-DAG:  %[[x_load:.*]] = fir.array_load %[[x]] : (!fir.box<!fir.array<?xf32>>) -> !fir.array<?xf32>
! CHECK-DAG:  %[[shape:.*]] = fir.shape %[[dim]]#1 : (index) -> !fir.shape<1>
! CHECK-DAG:  %[[temp_load:.*]] = fir.array_load %[[temp]](%[[shape]]) : (!fir.heap<!fir.array<?xf32>>, !fir.shape<1>) -> !fir.array<?xf32>
! CHECK:  %[[copyout:.*]] = fir.do_loop %[[i:.*]] = %{{.*}} to %{{.*}} step %{{.*}} iter_args(%[[res:.*]] = %[[x_load]]) -> (!fir.array<?xf32>) {
! CHECK:    %[[fetch:.*]] = fir.array_fetch %[[temp_load]], %[[i]] : (!fir.array<?xf32>, index) -> f32
! CHECK:    %[[update:.*]] = fir.array_update %[[res]], %[[fetch]], %[[i]] : (!fir.array<?xf32>, f32, index) -> !fir.array<?xf32>
! CHECK:    fir.result %[[update]] : !fir.array<?xf32>
! CHECK:  }
! CHECK:  fir.array_merge_store %[[x_load]], %[[copyout:.*]] to %[[x]] : !fir.array<?xf32>, !fir.array<?xf32>, !fir.box<!fir.array<?xf32>>

! CHECK: fir.freemem %[[temp]]

  call bar(x)
end subroutine

! Test that copy-in/copy-out does not trigger the re-evaluation of
! the designator expression.
! CHECK-LABEL: func @_QPeval_expr_only_once(
! CHECK-SAME: %[[x:.*]]: !fir.ref<!fir.array<200xf32>>{{.*}}) {
subroutine eval_expr_only_once(x)
  integer :: only_once
  real :: x(200)
! CHECK: fir.call @_QPonly_once()
! CHECK: %[[x_section:.*]] = fir.embox %[[x]](%{{.*}}) [%{{.*}}] : (!fir.ref<!fir.array<200xf32>>, !fir.shape<1>, !fir.slice<1>) -> !fir.box<!fir.array<?xf32>>
! CHECK: %[[temp:.*]] = fir.allocmem !fir.array<?xf32>
! CHECK-NOT: fir.call @_QPonly_once()
! CHECK:  fir.array_merge_store %{{.*}}, %{{.*}} to %[[temp]]
! CHECK-NOT: fir.call @_QPonly_once()

! CHECK:  %[[cast:.*]] = fir.convert %[[temp]] : (!fir.heap<!fir.array<?xf32>>) -> !fir.ref<!fir.array<?xf32>>
! CHECK:  fir.call @_QPbar(%[[cast]]) : (!fir.ref<!fir.array<?xf32>>) -> ()
  call bar(x(1:200:only_once()))

! CHECK-NOT: fir.call @_QPonly_once()
! CHECK:  fir.array_merge_store %{{.*}}, %{{.*}} to %[[x_section]]
! CHECK-NOT: fir.call @_QPonly_once()
! CHECK: fir.freemem %[[temp]]
end subroutine

! Test no copy-in/copy-out is generated for contiguous assumed shapes.
! CHECK-LABEL: func @_QPtest_contiguous(
! CHECK-SAME: %[[x:.*]]: !fir.box<!fir.array<?xf32>>
subroutine test_contiguous(x)
  real, contiguous :: x(:)
! CHECK: %[[addr:.*]] = fir.box_addr %[[x]] : (!fir.box<!fir.array<?xf32>>) -> !fir.ref<!fir.array<?xf32>>
! CHECK-NOT:  fir.array_merge_store
! CHECK: fir.call @_QPbar(%[[addr]]) : (!fir.ref<!fir.array<?xf32>>) -> ()
  call bar(x)
! CHECK-NOT:  fir.array_merge_store
! CHECK: return
end subroutine

! Test the parenthesis are preventing copy-out.
! CHECK: func @_QPtest_parenthesis(
! CHECK: %[[x:.*]]: !fir.box<!fir.array<?xf32>>{{.*}}) {
subroutine test_parenthesis(x)
  real :: x(:)
! CHECK: %[[dim:.*]]:3 = fir.box_dims %[[x]], %c0{{.*}} : (!fir.box<!fir.array<?xf32>>, index) -> (index, index, index)
! CHECK: %[[temp:.*]] = fir.allocmem !fir.array<?xf32>, %[[dim]]#1 {uniq_name = ".array.expr"}
! CHECK:  fir.array_merge_store %{{.*}}, %{{.*}} to %[[temp]]
! CHECK:  %[[cast:.*]] = fir.convert %[[temp]] : (!fir.heap<!fir.array<?xf32>>) -> !fir.ref<!fir.array<?xf32>>
! CHECK:  fir.call @_QPbar(%[[cast]]) : (!fir.ref<!fir.array<?xf32>>) -> ()
  call bar((x))
! CHECK-NOT:  fir.array_merge_store
! CHECK: fir.freemem %[[temp]]
! CHECK: return
end subroutine

! Test copy-in in is skipped for intent(out) arguments.
! CHECK: func @_QPtest_intent_out(
! CHECK: %[[x:.*]]: !fir.box<!fir.array<?xf32>>{{.*}}) {
subroutine test_intent_out(x)
  real :: x(:)
  interface
  subroutine bar_intent_out(x)
    real, intent(out) :: x(100)
  end subroutine
  end interface
! CHECK: %[[dim:.*]]:3 = fir.box_dims %[[x]], %c0{{.*}} : (!fir.box<!fir.array<?xf32>>, index) -> (index, index, index)
! CHECK: %[[temp:.*]] = fir.allocmem !fir.array<?xf32>, %[[dim]]#1
! CHECK-NOT:  fir.array_merge_store
! CHECK:  %[[cast:.*]] = fir.convert %[[temp]] : (!fir.heap<!fir.array<?xf32>>) -> !fir.ref<!fir.array<100xf32>>
! CHECK:  fir.call @_QPbar_intent_out(%[[cast]]) : (!fir.ref<!fir.array<100xf32>>) -> ()
  call bar_intent_out(x)
! CHECK:  fir.array_merge_store %{{.*}}, %{{.*}} to %[[x]]
! CHECK: fir.freemem %[[temp]]
! CHECK: return
end subroutine

! Test copy-out is skipped for intent(out) arguments.
! CHECK: func @_QPtest_intent_in(
! CHECK: %[[x:.*]]: !fir.box<!fir.array<?xf32>>{{.*}}) {
subroutine test_intent_in(x)
  real :: x(:)
  interface
  subroutine bar_intent_in(x)
    real, intent(in) :: x(100)
  end subroutine
  end interface
! CHECK: %[[dim:.*]]:3 = fir.box_dims %[[x]], %c0{{.*}} : (!fir.box<!fir.array<?xf32>>, index) -> (index, index, index)
! CHECK: %[[temp:.*]] = fir.allocmem !fir.array<?xf32>, %[[dim]]#1
! CHECK:  fir.array_merge_store %{{.*}}, %{{.*}} to %[[temp]]
! CHECK:  %[[cast:.*]] = fir.convert %[[temp]] : (!fir.heap<!fir.array<?xf32>>) -> !fir.ref<!fir.array<100xf32>>
! CHECK:  fir.call @_QPbar_intent_in(%[[cast]]) : (!fir.ref<!fir.array<100xf32>>) -> ()
  call bar_intent_in(x)
! CHECK-NOT:  fir.array_merge_store
! CHECK: fir.freemem %[[temp]]
! CHECK: return
end subroutine

! Test copy-in/copy-out is done for intent(inout)
! CHECK: func @_QPtest_intent_inout(
! CHECK: %[[x:.*]]: !fir.box<!fir.array<?xf32>>{{.*}}) {
subroutine test_intent_inout(x)
  real :: x(:)
  interface
  subroutine bar_intent_inout(x)
    real, intent(inout) :: x(100)
  end subroutine
  end interface
! CHECK: %[[dim:.*]]:3 = fir.box_dims %[[x]], %c0{{.*}} : (!fir.box<!fir.array<?xf32>>, index) -> (index, index, index)
! CHECK: %[[temp:.*]] = fir.allocmem !fir.array<?xf32>, %[[dim]]#1
! CHECK:  fir.array_merge_store %{{.*}}, %{{.*}} to %[[temp]]
! CHECK:  %[[cast:.*]] = fir.convert %[[temp]] : (!fir.heap<!fir.array<?xf32>>) -> !fir.ref<!fir.array<100xf32>>
! CHECK:  fir.call @_QPbar_intent_inout(%[[cast]]) : (!fir.ref<!fir.array<100xf32>>) -> ()
  call bar_intent_inout(x)
! CHECK:  fir.array_merge_store %{{.*}}, %{{.*}} to %[[x]]
! CHECK: fir.freemem %[[temp]]
! CHECK: return
end subroutine

! Test characters are handled correctly
! CHECK-LABEL: func @_QPtest_char(
! CHECK-SAME:    %[[VAL_0:.*]]: !fir.box<!fir.array<?x!fir.char<1,10>>>{{.*}}) {
subroutine test_char(x)
  ! CHECK: %[[VAL_1:.*]] = arith.constant 10 : index
  ! CHECK: %[[VAL_2:.*]] = arith.constant 0 : index
  ! CHECK: %[[VAL_3:.*]]:3 = fir.box_dims %[[VAL_0]], %[[VAL_2]] : (!fir.box<!fir.array<?x!fir.char<1,10>>>, index) -> (index, index, index)
  ! CHECK: %[[VAL_4:.*]] = fir.allocmem !fir.array<?x!fir.char<1,10>>, %[[VAL_3]]#1 {uniq_name = ".copyinout"}
  ! CHECK: %[[VAL_5:.*]] = fir.shape %[[VAL_3]]#1 : (index) -> !fir.shape<1>
  ! CHECK: %[[VAL_6:.*]] = fir.array_load %[[VAL_4]](%[[VAL_5]]) : (!fir.heap<!fir.array<?x!fir.char<1,10>>>, !fir.shape<1>) -> !fir.array<?x!fir.char<1,10>>
  ! CHECK: %[[VAL_7:.*]] = fir.array_load %[[VAL_0]] : (!fir.box<!fir.array<?x!fir.char<1,10>>>) -> !fir.array<?x!fir.char<1,10>>
  ! CHECK: %[[VAL_8:.*]] = arith.constant 1 : index
  ! CHECK: %[[VAL_9:.*]] = arith.constant 0 : index
  ! CHECK: %[[VAL_10:.*]] = arith.subi %[[VAL_3]]#1, %[[VAL_8]] : index
  ! CHECK: %[[VAL_11:.*]] = fir.do_loop %[[VAL_12:.*]] = %[[VAL_9]] to %[[VAL_10]] step %[[VAL_8]] unordered iter_args(%[[VAL_13:.*]] = %[[VAL_6]]) -> (!fir.array<?x!fir.char<1,10>>) {
  ! CHECK: %[[VAL_14:.*]] = fir.array_access %[[VAL_7]], %[[VAL_12]] : (!fir.array<?x!fir.char<1,10>>, index) -> !fir.ref<!fir.char<1,10>>
  ! CHECK: %[[VAL_15:.*]] = fir.array_access %[[VAL_13]], %[[VAL_12]] : (!fir.array<?x!fir.char<1,10>>, index) -> !fir.ref<!fir.char<1,10>>
  ! CHECK: %[[VAL_16:.*]] = arith.constant 10 : index
  ! CHECK: %[[VAL_17:.*]] = arith.constant 1 : i64
  ! CHECK: %[[VAL_18:.*]] = fir.convert %[[VAL_16]] : (index) -> i64
  ! CHECK: %[[VAL_19:.*]] = arith.muli %[[VAL_17]], %[[VAL_18]] : i64
  ! CHECK: %[[VAL_20:.*]] = arith.constant false
  ! CHECK: %[[VAL_21:.*]] = fir.convert %[[VAL_15]] : (!fir.ref<!fir.char<1,10>>) -> !fir.ref<i8>
  ! CHECK: %[[VAL_22:.*]] = fir.convert %[[VAL_14]] : (!fir.ref<!fir.char<1,10>>) -> !fir.ref<i8>
  ! CHECK: fir.call @llvm.memmove.p0i8.p0i8.i64(%[[VAL_21]], %[[VAL_22]], %[[VAL_19]], %[[VAL_20]]) : (!fir.ref<i8>, !fir.ref<i8>, i64, i1) -> ()
  ! CHECK: %[[VAL_23:.*]] = fir.array_amend %[[VAL_13]], %[[VAL_15]] : (!fir.array<?x!fir.char<1,10>>, !fir.ref<!fir.char<1,10>>) -> !fir.array<?x!fir.char<1,10>>
  ! CHECK: fir.result %[[VAL_23]] : !fir.array<?x!fir.char<1,10>>
  ! CHECK: }
  ! CHECK: fir.array_merge_store %[[VAL_6]], %[[VAL_24:.*]] to %[[VAL_4]] : !fir.array<?x!fir.char<1,10>>, !fir.array<?x!fir.char<1,10>>, !fir.heap<!fir.array<?x!fir.char<1,10>>>
  ! CHECK: %[[VAL_25:.*]] = fir.convert %[[VAL_4]] : (!fir.heap<!fir.array<?x!fir.char<1,10>>>) -> !fir.ref<!fir.char<1,?>>
  ! CHECK: %[[VAL_26:.*]] = fir.emboxchar %[[VAL_25]], %[[VAL_1]] : (!fir.ref<!fir.char<1,?>>, index) -> !fir.boxchar<1>
  ! CHECK: fir.call @_QPbar_char(%[[VAL_26]]) : (!fir.boxchar<1>) -> ()
  ! CHECK: %[[VAL_27:.*]] = fir.array_load %[[VAL_0]] : (!fir.box<!fir.array<?x!fir.char<1,10>>>) -> !fir.array<?x!fir.char<1,10>>
  ! CHECK: %[[VAL_28:.*]] = arith.constant 0 : index
  ! CHECK: %[[VAL_29:.*]]:3 = fir.box_dims %[[VAL_0]], %[[VAL_28]] : (!fir.box<!fir.array<?x!fir.char<1,10>>>, index) -> (index, index, index)
  ! CHECK: %[[VAL_30:.*]] = fir.shape %[[VAL_3]]#1 : (index) -> !fir.shape<1>
  ! CHECK: %[[VAL_31:.*]] = fir.array_load %[[VAL_4]](%[[VAL_30]]) : (!fir.heap<!fir.array<?x!fir.char<1,10>>>, !fir.shape<1>) -> !fir.array<?x!fir.char<1,10>>
  ! CHECK: %[[VAL_32:.*]] = arith.constant 1 : index
  ! CHECK: %[[VAL_33:.*]] = arith.constant 0 : index
  ! CHECK: %[[VAL_34:.*]] = arith.subi %[[VAL_29]]#1, %[[VAL_32]] : index
  ! CHECK: %[[VAL_35:.*]] = fir.do_loop %[[VAL_36:.*]] = %[[VAL_33]] to %[[VAL_34]] step %[[VAL_32]] unordered iter_args(%[[VAL_37:.*]] = %[[VAL_27]]) -> (!fir.array<?x!fir.char<1,10>>) {
  ! CHECK: %[[VAL_38:.*]] = fir.array_access %[[VAL_31]], %[[VAL_36]] : (!fir.array<?x!fir.char<1,10>>, index) -> !fir.ref<!fir.char<1,10>>
  ! CHECK: %[[VAL_39:.*]] = fir.array_access %[[VAL_37]], %[[VAL_36]] : (!fir.array<?x!fir.char<1,10>>, index) -> !fir.ref<!fir.char<1,10>>
  ! CHECK: %[[VAL_40:.*]] = arith.constant 10 : index
  ! CHECK: %[[VAL_41:.*]] = arith.constant 1 : i64
  ! CHECK: %[[VAL_42:.*]] = fir.convert %[[VAL_40]] : (index) -> i64
  ! CHECK: %[[VAL_43:.*]] = arith.muli %[[VAL_41]], %[[VAL_42]] : i64
  ! CHECK: %[[VAL_44:.*]] = arith.constant false
  ! CHECK: %[[VAL_45:.*]] = fir.convert %[[VAL_39]] : (!fir.ref<!fir.char<1,10>>) -> !fir.ref<i8>
  ! CHECK: %[[VAL_46:.*]] = fir.convert %[[VAL_38]] : (!fir.ref<!fir.char<1,10>>) -> !fir.ref<i8>
  ! CHECK: fir.call @llvm.memmove.p0i8.p0i8.i64(%[[VAL_45]], %[[VAL_46]], %[[VAL_43]], %[[VAL_44]]) : (!fir.ref<i8>, !fir.ref<i8>, i64, i1) -> ()
  ! CHECK: %[[VAL_47:.*]] = fir.array_amend %[[VAL_37]], %[[VAL_39]] : (!fir.array<?x!fir.char<1,10>>, !fir.ref<!fir.char<1,10>>) -> !fir.array<?x!fir.char<1,10>>
  ! CHECK: fir.result %[[VAL_47]] : !fir.array<?x!fir.char<1,10>>
  ! CHECK: }
  ! CHECK: fir.array_merge_store %[[VAL_27]], %[[VAL_48:.*]] to %[[VAL_0]] : !fir.array<?x!fir.char<1,10>>, !fir.array<?x!fir.char<1,10>>, !fir.box<!fir.array<?x!fir.char<1,10>>>
  ! CHECK: fir.freemem %[[VAL_4]]

  character(10) :: x(:)
  call bar_char(x)
  ! CHECK:         return
  ! CHECK:       }
end subroutine test_char

! CHECK-LABEL: func @_QPtest_scalar_substring_does_no_trigger_copy_inout
! CHECK-SAME: %[[arg0:.*]]: !fir.boxchar<1>
subroutine test_scalar_substring_does_no_trigger_copy_inout(c, i, j)
  character(*) :: c
  integer :: i, j
  ! CHECK: %[[unbox:.*]]:2 = fir.unboxchar %[[arg0]] : (!fir.boxchar<1>) -> (!fir.ref<!fir.char<1,?>>, index)
  ! CHECK: %[[c:.*]] = fir.convert %[[unbox]]#0 : (!fir.ref<!fir.char<1,?>>) -> !fir.ref<!fir.array<?x!fir.char<1>>>
  ! CHECK: %[[coor:.*]] = fir.coordinate_of %[[c]], %{{.*}} : (!fir.ref<!fir.array<?x!fir.char<1>>>, index) -> !fir.ref<!fir.char<1>>
  ! CHECK: %[[substr:.*]] = fir.convert %[[coor]] : (!fir.ref<!fir.char<1>>) -> !fir.ref<!fir.char<1,?>>
  ! CHECK: %[[boxchar:.*]] = fir.emboxchar %[[substr]], %{{.*}} : (!fir.ref<!fir.char<1,?>>, index) -> !fir.boxchar<1>
  ! CHECK: fir.call @_QPbar_char_2(%[[boxchar]]) : (!fir.boxchar<1>) -> ()
  call bar_char_2(c(i:j))
end subroutine

! CHECK-LABEL: func @_QPissue871(
! CHECK-SAME: %[[p:.*]]: !fir.ref<!fir.box<!fir.ptr<!fir.type<_QFissue871Tt{i:i32}>>>>{{.*}})
subroutine issue871(p)
  ! Test passing implicit derived from scalar pointer (no copy-in/out).
  type t
    integer :: i
  end type t
  type(t), pointer :: p
  ! CHECK: %[[box_load:.*]] = fir.load %[[p]]
  ! CHECK: %[[addr:.*]] = fir.box_addr %[[box_load]]
  ! CHECK: %[[cast:.*]] = fir.convert %[[addr]]
  ! CHECK: fir.call @_QPbar_derived(%[[cast]])
  call bar_derived(p)
end subroutine

! CHECK-LABEL: func @_QPissue871_array(
! CHECK-SAME: %[[p:.*]]: !fir.ref<!fir.box<!fir.ptr<!fir.array<?x!fir.type<_QFissue871_arrayTt{i:i32}>>>>>
subroutine issue871_array(p)
  ! Test passing implicit derived from contiguous pointer (no copy-in/out).
  type t
    integer :: i
  end type t
  type(t), pointer, contiguous :: p(:)
  ! CHECK: %[[box_load:.*]] = fir.load %[[p]]
  ! CHECK: %[[addr:.*]] = fir.box_addr %[[box_load]]
  ! CHECK: %[[cast:.*]] = fir.convert %[[addr]]
  ! CHECK: fir.call @_QPbar_derived_array(%[[cast]])
  call bar_derived_array(p)
end subroutine

! CHECK-LABEL: func @_QPwhole_components()
subroutine whole_components()
  ! Test no copy is made for whole components.
  type t
    integer :: i(100)
  end type
  ! CHECK: %[[a:.*]] = fir.alloca !fir.type<_QFwhole_componentsTt{i:!fir.array<100xi32>}>
  type(t) :: a
  ! CHECK: %[[field:.*]] = fir.field_index i, !fir.type<_QFwhole_componentsTt{i:!fir.array<100xi32>}>
  ! CHECK: %[[addr:.*]] = fir.coordinate_of %[[a]], %[[field]] : (!fir.ref<!fir.type<_QFwhole_componentsTt{i:!fir.array<100xi32>}>>, !fir.field) -> !fir.ref<!fir.array<100xi32>>
  ! CHECK: fir.call @_QPbar_integer(%[[addr]]) : (!fir.ref<!fir.array<100xi32>>) -> ()
  call bar_integer(a%i)
end subroutine

! CHECK-LABEL: func @_QPwhole_component_contiguous_pointer()
subroutine whole_component_contiguous_pointer()
  ! Test no copy is made for whole contiguous pointer components.
  type t
    integer, pointer, contiguous :: i(:)
  end type
  ! CHECK: %[[a:.*]] = fir.alloca !fir.type<_QFwhole_component_contiguous_pointerTt{i:!fir.box<!fir.ptr<!fir.array<?xi32>>>}>
  type(t) :: a
  ! CHECK: %[[field:.*]] = fir.field_index i, !fir.type<_QFwhole_component_contiguous_pointerTt{i:!fir.box<!fir.ptr<!fir.array<?xi32>>>}>
  ! CHECK: %[[coor:.*]] = fir.coordinate_of %[[a]], %[[field]] : (!fir.ref<!fir.type<_QFwhole_component_contiguous_pointerTt{i:!fir.box<!fir.ptr<!fir.array<?xi32>>>}>>, !fir.field) -> !fir.ref<!fir.box<!fir.ptr<!fir.array<?xi32>>>>
  ! CHECK: %[[box_load:.*]] = fir.load %[[coor]] : !fir.ref<!fir.box<!fir.ptr<!fir.array<?xi32>>>>
  ! CHECK: %[[addr:.*]] = fir.box_addr %[[box_load]] : (!fir.box<!fir.ptr<!fir.array<?xi32>>>) -> !fir.ptr<!fir.array<?xi32>>
  ! CHECK: %[[cast:.*]] = fir.convert %[[addr]] : (!fir.ptr<!fir.array<?xi32>>) -> !fir.ref<!fir.array<100xi32>>
  ! CHECK: fir.call @_QPbar_integer(%[[cast]]) : (!fir.ref<!fir.array<100xi32>>) -> ()
  call bar_integer(a%i)
end subroutine

! CHECK-LABEL: func @_QPwhole_component_contiguous_char_pointer()
subroutine whole_component_contiguous_char_pointer()
  ! Test no copy is made for whole contiguous character pointer components.
  type t
    character(:), pointer, contiguous :: i(:)
  end type
  ! CHECK: %[[a:.*]] = fir.alloca !fir.type<_QFwhole_component_contiguous_char_pointerTt{i:!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>}>
  type(t) :: a
  ! CHECK: %[[field:.*]] = fir.field_index i, !fir.type<_QFwhole_component_contiguous_char_pointerTt{i:!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>}>
  ! CHECK: %[[coor:.*]] = fir.coordinate_of %0, %1 : (!fir.ref<!fir.type<_QFwhole_component_contiguous_char_pointerTt{i:!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>}>>, !fir.field) -> !fir.ref<!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>>
  ! CHECK: %[[box_load:.*]] = fir.load %[[coor]] : !fir.ref<!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>>
  ! CHECK: %[[addr:.*]] = fir.box_addr %[[box_load]] : (!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>) -> !fir.ptr<!fir.array<?x!fir.char<1,?>>>
  ! CHECK: %[[len:.*]] = fir.box_elesize %[[box_load]] : (!fir.box<!fir.ptr<!fir.array<?x!fir.char<1,?>>>>) -> index
  ! CHECK: %[[cast:.*]] = fir.convert %[[addr]] : (!fir.ptr<!fir.array<?x!fir.char<1,?>>>) -> !fir.ref<!fir.char<1,?>>
  ! CHECK: %[[embox:.*]] = fir.emboxchar %[[cast]], %[[len]] : (!fir.ref<!fir.char<1,?>>, index) -> !fir.boxchar<1>
  ! CHECK: fir.call @_QPbar_char_3(%[[embox]]) : (!fir.boxchar<1>) -> ()
  call bar_char_3(a%i)
end subroutine

max_stars_repo_name: akashhansda/llvm-project

max_issues_repo_name: akashhansda/llvm-project

max_forks_repo_name: akashhansda/llvm-project


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program dart_to_ncommas

!----------------------------------------------------------------------
! purpose: interface between DART and the ncommas model
!
! method: Read DART state vector and overwrite values in a ncommas restart file.
!         If the DART state vector has an 'advance_to_time' present, a
!         file called ncommas_in.DART is created with a time_manager_nml namelist
!         appropriate to advance ncommas to the requested time.
!
!         The dart_to_ncommas_nml namelist setting for advance_time_present
!         determines whether or not the input file has an 'advance_to_time'.
!         Typically, only temporary files like 'assim_model_state_ic' have
!         an 'advance_to_time'.
!
! author: Tim Hoar 25 Jun 09, revised 12 July 2010
!----------------------------------------------------------------------

use        types_mod, only : r8
use    utilities_mod, only : initialize_utilities, finalize_utilities, &
                             find_namelist_in_file, check_namelist_read, &
                             logfileunit, open_file, close_file, &
                             error_handler, E_MSG
use  assim_model_mod, only : open_restart_read, aread_state_restart, close_restart
use time_manager_mod, only : time_type, print_time, print_date, operator(-), get_time
use        model_mod, only : static_init_model, sv_to_restart_file, &
                             get_model_size, get_base_time, get_ncommas_restart_filename

implicit none

! version controlled file description for error handling, do not edit
character(len=256), parameter :: source   = &
   "$URL$"
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate  = "$Date$"

!------------------------------------------------------------------
! The namelist variables
!------------------------------------------------------------------

character (len = 128) :: dart_to_ncommas_input_file = 'dart_restart'
logical               :: advance_time_present       = .false.

namelist /dart_to_ncommas_nml/ dart_to_ncommas_input_file, &
                           advance_time_present

!----------------------------------------------------------------------

character(len=256)    :: ncommas_restart_filename
integer               :: iunit, io, x_size, diff1, diff2
type(time_type)       :: model_time, adv_to_time, base_time
real(r8), allocatable :: statevector(:)
logical               :: verbose              = .FALSE.

!----------------------------------------------------------------------

call initialize_utilities(progname='dart_to_ncommas', output_flag=verbose)

!----------------------------------------------------------------------
! Call model_mod:static_init_model() which reads the ncommas namelists
! to set grid sizes, etc.
!----------------------------------------------------------------------

call static_init_model()

x_size = get_model_size()
allocate(statevector(x_size))

! Read the namelist to get the input filename.

call find_namelist_in_file("input.nml", "dart_to_ncommas_nml", iunit)
read(iunit, nml = dart_to_ncommas_nml, iostat = io)
call check_namelist_read(iunit, io, "dart_to_ncommas_nml")

call get_ncommas_restart_filename( ncommas_restart_filename )

write(*,*)
write(*,'(''dart_to_ncommas:converting DART file '',A, &
      &'' to ncommas restart file '',A)') &
     trim(dart_to_ncommas_input_file), trim(ncommas_restart_filename)

!----------------------------------------------------------------------
! Reads the valid time, the state, and the target time.
!----------------------------------------------------------------------

iunit = open_restart_read(dart_to_ncommas_input_file)

if ( advance_time_present ) then
   call aread_state_restart(model_time, statevector, iunit, adv_to_time)
else
   call aread_state_restart(model_time, statevector, iunit)
endif
call close_restart(iunit)

!----------------------------------------------------------------------
! update the current ncommas state vector
! Convey the amount of time to integrate the model ...
! time_manager_nml: stop_option, stop_count increments
!----------------------------------------------------------------------

call sv_to_restart_file(statevector, ncommas_restart_filename, model_time)

if ( advance_time_present ) then
   base_time = get_base_time(ncommas_restart_filename)
   call get_time((model_time  - base_time), diff1)
   call get_time((adv_to_time - base_time), diff2)
   iunit = open_file('times', action='write')
   write(iunit, '(I8, I8)') diff1, diff2
   call close_file(iunit)
endif

!----------------------------------------------------------------------
! Log what we think we're doing, and exit.
!----------------------------------------------------------------------

call print_date( model_time,'dart_to_ncommas:ncommas  model date')
call print_time( model_time,'dart_to_ncommas:DART model time')
call print_date( model_time,'dart_to_ncommas:ncommas  model date',logfileunit)
call print_time( model_time,'dart_to_ncommas:DART model time',logfileunit)

if ( advance_time_present ) then
call print_time(adv_to_time,'dart_to_ncommas:advance_to time')
call print_date(adv_to_time,'dart_to_ncommas:advance_to date')
call print_time(adv_to_time,'dart_to_ncommas:advance_to time',logfileunit)
call print_date(adv_to_time,'dart_to_ncommas:advance_to date',logfileunit)
endif

call error_handler(E_MSG,'dart_to_ncommas','Finished successfully.',source,revision,revdate)
call finalize_utilities()

end program dart_to_ncommas

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
C$Procedure      SPKS20 ( S/P Kernel, subset, type 20 )

      SUBROUTINE SPKS20 ( HANDLE, BADDR, EADDR, BEGIN, END )
      IMPLICIT NONE

C$ Abstract
C
C     Extract a subset of the data in a SPK segment of type 20
C     into a new segment.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     SPK
C
C$ Keywords
C
C     EPHEMERIS
C
C$ Declarations

      INTEGER               HANDLE
      INTEGER               BADDR
      INTEGER               EADDR
      DOUBLE PRECISION      BEGIN
      DOUBLE PRECISION      END

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     HANDLE     I   Handle of source segment.
C     BADDR      I   Beginning address of source segment.
C     EADDR      I   Ending address of source segment.
C     BEGIN      I   Beginning (initial epoch) of subset.
C     END        I   End (final epoch) of subset.
C
C$ Detailed_Input
C
C     HANDLE,
C     BADDR,
C     EADDR       are the file handle assigned to a SPK file, and the
C                 beginning and ending addresses of a segment within
C                 the file.  Together they determine a complete set of
C                 ephemeris data, from which a subset is to be
C                 extracted.
C
C     BEGIN,
C     END         are the initial and final epochs (ephemeris time)
C                 of the subset to be extracted.
C
C$ Detailed_Output
C
C     None. This routine writes data to the SPK file currently
C     open for write access.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)  Any errors that occur while reading data from the source SPK
C         file will be diagnosed by routines in the call tree of this
C         routine.
C
C     2)  Any errors that occur while writing data to the output SPK
C         file will be diagnosed by routines in the call tree of this
C         routine.
C
C$ Files
C
C     See argument HANDLE.
C
C$ Particulars
C
C     The exact structure of a segment of data type 20 is detailed in
C     the SPK Required Reading file.
C
C$ Examples
C
C     None.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     NAIF Document 168.0, "S- and P- Kernel (SPK) Specification and
C     User's Guide"
C
C$ Author_and_Institution
C
C     N.J. Bachman    (JPL)
C     R.E. Thurman    (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.0, 23-DEC-2013 (NJB) (RET)
C
C-&

C$ Index_Entries
C
C     subset type_20 spk segment
C
C-&

C
C     SPICELIB functions
C
      DOUBLE PRECISION      J2000
      DOUBLE PRECISION      SPD

      LOGICAL               RETURN

C
C     Local parameters
C
      INTEGER               BUFSIZ
      PARAMETER           ( BUFSIZ = 100 )

C
C     Local variables
C
      DOUBLE PRECISION      BTIME
      DOUBLE PRECISION      DATA   ( BUFSIZ )
      DOUBLE PRECISION      DSCALE
      DOUBLE PRECISION      TSCALE
      DOUBLE PRECISION      INITJD
      DOUBLE PRECISION      INITFR
      DOUBLE PRECISION      INTLEN
      DOUBLE PRECISION      INTRVL
      DOUBLE PRECISION      SUBBEG
      DOUBLE PRECISION      SUBIFR
      DOUBLE PRECISION      SUBIJD

      INTEGER               RECSIZ
      INTEGER               NREC
      INTEGER               FIRST
      INTEGER               LAST
      INTEGER               REMAIN
      INTEGER               ADDR
      INTEGER               MOVE

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      END IF

      CALL CHKIN ( 'SPKS20' )

C
C     The segment is made up of a number of logical records, each
C     having the same size, and covering the same length of time.
C
C     We can determine which record to return using the input epoch,
C     the integer and fractional parts of the initial time of the first
C     record's coverage interval, and the length of the interval
C     covered by each record. These constants are located at the end of
C     the segment, along with the size of each logical record and the
C     total number of records.
C
C     For convenience, we'll fetch the segment's distance and time
C     scales in the same call.
C
      CALL DAFGDA ( HANDLE, EADDR-6, EADDR, DATA )

      DSCALE =      DATA( 1 )
      TSCALE =      DATA( 2 )
      INITJD =      DATA( 3 )
      INITFR =      DATA( 4 )
      INTLEN =      DATA( 5 )
      RECSIZ = INT( DATA( 6 ) )
      NREC   = INT( DATA( 7 ) )

      BTIME  = ( ( INITJD - J2000() ) + INITFR ) * SPD()
      INTRVL = INTLEN * SPD()

      FIRST = INT( (BEGIN - BTIME) / INTRVL ) + 1
      FIRST = MAX( 1,  MIN( FIRST, NREC )  )

      LAST  = INT( (END   - BTIME) / INTRVL ) + 1
      LAST  = MAX( 1,  MIN( LAST,  NREC )  )

C
C     The number of records to be moved.
C
      NREC = LAST - FIRST + 1

C
C     We're going to move the data in chunks of BUFSIZ d.p. words.
C     Compute the number of words left to move, the address of the
C     beginning of the records to move, and the number to move this
C     time.
C
      REMAIN = NREC * RECSIZ
      ADDR   = BADDR + ( FIRST - 1 )*RECSIZ
      MOVE   = MIN( BUFSIZ, REMAIN )

      DO WHILE ( REMAIN .GT. 0 )

         CALL DAFGDA ( HANDLE, ADDR, ADDR + MOVE - 1, DATA )
         CALL DAFADA ( DATA, MOVE )
         REMAIN = REMAIN - MOVE
         ADDR   = ADDR   + MOVE
         MOVE   = MIN( BUFSIZ, REMAIN )

      END DO

C
C     That's all the records we have to move. But there are still seven
C     final numbers left to write:
C
C        1)  The distance scale (DSCALE).
C        2)  The time scale (TSCALE).
C        3)  The initial integer Julian date of the start time of the
C            first record.
C        4)  The fractional part of the state time of the first
C            record.
C        5)  The time interval length for each polynomial in days
C            (INTLEN).
C        6)  The record size (RECSIZ).
C        7)  The number of records (NREC).
C
C
C
C     Let SUBBEG be the subset begin time expressed as a TDB Julian
C     date.
C
      SUBBEG =  J2000()  +  (  BTIME + ( FIRST - 1 )*INTRVL  ) / SPD()

      SUBIJD =  AINT( SUBBEG )
      SUBIFR =  SUBBEG - SUBIJD

      DATA( 1 ) = DSCALE
      DATA( 2 ) = TSCALE
      DATA( 3 ) = SUBIJD
      DATA( 4 ) = SUBIFR
      DATA( 5 ) = INTLEN
      DATA( 6 ) = DBLE( RECSIZ )
      DATA( 7 ) = DBLE( NREC   )

      CALL DAFADA ( DATA, 7 )

      CALL CHKOUT ( 'SPKS20' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
Program test_to_str

    Use Kinds, Only: sp,dp
#   if EXTP
    Use Kinds, Only: qp
#   endif
    Use Kinds, Only: is,id,il
    Use FileIO_Utilities, Only: to_str
    Implicit None
    Real(sp) :: rs
    Real(dp) :: rd
#   if EXTP
    Real(qp) :: rq
#   endif
    Integer(is) :: js
    Integer(il) :: jl
    Integer(id) :: jd

    Write(*,'(A)') 'SP:'
    Write(*,'(A)') to_str(0._sp)
    Write(*,'(A)') to_str(TINY(rs))
    Write(*,'(A)') to_str(HUGE(rs))

    Write(*,'(A)') 'DP:'
    Write(*,'(A)') to_str(0._dp)
    Write(*,'(A)') to_str(TINY(rd))
    Write(*,'(A)') to_str(HUGE(rd))

#   if EXTP
    Write(*,'(A)') 'QP:'
    Write(*,'(A)') to_str(0._qp)
    Write(*,'(A)') to_str(TINY(rq))
    Write(*,'(A)') to_str(HUGE(rq))
#   endif

    Write(*,'(A)') 'I2:'
    Write(*,'(A)') to_str(0_is)
    Write(*,'(A)') to_str(-HUGE(js))
    Write(*,'(A)') to_str(HUGE(js))

    Write(*,'(A)') 'I4:'
    Write(*,'(A)') to_str(0_il)
    Write(*,'(A)') to_str(-HUGE(jl))
    Write(*,'(A)') to_str(HUGE(jl))

    Write(*,'(A)') 'I8:'
    Write(*,'(A)') to_str(0_id)
    Write(*,'(A)') to_str(-HUGE(jd))
    Write(*,'(A)') to_str(HUGE(jd))

End Program

max_stars_repo_name: wtdailey/commonTools

max_issues_repo_name: wtdailey/commonTools

max_forks_repo_name: wtdailey/commonTools


    Program dsyevd_example

!     DSYEVD Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_blas_damax_val, &
        nagf_file_print_matrix_real_gen
      Use lapack_interfaces, Only: dsyevd
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=dp), Parameter :: zero = 0.0_dp
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: r
      Integer :: i, ifail, info, k, lda, liwork, lwork, n
      Character (1) :: job, uplo
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), w(:), work(:)
      Integer, Allocatable :: iwork(:)
!     .. Executable Statements ..
      Write (nout, *) 'DSYEVD Example Program Results'
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n
      lda = n
      liwork = 5*n + 3
      lwork = 2*n*n + 6*n + 1
      Allocate (a(lda,n), w(n), work(lwork), iwork(liwork))

!     Read A from data file

      Read (nin, *) uplo
      If (uplo=='U') Then
        Read (nin, *)(a(i,i:n), i=1, n)
      Else If (uplo=='L') Then
        Read (nin, *)(a(i,1:i), i=1, n)
      End If

      Read (nin, *) job

!     Calculate all the eigenvalues and eigenvectors of A
      Call dsyevd(job, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)

      Write (nout, *)
      If (info>0) Then
        Write (nout, *) 'Failure to converge.'
      Else

!       Print eigenvalues and eigenvectors

        Write (nout, *) 'Eigenvalues'
        Write (nout, 100) w(1:n)
        Write (nout, *)
        Flush (nout)

!       Normalize the eigenvectors: largest element positive
        Do i = 1, n
          Call nagf_blas_damax_val(n, a(1,i), 1, k, r)
          If (a(k,i)<zero) Then
            a(1:n, i) = -a(1:n, i)
          End If
        End Do

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call nagf_file_print_matrix_real_gen('General', ' ', n, n, a, lda, &
          'Eigenvectors', ifail)

      End If

100   Format (3X, (8F8.4))
    End Program


text: Address(Golf View Circle) is a residential Culdesacs culdesac in El Macero. It branches off from South El Macero Drive.


module plusU_global

  use communication, only: comm_get_globalinfo

  implicit none
  private
  public :: PLUS_U_ON
  public :: read_Hubbard_parameters
  public :: prep_Hubbard_parameters
  public :: dm_mms_nla
  public :: U_eff
  public :: V_eff

  logical :: PLUS_U_ON = .false.
  complex(8),allocatable :: dm_mms_nla(:,:,:,:,:,:) ! density matrix
  real(8),allocatable :: U_eff(:,:,:)
  real(8),allocatable :: V_eff(:,:)

  integer,allocatable :: orb_info_Hubbard(:,:)
  real(8),allocatable :: U_Hubbard(:)
  real(8),allocatable :: J_Hubbard(:)
  logical :: read_parameters=.false.
  integer :: npid

contains

  subroutine read_Hubbard_parameters
    implicit none
    integer,parameter :: unit=100
    integer :: nlines,i,ngid,nprocs
    character(2) :: unit_of_parameters
    real(8),parameter :: au2eV=27.2116d0
    if ( read_parameters ) return
    call comm_get_globalinfo( ngid, npid, nprocs )
    call check_plusU_ON_OFF
    if ( .not.PLUS_U_ON ) return
    open(100,file="hubbard_param.dat",status="old")
    read(unit,*,END=9) ! the 1st line is assumed to be true/false switch which is already read
    read(unit,*,END=9) unit_of_parameters
    nlines=0
    do
      read(unit,*,END=9)
      nlines=nlines+1
    end do
9   continue
!
    allocate( U_Hubbard(nlines) ); U_Hubbard=0.0d0
    allocate( J_Hubbard(nlines) ); J_Hubbard=0.0d0
    allocate( orb_info_Hubbard(3,nlines) ); orb_info_Hubbard=0
!
    rewind unit
    read(unit,*)  ! skip
    read(unit,*)  ! skip
    if( npid == 0 )then
      write(*,'(1x,2x,2x,3a4,2a12,3x,"(unit is in ",a2,")")') "a","l","n","U","J",unit_of_parameters
    end if
    do i=1,nlines
      read(unit,*) orb_info_Hubbard(1:3,i), U_Hubbard(i), J_Hubbard(i)
      if( npid == 0 )then
        write(*,'(1x,i2,2x,3i4,2f12.5)') i, orb_info_Hubbard(1:3,i), U_Hubbard(i), J_Hubbard(i)
      end if
    end do
    close(unit)
    if( unit_of_parameters=="eV" .or. unit_of_parameters=="ev" )then
      U_Hubbard(1:nlines)=U_Hubbard(1:nlines)/au2eV
      J_Hubbard(1:nlines)=J_Hubbard(1:nlines)/au2eV
    end if
    read_parameters=.true.
  end subroutine read_Hubbard_parameters

  subroutine check_plusU_ON_OFF
    implicit none
    logical :: flag
    inquire( FILE='hubbard_param.dat', EXIST=flag )
    if ( flag ) then
      open(100,file='hubbard_param.dat',status='old')
      read(100,*) PLUS_U_ON
      close(100)
      if( npid == 0 )then
        write(*,*) "hubbard_param.dat is found: PLUS_U_ON=",PLUS_U_ON
        if ( PLUS_U_ON ) write(*,*) "+U calculation is performed."
      end if
    end if
  end subroutine check_plusU_ON_OFF

  subroutine prep_Hubbard_parameters( U_eff )
    implicit none
    real(8),intent(out) :: U_eff(:,0:,:)
    integer :: i,a,l,n
    U_eff=0.0d0
    do i=1,size(J_Hubbard)
      a=orb_info_Hubbard(1,i)
      l=orb_info_Hubbard(2,i)
      n=orb_info_Hubbard(3,i)
      U_eff(a,l,n) = U_Hubbard(i) - J_Hubbard(i)
    end do
  end subroutine prep_Hubbard_parameters

end module plusU_global