      module fmo_mod
      use mqc_general
      use mqc_algebra2
      use mqc_gaussian
      use iso_fortran_env

      contains

      function countNBasis(iAtom,atom2bfMap) Result(nBasis)
!
!     This function returns an (intrinsic) integer giving the number of basis
!     functions, <nBasis>, belonging to atomic center number <iAtom> according
!     to the atom --> basis function map <atom2bfMap).
!
!
      implicit none
      integer(kind=int64),intent(in)::iAtom
      type(MQC_Variable)::atom2bfMap
      integer(kind=int64)::nBasis
      integer(kind=int64)::bf1,bf2
!
      bf1 = INT(MQC_Variable_get_MQC(atom2bfMap,[iAtom]))
      bf2 = INT(MQC_Variable_get_MQC(atom2bfMap,[iAtom+1]))
      nBasis = bf2-bf1
!
      return
      end function countNBasis

      function build_atom2bfMap(nAtoms,bf2AtomMap)  &
        Result(atom2bfMap)
!
!     This function produces an atom --> basis function map. The output array is
!     dimensioned (nAtoms+1) long. Basis functions on atom center number i runs
!     from atom2bfMap(i) to atom2bfMap(i+1)-1.
!
!
      implicit none
      integer(kind=int64),intent(in)::nAtoms
      class(MQC_Variable)::bf2AtomMap
      type(MQC_Variable)::temp
      integer(kind=int64),dimension(:),allocatable::atom2bfMap
      integer(kind=int64),dimension(1)::tempArray
      integer(kind=int64)::i,currentAtom
!
!
!     Allocate memory for atom2bfMap and then fill it.
!
      allocate(atom2bfMap(nAtoms+1))
      atom2bfMap = 0
      atom2bfMap(nAtoms+1) = SIZE(bf2AtomMap)+1
      call bf2AtomMap%print(iout=6,header='bf2AtomMap')
      tempArray(1) = 1
      do i = 1,SIZE(bf2AtomMap)
        tempArray(1) = i
        temp = MQC_Variable_get_MQC(bf2AtomMap,tempArray)
        if(atom2bfMap(INT(temp)).eq.0) atom2bfMap(INT(temp)) = i
      endDo
!
!     In case there are atomic centers without any basis functions -- as there
!     will be in the fragment calculations -- we have to step through
!     <atom2bfMap> and replace zeros.
!
      do i = nAtoms,1,-1
        if(atom2bfMap(i).eq.0) atom2bfMap(i) = atom2bfMap(i+1)
      endDo
!
      return
      end function build_atom2bfMap

      function build_full2fragmentMap(iFragment,atom2fragment,  &
        atom2bfFull,bf2AtomFull,atom2bfFragment,bf2AtomFragment)  &
        Result(bfFull2FragmentMap)
!
!     This function produces a fragment --> full basis set map for a particular
!     fragment (number iFragment). To do thise work, the function needs the
!     fragment number being considered <iFragment>, the fragment list over atoms
!     <atom2fragment>, the full molecule's basis function --> atomic center
!     map <bf2AtomFull>, and the current fragment job's basis function -->
!     atomic center map <bf2AtomFragment>.
!
!
      implicit none
      integer(kind=int64),intent(in)::iFragment
      type(mqc_variable)::atom2fragment,atom2bfFull,bf2AtomFull
      type(mqc_variable)::atom2bfFragment,bf2AtomFragment
      type(mqc_variable)::bfFull2FragmentMap
!
      integer(kind=int64)::i,j,myAtom,myFragment
      integer(kind=int64),dimension(:,:),allocatable::bfFull2FragmentMapTemp
!
!     Do the map work.
!
      allocate(bfFull2FragmentMapTemp(SIZE(bf2AtomFull,1),2))
!hph      bfFull2FragmentMapTemp(:,:) = int(0.0)
      do i = 1,SIZE(bf2AtomFull,1)
        do j = 1,2
          bfFull2FragmentMapTemp(i,j) = 0
        endDo
      endDo
!      write(*,*)' Hrant: ',bfFull2FragmentMapTemp(1,1)
!      call mqc_print(6,bfFull2FragmentMapTemp,header='Hrant: Flag 1')
      write(*,*)
      write(*,*)
      write(*,*)' Hrant - Building frag->full map...'
      write(*,*)'         Fragment = ',iFragment
      write(*,*)
      do myAtom = 1,SIZE(atom2fragment)
        myFragment = INT(MQC_Variable_get_MQC(atom2fragment,[myAtom]))
        write(*,*)'         myAtom = ',myAtom
        write(*,*)'         myFrag = ',myFragment
        if(myFragment.eq.iFragment) then
          j = INT(MQC_Variable_get_MQC(atom2bfFull,[myAtom]))
          write(*,*)'           MATCH! j = ',j,  &
            ' nBasis = ',countNBasis(myAtom,atom2bfFull)
          do i = INT(MQC_Variable_get_MQC(atom2bfFragment,[myAtom])),  &
            INT(MQC_Variable_get_MQC(atom2bfFragment,[myAtom+1]))-1
            bfFull2FragmentMapTemp(j,1) = i
            bfFull2FragmentMapTemp(j,2) = iFragment
            j = j+1
          endDo
        endIf
      endDo
      bfFull2FragmentMap = bfFull2FragmentMapTemp
!
      return
      end function build_full2fragmentMap

      function build_blockFragmentMatrixAO(iFragment,bfFull2Fragment,  &
        fragmentMatrix) Result(blockMatrix)
!
!     This function pushes a matrix in a fragment AO basis (and sized NFragBasis
!     x NFragBasis) out to a block in a full molecule AO basis.
!
!
      implicit none
      integer(kind=int64),intent(in)::iFragment
      type(mqc_variable)::bfFull2Fragment,fragmentMatrix
      type(mqc_variable)::blockMatrix
      real(kind=real64),dimension(:,:),allocatable::blockMatrixTemp
!
      integer(kind=int64)::i,j,nBasisFull,iFrag,jFrag,iFragBF,jFragBF
      real(kind=real64)::temp
!
!     Do the work...
!
      nBasisFull = SIZE(bfFull2Fragment,1)
      allocate(blockMatrixTemp(nBasisFull,nBasisFull))
!hph      blockMatrixTemp = 0.0
      do i = 1,nBasisFull
        do j = 1,nBasisFull
          blockMatrixTemp(i,j) = dfloat(0)
        endDo
      endDo

      write(*,*)
      do i = 1,nBasisFull
        iFrag = INT(MQC_Variable_get_MQC(bfFull2Fragment,[i,2]))
        if(iFrag.eq.iFragment) then
          iFragBF = INT(MQC_Variable_get_MQC(bfFull2Fragment,[i,1]))
          do j = 1,nBasisFull
            jFrag = INT(MQC_Variable_get_MQC(bfFull2Fragment,[j,2]))
            if(jFrag.eq.iFragment) then
              jFragBF = INT(MQC_Variable_get_MQC(bfFull2Fragment,[j,1]))
              temp = float(MQC_Variable_get_MQC(fragmentMatrix, &
                          [iFragBF,jFragBF]))
              write(*,'(3x,''i='',I3,'' j='',I3,'' iFrag='',I3, &
                        '' jFrag='',I3,'' val='',f10.6)') i, &
                        j,iFrag,jFrag,temp
              blockMatrixTemp(i,j) = &
                float(MQC_Variable_get_MQC(fragmentMatrix, &
                [iFragBF,jFragBF]))
            endIf
          endDo
        endIf
      endDo
      write(*,*)
      write(*,*)
      blockMatrix = blockMatrixTemp
!
      return
      end function build_blockFragmentMatrixAO
      end module fmo_mod

!
!     *********************************************************************
!     *********************************************************************
!

        Program IPA_0001
        Use MQC_Gaussian
        Use MQC_Algebra2
        Use MQC_Files
        Use MQC_General
        Use iso_fortran_env

!       Written By: Samantha Bidwell
!       Code Version: 0001
!       Date: 2/25/20

!       Purpose
!               The program is developed to calculate the energy
!               contribution from the mixing of the MO's of three fragments.
!               This is done using the 1st, 2nd and 3rd order energy corrections
!               of perturbation theory.

!       General Information
!               - This code ising gaussian matrix files as its input
!               - The MQC package developed by the Hratchian group is
!                 used to analyze the matrix file.
!               - To compile the code run the associated makefile.

!       General Instructions for Using the Code
!               1. Run an optimization and frequency of the molecule you
!               would like to study with the keyword output=matrix to
!               create the matrix file.
!               2. Run individual jobs for each of fragments with the
!               remaining molecules as ghost atoms also using the
!               keyword output=matrix.
!               3. To complile the code...
!               4. To run the code (after it is compliled) use ./IPA_100
!               molecule.mat frag1.mat frag2.mat frag3.mat

!       Code Outline
!               -- Overall variable declaration
!               -- Read in the matrix files
!               -- Save the fock, overlap, and coefficent matricies and
!               the energy vectors
!               -- Run the basic_information subroutine to calculate the
!               HF_E, the IPM_E, the 2e_E, the 1e_E, and print the
!               number of atoms, basis functions and electrons for each
!               of the files.
!               -- Check the the sum of the fragment jobs is equivilant
!               to the overall molecule.
!               -- Form the composite matricies using the proper basis
!               function mapping technique.
!               -- Transform all of the matricies into the FMO basis.
!               -- Calculate the 1st, 2nd and 3rd order mixing for the
!               fragments
!               -- Print each of the energy corrections and emphasize
!               the orbitals with the largest contributions.

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

!       Variable Declarations
        Implicit None

        Character :: M, F1, F2, F3
        Integer :: iout=6,p,q,r,s,i
        Integer :: NAtom,NBasis,NEle
        Integer :: NAtomTotal,NBasisTotal,NEleTotal
        Integer :: M_NAtom,M_NBasis,M_NEle
        Integer :: F1_NAtom,F1_NBasis,F1_NEle
        Integer :: F2_NAtom,F2_NBasis,F2_NEle
        Integer :: F3_NAtom,F3_NBasis,F3_NEle
        Real :: N
        Type(MQC_Variable) :: temp,temp1,temp2,temp3
        Type(MQC_Variable) :: F_Comp,F_FMO_Comp,F_FMO_Molecule
        Type(MQC_Variable) :: C_Comp,E_Comp
        Type(MQC_Variable) :: First_Order_E,Second_Order_E,Third_Order_E
        Type(MQC_Variable) :: DeltaE,VonE,VVonE,V,VV,VVV,TermOne,TermTwo
        Type(MQC_Variable) :: DoubleDeltaE,EoneEtwo,term1minusterm2
        Type(MQC_Variable) :: M_F,M_S,M_D,M_C,M_H,M_E
        Type(MQC_Variable) :: F1_F,F1_S,F1_D,F1_C,F1_H,F1_E
        Type(MQC_Variable) :: F2_F,F2_S,F2_D,F2_C,F2_H,F2_E
        Type(MQC_Variable) :: F3_F,F3_S,F3_D,F3_C,F3_H,F3_E
        Type(MQC_Gaussian_Unformatted_Matrix_File) :: M_Mat
        Type(MQC_Gaussian_Unformatted_Matrix_File) :: F1_Mat
        Type(MQC_Gaussian_Unformatted_Matrix_File) :: F2_Mat
        Type(MQC_Gaussian_Unformatted_Matrix_File) :: F3_Mat

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

!       Read in each of the command line arguments.
!       Save the Fock, Overlap, Energy, Density and Coefficent matricies
!       from each of the matrix files using the Basic_In formation
!       subroutine.

        Write(*,*) 'STEP: Calling each of the command line arguments &
        and reading in all of the information'
        Write(*,*) 'Reading in Molecule information.'
        Call Get_Command_Argument(1,M)
        Call Basic_Information(M,M_Mat,M_NAtom,M_NBasis,M_NEle, &
             M_F,M,M_D,M_C,M_H,M_E)
        Write(*,*) 'Reading in Fragment_1 information.'
        Call Get_Command_Argument(1,F1)
        Call Basic_Information(F1,F1_Mat,F1_NAtom,F1_NBasis,F1_NEle, &
             F1_F,F1_S,F1_D,F1_C,F1_H,F1_E)
        Write(*,*) 'Reading in Fragment_2 information.'
        Call Get_Command_Argument(1,F2)
        Call Basic_Information(F2,F2_Mat,F2_NAtom,F2_NBasis,F2_NEle, &
             F2_F,F2_S,F2_D,F2_C,F2_H,F2_E)
        Write(*,*) 'Reading in Fragment_3 information.'
        Call Get_Command_Argument(1,F3)
        Call Basic_Information(F3,F3_Mat,F3_NAtom,F3_NBasis,F3_NEle, &
             F3_F,F3_S,F3_D,F3_C,F3_H,F3_E)

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

