text:
module Mod_sldsup_elmope_lin
  use Mod_CauchyElement
  use Mod_SUPSolids_lin
  use Mod_sld_BaseElmopeRoutines
  implicit none

  contains

  subroutine SetPointersElmopeSUP_lin
     use typre
     implicit none
     procedure() :: NULLSUB
     integer(ip) :: kfl_nonlinear,nelty

     call ConcatenateProcedures(ProcHook%AllocateArrays  ,AllocateBaseElmopeArrays)
     call ConcatenateProcedures(ProcHook%DeallocateArrays,DeallocateBaseElmopeArrays)

     call ConcatenateProcedures(ProcHook%AllocateArrays  ,AllocateAssemblyArrays)
     call ConcatenateProcedures(ProcHook%DeallocateArrays,DeallocateAssemblyArrays)

     call ConcatenateProcedures(ProcHook%AllocateArrays,AllocateBaseElmopeMatricesSUP)
     call ConcatenateProcedures(ProcHook%DeallocateArrays,DeallocateBaseElmopeMatricesSUP)

     call ConcatenateProcedures(ProcHook%AllocateArrays,AllocateBaseElmopeMatricesSUPLinear)
     call ConcatenateProcedures(ProcHook%DeallocateArrays,DeallocateBaseElmopeMatricesSUPLinear)

     call ConcatenateProcedures(ProcHook%AllocateArrays,AllocateBaseElmopeMatricesSUP_SGS)
     call ConcatenateProcedures(ProcHook%DeallocateArrays,DeallocateBaseElmopeMatricesSUP_SGS)

     call ConcatenateProcedures(ProcHook%AllocateArrays,AllocateSUPAssemblyMatrices)
     call ConcatenateProcedures(ProcHook%DeallocateArrays,DeallocateSUPAssemblyMatrices)

     call ConcatenateProcedures(ProcHook%ResetArrays,ResetSUPBaseArrays)
     call ConcatenateProcedures(ProcHook%ResetArrays,ResetSUPBaseArrays_SGS)

     call ConcatenateProcedures(ProcHook%ResetArrays,ResetAssemblyMatrices)

     call ConcatenateProcedures(ProcHook%Interpolates,InterpolateGpDisplacements)
     call ConcatenateProcedures(ProcHook%Interpolates,InterpolateGpSigma)
     call ConcatenateProcedures(ProcHook%Interpolates,InterpolateGpPress)

     call SetPointers_Linear

     if(a%kfl_linea == 0 ) then      !Standard RHS linearization
         call SetRHSPointers
     else if(a%kfl_linea == 2 ) then !Newton-Raphson linearization
         !call SetNRPointers
         write(*,*) "sldsup: NR linearization not ready, setting RHS"
         call SetRHSPointers
     endif

     ProcPointer%getTauParameters => getTauParameters_linear

     call ConcatenateProcedures(ProcHook%AssemblyLhs,AssembleStiffnessMatrixSUP)
     call ConcatenateProcedures(ProcHook%AssemblyLhs,AssembleStiffnessMatrixSUP_SGS)
     call ConcatenateProcedures(ProcHook%AssemblyRhs,AssembleRhsSUP)

     call ConcatenateProcedures(ProcHook%InGauss,InGaussResetExternal)

     !--------------------------------------------------------------------------
     if(sup_lin%kfl_timei==1) then
         !Dynamic allocations
         call ConcatenateProcedures(ProcHook%AllocateArrays  ,AllocateDynamicArray)
         call ConcatenateProcedures(ProcHook%DeallocateArrays,DeallocateDynamicArray)

         !Dynamic resets
         call ConcatenateProcedures(ProcHook%ResetArrays,ResetDynamicMatrix)

         !Dynamic elemental build
         call ConcatenateProcedures(ProcHook%DynamicMass,buildMassMat)

         !Dynamic build to LHS
         call ConcatenateProcedures(ProcHook%AssemblyLhs,AssembleMassMatrix)

         call ConcatenateProcedures(ProcHook%Initializations,AllocVelandAccel)
         call ConcatenateProcedures(ProcHook%Initializations,InitTimeIntegration)
         call ConcatenateProcedures(ProcHook%Gathers        ,GatherDynamic)
         call ConcatenateProcedures(ProcHook%Finalizations  ,DeallocVelandAccel)

         !Dynamic matrix routines
         call ConcatenateProcedures(ProcHook%Dynamic,ProcHook%DynamicForce)
         call ConcatenateProcedures(ProcHook%Dynamic,ProcHook%DynamicMass)

     endif

     !Non-linear elements
     call sup_lin%Mesh%IsNonLinearElement(ielem,kfl_nonlinear)
     if (kfl_nonlinear == 1) then
         call ConcatenateProcedures(ProcHook%InGaussElmats,ProcPointer%PostGauss)
         ProcPointer%PostGauss => NULLSUB
     end if

     call ConcatenateProcedures(ProcHook%ToLinearSystem,AssembleLinearSystem)
     !-------------------------------------------------------
     !If more than one element type, then pointers should be reset when ielty changes!
     call a%Mesh%GetNelty(nelty)
     if (nelty > 1) ProcHook%OnIeltyChange => OnIeltyChange

    !Mass matrix
    call ConcatenateProcedures(ProcHook%PostInterpolates,ProcHook%Dynamic)
    !K_mat and K_geo, and internal forces
    call ConcatenateProcedures(ProcHook%PostInterpolates,ProcHook%InGaussElmats)
    call ConcatenateProcedures(ProcHook%PostInterpolates,ProcHook%InGaussVec)

  end subroutine SetPointersElmopeSUP_lin

  subroutine SetPointers_Linear
     use typre
     implicit none

     !Set the Specific Pointers for this subroutine
     ProcPointer%PostGauss  => PostGaussElmats_linear
     ProcPointer%sld_elmrhu => sld_elmrhu

     call ConcatenateProcedures(ProcHook%PhysicalProp,PhysicalProp)

     call ConcatenateProcedures(ProcHook%Gathers,sigmaGather)
     call ConcatenateProcedures(ProcHook%Gathers,displacementGather)
     call ConcatenateProcedures(ProcHook%Gathers,pressGather)

  end subroutine SetPointers_Linear

  !----------------------------------------------------------
  !PostGauss Matrices
  subroutine PostGaussElmats_linear
     implicit none
     integer(ip) :: nd,tn

     call sup_lin%Mesh%GetNdime(nd)
     tn=(nd*(nd+1))/2

     call ProcHook%PhysicalProp

     call sup_elmVS(e,nd,tn,dvol,P,     elmsv)
     call sup_elmVP(e,nd,   dvol,       elmpv)
     call sup_elmEU(e,nd,tn,dvol,P,     elmut)
     call sup_elmES(e,nd,tn,dvol,D,     elmst)
     call sup_elmQU(e,nd,   dvol,       elmuq)
     call sup_elmQP(e,nd,   dvol,Dv_scalar,elmpq)

     !-------------------Tau_u matrix calculation----------
     call sup_elm_usgs_ES(e,nd,tn,dvol,elm_usgs_ES)
     call sup_elm_usgs_EP(e,nd,tn,dvol,elm_usgs_EP)
     call sup_elm_usgs_QS(e,nd,tn,dvol,elm_usgs_QS)
     call sup_elm_usgs_QP(e,nd,tn,dvol,elm_usgs_QP)

     !-------------------Tau_s matrix calculation----------
     call sup_elm_ssgs_VU(e,nd,tn,dvol,C_dev,elm_ssgs_VU)
     call sup_elm_ssgs_VS(e,nd,tn,dvol,P    ,elm_ssgs_VS)
     call sup_elm_ssgs_EU(e,nd,tn,dvol,P    ,elm_ssgs_EU)
     call sup_elm_ssgs_ES(e,nd,tn,dvol,D    ,elm_ssgs_ES)

     !-------------------Tau_p matrix calculation----------
     call sup_elm_psgs_VU(e,nd,   dvol,K,G      ,elm_psgs_VU)
     call sup_elm_psgs_VP(e,nd,   dvol          ,elm_psgs_VP)
     call sup_elm_psgs_QU(e,nd,   dvol          ,elm_psgs_QU)
     call sup_elm_psgs_QP(e,nd,   dvol,Dv_scalar,elm_psgs_QP)
  end subroutine

  subroutine PhysicalProp
      implicit none
      integer(ip) :: nd,tn

      call sup_lin%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2

      call sup_lin%SUPGetPhysicalParameters(nd,tn,densi,K,G,C_dev,D,D_dev,Dv_scalar,P)
  end subroutine

  !Gauss matrices
  subroutine InGaussElmats_linear
      implicit none

  end subroutine

  subroutine OnIeltyChange
     implicit none

     if (e%ielty /= ielty0) then
        call SetPointersElmopeSUP_lin
        ielty0 = e%ielty
     endif
  end subroutine

end module

!Solids_ELMOPE subroutine
subroutine sldsup_Elmope_lin(SldProblem)
   use Mod_sldsup_elmope_lin

   implicit none
   class(SUPSolidsProblem_lin), target :: SldProblem

   !Setting sld problem pointer in sld_BaseElmope
   up      => SldProblem
   sup_lin => SldProblem
   sup     => SldProblem
   up      => SldProblem
   a       => SldProblem

   call SetPointersGeneral
   call SetPointersPointForces
   call SetPointersElmopeSUP_lin

   call sld_elemLoop

end subroutine sldsup_Elmope_lin

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
PROGRAM memcpy_test

USE hip_fortran

IMPLICIT NONE

  REAL(8), ALLOCATABLE, TARGET :: a(:,:)
  TYPE(c_ptr) :: a_dev = c_null_ptr
  INTEGER(c_size_t) :: a_size

    ! Allocate and initialize host a
    ALLOCATE(a(0:10,0:10))
    a = 10.0D0
    a_size = SIZEOF(a)

    ! Allocate device a
    CALL hfMalloc(a_dev, a_size)

    ! Copy host memory to device memory
    CALL hipFortran(hipMemcpy(a_dev, c_loc(a), SIZEOF(a), hipMemcpyHostToDevice))

    ! Copy device memory to host memory
    CALL hipFortran(hipMemcpy(c_loc(a), a_dev, SIZEOF(a), hipMemcpyDeviceToHost))

    ! Clean up host and device memory
    CALL hfFree(a_dev)
    DEALLOCATE(a)

END PROGRAM memcpy_test

max_stars_repo_name: dmageeLANL/hip-fortran

max_issues_repo_name: dmageeLANL/hip-fortran

max_forks_repo_name: dmageeLANL/hip-fortran


text:
      program main
c
c ... array declarations.
c
      implicit double precision (a-h, o-z)
      dimension coef(120,5), rhs(100), u(100), wksp(600), ubar(1),
     a     rparm(30)
      integer jcoef(120,5), iwksp(300), iparm(30), p(100), ip(100)
      external cg, rs6
c
      ndim = 120
      mdim = 5
      nw = 600
      inw = 300
c
c ... generate coef, jcoef, and rhs.
c
      nx = 10
      ny = 10
      n = nx*ny
      h = 1.0d0/dble(nx + 1)
      maxnz = 5
      do i = 1,n
         coef(i,1) = 6.0d0
         coef(i,2) = -1.0d0
         coef(i,3) = -2.0d0
         coef(i,4) = -1.0d0
         coef(i,5) = -2.0d0
         jcoef(i,1) = i
         jcoef(i,2) = i + 1
         jcoef(i,3) = i + nx
         jcoef(i,4) = i - 1
         jcoef(i,5) = i - nx
         rhs(i) = 0.0d0
      enddo
      k = 0
      do j = 1,ny
         y = dble(j)*h
         do i = 1,nx
            x = dble(i)*h
            k = k + 1
            if (j .eq. 1) then
               rhs(k) = rhs(k) + 2.0d0
               coef(k,5) = 0.0d0
               jcoef(k,5) = 0
            endif
            if (j .eq. ny) then
               rhs(k) = rhs(k) + 2.0d0*(1.0d0 + x)
               coef(k,3) = 0.0d0
               jcoef(k,3) = 0
            endif
            if (i .eq. 1) then
               rhs(k) = rhs(k) + 1.0d0
               coef(k,4) = 0.0d0
               jcoef(k,4) = 0
            endif
            if (i .eq. nx) then
               rhs(k) = rhs(k) + 1.0d0 + y
               coef(k,2) = 0.0d0
               jcoef(k,2) = 0
            endif
         enddo
      enddo
      call dfault (iparm,rparm)