!       Check that the number of basis functions, atoms and electrons
!       are equal to their expected values.

        Write(*,*) 'STEP: Checking that the number of atoms in each &
        fragment input are equal:'

        NAtomTotal=F1_NAtom+F2_NAtom+F3_NAtom

        If (NAtomTotal.eq.3*M_NAtom) then
         Write(iout,*) 'Atoms are equal. Continuing...'
        Else
         Write(iout,*) 'Atoms not equal. Program terminating.'
         STOP
        EndIf

        NAtom=M_NAtom

        Write(*,*) 'STEP: Checking that the number of basisfunctions &
        in each fragment input are equal:'

        NBasisTotal=F1_NBasis+F2_NBasis+F3_NBasis

        If (NBasisTotal.eq.3*M_NBasis) then
         Write(iout,*) 'Basis functions are equal. Continuing...'
        Else
         Write(iout,*) 'Basis functions not equal. Program terminating.'
         STOP
        EndIf

        NBasis=M_NBasis

        Write(*,*) 'STEP: Checking that the number of electrons &
        in each fragment input are equal:'

        NEleTotal=F1_NEle+F2_NEle+F3_NEle

        If (NEleTotal.eq.3*M_NEle) then
         Write(iout,*) 'Electrons are equal. Continuing...'
        Else
         Write(iout,*) 'Electrons not equal. Program terminating.'
         STOP
        EndIf

        NEle=M_NEle

!       ---------------------------------------------------------------------------
!       ---------------------------------------------------------------------------
!       Form the composite matricies.
        Write(*,*) 'STEP: Forming the composite matricies.'
!       The composite matrices that are formed are: Fock_Composite,
!       Coefficent_Composite, Overlap_Composite and Energy_Composite.
!       ---------------------------------------------------------------------------
!       ---------------------------------------------------------------------------
!       Transform each of the matricies into the FMO basis.
        Write(*,*) 'STEP: Transforming each of the matricies into the &
        fragment molecular orbital (FMO) basis.'
!       These matricies are created to simplify the following energy
!       corrections.
        F_FMO_M =
        F_FMO_C =
        S_FMO_M =
        S_FMO_C =
!       ---------------------------------------------------------------------------
!       ---------------------------------------------------------------------------
!       First Order Energy Correction
        Write(*,*) 'STEP: Calculating the first order energy &
        correction.'
        Call First_Order_E%init(NBasisTotal,0)
        Do p=1,NBasis
          Call First_Order_E%put(F_FMO_Molecule%at(p,p)-F_FMO_Comp%at(p,p),p)
        EndDo
        Write(iout,*) ' '
        Write(iout,*) '**********************************************&
        ************************************************************'
        write(iout,*) 'The First Order Energy Correction'
        Write(iout,*) '**********************************************&
        ************************************************************'
        Call First_Order_E%print(iout,'First order:')
!       ---------------------------------------------------------------------------
!       ---------------------------------------------------------------------------
!       Second Order Energy Correction
        Write(*,*) 'STEP: Calculating the second order energy &
        correction.'
!       ---------------------------------------------------------------------------
!       ---------------------------------------------------------------------------
!       Third Order Energy Correction
        Write(*,*) 'STEP: Calculating the third order energy &
        correction.'
!       ---------------------------------------------------------------------------
!       ---------------------------------------------------------------------------
!       Print the energy corrections and emphasize the mixing that leads
!       to the largest contributions.
!       ---------------------------------------------------------------------------
!       ---------------------------------------------------------------------------
        End Program IPA_0001

!       ****************************************************************************
!       ****************************************************************************
        Subroutine Basic_Information(Commandline,Mat,NAtoms,NBasis, &
        Nele,F,S,D,C,H,E)

        Use MQC_Gaussian
        Use MQC_Algebra2
        Use MQC_Files
        Use MQC_General
        Use iso_fortran_env

!       This subroutine is written to define:
!               1. Restricted or Unrestricted?
!               2. Print natoms, nbasis and nelectrons
!               3. The hartree fock energy
!               4. The IPM energy
!               5. The 1 and 2 electron energy

        Implicit None
        Character(len=256) :: Commandline
        Integer :: iout=6, NAtoms, NBasis, Nele, i
        Type(MQC_Variable) :: HFE, IPME, OneEle, TwoEle
        Type(MQC_Variable) :: F, S, D, C, H, E
        Type(MQC_Variable) :: FD, HD, tempIPM
        Type(MQC_Gaussian_Unformatted_Matrix_File) :: Mat

        Call Mat%load(Commandline)
        Write(*,*) '  '
        Write(*,*) 'SUB-STEP: Loading information from the *.mat &
                    file', Commandline

!       Testing and printing if the job is restricted or unrestricted
        Write(iout,*) 'Is restricted?', Mat%isrestricted()
        Write(iout,*) 'Is unrestricted?', Mat%isunrestricted()

!       Printing the number of ele, atoms and basis from the mat file
!       and saving them.
        Nele=Mat%getval('NElectrons')
        Write(iout,*) 'NElectrons=', Nele
        NAtoms=Mat%getval('NAtoms')
        Write(iout,*) 'NAtoms=', NAtoms
        Nele=Mat%getval('NBasis')
        Write(iout,*) 'NBasis=', NBasis

!       Load and save all need matricies/vectors (F,S,D,C,H,E)
        Call Mat%getarray('alpha fock matrix', &
        mqcVarOut=F)
        Call Mat%getarray('overlap',mqcVarOut=S)
        Call Mat%getarray('alpha density matrix', &
        mqcVarOut=D)
        Call Mat%getarray('alpha mo coefficents', &
        mqcVarOut=C)
        Call Mat%getarray('core hamiltonian alpha', &
        mqcVarOut=H)
        Call Mat%getarray('alpha orbital energies', &
        mqcVarOut=E)

        Write(*,*) '  '
        Write(*,*) 'SUB-STEP: Calculating basic information using &
        information for the *.mat file', Commandline

!       Calculate the Hartree Fock Energy (HF=1/2(<FD>+1/2<HD>))
        Call MQC_Variable_initialize(HFE,'real')
        Call MQC_Variable_initialize(HD,'real')
        Call MQC_Variable_initialize(FD,'real')
        FD=contraction(F,D)
        HD=contraction(H,D)
        HFE=(1/2)*(float(FD)+((1/2)*float(HD)))
        Call HFE%print(iout,'The HF energy is: ')

!       Calculate the IPM Energy
        Call MQC_Variable_initialize(IPME,'real')
        Call MQC_Variable_initialize(tempIPM,'real')

        Do i=1,Nele/2
!hph           tempIPM=E%at(i)
           IPME=float(IPME)+(2*float(tempIPM))
        EndDo
        Call IPME%print(iout,'The IPM energy is: ')

!       Calculate the one-electron energy
        Call MQC_Variable_initialize(OneEle,'real')
        OneEle=contraction(F,D)-contraction(H,D)
        Call OneEle%print(iout,'The One Electron Energy is: ')

!       Calculate the two-electron energy
        Call MQC_Variable_initialize(TwoEle,'real')
        TwoEle=contraction(D,H)
        Call TwoEle%print(iout,'The Two Electron Energy is: ')

        write(*,*) '  '
        write(*,*) 'SUB-STEP: Finishing basic information subroutine.'
        End Subroutine Basic_Information


text:
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>!
!                                 !                                    !
!                                 !   Bojan Niceno                     !
!      Parameter definitions      !   Delft University of Technology   !
!         for all modules         !   Section Heat Transfer            !
!                                 !   niceno@duttwta.wt.tn.tudelft.nl  !
!                                 !                                    !
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>!
module allp_mod

  implicit none

  integer, parameter :: YES =  0
  integer, parameter :: NO  = -1

  integer, parameter ::          &
    MAXL     = 1000,             &
    MAXPRO   = 1024,             & ! max. n. of processors
    INFLOW   =    1,             & ! boundary condition
    WALL     =    2,             & ! boundary condition
    OUTFLOW  =    3,             & ! boundary condition
    SYMMETRY =    4,             & ! boundary condition
    CONVECT  =    5,             & ! boundary condition
    PRESSURE = 2312,             & ! boundary condition
    PERIODIC = 2313,             & ! boundary condition
    BUFFER   = 2311,             & ! boundary condition
                                   ! BUFFER CANNOT BE 11 AS THIS
                                   ! MAKES PROBLEM IF CASE HAS MORE THAN 10 BOUNDARY
                                   ! CONDITIONS. 11th BOUNDARY CONDITION IS CONFUSED
                                   ! FOR BUFFER IN CASE BUFFER = 11.
    WALLFL   =    6                ! boundary condition

  integer, parameter ::          &
    FLUID    =    7,             & ! material state: fluid
    SOLID    =    8,             & ! material state: solid
    CMN_FILE =    7                ! T-FlowS command file (T-FlowS.cmn)

  real, parameter ::                   &
    HUGE       = 1.e+30,               &
    TINY       = 1.e-30,               &
    PI         = 3.14159265359,        &
    ONE_THIRD  = 0.33333333333333333,  &
    TWO_THIRDS = 1.0 - ONE_THIRD

end module

max_stars_repo_name: mhadziabdic/T-FlowS

max_issues_repo_name: mhadziabdic/T-FlowS

max_forks_repo_name: mhadziabdic/T-FlowS


text:
! { dg-do compile }
! This checks the fix for PR33295 in which the A_type in initA was
! not promoted to module level and so not recognised as being the
! same as that emanating directly from module a.
!
! Contributed by Janus Weil <jaydub66@gmail.com>
!
module A
  type A_type
    real comp
  end type
end module A

module B
contains
  function initA()
    use A
    implicit none
    type(A_type):: initA
    initA%comp=1.0
  end function
end module B

program C
  use B
  use A
  implicit none
  type(A_type):: A_var
  A_var = initA()
end program C

! { dg-final { cleanup-modules "A B" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
c     nilsson model
      implicit real*8(a-h,o-z)
      parameter (km=100)
      dimension ham(km,km), eval(km), evec(km,km)
      dimension ial(km),ialam(km),iasig(km)
      dimension delta_list(60), fes_list(60, 6), temp_list(6)
      dimension spe(100,364),vk2(364)
      external secant, chem_pota, fes_calc
      dimension rkappa(0:11),rmu(0:11)
      common /param/ ei(364), temp, np
      data (rkappa(i),i=0,7) /8*0.05D0/
      data (rmu(i),i=0,5) /0.0D0,0.0D0,0.0D0,0.35D0,0.625D0,0.630D0 /
      data (rmu(i),i=6,7) /0.448D0,0.434D0 /
      data (temp_list(i), i=1,6) /0.0d0,1.0d0,2.0d0,3.0d0,4.0d0,5.0d0/
      ip=40
      in=44
      ra=dfloat(ip+in)
      icol=0
      do n=0,14,2  !all quantum numbers multiplied by 2
        do iom=1,n+1,2
          nbas=0
          do l=n,0,-4
            do lam=-l,l,2
              isig=iom-lam
              if(abs(isig).ne.1)cycle
              nbas=nbas+1
              ial(nbas)=l
              ialam(nbas)=lam
              iasig(nbas)=isig
            enddo
          enddo
          idef = 0
          do delta=-0.3d0, 0.301d0, 0.01d0

            fdel=((1.d0+(2.d0/3.d0)*delta)**2.d0
     &            *(1.d0-(4.d0/3.d0)*delta))**(-1.d0/6.d0)
            hw00=41.d0*ra**(-1.d0/3.d0)
            hw0=hw00*fdel
            c=-2.0D0*hw00*rkappa(n/2)
            d=rmu(n/2)*c/2.0D0
            do i=1,nbas
              l=ial(i)
              lam=ialam(i)
              isig=iasig(i)
              do j=i,nbas
                l1=ial(j)
                lam1=ialam(j)
                isig1=iasig(j)

                h00=0.d0
                hl2=0.d0
                hls=0.d0
                hr2=0.d0
                hy20=0.d0
                if(i.eq.j)then  !Diagonal matrix elements
                  h00=dfloat(n+3)/2.d0*hw0
                  hl2=dfloat(l*(l+2))/4.d0
                endif

                if(l1.eq.l)then
                  if((lam1.eq.lam+2).and.(isig1.eq.isig-2))then
                    hls=dsqrt(dfloat((l-lam)*(l+lam+2)))/4.d0
                  elseif((lam1.eq.lam).and.(isig1.eq.isig))then
                    hls=lam*isig/4.d0
                  elseif((lam1.eq.lam-2).and.(isig1.eq.isig+2))then
                    hls=dsqrt(dfloat((l+lam)*(l-lam+2)))/4.d0
                  endif
                  hr2=dfloat(n+3)/2.d0
                elseif(l1.eq.l-4)then
                  hr2=dsqrt(dfloat((n-l+4)*(n+l+2)))/2.d0
                elseif(l1.eq.l+4)then
                  hr2=dsqrt(dfloat((n-l)*(n+l+6)))/2.d0
                endif
                if((abs(hr2).gt.1e-5).and.(lam1.eq.lam))then
                  hy20=dsqrt(dfloat(l+1)/dfloat(l1+1))
     &                 *cgc(l,4,l1,lam,0)*cgc(l,4,l1,0,0)
                endif
                hdel=-delta*hw0*dsqrt(4.d0)/3.d0*hr2*hy20

                ham(i,j)=h00+hdel+c*hls+d*hl2
                ham(j,i)=ham(i,j)
              enddo
            enddo

            call diag(ham,km,nbas,eval,evec)
            idef=idef+1
            delta_list(idef) = delta
            do i=1,nbas
              spe(idef,icol+i)=eval(i)
            enddo
          enddo
          icol=icol+nbas
        enddo
      enddo

      open(unit=21, file='ef_delta.txt')
      open(unit=75, file="s_data.txt")
      do idef=1,60
         do itemp=1,6
          call sort(spe(idef,:),1000,icol)
          temp = (itemp-1)
          fes_list(idef, itemp) = fes_calc(ip, in, spe(idef,:),temp)
        enddo
      enddo
      close(21)
      open(unit=99, file="pes_data.txt")
      do itemp=1,6
         fe0 = fes_list(31, itemp)
         do idef=2,60
            fes_list(idef,itemp)=fes_list(idef,itemp)-fe0
c$$$            write(99, '(2f12.6)') delta_list(idef), fes_list(idef, itemp)
         enddo
      enddo
      write(99, '(A,6f5.2)') "Delta", temp_list(:)
      do idef=2,60
         write(99, '(7f12.6)') delta_list(idef), fes_list(idef, :)
      end do
      close(99)
      close(75)

c$$$      For neutrons

c$$$      do i =1,364
c$$$         ei(i) = spe(idef,i)
c$$$      end do
c$$$
c$$$      np=44
c$$$      ef0 = ei(np/2)
c$$$      write(*,*) ef0
c$$$      open(unit=45, file='hot_nuc.txt')
c$$$      do i=1,364
c$$$         if (ei(i) .eq. 0) exit
c$$$         fi = fermi_dist(ei(i), ef_root, temp)
c$$$         write(45,*) ei(i), fi
c$$$      end do
c$$$      close(45)
c$$$      open(unit=54, file='chempot_temp.txt')
c$$$      do temp = 0.05d0, 5.0d0, 0.05d0
c$$$         ef_root = secant(chem_pot, ef0)
c$$$         write(54, *) temp, ef_root
c$$$      end do
c$$$      close(54)
      end

      real*8 function fermi_dist(e, ef, T)
        implicit real*8(a-h,o-z)
        exp_arg = (e-ef)/T

        if (exp_arg .gt. 700.0d0) then
          fermi_dist =0.0d0
        else
          fermi_dist = 1.0d0/(1.0d0+dexp(exp_arg))
        end if
      end function

      real*8 function chem_pot(ef)
        implicit real*8(a-h,o-z)
        common /param/ ei(364), temp, np
        sumni=0.0d0
        do i=1,364
           if (ei(i) .lt. 1e-5) exit
           fni = fermi_dist(ei(i), ef, temp)
           sumni = sumni+fni
           if(fni .lt. 1e-8) exit
        end do
        chem_pot = dfloat(np) - 2.0d0*sumni
c$$$        write(*, *) temp, np, ef, 2.0d0*sumni
      end function

      real*8 function fes_calc(inp, inn, spec, t)
        implicit real*8(a-h,o-z)
        common /param/ ei(364), temp, np
        dimension spec(1000)
        external sumene, chem_pot, fermi_dist, secant
        sum = 0.0d0
        temp=t
        if (temp .eq. 0.0d0) then
          sum=sumene(inp,spec,1000)
          sum=sum+sumene(inn,spec,1000)
        else
c$$$           protons
           do i=1,364
              ei(i)=spec(i)
           end do
           np = inp
           ef0 = ei(np/2)
           efp = secant(chem_pot, ef0)
c$$$           Neutrons
           np = inn
           efn = secant(chem_pot, ef0)
           write(21, *) efp, efn
           do i=1,364
              fp = fermi_dist(ei(i), efp, temp)
              fn = fermi_dist(ei(i), efn, temp)
              if( (fp .lt. 1e-8) .or. (fn .lt. 1e-8)) exit
              if ((fp .ge. 1d-45) .and. (fp .lt. 1.0d0))
     &             sp = -(fp*dlog(fp) + (1-fp)*dlog(1-fp))
              if ((fn .ge. 1d-45) .and. (fn .lt. 1.0d0))
     &             sn = -(fn*dlog(fn) + (1-fn)*dlog(1-fn))
c$$$              write(75, '(i4,5f12.6)') i,ei(i),fp,fn,sp,sn
              sum=sum+ (fp+fn)*ei(i)-temp*(sp+sn)
           end do
        end if
        fes_calc = sum
        return
      end function

      function sumene(ip,spe,ndim)
      implicit real*8(a-h,o-z)
      dimension spe(ndim)
      sumene=0.d0
      do i=1,ip/2
        sumene=sumene+2*spe(i)
      end do
      if(mod(ip,2).ne.0)sumene=sumene+spe(ip/2+1)
      end

ext: f


text:
      program testcalc
      implicit none
      real*8 eepdf_tilde_calc,eepdf_tilde
      real*8 x,Q2,jkb
      real*8 res_y,res_x,res_grid
      real*8 factor
      integer id
      x=1.00000000d0
      Q2=10d0*10d0
      read *,id
      do while (x .gt. 1.97d-4)
      res_x = eepdf_tilde_calc(x,Q2,id,11,11)
      res_grid = eepdf_tilde(x,Q2,id,11,11)
      print *,x,res_grid,res_x
      x=x*0.99999d0
      if (x.lt.0.99d0)then
      x=x*0.999d0
      endif
      enddo
      end program


text:
C
C      $Id: nm04f.f,v 1.8 2010-03-15 22:49:24 haley Exp $
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                       C
C                Copyright (C)  1997                                    C
C        University Corporation for Atmospheric Research                C
C                All Rights Reserved                                    C
C                                                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  File:       nm04f.f
C
C  Author:     Mary Haley (taken from one of Fred Clare's examples)
C          National Center for Atmospheric Research
C          PO 3000, Boulder, Colorado
C
C  Date:       Mon Dec 22 13:53:37 MST 1997
C
C  Description: Simple 3D interpolation.
C
      external NhlFAppClass
      external NhlFNcgmWorkstationClass
      external NhlFPSWorkstationClass
      external NhlFPDFWorkstationClass
      external NhlFCairoPSPDFWorkstationClass
      external NhlFCairoImageWorkstationClass
      external NhlFCairoWindowWorkstationClass

      parameter(NUM=1000,NX=21,NY=21,NZ=21)

      real xi(NUM), yi(NUM), zi(NUM), u(NUM)

      real xo(NX), yo(NY), zo(NZ), output(NX,NY,NZ)
      real xmin,ymin,zmin,xmax,ymax,zmax
      data xmin,ymin,zmin,xmax,ymax,zmax / -2., -2., -2., 2., 2., 2./
      integer ier
      integer appid,wid,gkswid
      integer srlist, grlist
      integer i
      character*7  wks_type
C
C Default is to display output to an NCGM workstation.
C
      wks_type = "ncgm"
C
C Initialize the high level utility library
C
      call NhlFInitialize
C
C Create an application context. Set the app dir to the current directory
C so the application looks for a resource file in the working directory.
C In this example the resource file supplies the plot title only.
C
      call NhlFRLCreate(srlist,'setrl')
      call NhlFRLCreate(grlist,'getrl')

      call NhlFRLClear(srlist)
      call NhlFRLSetString(srlist,'appDefaultParent','True',ierr)
      call NhlFRLSetString(srlist,'appUsrDir','./',ierr)
      call NhlFCreate(appid,'nm04',NhlFAppClass,0,srlist,ierr)

      if (wks_type.eq."ncgm".or.wks_type.eq."NCGM") then
C
C Create an NCGM workstation.
C
         call NhlFRLClear(srlist)
         call NhlFRLSetString(srlist,'wkMetaName','./nm04f.ncgm',ierr)
         call NhlFCreate(wid,'nm04Work',
     +        NhlFNcgmWorkstationClass,0,srlist,ierr)
      else if (wks_type.eq."x11".or.wks_type.eq."X11") then
C
C Create an X11 workstation.
C
         call NhlFRLClear(srlist)
         call NhlFRLSetString(srlist,'wkPause','True',ierr)
         call NhlFCreate(wid,'nm04Work',
     +        NhlFCairoWindowWorkstationClass,
     +        0,srlist,ierr)
      else if (wks_type.eq."oldps".or.wks_type.eq."OLDPS") then
C
C Create an older-style PostScript workstation.
C
         call NhlFRLClear(srlist)
         call NhlFRLSetString(srlist,'wkPSFileName','./nm04f.ps',ierr)
         call NhlFCreate(wid,'nm04Work',
     +        NhlFPSWorkstationClass,0,srlist,ierr)
      else if (wks_type.eq."oldpdf".or.wks_type.eq."OLDPDF") then
C
C Create an older-style PDF workstation.
C
         call NhlFRLClear(srlist)
         call NhlFRLSetString(srlist,'wkPDFFileName','./nm04f.pdf',ierr)
         call NhlFCreate(wid,'nm04Work',
     +        NhlFPDFWorkstationClass,0,srlist,ierr)
      else if (wks_type.eq."pdf".or.wks_type.eq."PDF".or.
     +         wks_type.eq."ps".or.wks_type.eq."PS") then
C
C Create a cairo PS/PDF workstation.
C
         call NhlFRLClear(srlist)
         call NhlFRLSetString(srlist,'wkFileName','./nm04f',ierr)
         call NhlFRLSetString(srlist,'wkFormat',wks_type,ierr)
         call NhlFCreate(wid,'nm04Work',
     +        NhlFCairoPSPDFWorkstationClass,0,srlist,ierr)
      else if (wks_type.eq."png".or.wks_type.eq."PNG") then
C
C Create a cairo PNG workstation.
C
         call NhlFRLClear(srlist)
         call NhlFRLSetString(srlist,'wkFileName','./nm04f',ierr)
         call NhlFRLSetString(srlist,'wkFormat',wks_type,ierr)
         call NhlFCreate(wid,'nm04Work',
     +        NhlFCairoImageWorkstationClass,0,srlist,ierr)
      endif

      do 10 i=1,NUM
         xi(i) = xmin+(xmax-xmin)*drnm1()
         yi(i) = ymin+(ymax-ymin)*drnm1()
         zi(i) = zmin+(zmax-zmin)*drnm1()
         u(i) = xi(i)**2 + yi(i)**2 + zi(i)**2
 10   continue
C
C  Create the output grid.
C
      do 102 i=1,NX
         xo(i) = xmin+(real(i-1)/real(NX-1))*(xmax-xmin)
 102  continue

      do 103 j =1,NY
         yo(j)= ymin+(real(j-1)/real(NY-1))*(ymax-ymin)
 103  continue

      do 104 k=1,NZ
        zo(k) = zmin+(real(k-1)/real(NZ-1))*(zmax-zmin)
  104 continue
C
C  Interpolate.
C
      call dsgrid3s(num,xi,yi,zi,u,NX,NY,NZ,xo,yo,zo,output,ier)
C
C Get Workstation ID.
C
      call NhlFRLClear(grlist)
      call NhlFRLGetInteger(grlist,'wkGksWorkId',gkswid,ierr)
      call NhlFGetValues(wid,grlist,ierr)
C
C There's no HLU object for surface plots yet, so we need to call the
C LLUs to get a surface plot.
C
      call gacwk (gkswid)
      call tdez3d(NX,NY,NZ,xo,yo,zo,output,3.0,2.,-35.,65.,6)
      call gdawk (gkswid)
      call NhlFFrame(wid,ierr)
C
C NhlDestroy destroys the given id and all of its children.
C
      call NhlFDestroy(wid,ierr)
C
C Restores state.
C
      call NhlFClose
      stop
      end
      REAL FUNCTION DRNM1()
C
C  Portable random number generator.
C
      PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773,
     +           MOMDMP=2836)
C
      INTEGER HVLUE, LVLUE, TESTV, NEXTN
      SAVE    NEXTN
      DATA JSEED,IFRST/123456789,0/
C
      IF (IFRST .EQ. 0) THEN
        NEXTN = JSEED
        IFRST = 1
      ENDIF
C
      HVLUE = NEXTN / MOBYMP
      LVLUE = MOD(NEXTN, MOBYMP)
      TESTV = MPLIER*LVLUE - MOMDMP*HVLUE
      IF (TESTV .GT. 0) THEN
        NEXTN = TESTV
      ELSE
        NEXTN = TESTV + MODLUS
      ENDIF
      DRNM1 = REAL(NEXTN)/REAL(MODLUS)
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


module class_Graph
  use class_Stack
  use utility

  integer :: WHITE = 0, GRAY = 1, BLACK = 2

  type Vertex
    integer :: value, predecessor, color
    real :: distance
    logical :: is_leaf
  end type Vertex

  type, public :: Graph
    type(Vertex), private, allocatable, dimension(:) :: vertices
    class(Stack), private, allocatable, dimension(:) :: adj_list

    contains

      procedure :: init => graph_init
      procedure :: color_vertex => graph_color_vertex
      procedure :: get_color => graph_get_color
      procedure :: get_adj_list => graph_get_adj_list
      procedure :: set_predecessor => graph_set_predecessor
      procedure :: get_predecessor => graph_get_predecessor
  end type

contains

  subroutine graph_init(self, num_vertices, edges, num_edges)
    use, intrinsic :: ieee_arithmetic

    ! define the parameters
    class(Graph), intent(inout) :: self
    integer, intent(in) :: num_edges
    integer, intent(in) :: num_vertices
    integer, intent(in) :: edges(:, :) ! a 2xV array

    ! define the data variables
    integer :: i
    real :: inf

    ! define infinity
    inf = ieee_value(inf, ieee_positive_inf)

    ! allocate the vertices and the adjacency list
    allocate(self%vertices(num_vertices))
    allocate(self%adj_list(num_vertices))

    ! initialized the vertices
    do i = 1, num_vertices
      self%vertices(i)%value = i
      self%vertices(i)%distance = inf
      self%vertices(i)%predecessor = 0
      self%vertices(i)%color = WHITE
      self%vertices(i)%is_leaf = .false.
    enddo

    ! initialize the adjacency lists
    do i = 1, num_vertices
      call self%adj_list(i)%init()
    enddo

    ! create the adjacency list
    do i = 1, num_edges
      call self%adj_list(edges(i, 1))%push(edges(i, 2))
    enddo

  end subroutine graph_init

  subroutine graph_color_vertex(self, vertex_index, color)
    Class(Graph), intent(inout) :: self
    integer, intent(in) :: vertex_index, color

    self%vertices(vertex_index)%color = color

  end subroutine graph_color_vertex

  function graph_get_color(self, vertex_index) result(out_color)
    Class(Graph), intent(inout) :: self
    integer, intent(in) :: vertex_index
    integer :: out_color

    out_color = self%vertices(vertex_index)%color

  end function

  function graph_get_adj_list(self, vertex_index) result(out_stackiterator)
    Class(Graph), intent(inout) :: self
    integer, intent(in) :: vertex_index
    type(StackIterator) :: out_stackiterator

    call out_stackiterator%create(self%adj_list(vertex_index))

  end function graph_get_adj_list

  subroutine graph_set_predecessor(self, vertex_index, vertex_parent_index)
    Class(Graph), intent(inout) :: self
    integer, intent(in) :: vertex_index, vertex_parent_index

    self%vertices(vertex_index)%predecessor = vertex_parent_index

  end subroutine graph_set_predecessor

  function graph_get_predecessor(self, vertex_index) result(out_vertex)
    Class(Graph), intent(inout) :: self
    integer, intent(in) :: vertex_index
    integer :: out_vertex

    out_vertex = self%vertices(vertex_index)%predecessor

  end function graph_get_predecessor

end module class_Graph


text:
Lives at 2500 Loyola

i have a fuckin tight cat that just sleeps all day and beats up dogs.

Photographer extrodinare, you can check out his stuff at http://www.logicalexposure.com

if you want to model for me my aim is : JeffLongshadow

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!-------------------------------------------------------------------------------------------------------------------------------
! This module contains all the filters
MODULE Filters
!...............................................................................................................................

    IMPLICIT NONE

CONTAINS

    !-------------------------------------------------------------------------------------------------------------------------------
    ! Discrete time Low-Pass Filter
    REAL FUNCTION LPFilter(InputSignal, DT, CornerFreq, iStatus, reset, inst)
    !...............................................................................................................................

		IMPLICIT NONE

			! Inputs

		REAL(4), INTENT(IN)			:: InputSignal
		REAL(4), INTENT(IN)			:: DT						! time step [s]
		REAL(4), INTENT(IN)			:: CornerFreq				! corner frequency [rad/s]
		INTEGER(4), INTENT(IN)		:: iStatus					! A status flag set by the simulation as follows: 0 if this is the first call, 1 for all subsequent time steps, -1 if this is the final call at the end of the simulation.
		INTEGER(4), INTENT(INOUT)	:: inst						! Instance number. Every instance of this function needs to have an unique instance number to ensure instances don't influence each other.
		LOGICAL(4), INTENT(IN)		:: reset					! Reset the filter to the input signal

			! Local

        REAL(4), DIMENSION(99), SAVE	:: InputSignalLast		! Input signal the last time this filter was called. Supports 99 separate instances.
        REAL(4), DIMENSION(99), SAVE	:: OutputSignalLast	! Output signal the last time this filter was called. Supports 99 separate instances.

			! Initialization

		IF ((iStatus == 0) .OR. reset) THEN
			OutputSignalLast(inst) = InputSignal
			InputSignalLast(inst) = InputSignal
		ENDIF

			! Body

		LPFilter = (DT*CornerFreq*InputSignal + DT*CornerFreq*InputSignalLast(inst) - (DT*CornerFreq-2.0)*OutputSignalLast(inst))/(DT*CornerFreq+2.0)

			! Save signals for next time step

		InputSignalLast(inst)  = InputSignal
		OutputSignalLast(inst) = LPFilter
		inst = inst + 1

    END FUNCTION LPFilter
    !-------------------------------------------------------------------------------------------------------------------------------
    ! Discrete time second order Low-Pass Filter
    REAL FUNCTION SecLPFilter(InputSignal, DT, CornerFreq, Damp, iStatus, reset, inst)
    !...............................................................................................................................

        IMPLICIT NONE

			! Inputs

        REAL(4), INTENT(IN)   		:: InputSignal
		REAL(4), INTENT(IN) 		:: DT						! time step [s]
		REAL(4), INTENT(IN)    		:: CornerFreq				! corner frequency [rad/s]
		REAL(4), INTENT(IN)    		:: Damp						! Dampening constant
        INTEGER(4), INTENT(IN)		:: iStatus					! A status flag set by the simulation as follows: 0 if this is the first call, 1 for all subsequent time steps, -1 if this is the final call at the end of the simulation.
        INTEGER(4), INTENT(INOUT)	:: inst						! Instance number. Every instance of this function needs to have an unique instance number to ensure instances don't influence each other.
		LOGICAL(4), INTENT(IN)		:: reset					! Reset the filter to the input signal

			! Local

        REAL(4), DIMENSION(99), SAVE :: InputSignalLast1	! Input signal the last time this filter was called. Supports 99 separate instances.
		REAL(4), DIMENSION(99), SAVE :: InputSignalLast2	! Input signal the next to last time this filter was called. Supports 99 separate instances.
        REAL(4), DIMENSION(99), SAVE :: OutputSignalLast1	! Output signal the last time this filter was called. Supports 99 separate instances.
        REAL(4), DIMENSION(99), SAVE :: OutputSignalLast2	! Output signal the next to last time this filter was called. Supports 99 separate instances.

			! Initialization

        IF ((iStatus == 0) .OR. reset )  THEN
            OutputSignalLast1(inst)  = InputSignal
            OutputSignalLast2(inst)  = InputSignal
            InputSignalLast1(inst)   = InputSignal
            InputSignalLast2(inst)   = InputSignal
        ENDIF

			! Body

        SecLPFilter = 1/(4+4*DT*Damp*CornerFreq+DT**2*CornerFreq**2) * ( (8-2*DT**2*CornerFreq**2)*OutputSignalLast1(inst) &
						+ (-4+4*DT*Damp*CornerFreq-DT**2*CornerFreq**2)*OutputSignalLast2(inst) + (DT**2*CornerFreq**2)*InputSignal &
							+ (2*DT**2*CornerFreq**2)*InputSignalLast1(inst) + (DT**2*CornerFreq**2)*InputSignalLast2(inst) )

			! Save signals for next time step

        InputSignalLast2(inst)   = InputSignalLast1 (inst)
        InputSignalLast1(inst)   = InputSignal
        OutputSignalLast2(inst)  = OutputSignalLast1 (inst)
        OutputSignalLast1(inst)  = SecLPFilter
		inst = inst + 1

    END FUNCTION SecLPFilter
    !-------------------------------------------------------------------------------------------------------------------------------
    ! Discrete time High-Pass Filter
    REAL FUNCTION HPFilter( InputSignal, DT, CornerFreq, iStatus, reset, inst)
    !...............................................................................................................................

        IMPLICIT NONE

			! Inputs

        REAL(4), INTENT(IN)		:: InputSignal
		REAL(4), INTENT(IN)		:: DT						! time step [s]
		REAL(4), INTENT(IN)		:: CornerFreq				! corner frequency [rad/s]
        INTEGER, INTENT(IN)		:: iStatus					! A status flag set by the simulation as follows: 0 if this is the first call, 1 for all subsequent time steps, -1 if this is the final call at the end of the simulation.
        INTEGER, INTENT(INOUT)	:: inst						! Instance number. Every instance of this function needs to have an unique instance number to ensure instances don't influence each other.
		LOGICAL(4), INTENT(IN)	:: reset					! Reset the filter to the input signal

			! Local

		REAL(4)							:: K						! Constant gain
        REAL(4), DIMENSION(99), SAVE	:: InputSignalLast		! Input signal the last time this filter was called. Supports 99 separate instances.
        REAL(4), DIMENSION(99), SAVE	:: OutputSignalLast	! Output signal the last time this filter was called. Supports 99 separate instances.

			! Initialization

        IF ((iStatus == 0) .OR. reset)  THEN
            OutputSignalLast(inst) = InputSignal
            InputSignalLast(inst) = InputSignal
        ENDIF

        K = 2.0 / DT

			! Body

        HPFilter = K/(CornerFreq + K)*InputSignal - K/(CornerFreq + K)*InputSignalLast(inst) - (CornerFreq - K)/(CornerFreq + K)*OutputSignalLast(inst)

			! Save signals for next time step

        InputSignalLast(inst)   = InputSignal
        OutputSignalLast(inst)  = HPFilter
		inst = inst + 1

    END FUNCTION HPFilter
    !-------------------------------------------------------------------------------------------------------------------------------
    ! Discrete time inverted Notch Filter with descending slopes, G = CornerFreq*s/(Damp*s^2+CornerFreq*s+Damp*CornerFreq^2)
    REAL FUNCTION NotchFilterSlopes(InputSignal, DT, CornerFreq, Damp, iStatus, reset, inst)
    !...............................................................................................................................

        IMPLICIT NONE

			! Inputs

		REAL(4), INTENT(IN)     :: InputSignal
		REAL(4), INTENT(IN)     :: DT						! time step [s]
		REAL(4), INTENT(IN)     :: CornerFreq				! corner frequency [rad/s]
		REAL(4), INTENT(IN)     :: Damp						! Dampening constant
		INTEGER, INTENT(IN)		:: iStatus					! A status flag set by the simulation as follows: 0 if this is the first call, 1 for all subsequent time steps, -1 if this is the final call at the end of the simulation.
		INTEGER, INTENT(INOUT)	:: inst						! Instance number. Every instance of this function needs to have an unique instance number to ensure instances don't influence each other.
		LOGICAL(4), INTENT(IN)	:: reset					! Reset the filter to the input signal

			! Local

		REAL(4), DIMENSION(99), SAVE :: InputSignalLast1	! Input signal the last time this filter was called. Supports 99 separate instances.
		REAL(4), DIMENSION(99), SAVE :: InputSignalLast2	! Input signal the next to last time this filter was called. Supports 99 separate instances.
		REAL(4), DIMENSION(99), SAVE :: OutputSignalLast1	! Output signal the last time this filter was called. Supports 99 separate instances.
		REAL(4), DIMENSION(99), SAVE :: OutputSignalLast2	! Output signal the next to last time this filter was called. Supports 99 separate instances.

			! Initialization

		IF ((iStatus == 0) .OR. reset) THEN
			OutputSignalLast1(inst)  = InputSignal
			OutputSignalLast2(inst)  = InputSignal
			InputSignalLast1(inst)   = InputSignal
			InputSignalLast2(inst)   = InputSignal
		ENDIF

			! Body

		NotchFilterSlopes = 1.0/(4.0+2.0*DT*Damp*CornerFreq+DT**2.0*CornerFreq**2.0) * ( (8.0-2.0*DT**2.0*CornerFreq**2.0)*OutputSignalLast1(inst) &
						+ (-4.0+2.0*DT*Damp*CornerFreq-DT**2.0*CornerFreq**2.0)*OutputSignalLast2(inst) + &
							(2.0*DT*Damp*CornerFreq)*InputSignal + (-2.0*DT*Damp*CornerFreq)*InputSignalLast2(inst) )

			! Save signals for next time step

		InputSignalLast2(inst)   = InputSignalLast1(inst)
		InputSignalLast1(inst)   = InputSignal			!Save input signal for next time step
		OutputSignalLast2(inst)  = OutputSignalLast1(inst)		!Save input signal for next time step
		OutputSignalLast1(inst)  = NotchFilterSlopes
		inst = inst + 1

	END FUNCTION NotchFilterSlopes
	!-------------------------------------------------------------------------------------------------------------------------------
	! Discrete time Notch Filter, G = (s^2 + 2*omega*betaNum*s + omega^2)/(s^2 + 2*omega*betaDen*s + omega^2)
	REAL FUNCTION NotchFilter(InputSignal, DT, omega, betaNum, betaDen, iStatus, reset, inst)
	!...............................................................................................................................

		IMPLICIT NONE

			! Inputs

		REAL(4), INTENT(IN)     :: InputSignal
		REAL(4), INTENT(IN)     :: DT						! time step [s]
		REAL(4), INTENT(IN)     :: omega					! corner frequency [rad/s]
		REAL(4), INTENT(IN)     :: betaNum					! Dampening constant in numerator of filter transfer function
		REAL(4), INTENT(IN)     :: betaDen					! Dampening constant in denominator of filter transfer function
		INTEGER, INTENT(IN)		:: iStatus					! A status flag set by the simulation as follows: 0 if this is the first call, 1 for all subsequent time steps, -1 if this is the final call at the end of the simulation.
		INTEGER, INTENT(INOUT)	:: inst						! Instance number. Every instance of this function needs to have an unique instance number to ensure instances don't influence each other.
		LOGICAL(4), INTENT(IN)	:: reset					! Reset the filter to the input signal

			! Local
		REAL(4)							:: K, P1, P2, P3, P4, P5	! Constant gain
		REAL(4), DIMENSION(99), SAVE	:: InputSignalLast1			! Input signal the last time this filter was called. Supports 99 separate instances.
		REAL(4), DIMENSION(99), SAVE	:: InputSignalLast2			! Input signal the next to last time this filter was called. Supports 99 separate instances.
		REAL(4), DIMENSION(99), SAVE	:: OutputSignalLast1		! Output signal the last time this filter was called. Supports 99 separate instances.
		REAL(4), DIMENSION(99), SAVE	:: OutputSignalLast2		! Output signal the next to last time this filter was called. Supports 99 separate instances.

			! Initialization

		IF ((iStatus == 0) .OR. reset) THEN
			OutputSignalLast1(inst)  = InputSignal
			OutputSignalLast2(inst)  = InputSignal
			InputSignalLast1(inst)   = InputSignal
			InputSignalLast2(inst)   = InputSignal
		ENDIF

		K = 2/DT
		P1 = (K**2 + 2*omega*BetaNum*K + omega**2)/(K**2 + 2*omega*BetaDen*K + omega**2)
		P2 = (2*omega**2 - 2*K**2)  / (K**2 + 2*omega*BetaDen*K + omega**2);
		P3 = (K**2 - 2*omega*BetaNum*K + omega**2) / (K**2 + 2*omega*BetaDen*K + omega**2)
		P4 = (2*omega**2 - 2*K**2)  / (K**2 + 2*omega*BetaDen*K + omega**2)
		P5 = (K**2 - 2*omega*BetaDen*K + omega**2)/ (K**2 + 2*omega*BetaDen*K + omega**2)

			! Body

		NotchFilter = P1*InputSignal + P2*InputSignalLast1(inst) + P3*InputSignalLast2(inst) - P4*OutputSignalLast1(inst) - P5*OutputSignalLast2(inst)

			! Save signals for next time step

		InputSignalLast2(inst)   = InputSignalLast1(inst)
		InputSignalLast1(inst)   = InputSignal			!Save input signal for next time step
		OutputSignalLast2(inst)  = OutputSignalLast1(inst)		!Save input signal for next time step
		OutputSignalLast1(inst)  = NotchFilter
		inst = inst + 1

	END FUNCTION NotchFilter
	!-------------------------------------------------------------------------------------------------------------------------------
	END MODULE Filters


text:
C++**********************************************************************
C Set of subroutines
C to read Midas tables and spectra (1-D image)
C
C Contains:
C JLP_RDTABLE, JLP_RDSPEC, JLP_TRANMI
C
C JLP
C Version 15-07-94
C--**********************************************************************
C++**********************************************************************
C Subroutine JLP_RDTABLE
C To read MIDAS tables (the first two columns)
C Use MIDAS routines
C--**********************************************************************
	SUBROUTINE JLP_RDTABLE(XX,YY,NROW)
	REAL XX(*),YY(*)
	INTEGER ISTAT
	INTEGER NCOL,NROW,NSC1,JLP_INIT_MIDAS
	INTEGER MADRID(1),ITID,ACOL,AROW
	LOGICAL*4 NULL
	CHARACTER TABLE*40
	COMMON /JLP_MIDAS/JLP_INIT_MIDAS
	COMMON /VMR/ MADRID
	INCLUDE 'jlpsub:st_def.inc'
	INCLUDE 'jlpsub:st_dat.inc'

C Midas environment:
C Initializing MIDAS parameters:
	IF(JLP_INIT_MIDAS.NE.1)THEN
	  CALL JLP_SCSPRO('-1')
	  JLP_INIT_MIDAS = 1
	ENDIF

C Get the table name
	WRITE(6,*) ' Input table :==  '
	READ(5,10) TABLE
10	FORMAT(A)

C Read the table:
	CALL TBTOPN(TABLE,F_I_MODE,ITID,ISTAT)
	IF(ISTAT.NE.0) GO TO 1000
	CALL TBIGET(ITID,NCOL,NROW,NSC1,ACOL,AROW,ISTAT)
	IF(ISTAT.NE.0) GO TO 1000

C Displays the table description:
	WRITE(6,77) NCOL,NROW,NSC1
77	FORMAT(' NCOL=',I5,' NROW=',I5,' NSC1=',I5,
     1	/,' ACOL=',I5,' AROW=',I5)

C Read the first two columns:
	DO J=1,NROW
	  CALL TBERDR(ITID,J,1,XX(J),NULL,ISTAT)
	  IF(ISTAT.NE.0) GO TO 1000
	  CALL TBERDR(ITID,J,2,YY(J),NULL,ISTAT)
	  IF(ISTAT.NE.0) GO TO 1000

C Displays the output:
C	  WRITE(6,78) J,XX(J),YY(J)
C78	  FORMAT(' J,X,Y',X,I5,2X,F12.6,2X,F12.6)

	END DO

	RETURN

C When error:
1000 	WRITE(6,*) ' Fatal error reading MIDAS table'
	CALL JLP_END
	RETURN
	END

C++**********************************************************************
C Subroutine JLP_RDSPEC
C To read MIDAS spectra (1-D images)
C Uses MIDAS routines
C
C--**********************************************************************
	SUBROUTINE JLP_RDSPEC(XX,YY,NX)
	REAL*4 XX(*),YY(*)
	REAL*8 START(2),STEP(2)
	INTEGER*4 NPIX(2),PNTR_IMAGE
	INTEGER MADRID(1),JLP_INIT_MIDAS
	CHARACTER NAME*80,COMMENTS*80
	CHARACTER CUNIT*80,IDENT*72
	INCLUDE 'jlpsub:st_def.inc'
	INCLUDE 'jlpsub:st_dat.inc'
	COMMON /VMR/MADRID
	COMMON /JLP_MIDAS/JLP_INIT_MIDAS

C Initializing MIDAS parameters:
	IF(JLP_INIT_MIDAS.NE.1)THEN
	  CALL JLP_SCSPRO('-1')
	  JLP_INIT_MIDAS = 1
	ENDIF

C Opening the MIDAS file
5	WRITE(6,*) ' Input 1-D image ?'
	READ(5,10) NAME
10	FORMAT(A)
	CALL STIGET(NAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     1	1,NAXIS,NPIX,START,STEP,IDENT,CUNIT,
     1	PNTR_IMAGE,INUMB,ISTAT)
	IF(ISTAT.NE.0)THEN
	  PRINT *,' ERROR JLP_RDMIDAS/STIGET:  ISTAT=',ISTAT
	  PRINT *,' ERROR WHEN OPENING ',NAME
	  GOTO 5
	ENDIF

	NX=NPIX(1)
	COMMENTS=' '
	COMMENTS(1:72)=IDENT(1:72)
c
c Transfering the data:
c
	CALL JLP_TRANMI(YY,MADRID(PNTR_IMAGE),NX)

C Generating the wavelengths:
	DO I=1,NX
	   XX(I)=START(1)+(I-1)*STEP(1)
	END DO

	RETURN
	END
C----------------------------------------------------------------------
C Transfering the input MIDAS BDF file into the output array
C Called by JLP_WRMIDAS
C----------------------------------------------------------------------
	SUBROUTINE JLP_TRANMI(OUTPUT,INPUT,NX)
	REAL*4 INPUT(*),OUTPUT(*)

	  DO IX=1,NX
	    OUTPUT(IX)=INPUT(IX)
	  END DO

	RETURN
	END

ext: for


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE CALC_ELEM_STRESSES ( SIZE_ALLOCATED, NUM1, NUM_FEMAP_ROWS, WRITE_OGEL, WRITE_FEMAP )

! Calls routines that process the STRESS array calculated in subr ELEM_STRE_STRN_ARRAYS to obtain element specific stress values
! that will be written to the F06 file

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG
      USE IOUNT1, ONLY                :  WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE MODEL_STUF, ONLY            :  TYPE
      USE SUBR_BEGEND_LEVELS, ONLY    :  CALC_ELEM_STRESSES_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO

      USE CALC_ELEM_STRESSES_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'CALC_ELEM_STRESSES'
      CHARACTER( 1*BYTE), INTENT(IN)  :: WRITE_OGEL        ! If 'Y' then write data to array OGEL
      CHARACTER( 1*BYTE), INTENT(IN)  :: WRITE_FEMAP       ! If 'Y' then write data to array FEMAP arrays

      INTEGER(LONG), INTENT(IN)       :: SIZE_ALLOCATED    ! No. of rows allocated to array that will be written to in a subr
!                                                            called here so we can check that we don't try to write more rows
      INTEGER(LONG), INTENT(IN)       :: NUM_FEMAP_ROWS    ! Number of rows that will be written to FEMAP arrays
      INTEGER(LONG), INTENT(INOUT)    :: NUM1              ! Cum rows written to OGEL prior to running this subr
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = CALC_ELEM_STRESSES_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! Calculate STRESS(4-6) for elements that have nonzero STRESS(4-6)

      IF      ((TYPE(1:3) == 'BAR'  ) .OR. (TYPE(1:4) == 'BUSH') .OR. (TYPE(1:4) == 'ELAS') .OR. (TYPE(1:3) == 'ROD'  )) THEN
         CALL ONE_D_STRESS_OUTPUTS ( SIZE_ALLOCATED, NUM1, NUM_FEMAP_ROWS, WRITE_OGEL, WRITE_FEMAP )

      ELSE IF ((TYPE(1:5) == 'TRIA3') .OR. (TYPE(1:5) == 'QUAD4') .OR. (TYPE(1:5) == 'SHEAR') .OR. (TYPE(1:6) == 'USERIN')) THEN
         CALL SHELL_STRESS_OUTPUTS ( SIZE_ALLOCATED, NUM1, NUM_FEMAP_ROWS, WRITE_OGEL, WRITE_FEMAP )

      ELSE IF ((TYPE(1:4) == 'HEXA' ) .OR. (TYPE(1:5) == 'PENTA') .OR. (TYPE(1:5) == 'TETRA')) THEN
         CALL SOLID_STRESS_OUTPUTS ( SIZE_ALLOCATED, NUM1, NUM_FEMAP_ROWS, WRITE_OGEL, WRITE_FEMAP )

      ELSE

         FATAL_ERR = FATAL_ERR + 1
         WRITE(ERR,9203) SUBR_NAME, TYPE
         WRITE(F06,9203) SUBR_NAME, TYPE
         CALL OUTA_HERE ( 'Y' )

      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
! not used
 9200 FORMAT(' *ERROR  9200: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' ARRAY OGEL WAS ALLOCATED TO HAVE ',I12,' ROWS. ATTEMPT TO WRITE TO OGEL BEYOND THIS')

 9203 FORMAT(' *ERROR  9203: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' INCORRECT ELEMENT TYPE = "',A,'"')

! **********************************************************************************************************************************

      END SUBROUTINE CALC_ELEM_STRESSES


text:
      real function smtmax(x,ldx,nrow,ncol)
c--------------------------------------------------------------------
c        scan two dimensional floating point array x and returns
c  largest or smallest value.
c    entry points-
c      smtmax - returns maximum element of array x.
c      smtmin - returns minimum element of array x.
c
c  arguments-
c      x   - real nrow by ncol array to be scanned
c      ldx - first dimension of x in calling program.
c      nrow- number of rows in x
c      ncol- number of columns in x.
c
c  language - 1977 ansi standard fortran
c  written  - july 1981
c  author   - gary l. pavlis
c------------------------------------------------------------------
      integer nrow,ncol,ldx
      real x(ldx,ncol)
      smtmax = x(1,1)
      do 150 j=1,ncol
            do 100 i=1,nrow
                  smtmax = amax1(x(i,j),smtmax)
  100       continue
  150 continue
      return
      entry smtmin(x,ldx,nrow,ncol)
      smtmin = x(1,1)
      do 250 j=1,ncol
            do 200 i=1,nrow
                  smtmin = amin1(x(i,j),smtmin)
  200       continue
  250 continue
      return
      end

c $Id$


text:
C    Copyright(C) 2008 Sandia Corporation.  Under the terms of Contract
C    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
C    certain rights in this software
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C    * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C    * Redistributions in binary form must reproduce the above
C      copyright notice, this list of conditions and the following
C      disclaimer in the documentation and/or other materials provided
C      with the distribution.
C
C    * Neither the name of Sandia Corporation nor the names of its
C      contributors may be used to endorse or promote products derived
C      from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C
      SUBROUTINE SRCHC(IVEC,ILBIN,IUBIN,NUM,ICOD,LOC)
C
C     THIS SUBROUTINE SEARCHES ROW ILB THOURGH ROW IUB OF A
C     NUMERICALLY ORDERED CHARACTER COLUMN VECTOR FOR THE OCCURANCE
C     THE VALUE NUM.
C     IF NUM IS FOUND ICOD IS SET TO UNITY AND LOC IS THE ROW NUMBER
C     WHERE NUM RESIDES. IF NUM IS NOT FOUND ICOD IS ZERO AND LOC IS
C     THE ROW NUMBER WHERE NUM WOULD RESIDE IF IT WERE IN THE NUMER-
C     ICALLY ORDER LIST.   BOB LUST
C
C     THIS SUBROUTINE HAS BEEN CHANGED FROM BOB LUST'S VERSION
C     AND NOW ASSUMES THAT THERE IS NO MORE THAN ONE MATCH
C     IN THE ORDERED LIST 'IVEC'.  BILL MILLS-CURRAN  JAN. 1, 1983
C
C     IVEC    ORDERED CHARACTER LIST (SINGLE COLUMN)
C
C     ILBIN   LOW NUMBERED ROW OF SEARCH RANGE
C
C     IUBIN   HIGH NUMBERED ROW OF SEARCH RANGE
C
C     NUM     VALUE TO BE LOCATED IN IVEC
C
C     ICOD    RETURN CODE  0 = NO MATCH   1 = MATCH
C
C     LOC     LOCATION IN IVEC FOR NUM
C
      CHARACTER*(*) IVEC(1),NUM
C
      ILB = ILBIN
      IUB = IUBIN
      ICOD=0
      IF (IUB .LT. ILB) THEN
         LOC = 1
         RETURN
      END IF
C
C     CHECK TO SEE IF NUM IS AT EITHER END OF LIST
C
      IF(IVEC(ILB).GT.NUM)THEN
         LOC=ILB
         RETURN
      ELSE IF(IVEC(IUB).LT.NUM) THEN
         LOC=IUB+1
         RETURN
      END IF
C
C     NUM IS INTERNAL TO IVEC
C
  100 MID=(ILB+IUB)/2
      IF(MID.LE.ILB)GO TO 110
C
C     SEARCH RANGE IS MORE THAN 2
C
      IF(IVEC(MID).LT.NUM) THEN
C
C        UPPER PART OF LIST
C
         ILB=MID
         GO TO 100
      ELSE IF(IVEC(MID).GT.NUM) THEN
C
C        LOWER PART OF LIST
C
         IUB=MID
         GO TO 100
      ELSE
C
C        MATCH HAS OCCURED AT "MID"
C
         ICOD=1
         LOC=MID
         RETURN
      END IF
  110 CONTINUE
C
C     SEARCH RANGE IS 2 OR LESS
C
      IF(NUM.EQ.IVEC(ILB)) THEN
C
C        MATCH AT "ILB"
C
         ICOD=1
         LOC=ILB
         RETURN
      ELSE IF(NUM.EQ.IVEC(IUB)) THEN
C
C        MATCH AT "IUB"
C
         ICOD=1
         LOC=IUB
         RETURN
      ELSE
C
C        NO MATCH IN LIST.
C        LOCATION FOR NEW ENTRY IS "IUB".
C
         LOC=IUB
         RETURN
      END IF
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
!
! LBLRTM_Fhdr_netCDF_IO
!
! Module containing routine to read and write LBLRTM Fhdr objects as
! groups to a netCDF format file.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 19-Feb-2014
!                       paul.vandelst@noaa.gov
!

MODULE LBLRTM_Fhdr_netCDF_IO

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds        , ONLY: FP, IP, DP => Double, Long
  USE Message_Handler   , ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  USE String_Utility    , ONLY: StrClean
  USE LBLRTM_Parameters , ONLY: N_MOL => LBLRTM_MAX_N_MOLECULES
  USE LBLRTM_Fhdr_Define, ONLY: LBLRTM_Fhdr_type    , &
                                LBLRTM_Fhdr_IsValid , &
                                LBLRTM_Fhdr_SetValid, &
                                LBLRTM_Fhdr_Destroy , &
                                LBLRTM_Fhdr_Inspect
  USE netcdf
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  ! Procedures
  PUBLIC :: LBLRTM_Fhdr_netCDF_ReadGroup
  PUBLIC :: LBLRTM_Fhdr_netCDF_WriteGroup
  PUBLIC :: LBLRTM_Fhdr_netCDF_IOVersion

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  ! Default message string length
  INTEGER, PARAMETER :: ML = 1024
  ! Literal constants
  REAL(DP), PARAMETER :: ZERO = 0.0_DP
  REAL(DP), PARAMETER :: ONE  = 1.0_DP
  ! Extra parameters not in netCDF(?)
  INTEGER, PARAMETER :: MAX_N_GROUPS = 8096

  ! Global attribute names. Case sensitive
  CHARACTER(*), PARAMETER :: RELEASE_GATTNAME = 'Release'
  CHARACTER(*), PARAMETER :: VERSION_GATTNAME = 'Version'
  CHARACTER(*), PARAMETER :: TITLE_GATTNAME   = 'Title'
  CHARACTER(*), PARAMETER :: HISTORY_GATTNAME = 'History'
  CHARACTER(*), PARAMETER :: COMMENT_GATTNAME = 'Comment'

  ! Dimension names
  CHARACTER(*), PARAMETER :: MOLECULE_DIMNAME   = 'n_molecules'
  CHARACTER(*), PARAMETER :: ANCILLARY_DIMNAME  = 'n_ancillary'
  CHARACTER(*), PARAMETER :: UID_STRLEN_DIMNAME = 'uid_strlen'
  CHARACTER(*), PARAMETER :: SL_STRLEN_DIMNAME  = 'sl_strlen'

  ! Variable names
  CHARACTER(*), PARAMETER :: USER_ID_VARNAME          = 'User_ID'
  CHARACTER(*), PARAMETER :: COL_SCALE_FACTOR_VARNAME = 'Column_Scale_Factor'
  CHARACTER(*), PARAMETER :: AVG_LAYER_PRES_VARNAME   = 'Average_Layer_Pressure'
  CHARACTER(*), PARAMETER :: AVG_LAYER_TEMP_VARNAME   = 'Average_Layer_Temperature'
  CHARACTER(*), PARAMETER :: MOL_ID_VARNAME           = 'Molecule_Id'
  CHARACTER(*), PARAMETER :: MOL_COL_DENS_VARNAME     = 'Molecule_Column_Density'
  CHARACTER(*), PARAMETER :: BROAD_COL_DENS_VARNAME   = 'Broadening_Gas_Column_Density'
  CHARACTER(*), PARAMETER :: FREQ_INTERVAL_VARNAME    = 'Frequency_Interval'
  CHARACTER(*), PARAMETER :: BEGIN_FREQ_VARNAME       = 'Begin_Frequency'
  CHARACTER(*), PARAMETER :: END_FREQ_VARNAME         = 'End_Frequency'
  CHARACTER(*), PARAMETER :: BDRY_TEMP_VARNAME        = 'Boundary_Temperature'
  CHARACTER(*), PARAMETER :: BDRY_EMIS_VARNAME        = 'Boundary_Emissivity'
  CHARACTER(*), PARAMETER :: N_MOLECULES_VARNAME      = 'n_Molecules'
  CHARACTER(*), PARAMETER :: N_LAYER_VARNAME          = 'n_Layer'
  CHARACTER(*), PARAMETER :: OD_LAYER_FLAG_VARNAME    = 'OD_Layering_Control_Flag'
  CHARACTER(*), PARAMETER :: CALC_DATE_VARNAME        = 'Calculation_Date'
  CHARACTER(*), PARAMETER :: CALC_TIME_VARNAME        = 'Calculation_Time'
  CHARACTER(*), PARAMETER :: ANCILLARY_VARNAME        = 'ancillary'
  ! ...The run flags
  CHARACTER(*), PARAMETER :: HIRAC_VARNAME = 'hirac'
  CHARACTER(*), PARAMETER :: LBLF4_VARNAME = 'lblf4'
  CHARACTER(*), PARAMETER :: XSCNT_VARNAME = 'xscnt'
  CHARACTER(*), PARAMETER :: AERSL_VARNAME = 'aersl'
  CHARACTER(*), PARAMETER :: EMIT_VARNAME  = 'emit'
  CHARACTER(*), PARAMETER :: SCAN_VARNAME  = 'scan'
  CHARACTER(*), PARAMETER :: PLOT_VARNAME  = 'plot'
  CHARACTER(*), PARAMETER :: PATH_VARNAME  = 'path'
  CHARACTER(*), PARAMETER :: JRAD_VARNAME  = 'jrad'
  CHARACTER(*), PARAMETER :: TEST_VARNAME  = 'test'
  CHARACTER(*), PARAMETER :: MERGE_VARNAME = 'merge'
  CHARACTER(*), PARAMETER :: SCNID_VARNAME = 'scnid'
  CHARACTER(*), PARAMETER :: HWHM_VARNAME  = 'hwhm'
  CHARACTER(*), PARAMETER :: IDABS_VARNAME = 'idabs'
  CHARACTER(*), PARAMETER :: ATM_VARNAME   = 'atm'
  CHARACTER(*), PARAMETER :: LAYR1_VARNAME = 'layr1'
  CHARACTER(*), PARAMETER :: NLAYR_VARNAME = 'nlayr'

  ! Variable long name attribute
  CHARACTER(*), PARAMETER :: LONGNAME_ATTNAME = 'long_name'

  CHARACTER(*), PARAMETER :: USER_ID_LONGNAME          = 'User ID'
  CHARACTER(*), PARAMETER :: COL_SCALE_FACTOR_LONGNAME = 'Column Scale Factor'
  CHARACTER(*), PARAMETER :: AVG_LAYER_PRES_LONGNAME   = 'Average Layer Pressure'
  CHARACTER(*), PARAMETER :: AVG_LAYER_TEMP_LONGNAME   = 'Average Layer Temperature'
  CHARACTER(*), PARAMETER :: MOL_ID_LONGNAME           = 'Molecule Id'
  CHARACTER(*), PARAMETER :: MOL_COL_DENS_LONGNAME     = 'Molecule Column Density'
  CHARACTER(*), PARAMETER :: BROAD_COL_DENS_LONGNAME   = 'Broadening Gas Column Density'
  CHARACTER(*), PARAMETER :: FREQ_INTERVAL_LONGNAME    = 'Frequency Interval'
  CHARACTER(*), PARAMETER :: BEGIN_FREQ_LONGNAME       = 'Begin Frequency'
  CHARACTER(*), PARAMETER :: END_FREQ_LONGNAME         = 'End Frequency'
  CHARACTER(*), PARAMETER :: BDRY_TEMP_LONGNAME        = 'Boundary Temperature'
  CHARACTER(*), PARAMETER :: BDRY_EMIS_LONGNAME        = 'Boundary Emissivity'
  CHARACTER(*), PARAMETER :: N_MOLECULES_LONGNAME      = 'Number of molecules'
  CHARACTER(*), PARAMETER :: N_LAYER_LONGNAME          = 'Number of layer'
  CHARACTER(*), PARAMETER :: OD_LAYER_FLAG_LONGNAME    = 'OD Layering Control Flag'
  CHARACTER(*), PARAMETER :: CALC_DATE_LONGNAME        = 'Calculation Date'
  CHARACTER(*), PARAMETER :: CALC_TIME_LONGNAME        = 'Calculation Time'
  CHARACTER(*), PARAMETER :: ANCILLARY_LONGNAME        = 'ancillary'
  ! ...The run flags
  CHARACTER(*), PARAMETER :: HIRAC_LONGNAME = 'hirac run flag'
  CHARACTER(*), PARAMETER :: LBLF4_LONGNAME = 'lblf4 run flag'
  CHARACTER(*), PARAMETER :: XSCNT_LONGNAME = 'xscnt run flag'
  CHARACTER(*), PARAMETER :: AERSL_LONGNAME = 'aersl run flag'
  CHARACTER(*), PARAMETER :: EMIT_LONGNAME  = 'emit run flag'
  CHARACTER(*), PARAMETER :: SCAN_LONGNAME  = 'scan run flag'
  CHARACTER(*), PARAMETER :: PLOT_LONGNAME  = 'plot run flag'
  CHARACTER(*), PARAMETER :: PATH_LONGNAME  = 'path run flag'
  CHARACTER(*), PARAMETER :: JRAD_LONGNAME  = 'jrad run flag'
  CHARACTER(*), PARAMETER :: TEST_LONGNAME  = 'test run flag'
  CHARACTER(*), PARAMETER :: MERGE_LONGNAME = 'merge run flag'
  CHARACTER(*), PARAMETER :: SCNID_LONGNAME = 'scnid run flag'
  CHARACTER(*), PARAMETER :: HWHM_LONGNAME  = 'hwhm run flag'
  CHARACTER(*), PARAMETER :: IDABS_LONGNAME = 'idabs run flag'
  CHARACTER(*), PARAMETER :: ATM_LONGNAME   = 'atm run flag'
  CHARACTER(*), PARAMETER :: LAYR1_LONGNAME = 'layr1 run flag'
  CHARACTER(*), PARAMETER :: NLAYR_LONGNAME = 'nlayr run flag'

  ! Variable description attribute
  CHARACTER(*), PARAMETER :: DESCRIPTION_ATTNAME = 'description'

  CHARACTER(*), PARAMETER :: USER_ID_DESCRIPTION          = 'User Identification string'
  CHARACTER(*), PARAMETER :: COL_SCALE_FACTOR_DESCRIPTION = 'Column profile amount scaling factor'
  CHARACTER(*), PARAMETER :: AVG_LAYER_PRES_DESCRIPTION   = 'Average layer pressure'
  CHARACTER(*), PARAMETER :: AVG_LAYER_TEMP_DESCRIPTION   = 'Average layer temperature'
  CHARACTER(*), PARAMETER :: MOL_ID_DESCRIPTION           = 'Molecule identification string'
  CHARACTER(*), PARAMETER :: MOL_COL_DENS_DESCRIPTION     = 'Molecule column density'
  CHARACTER(*), PARAMETER :: BROAD_COL_DENS_DESCRIPTION   = 'Broadening gas column density'
  CHARACTER(*), PARAMETER :: FREQ_INTERVAL_DESCRIPTION    = 'Calculation frequency interval'
  CHARACTER(*), PARAMETER :: BEGIN_FREQ_DESCRIPTION       = 'Calculation begin frequency'
  CHARACTER(*), PARAMETER :: END_FREQ_DESCRIPTION         = 'Calculation end frequency'
  CHARACTER(*), PARAMETER :: BDRY_TEMP_DESCRIPTION        = 'Boundary temperature'
  CHARACTER(*), PARAMETER :: BDRY_EMIS_DESCRIPTION        = 'Boundary emissivity'
  CHARACTER(*), PARAMETER :: N_MOLECULES_DESCRIPTION      = 'Number of gaseous absorbers included used in calculation'
  CHARACTER(*), PARAMETER :: N_LAYER_DESCRIPTION          = 'Number of atmospheric layer'
  CHARACTER(*), PARAMETER :: OD_LAYER_FLAG_DESCRIPTION    = 'Optical depth layering control flag'
  CHARACTER(*), PARAMETER :: CALC_DATE_DESCRIPTION        = 'Calculation date'
  CHARACTER(*), PARAMETER :: CALC_TIME_DESCRIPTION        = 'Calculation time'
  CHARACTER(*), PARAMETER :: ANCILLARY_DESCRIPTION        = 'ancillary'
  ! ...The run flags
  CHARACTER(*), PARAMETER :: HIRAC_DESCRIPTION = 'LBLRTM control - hirac run flag'
  CHARACTER(*), PARAMETER :: LBLF4_DESCRIPTION = 'LBLRTM control - lblf4 run flag'
  CHARACTER(*), PARAMETER :: XSCNT_DESCRIPTION = 'LBLRTM control - xscnt run flag'
  CHARACTER(*), PARAMETER :: AERSL_DESCRIPTION = 'LBLRTM control - aersl run flag'
  CHARACTER(*), PARAMETER :: EMIT_DESCRIPTION  = 'LBLRTM control - emit run flag'
  CHARACTER(*), PARAMETER :: SCAN_DESCRIPTION  = 'LBLRTM control - scan run flag'
  CHARACTER(*), PARAMETER :: PLOT_DESCRIPTION  = 'LBLRTM control - plot run flag'
  CHARACTER(*), PARAMETER :: PATH_DESCRIPTION  = 'LBLRTM control - path run flag'
  CHARACTER(*), PARAMETER :: JRAD_DESCRIPTION  = 'LBLRTM control - jrad run flag'
  CHARACTER(*), PARAMETER :: TEST_DESCRIPTION  = 'LBLRTM control - test run flag'
  CHARACTER(*), PARAMETER :: MERGE_DESCRIPTION = 'LBLRTM control - merge run flag'
  CHARACTER(*), PARAMETER :: SCNID_DESCRIPTION = 'LBLRTM control - scnid run flag'
  CHARACTER(*), PARAMETER :: HWHM_DESCRIPTION  = 'LBLRTM control - hwhm run flag'
  CHARACTER(*), PARAMETER :: IDABS_DESCRIPTION = 'LBLRTM control - idabs run flag'
  CHARACTER(*), PARAMETER :: ATM_DESCRIPTION   = 'LBLRTM control - atm run flag'
  CHARACTER(*), PARAMETER :: LAYR1_DESCRIPTION = 'LBLRTM control - layr1 run flag'
  CHARACTER(*), PARAMETER :: NLAYR_DESCRIPTION = 'LBLRTM control - nlayr run flag'

  ! Variable units attribute.
  CHARACTER(*), PARAMETER :: UNITS_ATTNAME = 'units'

  CHARACTER(*), PARAMETER :: USER_ID_UNITS          = 'N/A'
  CHARACTER(*), PARAMETER :: COL_SCALE_FACTOR_UNITS = 'N/A'
  CHARACTER(*), PARAMETER :: AVG_LAYER_PRES_UNITS   = 'hPa'
  CHARACTER(*), PARAMETER :: AVG_LAYER_TEMP_UNITS   = 'K'
  CHARACTER(*), PARAMETER :: MOL_ID_UNITS           = 'N/A'
  CHARACTER(*), PARAMETER :: MOL_COL_DENS_UNITS     = 'mol/cm^2'
  CHARACTER(*), PARAMETER :: BROAD_COL_DENS_UNITS   = 'mol/cm^2'
  CHARACTER(*), PARAMETER :: FREQ_INTERVAL_UNITS    = 'cm^-1'
  CHARACTER(*), PARAMETER :: BEGIN_FREQ_UNITS       = 'cm^-1'
  CHARACTER(*), PARAMETER :: END_FREQ_UNITS         = 'cm^-1'
  CHARACTER(*), PARAMETER :: BDRY_TEMP_UNITS        = 'K'
  CHARACTER(*), PARAMETER :: BDRY_EMIS_UNITS        = 'dimensionless'
  CHARACTER(*), PARAMETER :: N_MOLECULES_UNITS      = 'N/A'
  CHARACTER(*), PARAMETER :: N_LAYER_UNITS          = 'N/A'
  CHARACTER(*), PARAMETER :: OD_LAYER_FLAG_UNITS    = 'N/A'
  CHARACTER(*), PARAMETER :: CALC_DATE_UNITS        = 'N/A'
  CHARACTER(*), PARAMETER :: CALC_TIME_UNITS        = 'N/A'
  CHARACTER(*), PARAMETER :: ANCILLARY_UNITS        = 'N/A'
  ! ...The run flags
  CHARACTER(*), PARAMETER :: HIRAC_UNITS = 'N/A'
  CHARACTER(*), PARAMETER :: LBLF4_UNITS = 'N/A'
  CHARACTER(*), PARAMETER :: XSCNT_UNITS = 'N/A'
  CHARACTER(*), PARAMETER :: AERSL_UNITS = 'N/A'
  CHARACTER(*), PARAMETER :: EMIT_UNITS  = 'N/A'
  CHARACTER(*), PARAMETER :: SCAN_UNITS  = 'N/A'
  CHARACTER(*), PARAMETER :: PLOT_UNITS  = 'N/A'
  CHARACTER(*), PARAMETER :: PATH_UNITS  = 'N/A'
  CHARACTER(*), PARAMETER :: JRAD_UNITS  = 'N/A'
  CHARACTER(*), PARAMETER :: TEST_UNITS  = 'N/A'
  CHARACTER(*), PARAMETER :: MERGE_UNITS = 'N/A'
  CHARACTER(*), PARAMETER :: SCNID_UNITS = 'N/A'
  CHARACTER(*), PARAMETER :: HWHM_UNITS  = 'N/A'
  CHARACTER(*), PARAMETER :: IDABS_UNITS = 'N/A'
  CHARACTER(*), PARAMETER :: ATM_UNITS   = 'N/A'
  CHARACTER(*), PARAMETER :: LAYR1_UNITS = 'N/A'
  CHARACTER(*), PARAMETER :: NLAYR_UNITS = 'N/A'

  ! Variable _FillValue attribute.
  CHARACTER(*), PARAMETER :: FILLVALUE_ATTNAME = '_FillValue'

  CHARACTER(*), PARAMETER :: USER_ID_FILLVALUE          = NF90_FILL_CHAR
  REAL(DP)    , PARAMETER :: COL_SCALE_FACTOR_FILLVALUE = 0.0_DP
  REAL(FP)    , PARAMETER :: AVG_LAYER_PRES_FILLVALUE   = 0.0_FP
  REAL(FP)    , PARAMETER :: AVG_LAYER_TEMP_FILLVALUE   = 0.0_FP
  CHARACTER(*), PARAMETER :: MOL_ID_FILLVALUE           = NF90_FILL_CHAR
  REAL(FP)    , PARAMETER :: MOL_COL_DENS_FILLVALUE     = 0.0_FP
  REAL(FP)    , PARAMETER :: BROAD_COL_DENS_FILLVALUE   = 0.0_FP
  REAL(FP)    , PARAMETER :: FREQ_INTERVAL_FILLVALUE    = 0.0_FP
  REAL(DP)    , PARAMETER :: BEGIN_FREQ_FILLVALUE       = 0.0_DP
  REAL(DP)    , PARAMETER :: END_FREQ_FILLVALUE         = 0.0_DP
  REAL(FP)    , PARAMETER :: BDRY_TEMP_FILLVALUE        = 0.0_FP
  REAL(FP)    , PARAMETER :: BDRY_EMIS_FILLVALUE        = 0.0_FP
  INTEGER(IP) , PARAMETER :: N_MOLECULES_FILLVALUE      = 0_IP
  INTEGER(IP) , PARAMETER :: N_LAYER_FILLVALUE          = 0_IP
  INTEGER(IP) , PARAMETER :: OD_LAYER_FLAG_FILLVALUE    = 0_IP
  CHARACTER(*), PARAMETER :: CALC_DATE_FILLVALUE        = NF90_FILL_CHAR
  CHARACTER(*), PARAMETER :: CALC_TIME_FILLVALUE        = NF90_FILL_CHAR
  CHARACTER(*), PARAMETER :: ANCILLARY_FILLVALUE        = NF90_FILL_CHAR
  ! ...The run flags
  INTEGER(IP), PARAMETER :: HIRAC_FILLVALUE = 0_IP
  INTEGER(IP), PARAMETER :: LBLF4_FILLVALUE = 0_IP
  INTEGER(IP), PARAMETER :: XSCNT_FILLVALUE = 0_IP
  INTEGER(IP), PARAMETER :: AERSL_FILLVALUE = 0_IP
  INTEGER(IP), PARAMETER :: EMIT_FILLVALUE  = 0_IP
  INTEGER(IP), PARAMETER :: SCAN_FILLVALUE  = 0_IP
  INTEGER(IP), PARAMETER :: PLOT_FILLVALUE  = 0_IP
  INTEGER(IP), PARAMETER :: PATH_FILLVALUE  = 0_IP
  INTEGER(IP), PARAMETER :: JRAD_FILLVALUE  = 0_IP
  INTEGER(IP), PARAMETER :: TEST_FILLVALUE  = 0_IP
  INTEGER(IP), PARAMETER :: MERGE_FILLVALUE = 0_IP
  REAL(FP)   , PARAMETER :: SCNID_FILLVALUE = 0.0_FP
  REAL(FP)   , PARAMETER :: HWHM_FILLVALUE  = 0.0_FP
  INTEGER(IP), PARAMETER :: IDABS_FILLVALUE = 0_IP
  INTEGER(IP), PARAMETER :: ATM_FILLVALUE   = 0_IP
  INTEGER(IP), PARAMETER :: LAYR1_FILLVALUE = 0_IP
  INTEGER(IP), PARAMETER :: NLAYR_FILLVALUE = 0_IP

  ! Variable netCDF datatypes
  INTEGER(Long), PARAMETER :: USER_ID_TYPE          = NF90_CHAR
  INTEGER(Long), PARAMETER :: COL_SCALE_FACTOR_TYPE = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: AVG_LAYER_PRES_TYPE   = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: AVG_LAYER_TEMP_TYPE   = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: MOL_ID_TYPE           = NF90_CHAR
  INTEGER(Long), PARAMETER :: MOL_COL_DENS_TYPE     = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: BROAD_COL_DENS_TYPE   = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: FREQ_INTERVAL_TYPE    = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: BEGIN_FREQ_TYPE       = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: END_FREQ_TYPE         = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: BDRY_TEMP_TYPE        = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: BDRY_EMIS_TYPE        = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: N_MOLECULES_TYPE      = NF90_INT
  INTEGER(Long), PARAMETER :: N_LAYER_TYPE          = NF90_INT
  INTEGER(Long), PARAMETER :: OD_LAYER_FLAG_TYPE    = NF90_INT
  INTEGER(Long), PARAMETER :: CALC_DATE_TYPE        = NF90_CHAR
  INTEGER(Long), PARAMETER :: CALC_TIME_TYPE        = NF90_CHAR
  INTEGER(Long), PARAMETER :: ANCILLARY_TYPE        = NF90_CHAR
  ! ...The run flags
  INTEGER(Long), PARAMETER :: HIRAC_TYPE = NF90_INT
  INTEGER(Long), PARAMETER :: LBLF4_TYPE = NF90_INT
  INTEGER(Long), PARAMETER :: XSCNT_TYPE = NF90_INT
  INTEGER(Long), PARAMETER :: AERSL_TYPE = NF90_INT
  INTEGER(Long), PARAMETER :: EMIT_TYPE  = NF90_INT
  INTEGER(Long), PARAMETER :: SCAN_TYPE  = NF90_INT
  INTEGER(Long), PARAMETER :: PLOT_TYPE  = NF90_INT
  INTEGER(Long), PARAMETER :: PATH_TYPE  = NF90_INT
  INTEGER(Long), PARAMETER :: JRAD_TYPE  = NF90_INT
  INTEGER(Long), PARAMETER :: TEST_TYPE  = NF90_INT
  INTEGER(Long), PARAMETER :: MERGE_TYPE = NF90_INT
  INTEGER(Long), PARAMETER :: SCNID_TYPE = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: HWHM_TYPE  = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: IDABS_TYPE = NF90_INT
  INTEGER(Long), PARAMETER :: ATM_TYPE   = NF90_INT
  INTEGER(Long), PARAMETER :: LAYR1_TYPE = NF90_INT
  INTEGER(Long), PARAMETER :: NLAYR_TYPE = NF90_INT

CONTAINS

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

  !----------------------------------------------------------
  ! Function to write an LBLRTM File Header object as a group
  !----------------------------------------------------------
  FUNCTION LBLRTM_Fhdr_netCDF_WriteGroup( &
    Fhdr     , &  ! Input
    FileId   , &  ! Input
    GroupName, &  ! Optional input
    Quiet    , &  ! Optional input
    Debug    ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    TYPE(LBLRTM_Fhdr_type), INTENT(IN) :: Fhdr
    INTEGER(Long),          INTENT(IN) :: FileId
    CHARACTER(*), OPTIONAL, INTENT(IN) :: GroupName
    LOGICAL,      OPTIONAL, INTENT(IN) :: Quiet
    LOGICAL,      OPTIONAL, INTENT(IN) :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'LBLRTM_Fhdr_netCDF_IO::WriteGroup'
    ! Local variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: group_name
    LOGICAL :: noisy
    LOGICAL :: debug_output
    INTEGER(Long) :: nf90_stat
    INTEGER(Long) :: groupid
    INTEGER(Long) :: n_mol_dimid, n_ancillary_dimid
    INTEGER(Long) :: uid_strlen_dimid, sl_strlen_dimid

    ! Setup
    err_stat = SUCCESS
    ! ...Check structure
    IF ( .NOT. (LBLRTM_Fhdr_IsValid( Fhdr )) ) THEN
      msg = 'LBLRTM Fhdr object is invalid. Nothing to do!'
      CALL Write_CleanUp(); RETURN
    END IF
    ! ...Check GroupName argument, defining default.
    group_name = 'Fhdr'
    IF ( PRESENT(GroupName) ) THEN
      group_name = ADJUSTL(GroupName)
    END IF
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Set debug option
    debug_output = .FALSE.
    IF ( PRESENT(debug) ) debug_output = debug
    IF ( debug_output ) THEN
      CALL Display_Message(ROUTINE_NAME,'Entering...',INFORMATION)
      noisy = .TRUE.
    END IF

    ! Create a new group for the file header data
    nf90_stat = NF90_DEF_GRP( &
      fileid, &
      group_name, &
      groupid )
    IF ( nf90_stat /= NF90_NOERR ) THEN
      msg = 'Error creating '//TRIM(group_name)//' group - '//&
            ' - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Define the dimensions for the group
    err_stat = DefineDimensions( &
      Fhdr             , &
      groupid          , &
      n_mol_dimid      , &
      n_ancillary_dimid, &
      uid_strlen_dimid , &
      sl_strlen_dimid    )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error defining dimensions for the '//TRIM(group_name)//&
            ' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Define the variables for the group
    err_stat = DefineVariables( &
      groupid          , &
      n_mol_dimid      , &
      n_ancillary_dimid, &
      uid_strlen_dimid , &
      sl_strlen_dimid    )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error defining variables for the '//TRIM(group_name)//&
            ' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Take netCDF file out of define mode
    nf90_stat = NF90_ENDDEF( fileid )
    IF ( nf90_stat /= NF90_NOERR ) THEN
      msg = 'Error taking file out of define mode to write the '//&
            TRIM(group_name)//' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the variables for the group
    err_stat = WriteVariables( Fhdr, groupid )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error writing variables for the '//TRIM(group_name)//&
            ' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Put netCDF file back into define mode
    nf90_stat = NF90_REDEF( fileid )
    IF ( nf90_stat /= NF90_NOERR ) THEN
      msg = 'Error putting file back into define mode after writing the '//&
            TRIM(group_name)//' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE Write_CleanUp()
      nf90_stat = NF90_CLOSE( fileid )
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE Write_CleanUp

  END FUNCTION LBLRTM_Fhdr_netCDF_WriteGroup

  !---------------------------------------------------------
  ! Function to read an LBLRTM File Header object as a group
  !---------------------------------------------------------
  FUNCTION LBLRTM_Fhdr_netCDF_ReadGroup( &
    Fhdr     , &  ! Output
    FileId   , &  ! Input
    GroupName, &  ! Optional input
    Quiet    , &  ! Optional input
    Debug    ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    TYPE(LBLRTM_Fhdr_type), INTENT(OUT) :: Fhdr
    INTEGER(Long),          INTENT(IN)  :: FileId
    CHARACTER(*), OPTIONAL, INTENT(IN)  :: GroupName
    LOGICAL,      OPTIONAL, INTENT(IN)  :: Quiet
    LOGICAL,      OPTIONAL, INTENT(IN)  :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'LBLRTM_Fhdr_netCDF_IO::ReadGroup'
    ! Local variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: group_name
    LOGICAL :: noisy
    LOGICAL :: debug_output
    INTEGER(Long) :: nf90_stat
    INTEGER(Long) :: groupid

    ! Setup
    err_stat = SUCCESS
    ! ...Check GroupName argument, defining default.
    group_name = 'Fhdr'
    IF ( PRESENT(GroupName) ) THEN
      group_name = ADJUSTL(GroupName)
    END IF
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Set debug option
    debug_output = .FALSE.
    IF ( PRESENT(debug) ) debug_output = debug
    IF ( debug_output ) THEN
      CALL Display_Message(ROUTINE_NAME,'Entering...',INFORMATION)
      noisy = .TRUE.
    END IF

    ! Get the group id
    nf90_stat = NF90_INQ_GRP_NCID(fileid, group_name, groupid)
    IF ( nf90_stat /= NF90_NOERR ) THEN
      msg = 'Error inquiring '//TRIM(group_name)//' group for its group id - '//&
            TRIM(NF90_STRERROR( nf90_stat ))
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read the variables for the group
    err_stat = ReadVariables( Fhdr, groupid )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading variables for the '//TRIM(group_name)//&
            ' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Read_Cleanup(); RETURN
    END IF

    ! Tag object as valid
    CALL LBLRTM_Fhdr_SetValid(Fhdr)
    IF ( debug_output ) CALL LBLRTM_Fhdr_Inspect(fhdr)

  CONTAINS

    SUBROUTINE Read_CleanUp()
      CALL LBLRTM_Fhdr_Destroy(Fhdr)
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE Read_CleanUp

  END FUNCTION LBLRTM_Fhdr_netCDF_ReadGroup

  !------------------------------------------------
  ! Subroutine to return module version information
  !------------------------------------------------
  SUBROUTINE LBLRTM_Fhdr_netCDF_IOVersion( Id )
    CHARACTER(*), INTENT(OUT) :: Id
    Id = MODULE_VERSION_ID
  END SUBROUTINE LBLRTM_Fhdr_netCDF_IOVersion

!################################################################################
!################################################################################
!##                                                                            ##
!##                          ## PRIVATE MODULE ROUTINES ##                     ##
!##                                                                            ##
!################################################################################
!################################################################################

  INCLUDE 'LBLRTM_Fhdr_netCDF_IO.inc'

END MODULE LBLRTM_Fhdr_netCDF_IO

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
R.I.P.

Created by Users/ScottBiggart and Tim Carson (he XhasX had a http://www.timxdrama.com site with erotic pictures) during the summer of 2002, DavisCoeds is a social profile site.  Like many sites created in the image of http://www.friendster.com FriendsterEDIT: Friendster and DavisCoeds.com went public within a month of each other, http://www.daviscoeds.com DavisCoeds allows for a profile (for listing interests and other things) and a set of pictures for each member.  Unlike http://www.friendster.com Friendster or Facebook, popular socialnetworking sites, the aim of Coeds isnt to weave an electronic social network  hence there is no friends list.

DavisCoeds had seen two major updates since inception.  The site was originally written in ASP with a flatfile backend. The first major update was a rewrite in PHP using a MySQL backend.  The last update, which took place February 5th, 2006, featured layout changes (primarily to user profiles) as well as a rewritten PHP backend.  Daviscoeds is in the know about the Wiki, and allows you to perform a search based on a Daviswiki username!

Davis Coeds appears to have gone offline as of August 2007 after the domain name was swiped out from under Users/ScottBiggart.

As of October 26th, 2011, the domain has been returned....

One could say that the purpose of http://www.daviscoeds.com DavisCoeds is just to put up some pictures and meet people.  Most use it to meet people; a few use it to try and hook up with girls in the creepiest way possible; Users/ScottBiggart Scott uses it as a pool of Freshman girls to make out with  I think hes made out with... 90% of the girls on that site and god only knows how many guys.  Users/PhilipNeustrom

Another person could say, that if you arent white you wont get any on that website. An AfroAmerican buddy and myself being Indian have gotten nothing from that worthless website. Most people on daviscoeds are worthless. Its cuz we are brown and black. Racists!  Users/ZakirKhan

DavisCoeds is a meeting person website. Not a dating website. I do not believe that the site is racist. Scott Biggart is most certainly not. There is a diverse group of members on the site and if you go into the site thinking that you are going to get some then that may be the very reason why you are upset with the site. It has nothing to do will ethnicity and everything to do with personality.  Users/RobRoy

Yeah, Ive been hit on by some really creepy, really horny guys that had the social skills of a turnip...  Users/MarieHuynh

You should learn to take a joke Rob. I still think its racist. Unfortunately, some people will never see the humor in racism! Therefore I implore you to go see a play done by my friends the play is called Ngger Wtback Chnk. If you find that statement in itself offensive, then see the play before you judge it.http://www.dlartists.com/pdf/SP_032104_LATimes.pdf  Users/ZakirKhan

20070415 23:15:28 nbsp Saying that the creeps on there have the social skills of a turnip is an insult to turnips everywhere. Users/AMcKean

20070509 13:55:59 nbsp You dont have to be white to get any from Davis Coeds. Youll get plenty attention if youre an Asian female. Users/ThuHo

20070821 18:30:58 nbsp I have room on my server to host DavisCoeds.  MYSQL 4.1 and PHP5, Id host it for free. Users/WilliamBeamish

20080328 20:51:19 nbsp Its too bad too many people had to abuse the site or attempt to use it to get laid.  It was a spiffy way of meeting new people.  Ive made quite a few friends cause of that website that I still talk to.  Its a shame it had to go.  And in reference to above comments, Im Indian and had no issues meeting people.  But then again, my primary use of the website was not to get any from someone I met online (thats a bit too creepycraigslistsounding to me). Users/SunjeetBaadkar


text:
C
C	$Id: ffphol.f,v 1.4 2008-07-27 12:23:42 haley Exp $
C
C			     Copyright (C)  1997
C	     University Corporation for Atmospheric Research
C			     All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C

      SUBROUTINE FFPHOL(IER)
C
C  Process holes.
C
      include 'fntcom.h'
C
      DIMENSION IVALS(10)
C
C  Put out the begin hole flag.
C
      CALL FFPPKT(BITPNT,BEGINH,0,0)
      BITPNT = BITPNT+PKWID
   40 CONTINUE
      CALL CFRDLN(IDUM1,IDUM2,IER)
      IF (LINE(1:3) .NE. 'END') THEN
C
C  Has to be coordinate data or Bezier control points to get here.
C
        CALL FFTKIN(LINE,NUMI,IVALS,IER)
        IF (IER .NE. 0) GO TO 30
        IF (NUMI .EQ. 2) THEN
C
C  Coordinate data.
C
          IVALS(1) = IVALS(1)+XBIAS
          IVALS(2) = IVALS(2)+YBIAS
          CALL FFPPKT(BITPNT,COORD,IVALS(1),IVALS(2))
          BITPNT = BITPNT+PKWID
          GO TO 40
        ELSE IF (NUMI .EQ. 6) THEN
C
C  Bezier control points.
C

          IVALS(1) = IVALS(1)+XBIAS
          IVALS(2) = IVALS(2)+YBIAS
          CALL FFPPKT(BITPNT,BEZIER,IVALS(1),IVALS(2))
          BITPNT = BITPNT+PKWID
          DO 60 J=2,3
            IX = 2*J-1
            IY = IX+1
            IVALS(IX) = IVALS(IX)+XBIAS
            IVALS(IY) = IVALS(IY)+YBIAS
            CALL FFPPKT(BITPNT,COORD,IVALS(IX),IVALS(IY))
            BITPNT = BITPNT+PKWID
   60     CONTINUE
          GO TO 40
        ELSE
          IER = EORD
          GO TO 30
        ENDIF
      ELSE IF (LINE(1:10) .EQ. 'END_HOLE') THEN
C
C  Put out end hole flag.
C
        CALL FFPPKT(BITPNT,ENDH,0,0)
        BITPNT = BITPNT+PKWID
        RETURN
      ELSE
        IER = EORD
        GO TO 30
      ENDIF
C
   30 CONTINUE
C
      RETURN
      END


text:
module errs

use, intrinsic:: iso_fortran_env, only: error_unit
use cinter, only: endwin

implicit none (type, external)

contains

subroutine err(msg)

character(*),intent(in) :: msg

call endwin()
write(error_unit,*) msg
error stop

end subroutine err

end module

max_stars_repo_name: fortran-gaming/tetran

max_issues_repo_name: fortran-gaming/blocktran

max_forks_repo_name: fortran-gaming/blocktran


text:
      program w17

C     Check refined transformer for while loop with anded condition

      integer i

      i = 0

      if(.TRUE.) then
         do while(i.gt.-5.and.i.lt.5)
            if(y.gt.0.) then
               i = i + 1
            else
               i = i - 1
            endif
         enddo
      endif

      print *, i

      end