c
c ... now, reset some default values.
c
      iparm(3) = 3
      iparm(12) = 1
      iparm(14) = 1
      iparm(23) = 0
c
c ... generate an initial guess for u and call nspcg.
c
      call vfill (n,u,0.0d0)
c
      call redblk (ndim,n,maxnz,coef,jcoef,p,ip,1,iwksp,ier)
      call nspcg (rs6,cg,ndim,mdim,n,maxnz,coef,jcoef,p,ip,
     a            u,ubar,rhs,wksp,iwksp,nw,inw,iparm,rparm,ier)
      stop
      end

ext: f


text:
      SUBROUTINE W3FA03V(PRESS,HEIGHT,TEMP,THETA,N)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    W3FA03V     COMPUTE STANDARD HEIGHT, TEMP, AND POT TEMP
C   PRGMMR: KEYSER           ORG: W/NMC22    DATE: 92-06-29
C
C ABSTRACT: COMPUTES THE STANDARD HEIGHT, TEMPERATURE, AND POTENTIAL
C   TEMPERATURE GIVEN THE PRESSURE IN MILLIBARS ( > 8.68 MB ).  FOR
C   HEIGHT AND TEMPERATURE THE RESULTS DUPLICATE THE VALUES IN THE
C   U.S. STANDARD ATMOSPHERE (L962), WHICH IS THE ICAO STANDARD
C   ATMOSPHERE TO 54.7487 MB (20 KM) AND THE PROPOSED EXTENSION TO
C   8.68 MB (32 KM).  FOR POTENTIAL TEMPERATURE A VALUE OF  2/7  IS
C   USED FOR  RD/CP.
C
C PROGRAM HISTORY LOG:
C   74-06-01  J. MCDONELL W345     -- ORIGINAL AUTHOR
C   84-06-01  R.E.JONES W342       -- CHANGE TO IBM VS FORTRAN
C   92-06-29  D. A. KEYSER W/NMC22 -- CONVERT TO CRAY CFT77 FORTRAN
C   94-09-13  R.E.JONES            -- VECTORIZED VERSION TO DO ARRAY
C                                     INSTEAD OF ONE WORD
C
C USAGE:    CALL W3FA03V(PRESS,HEIGHT,TEMP,THETA,N)
C   INPUT ARGUMENT LIST:
C     PRESS    - PRESSURE ARRAY IN MILLIBARS
C
C   OUTPUT ARGUMENT LIST:
C     HEIGHT   - HEIGHT ARRAY IN METERS
C     TEMP     - TEMPERATURE ARRAY IN DEGREES KELVIN
C     THETA    - POTENTIAL TEMPERATURE ARRAY IN DEGREES KELVIN
C     N        - NUMBER OF POINTS IN ARRAY PRESS
C
C   SUBPROGRAMS CALLED:
C     LIBRARY:
C       CRAY     - ALOG
C
C REMARKS: NOT VALID FOR PRESSURES LESS THAN 8.68 MILLIBARS, DECLARE
C   ALL PARAMETERS AS TYPE REAL.
C
C WARNING: HEIGHT, TEMP, THETA ARE NOW ALL ARRAYS, YOU MUST
C          HAVE ARRAYS OF SIZE N OR YOU WILL WIPE OUT MEMORY.
C
C ATTRIBUTES:
C   LANGUAGE: CRAY CFT77 FORTRAN
C   MACHINE:  CRAY C916/128, Y-MP8/864, Y-MP EL92/256
C
C$$$
C
      REAL  M0
      REAL  HEIGHT(*)
      REAL  PRESS(*)
      REAL  TEMP(*)
      REAL  THETA(*)
C
      SAVE
C
      DATA  G/9.80665/,RSTAR/8314.32/,M0/28.9644/,PISO/54.7487/,
     $ ZISO/20000./,SALP/-.0010/,PZERO/1013.25/,T0/288.15/,ALP/.0065/,
     $ PTROP/226.321/,TSTR/216.65/
C
      ROVCP = 2.0/7.0
      R     = RSTAR/M0
      ROVG  = R/G
      FKT   = ROVG * TSTR
      AR    = ALP  * ROVG
      PP0   = PZERO**AR
      AR1   = SALP * ROVG
      PP01  = PISO**AR1
C
      DO J = 1,N
      IF (PRESS(J).LT.PISO) THEN
C
C     COMPUTE LAPSE RATE = -.0010 CASES
C
        HEIGHT(J) = ((TSTR/(PP01 * SALP )) * (PP01-(PRESS(J) ** AR1)))
     &              + ZISO
        TEMP(J)   = TSTR - ((HEIGHT(J) - ZISO) * SALP)
C
      ELSE IF (PRESS(J).GT.PTROP) THEN
C
        HEIGHT(J) = (T0/(PP0 * ALP)) * (PP0 - (PRESS(J) ** AR))
        TEMP(J)   = T0 - (HEIGHT(J) * ALP)
C
      ELSE
C
C     COMPUTE ISOTHERMAL CASES
C
        HEIGHT(J) = 11000.0 + (FKT * ALOG(PTROP/PRESS(J)))
        TEMP(J)   = TSTR
C
      END IF
        THETA(J)  = TEMP(J) * ((1000./PRESS(J))**ROVCP)
      END DO
C
      RETURN
      END

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack

max_forks_repo_name: GMAO-SI-Team/jedi-stack


text:
subroutine rdusrp(lunmd     ,lundia    ,error     ,mdfrec    ,nrrec     , &
                & gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: rdusrp.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/input/rdusrp.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: - Reads the user defined processes from the
!                MD-file: NPRCUS and nrpcus times PRCUSR & NPRINP
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer                        , pointer :: itis
    integer                        , pointer :: nprocs
    integer      , dimension(:, :) , pointer :: nprdim
    integer      , dimension(:)    , pointer :: nread
    integer                        , pointer :: nprcus
    integer      , dimension(:, :) , pointer :: nprinp
    integer      , dimension(:)    , pointer :: icousr
    real(fp)     , dimension(:)    , pointer :: rcousr
    character*256, dimension(:)    , pointer :: filusr
    character*20 , dimension(:)    , pointer :: procs
    character*20 , dimension(:)    , pointer :: prcusr
    character*20 , dimension(:)    , pointer :: ccousr
!
! Global variables
!
    integer                    :: lundia !  Description and declaration in inout.igs
    integer                    :: lunmd  !  Description and declaration in inout.igs
    integer                    :: nrrec  !!  Pointer to the record number in the MD-file
    logical      , intent(out) :: error  !!  Flag=TRUE if an error is encountered
    character(*)               :: mdfrec !!  Standard rec. length in MD-file (300)
!
!
! Local variables
!
    integer                     :: idef   ! Default value for integer
    integer                     :: k
    integer                     :: lenc   ! Help var. (length of var. chulp to be looked for in the MD-file)
    integer                     :: lkw
    integer                     :: n
    integer                     :: nk
    integer                     :: nlook  ! Help var.: nr. of data to look for in the MD-file
    integer                     :: ntrec  ! Help. var to keep track of NRREC
    integer      , dimension(4) :: isom
    integer      , dimension(4) :: ival   ! Help array for reading
    logical                     :: defaul ! Flag set to YES if default value may be applied in case var. read is empty (ier <= 0, or nrread < nlook)
    logical                     :: found
    logical                     :: lerror ! Flag=TRUE if a local error is encountered
    logical                     :: newkw  ! Logical var. specifying whether a new recnam should be read from the MD-file or just new data in the continuation line
    logical                     :: nodef  ! Flag set to YES if default value may NOT be applied in case var. read is empty (ier <= 0, or nrread < nlook)
    character(20)               :: chdef  ! Default file name (blank)
    character(20)               :: chhlp  ! Help variable for reading
    character(3)                :: errtxt ! errormessage text
    character(6)                :: keyw   ! Name of record to look for in the MD-file (usually KEYWRD or RECNAM)
!
!
!! executable statements -------------------------------------------------------
!
    !
    !
    nprocs  => gdp%gdusrpar%nprocs
    nprdim  => gdp%gdusrpar%nprdim
    nread   => gdp%gdusrpar%nread
    nprcus  => gdp%gdusrpar%nprcus
    nprinp  => gdp%gdusrpar%nprinp
    icousr  => gdp%gdusrpar%icousr
    rcousr  => gdp%gdusrpar%rcousr
    filusr  => gdp%gdusrpar%filusr
    procs   => gdp%gdusrpar%procs
    prcusr  => gdp%gdusrpar%prcusr
    ccousr  => gdp%gdusrpar%ccousr
    itis    => gdp%gdrdpara%itis
    !
    do n = 1, nprocs
       nread(n) = 0
    enddo
    !
    do n = 1, mxusrp
       prcusr(n) = ' '
       filusr(n) = ' '
       rcousr(n) = 0.0
       icousr(n) = 0
       ccousr(n) = ' '
       do k = 1, 4
          nprinp(k, n) = 0
       enddo
    enddo
    !
    nprcus = 0
    !
    !-----initialize local parameters
    !
    lerror = .false.
    newkw = .true.
    defaul = .true.
    nodef = .not.defaul
    found = .true.
    idef = 0
    chdef = ' '
    !
    !-----Test value of nprocs
    !
    if (nprocs==0) goto 9999
    !
    !-----locate 'Nprcus' record for number of user defined input processes
    !
    keyw = 'Nprcus'
    ntrec = nrrec
    lkw = 6
    call search(lunmd     ,lerror    ,newkw     ,nrrec     ,found     , &
              & ntrec     ,mdfrec    ,itis      ,keyw      ,lkw       , &
              & 'NO'      )
    !
    !-----found ? then read NPRCUS
    !
    if (found) then
       nlook = 1
       call read2i(lunmd     ,lerror    ,keyw      ,newkw     ,nlook     , &
                 & mdfrec    ,ival      ,idef      ,defaul    ,nrrec     , &
                 & ntrec     ,lundia    ,gdp       )
       !
       !
       !--------reading error?
       !
       if (lerror) then
          lerror = .false.
          ival(1) = idef
       endif
       nprcus = ival(1)
       !
       !--------test NPRCUS inside dimension MXUSRP
       !
       if (nprcus>nprocs) then
          write (errtxt(:2), '(i2)') nprocs
          call prterr(lundia    ,'V211'    ,errtxt(:2))
          !
          error = .true.
          lerror = .true.
          goto 9999
       endif
       !
       !--------locate user defined processes 'Prcusr' record
       !        for name of user defined process
       !        Note: be aware that the name 'Prcusr' is case-sensitive !!
       !
       do k = 1, 4
          isom(k) = 0
       enddo
       !
       ntrec = nrrec
       do n = 1, nprcus
          keyw = 'Prcusr'
          lenc = 20
          nlook = 1
          call read2c(lunmd     ,lerror    ,keyw      ,newkw     ,nlook     , &
                    & mdfrec    ,chhlp     ,chdef     ,lenc      ,nrrec     , &
                    & ntrec     ,lundia    ,gdp       )
          !
          !
          !-----------reading error?
          !
          if (lerror) then
             lerror = .false.
             chhlp = chdef
          endif
          prcusr(n) = chhlp
          !
          !-----------test contents of PRCUSR(N) Should not be blank !!
          !
          if (prcusr(n)==chdef) then
             call prterr(lundia    ,'V212'    ,' '       )
             !
             error = .true.
             lerror = .true.
             exit
          endif
          !
          !-----------test contents of PRCUSR(N) should be one of the pre-
          !           defined processes
          !
          nk = 0
          do k = 1, nprocs
             if (prcusr(n)==procs(k)) then
                if (nread(k)/=0) then
                   call prterr(lundia    ,'V213'    ,prcusr(n) )
                   !
                   error = .true.
                   lerror = .true.
                   goto 9999
                endif
                nk = k
                exit
             endif
          enddo
          !
          if (nk==0) then
             call prterr(lundia    ,'V202'    ,prcusr(n) )
             !
             error = .true.
             lerror = .true.
             exit
          endif
          nread(nk) = n
          !
          !-----------locate user defined processes 'Prcusr' record
          !           for name of user defined process
          !
          keyw = 'Nprinp'
          nlook = 4
          call read2i(lunmd     ,lerror    ,keyw      ,newkw     ,nlook     , &
                    & mdfrec    ,ival      ,idef      ,nodef     ,nrrec     , &
                    & ntrec     ,lundia    ,gdp       )
          !
          !
          !-----------reading error?
          !
          if (lerror) then
             lerror = .false.
             error = .true.
             do k = 1, 4
                ival(k) = idef
             enddo
          endif
          !
          !-----------Define NPRINP and compare contents NPRINP with NPRDIM
          !           test inside dimension MXUSRP
          !
          do k = 1, 4
             nprinp(k, n) = ival(k)
             !
             if (nprinp(k, n)/=nprdim(k, nk)) then
                call prterr(lundia    ,'V214'    ,prcusr(n) )
                !
                error = .true.
                lerror = .true.
                goto 9999
             endif
             !
             isom(k) = isom(k) + ival(k)
             if (isom(k)>mxusrp) then
                write (errtxt(:2), '(i2)') mxusrp
                call prterr(lundia    ,'V215'    ,errtxt(:2))
                !
                error = .true.
                lerror = .true.
                goto 9999
             endif
          !
          enddo
       !
       enddo
    endif
    !
    !
    !-----Reset read input to avoid inconstant reading later on
    !
 9999 continue
    if (lerror) then
       call prterr(lundia    ,'V201'    ,' '       )
       !
       nprcus = 0
       do n = 1, mxusrp
          prcusr(n) = chdef
          do k = 1, 4
             nprinp(k, n) = idef
          enddo
       enddo
    endif
end subroutine rdusrp

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


! { dg-do run }
!
! Run-time test for EVENT_TYPE
!
use iso_fortran_env, only: event_type
implicit none

type(event_type), save :: var[*]
integer :: count, stat

count = -42
call event_query (var, count)
if (count /= 0) call abort()

stat = 99
event post (var, stat=stat)
if (stat /= 0) call abort()
call event_query(var, count, stat=stat)
if (count /= 1 .or. stat /= 0) call abort()

stat = 99
event post (var[this_image()])
call event_query(var, count)
if (count /= 2) call abort()

stat = 99
event wait (var)
call event_query(var, count)
if (count /= 1) call abort()

stat = 99
event post (var)
call event_query(var, count)
if (count /= 2) call abort()

stat = 99
event post (var)
call event_query(var, count)
if (count /= 3) call abort()

stat = 99
event wait (var, until_count=2)
call event_query(var, count)
if (count /= 1) call abort()

stat = 99
event wait (var, stat=stat, until_count=1)
if (stat /= 0) call abort()
call event_query(event=var, stat=stat, count=count)
if (count /= 0 .or. stat /= 0) call abort()
end


text:
      subroutine pol_massless(p,i,f)
      implicit none
      include 'constants.f'
      include 'swapxz.f'
c     massless vector polarization  subroutine
c     mcfm style momentum assignment
      integer i,pol
      double complex p(4),f(4)
      double precision p0,px,py,pz,pv,ct,st,cphi,sphi
      logical,save::first
      data first/.true./

      if (first) then
      write(6,*) 'pol_massless:swapxz=',swapxz
      first=.false.
      if ((i .ne. 1) .and. (i .ne. -1)) then
      write(6,*) 'pol_massless:pol out of bounds, i= ',i
      stop
      endif
      endif

      if (swapxz) then
      p0=dreal(p(4))
      px=+dreal(p(3))
      py=-dreal(p(2))
      pz=+dreal(p(1))
      else
      p0=dreal(p(4))
      px=dreal(p(1))
      py=dreal(p(2))
      pz=dreal(p(3))
      endif
      pv=dsqrt(dabs(p0**2))
      ct=pz/pv
      st=dsqrt(dabs(1d0-ct**2))

      if (st.lt.1d-8) then
        cphi=1d0
        sphi=0d0
      else
        cphi= px/pv/st
        sphi= py/pv/st
      endif

c     the following ifstatement distinguishes between
c     positive and negative energies
      if ( p0.gt.0d0) then
        pol=i
      else
        pol=-i
      endif

      if (swapxz) then
      f(4)=dcmplx(0d0,0d0)
      f(3)=+dcmplx(ct*cphi/rt2,+pol*sphi/rt2)
      f(2)=-dcmplx(ct*sphi/rt2,-pol*cphi/rt2)
      f(1)=+dcmplx(-st/rt2,0d0)
      else
      f(4)=dcmplx(0d0,0d0)
      f(1)=+dcmplx(ct*cphi/rt2,+pol*sphi/rt2)
      f(2)=+dcmplx(ct*sphi/rt2,-pol*cphi/rt2)
      f(3)=+dcmplx(-st/rt2,0d0)
      endif
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
c ---------------------------------------------------------------
      SUBROUTINE DHYDRO(P,TKV,ZSFC,NLVL,ZH,IER)
      IMPLICIT NONE

c NCL: zh = hydro (p,tkv,zsfc)

c use the hydrostatic eqn to get geopotential height
c .   it is assumed that p(1) and tkv(1) contain the
c .   surface quantities corresponding to zsfc.

c .   missing data not allowed
C input

C no. levels; error code
      INTEGER NLVL,IER
C pres (Pa)
C temp (K) at each "p"
C sfc geopotential height (gpm)
      DOUBLE PRECISION P(NLVL),TKV(NLVL),ZSFC
C output

C calculated geopotential (gpm)
      DOUBLE PRECISION ZH(NLVL)
C local

      INTEGER NL
      DOUBLE PRECISION G,RDAIR,RDAG,TVBAR

c     --------
      IER = 0
      IF (NLVL.LT.1) THEN
C unrealistic nlvl value
          IER = -1
          RETURN
      END IF

C gravity at 45 deg lat used by the WMO
      G = 9.80665D0
C gas const dry air (j/{kg-k})
      RDAIR = 287.04D0
      RDAG = RDAIR/G

c calculate geopotential height if initial z available [hydrostatic eq]

      ZH(1) = ZSFC
      DO NL = 2,NLVL
C same as [ln(p(nl)+ln(p(nl-1)]
          TVBAR = (TKV(NL)*DLOG(P(NL))+TKV(NL-1)*DLOG(P(NL-1)))/
     +            DLOG(P(NL)*P(NL-1))
          ZH(NL) = ZH(NL-1) + RDAG*TVBAR*DLOG(P(NL-1)/P(NL))
      END DO

      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      SUBROUTINE ks2d2s(x1,y1,n1,x2,y2,n2,d,prob)
      INTEGER n1,n2
      REAL d,prob,x1(n1),x2(n2),y1(n1),y2(n2)
CU    USES pearsn,probks,quadct
      INTEGER j
      REAL d1,d2,dum,dumm,fa,fb,fc,fd,ga,gb,gc,gd,r1,r2,rr,sqen,probks
      d1=0.0
      do 11 j=1,n1
        call quadct(x1(j),y1(j),x1,y1,n1,fa,fb,fc,fd)
        call quadct(x1(j),y1(j),x2,y2,n2,ga,gb,gc,gd)
        d1=max(d1,abs(fa-ga),abs(fb-gb),abs(fc-gc),abs(fd-gd))
11    continue
      d2=0.0
      do 12 j=1,n2
        call quadct(x2(j),y2(j),x1,y1,n1,fa,fb,fc,fd)
        call quadct(x2(j),y2(j),x2,y2,n2,ga,gb,gc,gd)
        d2=max(d2,abs(fa-ga),abs(fb-gb),abs(fc-gc),abs(fd-gd))
12    continue
      d=0.5*(d1+d2)
      sqen=sqrt(float(n1)*float(n2)/float(n1+n2))
      call pearsn(x1,y1,n1,r1,dum,dumm)
      call pearsn(x2,y2,n2,r2,dum,dumm)
      rr=sqrt(1.0-0.5*(r1**2+r2**2))
      prob=probks(d*sqen/(1.0+rr*(0.25-0.75/sqen)))
      return
      END

ext: for

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


text:
	SUBROUTINE GGETR  ( data, nw, iret )
C************************************************************************
C* GGETR								*
C*									*
C* This subroutine gets part or all of a message.  This subroutine	*
C* waits for a message to be received.					*
C*									*
C* NOTE: NW is an INPUT parameter and DATA is an OUTPUT parameter.	*
C* Most GEMPAK subroutines have the input parameters before the		*
C* output parameters.							*
C*									*
C* This subroutine was created to avoid conflict of return value type	*
C* on linux platforms							*
C*									*
C* GGETR ( DATA, NW, IRET )						*
C*									*
C* Output parameters:							*
C*	DATA (NW)	REAL		Message part			*
C*									*
C* Input parameters:							*
C*	NW		INTEGER		Length of message part		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* m.gamazaychikov/SAIC	01/04	CREATED					*
C************************************************************************
        INCLUDE         'ERROR.PRM'
C*
	REAL	data(*)
C------------------------------------------------------------------------
	iret = NORMAL
	CALL GGET(data, nw, iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


!123456
      real, allocatable :: eigup(:,:),eigdn(:,:),feigup(:,:)
      real, allocatable :: feigdn(:,:),xrk(:,:), xk(:,:), kval(:)
      integer :: nband, nks, nline, istart, nrem, check, line, next_line
      integer :: nion
      real :: k, kpt, dis_k
      real :: fermi,spin
      real :: b11, b12, b13, b21, b22, b23, b31, b32, b33
      real :: c11, c12, c13, c21, c22, c23, c31, c32, c33
      real :: a11, a12, a13, a21, a22, a23, a31, a32, a33
      real, allocatable :: ts(:,:,:), tp(:,:,:), td(:,:,:), tt(:,:,:)
      real, allocatable :: xs(:,:,:), xp(:,:,:), xd(:,:,:), xt(:,:,:)
      real, allocatable :: ys(:,:,:), yp(:,:,:), yd(:,:,:), yt(:,:,:)
      real, allocatable :: zs(:,:,:), zp(:,:,:), zd(:,:,:), zt(:,:,:)
      real, allocatable :: tts(:,:), ttp(:,:), ttd(:,:), ttt(:,:)
      real, allocatable :: txs(:,:), txp(:,:), txd(:,:), txt(:,:)
      real, allocatable :: tys(:,:), typ(:,:), tyd(:,:), tyt(:,:)
      real, allocatable :: tzs(:,:), tzp(:,:), tzd(:,:), tzt(:,:)
      character (len=20) :: tot
      real, allocatable :: xtot(:,:), ytot(:,:), ztot(:,:)
      real, allocatable :: fspinup(:,:),fspindown(:,:)
      real :: Pi
!
!   please remember the number of writting field in the format should be larger than the number of bands
!123456
      OPEN(1,file='band.dat')
      read(1,*) fermi
      write(6,*) "spin unpolar(type '0') or spin polar(type '1')"
      write(6,*) "fermi level"
      write(6,*) "b vector"

      open(17,file='CONTCAR')
      read(17,*)
      read(17,*)
      read(17,*) a11, a12, a13
      read(17,*) a21, a22, a23
      read(17,*) a31, a32, a33

      c11 = a22*a33 - a23*a32
      c12 = a23*a31 - a21*a33
      c13 = a21*a32 - a22*a31
      c21 = a32*a13 - a33*a12
      c22 = a33*a11 - a31*a13
      c23 = a31*a12 - a32*a11
      c31 = a12*a23 - a13*a22
      c32 = a13*a21 - a11*a23
      c33 = a11*a22 - a12*a21

      Pi = Atan(1.0)*4

      b11 = 2*Pi*c11/(a11*c11 + a12*c12 + a13*c13)
      b12 = 2*Pi*c12/(a11*c11 + a12*c12 + a13*c13)
      b13 = 2*Pi*c13/(a11*c11 + a12*c12 + a13*c13)
      b21 = 2*Pi*c21/(a21*c21 + a22*c22 + a23*c23)
      b22 = 2*Pi*c22/(a21*c21 + a22*c22 + a23*c23)
      b23 = 2*Pi*c23/(a21*c21 + a22*c22 + a23*c23)
      b31 = 2*Pi*c31/(a31*c31 + a32*c32 + a33*c33)
      b32 = 2*Pi*c32/(a31*c31 + a32*c32 + a33*c33)
      b33 = 2*Pi*c33/(a31*c31 + a32*c32 + a33*c33)

      open(15,file='EIGENVAL')
      read(15,*) nion, c2, c3, spin

       If(spin.eq.1) then
      do line=1,4
        read(15,*)
        next_line=next_line+line
      enddo
      read(15,*)Nelectron,nk,nband
      allocate (feigup(nk,nband),eigup(nk,nband))
      allocate (xrk(nk,3),xk(nk,3),kval(nk))

      do i =1,nk
         read(15,*)
         read(15,*)xrk(i,1),xrk(i,2),xrk(i,3),tmp123
         do j = 1,nband
            read(15,*)itmp,eigup(i,j)
         enddo
      enddo

      xk(:,1) = b11*xrk(:,1) + b21*xrk(:,2) + b31*xrk(:,3)
      xk(:,2) = b12*xrk(:,1) + b22*xrk(:,2) + b32*xrk(:,3)
      xk(:,3) = b13*xrk(:,1) + b23*xrk(:,2) + b33*xrk(:,3)

      k = 0.0
      kpt = 0.0
      kval(1) = 0.0
      do i=2,nk
         k=sqrt((xk(i,1)-xk(i-1,1))**2 + (xk(i,2)-xk(i-1,2))**2 + &
                 (xk(i,3)-xk(i-1,3))**2)
         kpt = kpt + k
         kval(i) = kpt
      enddo

      feigup(:,:) = eigup(:,:) - fermi

      open(30,file='PROCAR')
      allocate (ts(nk,nband,nion), tp(nk,nband,nion), td(nk,nband,nion))
      allocate (tt(nk,nband,nion), xs(nk,nband,nion), xp(nk,nband,nion))
      allocate (xd(nk,nband,nion), xt(nk,nband,nion), ys(nk,nband,nion))
      allocate (yp(nk,nband,nion), yd(nk,nband,nion), yt(nk,nband,nion))
      allocate (zs(nk,nband,nion), zp(nk,nband,nion), zd(nk,nband,nion))
      allocate (zt(nk,nband,nion))
      allocate (tts(nk,nband), ttp(nk,nband), ttd(nk,nband))
      allocate (ttt(nk,nband), txs(nk,nband), txp(nk,nband))
      allocate (txd(nk,nband), txt(nk,nband), tys(nk,nband))
      allocate (typ(nk,nband), tyd(nk,nband), tyt(nk,nband))
      allocate (tzs(nk,nband), tzp(nk,nband), tzd(nk,nband))
      allocate (tzt(nk,nband))

      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)

      do i = 1,nk-1
        do j = 1,nband
          do k = 1,nion
            read(30,*)n,ts(i,j,k),tp(i,j,k),td(i,j,k),tt(i,j,k)
          enddo

          read(30,*)tot,tts(i,j),ttp(i,j),ttd(i,j),ttt(i,j)

          do k = 1,nion
            read(30,*)n,xs(i,j,k),xp(i,j,k),xd(i,j,k),xt(i,j,k)
          enddo

          read(30,*)tot,txs(i,j),txp(i,j),txd(i,j),txt(i,j)

          do k = 1,nion
            read(30,*)n,ys(i,j,k),yp(i,j,k),yd(i,j,k),yt(i,j,k)
          enddo

          read(30,*)tot,tys(i,j),typ(i,j),tyd(i,j),tyt(i,j)

          do k = 1,nion
            read(30,*)n,zs(i,j,k),zp(i,j,k),zd(i,j,k),zt(i,j,k)

          enddo

          read(30,*)tot,tzs(i,j),tzp(i,j),tzd(i,j),tzt(i,j)

          read(30,*)
          read(30,*)
          read(30,*)
          read(30,*)
      enddo

          read(30,*)
          read(30,*)
          read(30,*)
      enddo

       i = nk
        do j = 1,nband-1
          do k = 1,nion
            read(30,*)n,ts(i,j,k),tp(i,j,k),td(i,j,k),tt(i,j,k)
          enddo

          read(30,*)tot,tts(i,j),ttp(i,j),ttd(i,j),ttt(i,j)

          do k = 1,nion
            read(30,*)n,xs(i,j,k),xp(i,j,k),xd(i,j,k),xt(i,j,k)
          enddo

          read(30,*)tot,txs(i,j),txp(i,j),txd(i,j),txt(i,j)

          do k = 1,nion
            read(30,*)n,ys(i,j,k),yp(i,j,k),yd(i,j,k),yt(i,j,k)
          enddo

          read(30,*)tot,tys(i,j),typ(i,j),tyd(i,j),tyt(i,j)

          do k = 1,nion
            read(30,*)n,zs(i,j,k),zp(i,j,k),zd(i,j,k),zt(i,j,k)

          enddo

          read(30,*)tot,tzs(i,j),tzp(i,j),tzd(i,j),tzt(i,j)

          read(30,*)
          read(30,*)
          read(30,*)
          read(30,*)
        enddo

      i = nk
      j = nband
          do k = 1,nion
            read(30,*)n,ts(i,j,k),tp(i,j,k),td(i,j,k),tt(i,j,k)
          enddo

          read(30,*)tot,tts(i,j),ttp(i,j),ttd(i,j),ttt(i,j)

          do k = 1,nion
            read(30,*)n,xs(i,j,k),xp(i,j,k),xd(i,j,k),xt(i,j,k)
          enddo

          read(30,*)tot,txs(i,j),txp(i,j),txd(i,j),txt(i,j)

          do k = 1,nion
            read(30,*)n,ys(i,j,k),yp(i,j,k),yd(i,j,k),yt(i,j,k)
          enddo

          read(30,*)tot,tys(i,j),typ(i,j),tyd(i,j),tyt(i,j)

          do k = 1,nion
            read(30,*)n,zs(i,j,k),zp(i,j,k),zd(i,j,k),zt(i,j,k)

          enddo

          read(30,*)tot,tzs(i,j),tzp(i,j),tzd(i,j),tzt(i,j)

      open(24,file='band.txt')
 118  format(f14.7,3x,f14.7)
        write(24,*) '#kval, total'
      do j = 1, nband
      do i = 1, nk
          write(24,118)kval(i), feigup(i,j)
      enddo
          write(24,*)
      enddo

        allocate(xtot(nk,nband),ytot(nk,nband),ztot(nk,nband))

        do i = 1, nk
          do j = 1, nband
        if(txt(i,j)**2.GT.tyt(i,j)**2.and.txt(i,j)**2.GT.tzt(i,j)**2)  then
              xtot(i,j) = txt(i,j)/ttt(i,j)
            endif

         if(tyt(i,j)**2.GT.tzt(i,j)**2.and.tyt(i,j)**2.GT.txt(i,j)**2) then
              ytot(i,j) = tyt(i,j)/ttt(i,j)
            endif

         if(tzt(i,j)**2.GT.txt(i,j)**2.and.tzt(i,j)**2.GT.tyt(i,j)**2) then
              ztot(i,j) = tzt(i,j)/ttt(i,j)
            endif
          enddo
        enddo

        allocate(fspinup(nk,nband), fspindown(nk,nband))

        do i = 1, nk
          do j = 1, nband
            if(xtot(i,j) .GT. 0) then
              fspinup(i,j) = feigup(i,j)
             endif

            if(xtot(i,j) .LT. 0 ) then
              fspindown(i,j) = feigup(i,j)
            endif

            if(ytot(i,j) .GT. 0) then
              fspinup(i,j) = feigup(i,j)
           endif

            if(ytot(i,j) .LT. 0) then
              fspindown(i,j) = feigup(i,j)
            endif

            if(ztot(i,j) .GT. 0) then
              fspinup(i,j) = feigup(i,j)
            endif

            if(ztot(i,j) .LT. 0) then
              fspindown(i,j) = feigup(i,j)
            endif

          enddo
        enddo

        do i = 1,nk
        do j = 1, nband
          if(fspinup(i,j).eq.0) then
        fspinup(i,j) = 1000000
          endif
          if(fspindown(i,j).eq.0) then
        fspindown(i,j) = 1000000
         endif
        enddo
        enddo

     open(40,file='band2.txt')
 140  format(f14.7,3x,f14.7,3x,f14.7,3x,f14.7)
        write(40,*) '#kval, total'
      do j = 1, nband
      do i = 1, nk
        write(40,140)kval(i), feigup(i,j), fspinup(i,j), fspindown(i,j)
      enddo
          write(40,*)
      enddo

       else

      read(15,*)
      read(15,*)
      read(15,*)
      read(15,*)
      read(15,*)Nelectron,nk,nband
      allocate (feigup(nk,nband),feigdn(nk,nband),eigup(nk,nband))
      allocate (eigdn(nk,nband),xk(nk,3),xrk(nk,3))

      do i =1,nk
         read(15,*)
         read(15,*)xrk(i,1),xrk(i,2),xrk(i,3),tmp123
         do j = 1,nband
            read(15,*)itmp,eigup(i,j),eigdn(i,j)
         enddo
      enddo

      xk(:,1) = b11*xrk(:,1) + b21*xrk(:,2) + b31*xrk(:,3)
      xk(:,2) = b12*xrk(:,1) + b22*xrk(:,2) + b32*xrk(:,3)
      xk(:,3) = b13*xrk(:,1) + b23*xrk(:,2) + b33*xrk(:,3)

      k = 0.0
      kpt = 0.0
      kval(1) = 0.0
      do i=2,nks
         k=sqrt((xk(i,1)-xk(i-1,1))**2 + (xk(i,2)-xk(i-1,2))**2 + &
                (xk(i,3)-xk(i-1,3))**2)
         kpt = kpt + k
         kval(i) = kpt
      enddo

      feigup(:,:) = eigup(:,:) - fermi
      feigdn(:,:) = eigdn(:,:) - fermi

      open(30,file='PROCAR')
      allocate (ts(nk,nband,nion), tp(nk,nband,nion), td(nk,nband,nion))
      allocate (tt(nk,nband,nion), xs(nk,nband,nion), xp(nk,nband,nion))
      allocate (xd(nk,nband,nion), xt(nk,nband,nion), ys(nk,nband,nion))
      allocate (yp(nk,nband,nion), yd(nk,nband,nion), yt(nk,nband,nion))
      allocate (zs(nk,nband,nion), zp(nk,nband,nion), zd(nk,nband,nion))
      allocate (zt(nk,nband,nion))
      allocate (tts(nk,nband), ttp(nk,nband), ttd(nk,nband))
      allocate (ttt(nk,nband), txs(nk,nband), txp(nk,nband))
      allocate (txd(nk,nband), txt(nk,nband), tys(nk,nband))
      allocate (typ(nk,nband), tyd(nk,nband), tyt(nk,nband))
      allocate (tzs(nk,nband), tzp(nk,nband), tzd(nk,nband))
      allocate (tzt(nk,nband))

      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)

      do i = 1,nk-1
        do j = 1,nband
          do k = 1,nion
            read(30,*)n,ts(i,j,k),tp(i,j,k),td(i,j,k),tt(i,j,k)
          enddo

          read(30,*)tot,tts(i,j),ttp(i,j),ttd(i,j),ttt(i,j)

          do k = 1,nion
            read(30,*)n,xs(i,j,k),xp(i,j,k),xd(i,j,k),xt(i,j,k)
          enddo

          read(30,*)tot,txs(i,j),txp(i,j),txd(i,j),txt(i,j)

          do k = 1,nion
            read(30,*)n,ys(i,j,k),yp(i,j,k),yd(i,j,k),yt(i,j,k)
          enddo

          read(30,*)tot,tys(i,j),typ(i,j),tyd(i,j),tyt(i,j)

          do k = 1,nion
            read(30,*)n,zs(i,j,k),zp(i,j,k),zd(i,j,k),zt(i,j,k)

          enddo

          read(30,*)tot,tzs(i,j),tzp(i,j),tzd(i,j),tzt(i,j)

          read(30,*)
          read(30,*)
          read(30,*)
          read(30,*)
      enddo

          read(30,*)
          read(30,*)
          read(30,*)
      enddo

       i = nk
        do j = 1,nband-1
          do k = 1,nion
            read(30,*)n,ts(i,j,k),tp(i,j,k),td(i,j,k),tt(i,j,k)
          enddo

          read(30,*)tot,tts(i,j),ttp(i,j),ttd(i,j),ttt(i,j)

          do k = 1,nion
            read(30,*)n,xs(i,j,k),xp(i,j,k),xd(i,j,k),xt(i,j,k)
          enddo

          read(30,*)tot,txs(i,j),txp(i,j),txd(i,j),txt(i,j)

          do k = 1,nion
            read(30,*)n,ys(i,j,k),yp(i,j,k),yd(i,j,k),yt(i,j,k)
          enddo

          read(30,*)tot,tys(i,j),typ(i,j),tyd(i,j),tyt(i,j)

          do k = 1,nion
            read(30,*)n,zs(i,j,k),zp(i,j,k),zd(i,j,k),zt(i,j,k)

          enddo

          read(30,*)tot,tzs(i,j),tzp(i,j),tzd(i,j),tzt(i,j)

          read(30,*)
          read(30,*)
          read(30,*)
          read(30,*)
        enddo

      i = nk
      j = nband
          do k = 1,nion
            read(30,*)n,ts(i,j,k),tp(i,j,k),td(i,j,k),tt(i,j,k)
          enddo

          read(30,*)tot,tts(i,j),ttp(i,j),ttd(i,j),ttt(i,j)

          do k = 1,nion
            read(30,*)n,xs(i,j,k),xp(i,j,k),xd(i,j,k),xt(i,j,k)
          enddo

          read(30,*)tot,txs(i,j),txp(i,j),txd(i,j),txt(i,j)

          do k = 1,nion
            read(30,*)n,ys(i,j,k),yp(i,j,k),yd(i,j,k),yt(i,j,k)
          enddo

          read(30,*)tot,tys(i,j),typ(i,j),tyd(i,j),tyt(i,j)

          do k = 1,nion
            read(30,*)n,zs(i,j,k),zp(i,j,k),zd(i,j,k),zt(i,j,k)

          enddo

          read(30,*)tot,tzs(i,j),tzp(i,j),tzd(i,j),tzt(i,j)

        write(24,*) '#kval, up, down'
      do j = 1, nband
      do i = 1, nk
          write(24,118)kval(i),feigup(i,j),feigdn(i,j)
      enddo
          write(24,*)
      enddo

      endif

      stop

      end


text:
C    Copyright(C) 2014-2017 National Technology & Engineering Solutions of
C    Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    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 NTESS 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

C $Id: rgdstr.f,v 1.2 1998/07/14 18:19:58 gdsjaar Exp $
C $Log: rgdstr.f,v $
C Revision 1.2  1998/07/14 18:19:58  gdsjaar
C Removed unused variables, cleaned up a little.
C
C Changed BLUE labels to GREEN to help visibility on black background
C (indirectly requested by a couple users)
C
C Revision 1.1.1.1  1990/11/30 11:14:57  gdsjaar
C FASTQ Version 2.0X
C
c Revision 1.1  90/11/30  11:14:56  gdsjaar
c Initial revision
c
C
CC* FILE: [.MAIN]RGDSTR.FOR
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE RGDSTR (NPNODE, NPELEM, KKK, NNXK, XN, YN, NXK)
C************************************************************************
C
C  SUBROUTINE RGDSTR = CALCULATES A REGION DISTORTION MEASURE
C
C***********************************************************************
C
      DIMENSION XN (NPNODE), YN (NPNODE), NXK (NNXK, NPELEM)
C
C  CALCULATE THE ELEMENT DISTORTION
C
      N1 = NXK (1,1)
      N2 = NXK (2,1)
      N3 = NXK (3,1)
      N4 = NXK (4,1)
      CALL DSTORT (XN (N1), XN (N2), XN (N3), XN (N4),
     &   YN (N1), YN (N2), YN (N3), YN (N4), VALUE)
      VMIN = VALUE
      VMAX = VALUE
      SUM = VALUE
      KMIN = 1
      KMAX = 1
      DO 100 I = 1, KKK
         N1 = NXK (1,I)
         N2 = NXK (2,I)
         N3 = NXK (3,I)
         N4 = NXK (4,I)
         CALL DSTORT (XN (N1), XN (N2), XN (N3), XN (N4),
     &      YN (N1), YN (N2), YN (N3), YN (N4), VALUE)
         IF (VMIN .GT. VALUE) THEN
            VMIN = VALUE
            KMIN = I
         ELSE IF (VMAX .LT. VALUE) THEN
            VMAX = VALUE
            KMAX = I
         ENDIF
         SUM = SUM + VALUE
  100 CONTINUE
C
C  PRINT OUT THE RESULTS
C
      SUM = SUM / FLOAT (KKK)
      WRITE (*, 10000) VMIN, KMIN, VMAX, KMAX, SUM
C
      RETURN
C
10000 FORMAT ('  THE MINIMUM DISTORTION IS: ',G14.7,' IN ELEMENT: ',I10,
     &   /,
     &   '  THE MAXIMUM DISTORTION IS: ',G14.7,' IN ELEMENT: ',I10, /,
     &   '  THE AVERAGE DISTORTION IS: ',G14.7)
C
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


! Copyright (c) 2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

module m1
  !ERROR: Logical constant '.true.' may not be used as a defined operator
  interface operator(.TRUE.)
  end interface
  !ERROR: Logical constant '.false.' may not be used as a defined operator
  generic :: operator(.false.) => bar
end

module m2
  interface operator(+)
    module procedure foo
  end interface
  interface operator(.foo.)
    module procedure foo
  end interface
  interface operator(.ge.)
    module procedure bar
  end interface
contains
  integer function foo(x, y)
    logical, intent(in) :: x, y
    foo = 0
  end
  logical function bar(x, y)
    complex, intent(in) :: x, y
    bar = .false.
  end
end

!ERROR: Intrinsic operator '.le.' may not be used as a defined operator
use m2, only: operator(.le.) => operator(.ge.)
!ERROR: Intrinsic operator '.not.' may not be used as a defined operator
use m2, only: operator(.not.) => operator(.foo.)
!ERROR: Logical constant '.true.' may not be used as a defined operator
use m2, only: operator(.true.) => operator(.foo.)
end


text:
/* This is a multi-line
   comment. The lexer discards it. */
T = Nat->Nat;
lambda f:T. lambda x:Nat. f (f x);

max_stars_repo_name: brunoflores/fulluntyped

max_issues_repo_name: brunoflores/fulluntyped

max_forks_repo_name: brunoflores/fulluntyped


!*==zhptrd.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b ZHPTRD
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download ZHPTRD + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          UPLO
!       INTEGER            INFO, N
!       ..
!       .. Array Arguments ..
!       DOUBLE PRECISION   D( * ), E( * )
!       COMPLEX*16         AP( * ), TAU( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> ZHPTRD reduces a complex Hermitian matrix A stored in packed form to
!> real symmetric tridiagonal form T by a unitary similarity
!> transformation: Q**H * A * Q = T.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] UPLO
!> \verbatim
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!> \endverbatim
!>
!> \param[in,out] AP
!> \verbatim
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>          On exit, if UPLO = 'U', the diagonal and first superdiagonal
!>          of A are overwritten by the corresponding elements of the
!>          tridiagonal matrix T, and the elements above the first
!>          superdiagonal, with the array TAU, represent the unitary
!>          matrix Q as a product of elementary reflectors; if UPLO
!>          = 'L', the diagonal and first subdiagonal of A are over-
!>          written by the corresponding elements of the tridiagonal
!>          matrix T, and the elements below the first subdiagonal, with
!>          the array TAU, represent the unitary matrix Q as a product
!>          of elementary reflectors. See Further Details.
!> \endverbatim
!>
!> \param[out] D
!> \verbatim
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The diagonal elements of the tridiagonal matrix T:
!>          D(i) = A(i,i).
!> \endverbatim
!>
!> \param[out] E
!> \verbatim
!>          E is DOUBLE PRECISION array, dimension (N-1)
!>          The off-diagonal elements of the tridiagonal matrix T:
!>          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
!> \endverbatim
!>
!> \param[out] TAU
!> \verbatim
!>          TAU is COMPLEX*16 array, dimension (N-1)
!>          The scalar factors of the elementary reflectors (see Further
!>          Details).
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup complex16OTHERcomputational
!
!> \par Further Details:
!  =====================
!>
!> \verbatim
!>
!>  If UPLO = 'U', the matrix Q is represented as a product of elementary
!>  reflectors
!>
!>     Q = H(n-1) . . . H(2) H(1).
!>
!>  Each H(i) has the form
!>
!>     H(i) = I - tau * v * v**H
!>
!>  where tau is a complex scalar, and v is a complex vector with
!>  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
!>  overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
!>
!>  If UPLO = 'L', the matrix Q is represented as a product of elementary
!>  reflectors
!>
!>     Q = H(1) H(2) . . . H(n-1).
!>
!>  Each H(i) has the form
!>
!>     H(i) = I - tau * v * v**H
!>
!>  where tau is a complex scalar, and v is a complex vector with
!>  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
!>  overwriting A(i+2:n,i), and tau is stored in TAU(i).
!> \endverbatim
!>
!  =====================================================================
      SUBROUTINE ZHPTRD(Uplo,N,Ap,D,E,Tau,Info)
      IMPLICIT NONE
!*--ZHPTRD155
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Uplo
      INTEGER Info , N
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION D(*) , E(*)
      COMPLEX*16 Ap(*) , Tau(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      COMPLEX*16 ONE , ZERO , HALF
      PARAMETER (ONE=(1.0D+0,0.0D+0),ZERO=(0.0D+0,0.0D+0),              &
     &           HALF=(0.5D+0,0.0D+0))
!     ..
!     .. Local Scalars ..
      LOGICAL upper
      INTEGER i , i1 , i1i1 , ii
      COMPLEX*16 alpha , taui
!     ..
!     .. External Subroutines ..
      EXTERNAL XERBLA , ZAXPY , ZHPMV , ZHPR2 , ZLARFG
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      COMPLEX*16 ZDOTC
      EXTERNAL LSAME , ZDOTC
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC DBLE
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters
!
      Info = 0
      upper = LSAME(Uplo,'U')
      IF ( .NOT.upper .AND. .NOT.LSAME(Uplo,'L') ) THEN
         Info = -1
      ELSEIF ( N<0 ) THEN
         Info = -2
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('ZHPTRD',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible
!
      IF ( N<=0 ) RETURN
!
      IF ( upper ) THEN
!
!        Reduce the upper triangle of A.
!        I1 is the index in AP of A(1,I+1).
!
         i1 = N*(N-1)/2 + 1
         Ap(i1+N-1) = DBLE(Ap(i1+N-1))
         DO i = N - 1 , 1 , -1
!
!           Generate elementary reflector H(i) = I - tau * v * v**H
!           to annihilate A(1:i-1,i+1)
!
            alpha = Ap(i1+i-1)
            CALL ZLARFG(i,alpha,Ap(i1),1,taui)
            E(i) = alpha
!
            IF ( taui/=ZERO ) THEN
!
!              Apply H(i) from both sides to A(1:i,1:i)
!
               Ap(i1+i-1) = ONE
!
!              Compute  y := tau * A * v  storing y in TAU(1:i)
!
               CALL ZHPMV(Uplo,i,taui,Ap,Ap(i1),1,ZERO,Tau,1)
!
!              Compute  w := y - 1/2 * tau * (y**H *v) * v
!
               alpha = -HALF*taui*ZDOTC(i,Tau,1,Ap(i1),1)
               CALL ZAXPY(i,alpha,Ap(i1),1,Tau,1)
!
!              Apply the transformation as a rank-2 update:
!                 A := A - v * w**H - w * v**H
!
               CALL ZHPR2(Uplo,i,-ONE,Ap(i1),1,Tau,1,Ap)
!
            ENDIF
            Ap(i1+i-1) = E(i)
            D(i+1) = Ap(i1+i)
            Tau(i) = taui
            i1 = i1 - i
         ENDDO
         D(1) = Ap(1)
      ELSE
!
!        Reduce the lower triangle of A. II is the index in AP of
!        A(i,i) and I1I1 is the index of A(i+1,i+1).
!
         ii = 1
         Ap(1) = DBLE(Ap(1))
         DO i = 1 , N - 1
            i1i1 = ii + N - i + 1
!
!           Generate elementary reflector H(i) = I - tau * v * v**H
!           to annihilate A(i+2:n,i)
!
            alpha = Ap(ii+1)
            CALL ZLARFG(N-i,alpha,Ap(ii+2),1,taui)
            E(i) = alpha
!
            IF ( taui/=ZERO ) THEN
!
!              Apply H(i) from both sides to A(i+1:n,i+1:n)
!
               Ap(ii+1) = ONE
!
!              Compute  y := tau * A * v  storing y in TAU(i:n-1)
!
               CALL ZHPMV(Uplo,N-i,taui,Ap(i1i1),Ap(ii+1),1,ZERO,Tau(i),&
     &                    1)
!
!              Compute  w := y - 1/2 * tau * (y**H *v) * v
!
               alpha = -HALF*taui*ZDOTC(N-i,Tau(i),1,Ap(ii+1),1)
               CALL ZAXPY(N-i,alpha,Ap(ii+1),1,Tau(i),1)
!
!              Apply the transformation as a rank-2 update:
!                 A := A - v * w**H - w * v**H
!
               CALL ZHPR2(Uplo,N-i,-ONE,Ap(ii+1),1,Tau(i),1,Ap(i1i1))
!
            ENDIF
            Ap(ii+1) = E(i)
            D(i) = Ap(ii)
            Tau(i) = taui
            ii = i1i1
         ENDDO
         D(N) = Ap(ii)
      ENDIF
!
!
!     End of ZHPTRD
!
      END SUBROUTINE ZHPTRD


text:
!
! NAME:
!       MW_TauProfile
!
! PURPOSE:
!       Program to general microwave TauProfile datafiles based
!       on oSRF inputs.
!
! SIDE EFFECTS:
!       If the output file already exists, it is overwritten.
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 03-Sep-2013
!                       paul.vandelst@noaa.gov
!

PROGRAM MW_TauProfile

  ! ------------------
  ! Environment set up
  ! ------------------
  ! Module usage
  USE Type_Kinds                , ONLY: fp
  USE Message_Handler           , ONLY: SUCCESS, FAILURE, WARNING, &
                                        Display_Message, Program_Message
  USE Profile_Utility           , ONLY: ID_H2O    , &
                                        PPMV_UNITS, &
                                        ND_UNITS  , &
                                        MR_UNITS  , &
                                        MD_UNITS  , &
                                        PP_UNITS  , &
                                        PPMV_to_PP, &
                                        ND_to_PP  , &
                                        MR_to_PP  , &
                                        MD_to_PP
  USE Timing_Utility            , ONLY: Timing_type , &
                                        Timing_Begin, &
                                        Timing_End  , &
                                        Timing_Display
  USE Spectral_Units_Conversion , ONLY: inverse_cm_to_GHz
  USE SensorInfo_Define
  USE SensorInfo_LinkedList
  USE SensorInfo_IO
  USE AtmProfile_Define
  USE AtmProfile_netCDF_IO
  USE TauProfile_Define
  USE TauProfile_netCDF_IO
  USE PtrArr_Define             , ONLY: PtrArr_type, &
                                        PtrArr_FromVector, &
                                        PtrArr_ToVector
  USE oSRF_Define               , ONLY: oSRF_type, &
                                        oSRF_Destroy , &
                                        oSRF_GetValue, &
                                        oSRF_Convolve
  USE oSRF_File_Define          , ONLY: oSRF_File_type, &
                                        oSRF_File_GetValue, &
                                        oSRF_File_GetFrom , &
                                        oSRF_File_Read    , &
                                        oSRF_File_Destroy
  USE MWLBL_Transmittance
  USE CRTM_GeometryInfo
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------
  CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'MW_TauProfile'
  CHARACTER(*), PARAMETER :: PROGRAM_VERSION_ID = &
  ! Keyword argument set value
  INTEGER, PARAMETER :: SET = 1
  ! direction names
  INTEGER, PARAMETER :: N_DIRECTIONS = 2
  CHARACTER(*), PARAMETER :: DIRECTION_NAME(N_DIRECTIONS) = [ 'upwelling  ', &
                                                              'downwelling'  ]

  ! Model names
  INTEGER, PARAMETER :: N_MODELS = 2
  INTEGER, PARAMETER :: MODEL_LIEBE      = 1
  INTEGER, PARAMETER :: MODEL_ROSENKRANZ = 2
  INTEGER, PARAMETER :: MODEL_INDEX(N_MODELS) = [ MODEL_LIEBE, &
                                                  MODEL_ROSENKRANZ ]
  CHARACTER(*), PARAMETER :: MODEL_NAME(N_MODELS) = [ 'Liebe89/93  ', &
                                                      'Rosenkranz03'  ]
  ! Define the secant of the zenith angles to be used
  INTEGER, PARAMETER :: N_ANGLES = 7
  REAL(fp), PARAMETER :: ANGLE_SECANT(N_ANGLES) = &
    [ 1.00_fp, 1.25_fp, 1.50_fp, &
      1.75_fp, 2.00_fp, 2.25_fp, &
      3.00_fp ]
  ! The molecular set IDs.
  !   1 == WLO; Wet lines only
  !  10 == ALL; All absorbers with continua
  !  12 == WET; Wet (and the wet continua)
  !  13 == DRY; Dry (and the dry continua)
  !  15 == WCO; Wet continua only
  ! 101 == EffWLO; Effective wet lines (wet/wco)
  ! 113 == EffDRY; Effective dry (all/wet)
  INTEGER, PARAMETER :: N_MOLECULE_SETS = 7
  INTEGER, PARAMETER :: MOLECULE_SET_LIST( N_MOLECULE_SETS ) = &
    [ 1, 10, 12, 13, 15, 101, 113 ]
  CHARACTER(*), PARAMETER :: MOLECULE_SET_NAME( N_MOLECULE_SETS ) = &
    [ 'wlo   ', 'all   ', 'wet   ', 'dry   ', 'wco   ', 'effwlo', 'effdry' ]
  INTEGER, PARAMETER :: WLO_IDX = 1
  INTEGER, PARAMETER :: ALL_IDX = 2
  INTEGER, PARAMETER :: WET_IDX = 3
  INTEGER, PARAMETER :: DRY_IDX = 4
  INTEGER, PARAMETER :: WCO_IDX = 5
  INTEGER, PARAMETER :: EFFECTIVE_WLO_IDX = 6
  INTEGER, PARAMETER :: EFFECTIVE_DRY_IDX = 7
  ! Numeric literal
  REAL(fp), PARAMETER :: ZERO = 0.0_fp
  REAL(fp), PARAMETER :: ONE  = 1.0_fp
  REAL(fp), PARAMETER :: TOLERANCE = EPSILON(ONE)

  ! ---------
  ! Variables
  ! ---------
  CHARACTER(256) :: msg
  INTEGER :: err_stat
  INTEGER :: alloc_stat
  INTEGER :: IO_Status
  INTEGER :: direction
  LOGICAL :: downwelling
  INTEGER :: Model
  LOGICAL :: rosenkranz
  CHARACTER(256)             :: sinfo_filename
  TYPE(SensorInfo_type)      :: sinfo
  TYPE(Sensorinfo_List_type) :: sinfo_list
  CHARACTER(256)        :: atmprofile_filename
  TYPE(AtmProfile_type), ALLOCATABLE :: atmprofile(:)
  CHARACTER(256)        :: tauprofile_filename
  TYPE(TauProfile_type) :: tauprofile
  CHARACTER(256)        :: osrf_filename
  TYPE(oSRF_File_type)  :: osrf_file
  TYPE(oSRF_type)       :: osrf
  CHARACTER(256)        :: osrf_history, osrf_comment
  CHARACTER(256)        :: Comment
  INTEGER :: i ! N_ANGLES
  INTEGER :: j ! N_MOLECULE_SETS
  INTEGER :: k, n_layers
  INTEGER :: l, n_channels
  INTEGER :: n_frequencies
  INTEGER :: m, n_profiles
  INTEGER :: n, n_sensors
  INTEGER :: n_bands
  INTEGER, ALLOCATABLE :: n_points(:)
  CHARACTER(256) :: profile_set = ''
  CHARACTER(512) :: mwlbl_version_id
  INTEGER :: j_idx(1)
  INTEGER :: index_h2o
  REAL(fp), ALLOCATABLE :: h2o_pressure(:)
  REAL(fp), ALLOCATABLE :: frequency(:), wavenumber(:)
  REAL(fp), ALLOCATABLE, TARGET  :: tauall(:,:,:)     ! L x K x I
  REAL(fp), ALLOCATABLE, TARGET  :: tauwlo(:,:,:)     ! L x K x I
  REAL(fp), ALLOCATABLE, TARGET  :: tauwco(:,:,:)     ! L x K x I
  REAL(fp), ALLOCATABLE, TARGET  :: tauwet(:,:,:)     ! L x K x I
  REAL(fp), ALLOCATABLE, TARGET  :: taudry(:,:,:)     ! L x K x I
  REAL(fp),              POINTER :: tau(:,:,:) => NULL()
  REAL(fp) :: integrated_srf
  REAL(fp) :: zenith_radian
  TYPE(Timing_type) :: sensor_timing, total_timing
  TYPE(PtrArr_type), ALLOCATABLE :: passband_frequency(:)
  TYPE(PtrArr_type), ALLOCATABLE :: tau_bands(:)

  ! Program header
  CALL Program_Message( PROGRAM_NAME, &
                        'Program to compute transmittance profiles for user-defined '//&
                        'microwave instruments.', &
                        '$Revision$' )

  ! Get the calculation direction
  WRITE( *,FMT='(/5x,"Select atmospheric path")' )
  DO i = 1, N_DIRECTIONS
    WRITE( *,FMT='(10x,i1,") ",a)' ) i, DIRECTION_NAME(i)
  END DO
  WRITE( *,FMT='(5x,"Enter choice: ")',ADVANCE='NO' )
  READ( *,FMT='(i5)',IOSTAT=IO_Status ) direction
  IF ( IO_Status /= 0 ) THEN
    msg = 'Invalid ATMOSPHERIC PATH identifier input.'
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF
  IF ( direction < 1 .OR. direction > N_DIRECTIONS ) THEN
    msg = 'Invalid ATMOSPERIC PATH identifier value.'
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  ENDIF
  downwelling = .FALSE.
  IF ( TRIM(DIRECTION_NAME(direction)) == 'downwelling' ) THEN
    downwelling = .TRUE.
  END IF
  ! ...The absorption model
  model      = 2       ! Index
  rosenkranz = .TRUE.  ! Flag

  ! Read the SensorInfo file
  WRITE( *,FMT='(/5x,"Enter the SensorInfo filename: ")',ADVANCE='NO' )
  READ( *,FMT='(a)' ) sinfo_filename
  sinfo_filename = ADJUSTL(sinfo_filename)
  err_stat = Read_SensorInfo( sinfo_filename, sinfo_list, Quiet = SET )
  IF ( err_stat /= SUCCESS ) THEN
    msg = 'Error reading SensorInfo file '//TRIM(sinfo_filename)
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF
  ! ...Count the number of sensors
  n_sensors = Count_SensorInfo_Nodes( sinfo_list )
  IF ( n_sensors < 1 ) THEN
    msg = 'SensorInfo list is empty.'
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF
  WRITE( *,'(7x,"Number of sensor entries read: ",i0)' ) n_sensors

  ! Read the AtmProfile file
  WRITE( *,FMT='(/5x,"Enter the netCDF AtmProfile filename: ")',ADVANCE='NO' )
  READ( *,'(a)' ) atmprofile_filename
  atmprofile_filename = ADJUSTL(atmprofile_filename)
  ! ...Inquire the file for dimensions
  err_stat = AtmProfile_netCDF_InquireFile( atmprofile_filename, n_Profiles = n_profiles )
  IF ( err_stat /= SUCCESS ) THEN
    msg = 'Error inquiring AtmProfile input file '//TRIM(atmprofile_filename)
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF
  ! ...Allocate the AtmProfile array
  ALLOCATE( atmprofile( n_profiles ), STAT=alloc_stat )
  IF ( alloc_stat /= 0 ) THEN
    msg = 'Error allocating AtmProfile array'
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF
  ! ...Read the data
  err_stat = AtmProfile_netCDF_ReadFile( &
    atmprofile, &
    atmprofile_filename, &
    Profile_Set = profile_set, &
    Quiet = .TRUE. )
  IF ( err_stat /= SUCCESS ) THEN
    msg = 'Error reading AtmProfile file '//TRIM(atmprofile_filename)
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF
  WRITE( *,'(7x,"Number of profiles read: ",i0)' ) n_profiles
  ! ...Assign layer dimension to local variable
  n_layers = atmprofile(1)%n_Layers

  ! Convert water vapour to partial pressure
  H2O_Convert_Loop: DO m = 1, n_profiles

    ! Determine the H2O index in the absorber array
    n = COUNT( atmprofile(m)%Absorber_ID == ID_H2O )
    IF ( n == 0 ) THEN
      WRITE( msg,'("No H2O in absorber set for AtmProfile #",i0)' ) m
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF
    IF ( n > 1 ) THEN
      WRITE( msg,'("More than one H2O identifier in absorber set for AtmProfile #",i0)' ) m
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF
    j_idx = PACK((/(j,j=1,atmprofile(m)%n_Absorbers)/), atmprofile(m)%Absorber_ID == ID_H2O )
    index_h2o = j_idx(1)

    ! Allocate the H2O partial pressure array for units conversion
    ALLOCATE( h2o_pressure( atmprofile(m)%n_Layers ), STAT=alloc_stat )
    IF ( alloc_stat /= 0 ) THEN
      WRITE( msg,'("Error allocating H2O partial pressure array for AtmProfile #",i0)' ) m
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF

    ! Convert the water vapor amounts
    SELECT CASE ( atmprofile(m)%Absorber_Units_ID( index_h2o ) )
      ! Convert from ppmv
      CASE ( PPMV_UNITS )
        CALL PPMV_to_PP( atmprofile(m)%Layer_Pressure, &
                         atmprofile(m)%Layer_Absorber(:,index_h2o), &
                         h2o_pressure )
      ! Convert from number density
      CASE ( ND_UNITS )
        CALL ND_to_PP( atmprofile(m)%Layer_Absorber(:,index_h2o), &
                       atmprofile(m)%Layer_Temperature, &
                       h2o_pressure )
      ! Convert from mixing ratio
      CASE ( MR_UNITS )
        CALL MR_to_PP( atmprofile(m)%Layer_Pressure, &
                       atmprofile(m)%Layer_Absorber(:,index_h2o), &
                       h2o_pressure )
      ! Convert from mass density
      CASE ( MD_UNITS )
        CALL MD_to_PP( atmprofile(m)%Layer_Absorber(:,index_h2o), &
                       atmprofile(m)%Layer_Temperature, &
                       h2o_pressure )
      ! Partial pressure input
      CASE ( PP_UNITS )
        h2o_pressure = atmprofile(m)%Layer_Absorber(:,index_h2o)
      ! Any other input
      CASE DEFAULT
        WRITE( msg,'("Unrecognised water vapour units for AtmProfile #",i0)' ) m
        CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END SELECT

    ! Check the result
    IF ( ANY(h2o_pressure < ZERO) ) THEN
      WRITE( msg,'("Error converting water vapor units to hPa for AtmProfile #",i0)' ) m
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF

    ! Save the partial pressure in the structure array
    atmprofile(m)%Layer_Absorber(:,index_h2o) = h2o_pressure
    ! ...Update the absorber units ID
    atmprofile(m)%Absorber_Units_ID(index_h2o) = PP_UNITS
    ! ...Deallocate the H2O partial pressure array
    DEALLOCATE( h2o_pressure, STAT=alloc_stat )
    IF ( alloc_stat /= 0 ) THEN
      WRITE( msg,'("Error deallocating H2O partial pressure array for AtmProfile #",i0)' ) m
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF

  END DO H2O_Convert_Loop

  ! Get the LBL module id
  CALL MWLBL_Transmittance_Version(mwlbl_version_id)

  ! Begin the total timer
  CALL Timing_Begin(total_timing)

  ! Begin the loop over sensors
  Sensor_Loop: DO n = 1, n_sensors

    ! Begin the sensor timer
    CALL Timing_Begin(sensor_timing)

    ! Get the current sensor
    err_stat = GetFrom_Sensorinfo_List( sinfo_list, n, sinfo )
    IF ( err_stat /= SUCCESS ) THEN
      err_stat = FAILURE
      WRITE( msg,'("Error retrieving SensorInfo data for sensor # ",i0)' ) n
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF

    ! Only process microwave sensors
    IF ( sinfo%Sensor_Type /= MICROWAVE_SENSOR ) CYCLE Sensor_Loop

    ! Output an info message
    WRITE( *,'(//2x,"Calculating ",a," transmittances for ",a,/)' ) &
             TRIM(DIRECTION_NAME(direction)), TRIM(sinfo%Sensor_id)

    ! Construct sensor filenames
    osrf_filename = TRIM(sinfo%Sensor_Id)//'.osrf.nc'

    ! Read the spectral response function data file
    err_stat = oSRF_File_Read( osrf_file, osrf_filename )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading oSRF data from '//TRIM(osrf_filename)
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF
    ! ...Retrieve the oSRF file attributes
    err_stat = oSRF_File_GetValue( &
      osrf_file                , &
      n_Channels = n_channels  , &
      History    = osrf_history, &
      Comment    = osrf_comment  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error retrieving attributes from'//TRIM(osrf_filename)
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF

    ! Allocate the TauProfile structure
    err_stat = Allocate_TauProfile( n_layers       , &
                                    n_channels     , &
                                    N_ANGLES       , &
                                    n_profiles     , &
                                    N_MOLECULE_SETS, &
                                    tauprofile       )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error allocating TauProfile structure for '//TRIM(sinfo%Sensor_Id)
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF

    ! Begin loop over channels
    Channel_Loop: DO l = 1, n_channels

      WRITE( *,'(/4x,"Processing channel #",i0,"....")' ) sinfo%Sensor_Channel(l)

      ! Get the current oSRF object from the file container
      err_stat = oSRF_File_GetFrom( &
        osrf_file, &
        osrf     , &
        pos = l )
      IF ( err_stat /= SUCCESS ) THEN
        WRITE( msg,'("Error retrieving oSRF #",i0," from ",a)') l, TRIM(osrf_filename)
        CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
      END IF

      ! Get the band and frequency dimensions, and the SRF integral
      err_stat = oSRF_GetValue( &
        osrf, &
        n_Bands            = n_bands, &
        n_Points           = n_points, &
        Integral           = integrated_srf, &
        Passband_Frequency = passband_frequency )
      IF ( err_stat /= SUCCESS ) THEN
        WRITE( msg,'("Error retrieving oSRF #",i0," attributes")') l
        CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
      END IF
      ! ...The total number of points for all bands
      n_frequencies = SUM(n_points)

      ! Allocate the LBL and PtrArr transmittance arrays
      ALLOCATE( tauall( n_frequencies, n_layers, N_ANGLES ), &
                tauwlo( n_frequencies, n_layers, N_ANGLES ), &
                tauwco( n_frequencies, n_layers, N_ANGLES ), &
                tauwet( n_frequencies, n_layers, N_ANGLES ), &
                taudry( n_frequencies, n_layers, N_ANGLES ), &
                tau_bands( n_bands ), &
                STAT=alloc_stat )
      IF ( alloc_stat /= 0 ) THEN
        WRITE( msg,'("Error allocating frequency and transmittance arrays for ",a," oSRF #",i0)') &
                   TRIM(sinfo%Sensor_Id), l
        CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
      END IF

      ! Transfer frequencies to contiguous vector
      err_stat = PtrArr_ToVector( passband_frequency, wavenumber )
      IF ( err_stat /= SUCCESS ) THEN
        WRITE( msg,'("Error transferring passband frequencies to array for oSRF #",i0)') l
        CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
      END IF
      frequency = inverse_cm_to_GHz(wavenumber)

      ! Begin profile loop
      Profile_Loop: DO m = 1, n_profiles

        WRITE( *,'(6x,"Processing profile #",i0,"....")' ) m
        n_layers = atmprofile(m)%n_Layers

        ! Call the LBL transmittance function
        err_stat = MWLBL_Compute_Tau( atmprofile(m)%Layer_Pressure,              &  ! Input
                                      atmprofile(m)%Layer_Temperature,           &  ! Input
                                      atmprofile(m)%Layer_Absorber(:,index_h2o), &  ! Input
                                      atmprofile(m)%Layer_Delta_Z,               &  ! Input
                                      ANGLE_SECANT,                              &  ! Input
                                      frequency,                                 &  ! Input
                                      tauall,                                    &  ! Output
                                      tauwlo,                                    &  ! Output
                                      tauwco,                                    &  ! Output
                                      tauwet,                                    &  ! Output
                                      taudry,                                    &  ! Output
                                      Downwelling=downwelling,                   &  ! Optional input
                                      Rosenkranz =rosenkranz,                    &  ! Optional input
                                      Quiet      =.TRUE.                         )  ! Optional input
        IF ( err_stat /= SUCCESS ) THEN
          WRITE( msg,'("Error computing MW LBL transmittance for profile #",i0,&
                      &", channel ",i0," of ",a)' ) &
                      m, sinfo%Sensor_Channel(l), TRIM(sinfo%Sensor_Id)
          CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
        END IF

        ! Average the data over the frequencies, but not
        ! for the effective transmittances
        DO j = 1, N_MOLECULE_SETS - 2  ! Last two are effective WLO and DRY

          ! Use a pointer to the particular transmittances
          SELECT CASE (j)
            CASE (WLO_IDX)
              tau => tauwlo
            CASE (ALL_IDX)
              tau => tauall
            CASE (WET_IDX)
              tau => tauwet
            CASE (DRY_IDX)
              tau => taudry
            CASE (WCO_IDX)
              tau => tauwco
            CASE DEFAULT
              WRITE( msg,'("Invalid molecule list index found at profile #",i0,&
                          &", channel ",i0," of ",a)' ) &
                          m, sinfo%Sensor_Channel(l), TRIM(sinfo%Sensor_Id)
              CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
          END SELECT

          ! Loop over the other transmittance dimensions.
          DO i = 1, N_ANGLES
            DO k = 1, n_layers

              ! Construct a PtrArr of the transmittance for convolution
              err_stat = PtrArr_FromVector( tau_bands, n_points, tau(:,k,i) )
              IF ( err_stat /= SUCCESS ) THEN
                WRITE( msg,'("Error extracting transmittance vector for layer# ",i0,&
                            &", angle# ",i0,", and molecule set ",a)' ) &
                            k, i, MOLECULE_SET_NAME(j)
                CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
              END IF

              ! Convolve the LBL transmittance with the SRF data
              err_stat = oSRF_Convolve( osrf, tau_bands, tauprofile%Tau(k,l,i,m,j) )
              IF ( err_stat /= SUCCESS ) THEN
                WRITE( msg,'("Error convolving transmittance for layer# ",i0,&
                            &", angle# ",i0,", and molecule set ",a)' ) &
                            k, i, MOLECULE_SET_NAME(j)
                CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
              END IF

            END DO  ! k : 1, n_layers
          END DO  ! i : 1, N_ANGLES
        END DO  ! j : 1, N_MOLECULE_SETS-2

        ! Nullify the transmittance array pointer
        NULLIFY( Tau )

        ! Compute the EFFECTIVE transmittances.
        DO i = 1, N_ANGLES
          DO k = 1, n_layers
            ! The EFFECTIVE WLO transmittances
            IF( tauprofile%Tau( k,l,i,m,WCO_IDX ) > TOLERANCE ) THEN
              tauprofile%Tau( k,l,i,m,EFFECTIVE_WLO_IDX ) = tauprofile%Tau( k,l,i,m,WET_IDX ) / &
                                                            tauprofile%Tau( k,l,i,m,WCO_IDX )
            ELSE
              tauprofile%Tau( k,l,i,m,EFFECTIVE_WLO_IDX ) = -ONE
            END IF

            ! The EFFECTIVE DRY transmittances
            IF( tauprofile%Tau( k,l,i,m,WET_IDX ) > TOLERANCE ) THEN
              tauprofile%Tau( k,l,i,m,EFFECTIVE_DRY_IDX ) = tauprofile%Tau( k,l,i,m,ALL_IDX ) / &
                                                            tauprofile%Tau( k,l,i,m,WET_IDX )

            ELSE
              tauprofile%Tau( k,l,i,m,EFFECTIVE_DRY_IDX ) = -ONE
            END IF

          END DO  ! k : 1, n_layers
        END DO  ! i : 1, N_ANGLES

        ! Compute the geometric angle profiles
        DO i = 1, N_ANGLES
          zenith_radian = ACOS(ONE/ANGLE_SECANT(i))
          CALL SACONV( zenith_radian, &
                       atmprofile(m)%Level_Altitude(1:atmprofile(m)%n_Layers), &
                       tauprofile%Geometric_Angle(:,i,m) )
        END DO  ! i : 1, N_ANGLES

      END DO Profile_Loop

      ! Deallocate the LBL arrays
      DEALLOCATE( wavenumber, frequency, &
                  tauall, tauwlo, tauwco, tauwet, taudry, &
                  tau_bands, &
                  STAT=alloc_stat )
      IF ( alloc_stat /= 0 ) THEN
        WRITE( msg,'("Error deallocating frequency and transmittance arrays for ",a," oSRF #",i0)') &
                   TRIM(sinfo%Sensor_Id), l
        CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
      END IF

    END DO Channel_Loop

    ! Output the TauProfile data
    WRITE( Comment,'("Number of points used per channel in LBL calculation: ",i0,&
                    &". Absorption model used: ",a)' ) &
                    n_frequencies, MODEL_NAME(model)
    ! ...Create the output data file
    tauprofile_filename = TRIM(DIRECTION_NAME(direction))//'.'//&
                          TRIM(sinfo%Sensor_Id)//&
                          '.TauProfile.nc'
    err_stat = Create_TauProfile_netCDF( tauprofile_filename, &
                                         atmprofile(1)%Level_Pressure, &
                                         sinfo%Sensor_Channel, &
                                         ANGLE_SECANT, &
                                         (/(m, m=1,n_Profiles)/), &
                                         MOLECULE_SET_LIST, &
                                         Release         =tauprofile%Release, &
                                         Version         =tauprofile%Version, &
                                         Sensor_ID       =sinfo%Sensor_ID, &
                                         WMO_Satellite_ID=sinfo%WMO_Satellite_ID, &
                                         WMO_Sensor_ID   =sinfo%WMO_Sensor_ID, &
                                         ID_Tag = TRIM(profile_set), &
                                         Title = TRIM(sinfo%Sensor_Name)//' '//&
                                                 TRIM(DIRECTION_NAME(direction))//&
                                                 ' transmittances for '//&
                                                 TRIM(sinfo%Satellite_Name), &
                                         History = PROGRAM_VERSION_ID//'; '//&
                                                   TRIM(mwlbl_version_id),  &
                                         Comment = TRIM(Comment) )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error creating netCDF TauProfile file '//TRIM(tauprofile_filename)
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF
    ! ...Write the data to file
    err_stat = Write_TauProfile_netCDF( tauprofile_filename, tauprofile, profile_angle=1 )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error writing TauProfile structure for '//&
            TRIM(sinfo%Sensor_Id)//' to '//TRIM(tauprofile_filename)
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF

    ! Display the sensor timing results
    CALL Timing_End(sensor_timing)
    WRITE( *,'(7x,"Timing for ",a,":")') TRIM(sinfo%Sensor_id)
    CALL Timing_Display(sensor_timing)

    ! Clean up sensor dependent data structures
    ! ...TauProfile
    err_stat = Destroy_TauProfile( tauprofile )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error destroying TauProfile structure for '//TRIM(sinfo%Sensor_Id)
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF
    ! ...oSRF
    CALL oSRF_File_Destroy( osrf_file )
    ! ...SensorInfo
    err_stat = Destroy_SensorInfo( sinfo )
    IF ( err_stat /= SUCCESS ) THEN
      WRITE( msg,'("Error destroying SensorInfo structure for sensor #",i0)' ) n
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF

  END DO Sensor_Loop

  ! Display the total timing results
  CALL Timing_End(total_timing)
  WRITE( *,'(/5x,"Total time for all sensors:")')
  CALL Timing_Display(total_timing)

  ! Clean up sensor independent data structures
  ! ...The AtmProfile structure array
  DEALLOCATE( atmprofile, STAT=alloc_stat )
  IF ( alloc_stat /= 0 ) THEN
    msg = 'Error deallocating AtmProfile structure array'
    CALL Display_Message( PROGRAM_NAME, msg, WARNING )
  END IF
  ! ...The linked list
  err_stat = Destroy_Sensorinfo_List( sinfo_list, Quiet=1 )
  IF ( err_stat /= SUCCESS ) THEN
    msg = 'Error destroying SensorInfo linked list.'
    CALL Display_Message( PROGRAM_NAME, msg, WARNING )
  END IF

END PROGRAM MW_TauProfile

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
! Extends the original JSON module for usability
module json_xtnsn_mod

    use json_mod
    implicit none

    logical :: json_found

    interface json_xtnsn_get
        module procedure :: json_xtnsn_value_get_real, json_xtnsn_file_get_real
        module procedure :: json_xtnsn_value_get_integer, json_xtnsn_file_get_integer
        module procedure :: json_xtnsn_value_get_string, json_xtnsn_file_get_string
        module procedure :: json_xtnsn_value_get_logical, json_xtnsn_file_get_logical
    end interface json_xtnsn_get

contains

    subroutine json_xtnsn_value_get_real(json, name, value, default_value)
        implicit none
        type(json_value),intent(in),pointer :: json
        character(len=*), intent(in) :: name
        real, intent(out) :: value
        real, intent(in), optional :: default_value

        call json_get(json, name, value, json_found)
        if(json_failed() .or. (.not. json_found)) then
            if (present(default_value)) then
                value = default_value
                call json_clear_exceptions()
            else
                write(*,*) 'Error: Unable to read required value: ',name
                STOP
            end if
        end if
    end subroutine json_xtnsn_value_get_real

    subroutine json_xtnsn_value_get_integer(json, name, value, default_value)
        implicit none
        type(json_value),intent(in),pointer :: json
        character(len=*), intent(in) :: name
        integer, intent(out) :: value
        integer, intent(in), optional :: default_value

        call json_get(json, name, value, json_found)
        if((.not.json_found) .or. json_failed()) then
            if (present(default_value)) then
                value = default_value
                call json_clear_exceptions()
            else
                write(*,*) 'Error: Unable to read required value: ', name
                STOP
            end if
        end if

    end subroutine json_xtnsn_value_get_integer

    subroutine json_xtnsn_value_get_string(json, name, value, default_value)
        implicit none
        type(json_value), intent(in), pointer :: json
        character(len=*), intent(in) :: name
        character(:), allocatable, intent(out) :: value
        character(len=*), intent(in), optional :: default_value

        call json_get(json, name, value, json_found)
        if((.not.json_found) .or. json_failed()) then
            if (present(default_value)) then
                value = default_value
                call json_clear_exceptions()
            else
                write(*,*) 'Error: Unable to read required value: ', name
                STOP
            end if
        end if

    end subroutine json_xtnsn_value_get_string

    subroutine json_xtnsn_value_get_logical(json, name, value, default_value)
        implicit none
        type(json_value), intent(in), pointer :: json
        character(len=*), intent(in) :: name
        logical, intent(out) :: value
        logical, intent(in), optional :: default_value

        call json_get(json, name, value, json_found)
        if((.not.json_found) .or. json_failed()) then
            if (present(default_value)) then
                value = default_value
                call json_clear_exceptions()
            else
                write(*,*) 'Error: Unable to read required value: ', name
                STOP
            end if
        end if

    end subroutine json_xtnsn_value_get_logical

    subroutine json_xtnsn_file_get_real(json, name, value, default_value)
        implicit none
        type(json_file) :: json
        character(len=*), intent(in) :: name
        real, intent(out) :: value
        real, intent(in), optional :: default_value

        call json%get(name, value)
        if(json_failed()) then
            if (present(default_value)) then
                value = default_value
                call json_clear_exceptions()
            else
                write(*,*) 'Error: Unable to read required value: ', trim(name)
                STOP
            end if
        end if

    end subroutine json_xtnsn_file_get_real

    subroutine json_xtnsn_file_get_integer(json, name, value, default_value)
        implicit none
        type(json_file) :: json
        character(len=*), intent(in) :: name
        integer, intent(out) :: value
        integer, intent(in), optional :: default_value

        call json%get(name, value)
        if(json_failed()) then
            if (present(default_value)) then
                value = default_value
                call json_clear_exceptions()
            else
                write(*,*) 'Error: Unable to read required value: ',name
                STOP
            end if
        end if
    end subroutine json_xtnsn_file_get_integer

    subroutine json_xtnsn_file_get_string(json, name, value)
        implicit none
        type(json_file) :: json
        character(len=*), intent(in) :: name
        character(:), allocatable, intent(out) :: value

        call json%get(name, value)
        if(json_failed()) then
            write(*,*) 'Error: Unable to read required value: ',name
            STOP
        end if

        value = trim(value)
    end subroutine json_xtnsn_file_get_string

    subroutine json_xtnsn_file_get_logical(json, name, value)
        implicit none
        type(json_file) :: json
        character(len=*), intent(in) :: name
        logical, intent(out) :: value

        call json%get(name, value)
        if(json_failed()) then
            write(*,*) 'Error: Unable to read required value: ',name
            STOP
        end if

    end subroutine json_xtnsn_file_get_logical

    subroutine json_check()
        if(json_failed()) then
            call print_json_error_message()
            STOP
        end if
    end subroutine json_check

    subroutine print_json_error_message()
        implicit none
        character(len=:),allocatable :: error_msg
        logical :: status_ok

        !get error message:
        call json_check_for_errors(status_ok, error_msg)

        !print it if there is one:
        if (.not. status_ok) then
            write(*,'(A)') error_msg
            deallocate(error_msg)
            call json_clear_exceptions()
        end if

    end subroutine print_json_error_message

end module json_xtnsn_mod

max_stars_repo_name: usuaero/MFTran

max_issues_repo_name: usuaero/MFTran

max_forks_repo_name: usuaero/MFTran


text:
! Ensure argument -flarge-sizes works as expected.
! TODO: Add checks when actual codegen is possible.

!--------------------------
! FLANG DRIVER (flang-new)
!--------------------------
! RUN: rm -rf %t/dir-flang-new  && mkdir -p %t/dir-flang-new && %flang -fsyntax-only -module-dir %t/dir-flang-new %s  2>&1
! RUN: cat %t/dir-flang-new/m.mod | FileCheck %s --check-prefix=NOLARGE
! RUN: rm -rf %t/dir-flang-new  && mkdir -p %t/dir-flang-new && %flang -fsyntax-only -flarge-sizes -module-dir %t/dir-flang-new %s  2>&1
! RUN: cat %t/dir-flang-new/m.mod | FileCheck %s --check-prefix=LARGE

!-----------------------------------------
! FRONTEND FLANG DRIVER (flang-new -fc1)
!-----------------------------------------
! RUN: rm -rf %t/dir-flang-new  && mkdir -p %t/dir-flang-new && %flang_fc1 -fsyntax-only -module-dir %t/dir-flang-new %s  2>&1
! RUN: cat %t/dir-flang-new/m.mod | FileCheck %s --check-prefix=NOLARGE
! RUN: rm -rf %t/dir-flang-new  && mkdir -p %t/dir-flang-new && %flang_fc1 -fsyntax-only -flarge-sizes -module-dir %t/dir-flang-new %s  2>&1
! RUN: cat %t/dir-flang-new/m.mod | FileCheck %s --check-prefix=LARGE

!-----------------------------------------
! EXPECTED OUTPUT WITHOUT -flarge-sizes
!-----------------------------------------
! NOLARGE: real(4)::z(1_8:10_8)
! NOLARGE-NEXT: integer(4),parameter::size_kind=4_4

!-----------------------------------------
! EXPECTED OUTPUT FOR -flarge-sizes
!-----------------------------------------
! LARGE: real(4)::z(1_8:10_8)
! LARGE-NEXT: integer(4),parameter::size_kind=8_4

module m
  implicit none
  real :: z(10)
  integer, parameter :: size_kind = kind(ubound(z, 1)) !-flarge-sizes
end

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm