text:
        SUBROUTINE  LS_NDDF ( lsfrpt, ipt, iret )
C************************************************************************
C* LS_NDDF                                                              *
C*                                                                      *
C* This subroutine decodes the Nddff group.  This group contains the    *
C* total cloud cover and wind direction and speed parameters.           *
C*                                                                      *
C* LS_NDDF  ( LSFRPT, IPT, IRET )                                       *
C*                                                                      *
C* Input parameters:                                                    *
C*      LSFRPT         	CHAR*           Report array                    *
C*					                                *
C* Input and Output parameters:                                         *
C*	IPT		INTEGER 	On input, points to start of    *
C*					Nddff group; on output, to end  *
C*					                                *
C* Output parameters:                                                   *
C*      RIVALS(IRCFRT) 	REAL            Total cloud cover in oktas      *
C*      RIVALS(IRDRCT)  REAL            Wind direction in degrees       *
C*      RIVALS(IRSPED) 	REAL            Wind speed in meters/second     *
C*      IRET            INTEGER         Return code                     *
C*                                        0 = normal return             *
C*                                        1 = problems                  *
C*                                                                      *
C**								        *
C* Log:							                *
C* R. Hollern/NCEP      6/96                                            *
C* R. Hollern/NCEP      8/96    Set wind direction to missing, if code  *
C*                              figure is 99                            *
C* R. Hollern/NCEP      1/98    Changes based on MA_NDDF                *
C* A. Hardy/GSC         1/98    Reordered calling sequence              *
C* R. Hollern/NCEP      2/98    Removed check that Nddff group exists   *
C************************************************************************
        INCLUDE		'GEMPRM.PRM'
        INCLUDE  	'lscmn.cmn'
C
        CHARACTER*(*)    lsfrpt
C
        CHARACTER        fld2*2, fld1*1
C------------------------------------------------------------------------
        iret   = 0
        wdrct  = RMISSD
	ival   = IMISSD
C
C*      There is at most one space before start of group.
C
        IF ( lsfrpt ( ipt:ipt ) .eq. ' ' ) ipt = ipt + 1
C
C*      Get total cloud cover.  WMO Table 2700 value.
C
        IF ( lsfrpt (ipt:ipt) .ne. '/' ) THEN
            fld1 = lsfrpt (ipt:ipt)
            CALL  ST_INTG ( fld1, ival, ier )
            IF ( ier .eq. 0 ) rivals ( ircfrt ) = FLOAT ( ival )
        END IF
C
        ipt = ipt + 1
C
C*      Get wind direction in degrees.
C
        IF ( lsfrpt ( ipt:ipt ) .ne. '/' .and.
     +       lsfrpt ( ipt:ipt+1 ) .ne. '/' ) THEN
            fld2 = lsfrpt ( ipt:ipt+1 )
            CALL  ST_INTG ( fld2, ival, ier )
            IF ( ier .eq. 0 ) THEN
                IF ( ival .ge. 0 .and. ival .le. 36 ) THEN
                    wdrct = 10. * FLOAT ( ival )
                  ELSE IF ( ival .eq. 99 ) THEN
C
C*		    Direction is variable, or all directions (WMO Code
C*		    Table 0877).
C
		    wdrct = 99.
                END IF
            END IF
        END IF
C
C*      Get wind speed in meters/second.
C
        ipt    = ipt + 2
        iparam = 0
        CALL  LS_WSPD  ( lsfrpt, iparam, ipt, iret )
C
C*	If direction is variable but speed is 0, set direction missing.
C
	IF ( ival .eq. 99 .and. rivals ( irsped ) .eq. 0.0 )
     +       wdrct = RMISSD
C
	rivals ( irdrct ) = wdrct
C*
        RETURN
        END


text:
*----------------------------------------------------------------------*
!>      multiplies two me_lists
!
!!      @param[in] me1 -
!!      @param[in] me2 -
!!      @param[inout] buf1,buf2 buffer for the
!!      @param nincore without function
!!      @param lbuf lenght of the buffers
*----------------------------------------------------------------------*
      function me_ddot(me1,me2, buf1, buf2, nincore, lbuf)
*----------------------------------------------------------------------*
      implicit none
      include 'stdunit.h'
      include 'mdef_me_list.h'
      include 'def_file_array.h'
      include 'ioparam.h'

      integer, parameter::
     &     ntest = 00
      character(len=*),parameter::
     &     i_am="me_ddot"

      real(8)::
     &     me_ddot

      type(me_list), intent(in)::
     &     me1,me2
      real(8),intent(inout)::
     &     buf1(*),buf2(*)

      integer,intent(in)::
     &     lbuf,nincore

      integer::
     &     llist, irec1, irec2

      real(8),external ::
     &     da_ddot

      if (ntest.ge.100) then
         call write_title(lulog,wst_dbg_subr,i_am)
      end if

      if (me1%len_op.ne.me2%len_op)
     &     call quit(0,i_am,"list have to be of identical length")
      llist=me1%len_op

      irec1=me1%fhand%current_record
      irec2=me2%fhand%current_record
      me_ddot=da_ddot(me1%fhand, irec1, me2%fhand, irec2,
     &     llist,buf1,buf2,lbuf)
      end function

ext: f

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
	SUBROUTINE G2T_CHECK ( ktype, kwave, kwind, wok, sok, iret )
C************************************************************************
C* G2T_CHECK								*
C*									*
C* This subroutine checks if the range of wind/wave heights lies within	*
C* guideline.								*
C*									*
C* G2T_CHECK ( KTYPE, KWAVE, KWIND, WOK, SOK, IRET )			*
C*									*
C* Input parameters:							*
C*	KTYPE		INTEGER		Grid parameter type		*
C*					 0 = wave and wind		*
C*					 1 = wave			*
C*					 2 = wind			*
C*	KWAVE (2)	INTEGER		Wave min/max			*
C*	KWIND (2)	INTEGER		Wind min/max			*
C*									*
C* Output Parameters:							*
C*	WOK		LOGICAL		Wave range flag			*
C*	SOK		LOGICAL		Wind range flag			*
C*	IRET		INTEGER		Return code			*
C*					 0 = normal return		*
C**									*
C* Log:									*
C* T. Lee/SAIC		11/06						*
C* T. Lee/SAIC		11/06	Handled single values			*
C* T. Lee/SAIC		07/07	Added NCNTER check			*
C* T. Lee/SAIC		11/07	Table-driven range values		*
C************************************************************************
	INCLUDE		'goftxt.cmn'
	INTEGER		kwave (2), kwind (2)
	LOGICAL		wok, sok, contin
C-----------------------------------------------------------------------
	iret = 0
C
	wok = .true.
	sok = .true.
	kdwave = ABS ( kwave ( 2 ) - kwave ( 1 ) )
	kdwind = ABS ( kwind ( 2 ) - kwind ( 1 ) )
C
C*	Wave heights.
C
	IF  ( ktype .le. 1 )  THEN
	    kmax  = AMAX0 ( kwave ( 1 ), kwave ( 2 ) )
	    IF ( ncnter .eq.  1 )  THEN
		ki = 1
		contin = .true.
		DO WHILE ( contin )
		    IF ( kmax .ge. mnwgap ( ki, 1 ) .and.
     +			 kmax .lt. mnwgap ( ki, 2 ) )  THEN
			IF ( kdwave .gt. mnwgap (ki, 3) )  wok = .false.
			ki = ki + 1
		      ELSE
			ki = ki + 1
		    END IF
		    IF ( ki .gt. llwgap ) contin = .false.
		END DO
	      ELSE IF ( ncnter .eq. 2 )  THEN
		kmin = AMIN0 ( kwave ( 1 ), kwave ( 2 ) )
		kdif = kmax - kmin
		IF  ( kdif .gt. NINT ( kmin / 2. ) ) wok = .false.
	    END IF
	END IF
C
C*	Wind speed.
C
	IF ( ktype .eq. 0 .or. ktype .eq. 2 )  THEN
	    kmax  = AMAX0 ( kwind ( 1 ), kwind ( 2 ) )
	    kj = 1
	    contin = .true.
	    DO WHILE ( contin )
		IF ( kmax .gt. mnsgap ( kj, 1 ) .and.
     +		     kmax .le. mnsgap ( kj, 2 ) )  THEN
		    IF ( kdwind .gt. mnsgap ( kj, 3 ) )  sok = .false.
		    kj = kj + 1
		  ELSE
		    kj = kj + 1
		END IF
		IF ( kj .gt. llsgap ) contin = .false.
	    END DO
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
module m1

interface
	module subroutine  add(x, y, r)
		real, intent(in) :: x
		real, intent(in) :: y
		real, intent(out) :: r
	end subroutine
end interface

end module

max_stars_repo_name: nicolasessisbreton/fython

max_issues_repo_name: nicolasessisbreton/fython

max_forks_repo_name: nicolasessisbreton/fython


text:
* Example with array reference in the condition of test
      PROGRAM TEST
      PARAMETER (N=10)
      INTEGER T(N),X

      PRINT *, 'THIS IS A TEST'

      IF (T(5).LE.10) THEN
         T(X) = 10
      ELSE
         T(5) = 10
      ENDIF

      PRINT *,T
      END


! ::: -----------------------------------------------------------
! ::: This routine will tag high error cells based on overdensity
! :::
! ::: INPUTS/OUTPUTS:
! :::
! ::: tag      <=  integer tag array
! ::: lo,hi     => index extent of tag array
! ::: set       => integer value to tag cell for refinement
! ::: clear     => integer value to untag cell
! ::: den       => density array
! ::: nc        => number of components in density array
! ::: domlo,hi  => index extent of problem domain
! ::: delta     => cell spacing
! ::: xlo       => physical location of lower left hand
! :::              corner of tag array
! ::: problo    => phys loc of lower left corner of prob domain
! ::: time      => problem evolution time
! ::: level     => refinement level of this array
! ::: -----------------------------------------------------------
      subroutine tag_overdensity(tag,tagl1,tagl2,tagl3,tagh1,tagh2,tagh3, &
                                 set,clear, &
                                 den,denl1,denl2,denl3,denh1,denh2,denh3, &
                                 lo,hi,nc,domlo,domhi,delta,level,avg_den)

      use amrex_fort_module, only : rt => amrex_real
      use probdata_module
      implicit none

      integer  :: set, clear, nc, level
      integer  :: tagl1,tagl2,tagl3,tagh1,tagh2,tagh3
      integer  :: denl1,denl2,denl3,denh1,denh2,denh3
      integer  :: lo(3), hi(3), domlo(3), domhi(3)
      integer  :: tag(tagl1:tagh1,tagl2:tagh2,tagl3:tagh3)
      real(rt) :: den(denl1:denh1,denl2:denh2,denl3:denh3,nc)
      real(rt) :: delta(3), avg_den

      integer  :: i,j,k
      real(rt) :: m_cell, V_cell

      V_cell = delta(1)*delta(2)*delta(3)
      do k = lo(3), hi(3)
         do j = lo(2), hi(2)
            do i = lo(1), hi(1)
               m_cell = den(i,j,k,1)*V_cell
               if ( m_cell .gt. 3.5d9 ) then
                  tag(i,j,k) = set
               endif
            enddo
         enddo
      enddo

      end subroutine tag_overdensity


text:
!
! Fundamental_Constants
!
! Module containing various fundamental mathematical and physical constants.
!
! The fundamental constants and equations used are taken from the NIST
! Reference on Constants, Units, and Uncertainty website:
!
!   http://physics.nist.gov/cuu/Constants/
!
! See also:
!
!   Mohr, P.J. and B.N. Taylor, "CODATA recommended values of the
!     fundamental physical constants: 1998", Reviews of Modern Physics,
!     Vol.72, No.2, 2000.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 02-May-2000
!                       paul.vandelst@noaa.gov
!

MODULE Fundamental_Constants

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Modules used
  USE Type_Kinds, ONLY: fp
  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------------
  ! Default visibility
  ! ------------------
  PRIVATE
  ! Procedures

  ! -----------------
  ! Module parameters
  ! -----------------
  ! Numeric literals
  REAL(fp), PARAMETER :: ZERO = 0.0_fp
  REAL(fp), PARAMETER :: ONE  = 1.0_fp
  REAL(fp), PARAMETER :: TWO  = 2.0_fp

  !#----------------------------------------------------------------------------#
  !#                -- IRRATIONAL NUMBERS AND ASSOCIATED BITS --                #
  !#----------------------------------------------------------------------------#

  ! PI
  REAL(fp), PARAMETER, PUBLIC :: PI             = 3.141592653589793238462643_fp
  REAL(fp), PARAMETER, PUBLIC :: PI_RECIPROCAL  = 0.318309886183790671537767_fp
  REAL(fp), PARAMETER, PUBLIC :: PI_SQUARED     = 9.869604401089358618834491_fp
  REAL(fp), PARAMETER, PUBLIC :: PI_SQUARE_ROOT = 1.772453850905516027298167_fp
  REAL(fp), PARAMETER, PUBLIC :: PI_LN          = 1.144729885849400174143427_fp
  REAL(fp), PARAMETER, PUBLIC :: PI_LOG10       = 0.497149872694133854351268_fp

  ! E
  REAL(fp), PARAMETER, PUBLIC :: E              = 2.718281828459045235360287_fp
  REAL(fp), PARAMETER, PUBLIC :: E_RECIPROCAL   = 0.367879441171442321595523_fp
  REAL(fp), PARAMETER, PUBLIC :: E_SQUARED      = 7.389056098930650227230427_fp
  REAL(fp), PARAMETER, PUBLIC :: E_LOG10        = 0.434294481903251827651129_fp

  ! Other transcendentals
  REAL(fp), PARAMETER, PUBLIC :: LN2            = 0.693147180559945309417232_fp

  !#----------------------------------------------------------------------------#
  !#                            -- UNIVERAL CONSTANTS --                        #
  !#----------------------------------------------------------------------------#

  ! ----------------------------------------------
  ! Speed of light
  ! Symbol:c,  Units:m/s,  Rel.Uncert.(ppm): exact
  ! ----------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: SPEED_OF_LIGHT = 2.99792458e+08_fp

  ! --------------------------------------------------
  ! Permeability of vacuum
  ! Symbol:mu0,  Units:N/A^2,  Rel.Uncert.(ppm): exact
  ! --------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: PERMEABILITY = PI * 4.0e-07_fp

  ! -----------------------------------------------------
  ! Permittivity of vacuum
  ! Symbol:epsilon0,  Units:F/m,  Rel.Uncert.(ppm): exact
  ! -----------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: PERMITTIVITY =                ONE                  / &
  !                                             ------------------------------------
                                                ( PERMEABILITY * SPEED_OF_LIGHT**2 )

  ! ---------------------------------------------
  ! Planck constant
  ! Symbol:h,  Units:Js,  Rel.Uncert.(ppm): 0.078
  ! ---------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: PLANCK_CONSTANT = 6.62606876e-34_fp

  ! ----------------------------------------------------
  ! Gravitational constant
  ! Symbol:G,  Units:m^3/kg/s^2,  Rel.Uncert.(ppm): 1500
  ! ----------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: GRAVITATIONAL_CONSTANT = 6.673e-11_fp

  !#----------------------------------------------------------------------------#
  !#                          -- CONVERSION FACTORS --                          #
  !#----------------------------------------------------------------------------#

  ! ---------------------------------------------
  ! Electron volt
  ! Symbol:eV,  Units:J,  Rel.Uncert.(ppm): 0.039
  ! ---------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: ELECTRON_VOLT = 1.602176462e-19_fp

  ! ---------------------------------------------
  ! Unified atomic mass unit
  ! Symbol:u,  Units:kg,  Rel.Uncert.(ppm): 0.079
  ! ---------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: UNIFIED_ATOMIC_MASS_UNIT = 1.66053873e-27_fp

  ! ----------------------------------------------
  ! Standard atmosphere
  ! Symbol:P0,  Units:Pa,  Rel.Uncert.(ppm): exact
  ! ----------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: STANDARD_ATMOSPHERE = 101325.0_fp

  ! ----------------------------------------------------------------------
  ! Standard temperature
  ! Symbol:T0,  Units:Kelvin,  Rel.Uncert.(ppm): exact
  !
  ! Note that the unit of thermodynamic temperature, the Kelvin, is the
  ! fraction 1/273.16 of the thermodynamic temperature of the triple point
  ! of water. The standard temperature is the ice point of water, NOT the
  ! triple point, hence the 0.01K difference.
  ! ----------------------------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: STANDARD_TEMPERATURE = 273.15_fp

  ! ------------------------------------------------
  ! Standard gravity
  ! Symbol:g,  Units:m/s^2,  Rel.Uncert.(ppm): exact
  ! ------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: STANDARD_GRAVITY = 9.80665_fp

  !#----------------------------------------------------------------------------#
  !#                        -- PHYSICOCHEMICAL CONSTANTS --                     #
  !#----------------------------------------------------------------------------#

  ! -----------------------------------------------------
  ! Avogadro constant
  ! Symbol:N(A),  Units:mole^-1,  Rel.Uncert.(ppm): 0.079
  ! -----------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: AVOGADRO_CONSTANT = 6.02214199e+23_fp

  ! -------------------------------------------------
  ! Molar gas constant
  ! Symbol:R,  Units:J/mole/K,  Rel.Uncert.(ppm): 1.7
  ! -------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: MOLAR_GAS_CONSTANT = 8.314472_fp

  ! --------------------------------------------
  ! Boltzmann constant
  ! Symbol:k,  Units:J/K,  Rel.Uncert.(ppm): 1.7
  !
  !         R
  !   k = ------
  !        N(A)
  !
  !     = 1.3806503(24)e-23
  !
  ! --------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: BOLTZMANN_CONSTANT = MOLAR_GAS_CONSTANT / &
  !                                                   ------------------
                                                       AVOGADRO_CONSTANT

  ! ------------------------------------------------------
  ! Stefan-Boltzmann constant
  ! Symbol:sigma,  Units:W/m^2/K^4,  Rel.Uncert.(ppm): 7.0
  !
  !             PI^2
  !             ----.k^4
  !              60                     h
  !   sigma = ------------   ( hbar = ----- )
  !            hbar^3.c^2              2PI
  !
  !         = 5.670400(40)e-08
  !
  ! I just placed the value here due to the mathematical
  ! gymnastics required to calculate it directly.
  ! ------------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: STEFAN_BOLTZMANN_CONSTANT = 5.670400e-08_fp

  ! -------------------------------------------------------
  ! First Planck function constant
  ! Symbol:c1,  Units:W.m^2.sr^-1,  Rel.Uncert.(ppm): 0.078
  !
  !   c1 = 2.h.c^2
  !
  !      = 1.191042722(93)e-16
  !
  ! -------------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: C_1 = TWO * PLANCK_CONSTANT * SPEED_OF_LIGHT**2

  ! ---------------------------------------------
  ! Second Planck function constant
  ! Symbol:c2,  Units:K.m,  Rel.Uncert.(ppm): 1.7
  !
  !         h.c
  !   c2 = -----
  !          k
  !
  !      = 1.4387752(25)e-02
  !
  ! ---------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: C_2 = PLANCK_CONSTANT * SPEED_OF_LIGHT / &
  !                                    ----------------------------------
                                               BOLTZMANN_CONSTANT

  ! -----------------------------------------------------------------
  ! Molar volume of an ideal gas at standard temperature and pressure
  ! Symbol:Vm,  Units:m^3/mol,  Rel.Uncert.(ppm): 1.7
  !
  !         R.T0
  !   Vm = ------
  !          P0
  !
  !      = 2.2413996(39)e-02
  !
  ! -----------------------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: STP_MOLAR_VOLUME = ( MOLAR_GAS_CONSTANT * STANDARD_TEMPERATURE ) / &
  !                                                 ---------------------------------------------
                                                                 STANDARD_ATMOSPHERE

  ! ------------------------------------------------------------------
  ! Loschmidt constant: The number density of one mole of an ideal gas
  ! at standard temperature and pressure
  ! Symbol:n0,  Units:m^-3,  Rel.Uncert.(ppm): 1.7
  !
  !         N(A).P0
  !   n0 = ---------
  !          R.T0
  !
  !         N(A)
  !      = ------     .....(1)
  !          Vm
  !
  !      = 2.6867775(47)e+25
  !
  ! Alternatively, using the ideal gas law directly, we know,
  !
  !   P.V = n.k.T     .....(2)
  !
  ! For V = 1m^3 (unit volume), and P = P0, T = T0, then eqn.(2)
  ! becomes,
  !
  !   P0 = n0.k.T0
  !
  ! which rearranges to
  !
  !          P0
  !   n0 = ------     .....(3)
  !         k.T0
  !
  ! Equation (1) rather than eqn(3) is used here.
  ! ------------------------------------------------------------------
  REAL(fp), PARAMETER, PUBLIC :: LOSCHMIDT_CONSTANT = AVOGADRO_CONSTANT / &
  !                                                   -----------------
                                                      STP_MOLAR_VOLUME
END MODULE Fundamental_Constants

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
!RUN: %flang %s -g -S -emit-llvm -o - | FileCheck %s

!Ensure that for an allocatable variable, we're taking the type of
!allocatable variable as DW_TAG_pointer_type.
!CHECK: call void @llvm.dbg.declare(metadata double** %{{.*}}, metadata ![[DILocalVariable:[0-9]+]], metadata !DIExpression())
!CHECK: ![[DILocalVariable]] = !DILocalVariable(name: "alcvar"
!CHECK-SAME: type: ![[PTRTYPE:[0-9]+]]
!CHECK: ![[PTRTYPE]] = !DIDerivedType(tag: DW_TAG_pointer_type, baseType: ![[TYPE:[0-9]+]]
!CHECK: ![[TYPE]] = !DIBasicType(name: "double precision",{{.*}}

program main
  real(kind=8), allocatable :: alcvar
  allocate(alcvar)
  alcvar = 7.7
  print *, alcvar
end program main

max_stars_repo_name: vigbalu/flang

max_issues_repo_name: vigbalu/flang

max_forks_repo_name: vigbalu/flang


text:
subroutine poisson_fort_header_ascii(matname1,matname2, time,mfields,maux,ngrids)
    implicit none

    integer iframe,mfields,ngrids, maux

    character*11 matname1
    character*11 matname2
    double precision time
    integer matunit1, matunit2
    integer mfieldsp2

    matunit1 = 10
    matunit2 = 15

    open(unit=matunit2,file=matname2)

    mfieldsp2 = mfields + 2  !! include soln, error
    write(matunit2,1000) time,mfieldsp2,ngrids,maux,2
 1000 format(e30.20,'    time', /,         &
           i5,'                 mfields'/,    &
           i5,'                 ngrids'/,  &
           i5,'                 num_aux'/, &
           i5,'                 num_dim')

    close(matunit2)

    open(unit=matunit1,file=matname1,status='replace')
    close(matunit1)

end subroutine poisson_fort_header_ascii

subroutine poisson_fort_output_ascii(matname1, &
         mx,my,mfields,mbc, xlower,ylower, dx,dy,  &
         rhs,soln,error,patch_num,level,blockno,mpirank)

    implicit none

    character(len=11) matname1
    integer mfields,mbc,mx,my
    integer patch_num
    integer level, blockno, mpirank
    double precision xlower, ylower,dx,dy

    double precision rhs(1-mbc:mx+mbc,1-mbc:my+mbc,mfields)
    double precision error(1-mbc:mx+mbc,1-mbc:my+mbc,mfields)
    double precision soln(1-mbc:mx+mbc,1-mbc:my+mbc,mfields)

    integer matunit1
    integer i,j,mq

    matunit1 = 10
    open(matunit1,file=matname1,position='append');

    call fclaw2d_clawpatch46_fort_write_grid_header(matunit1, &
           mx,my,xlower,ylower, dx,dy,patch_num,level, &
           blockno,mpirank)

    if (mfields .gt. 5) then
        write(6,'(A,A,A,I5,A)')     &
             'Warning (fclaw2d_fort_write_grid_header.f) ',  &
             ': meqn > 5; change format statement 120.',     &
             '(meqn = ',mfields,')'
        stop
    endif

    do j = 1,my
        do i = 1,mx
            do mq = 1,mfields
                if (abs(rhs(i,j,mq)) .lt. 1d-99) then
                    rhs(i,j,mq) = 0.d0
                elseif (abs(rhs(i,j,mq)) .gt. 1d99) then
                    rhs(i,j,mq) = 1d99
                endif
            end do
            if (abs(error(i,j,1)) .lt. 1d-99) then
                error(i,j,1) = 0.d0
            elseif (abs(error(i,j,1)) .gt. 1d99) then
                error(i,j,1) = 1d99
            endif
            if (abs(soln(i,j,1)) .lt. 1d-99) then
                soln(i,j,1) = 0.d0
            elseif (abs(soln(i,j,1)) .gt. 1d99) then
                soln(i,j,1) = 1d99
            endif
            write(matunit1,120) (rhs(i,j,mq),mq=1,mfields), soln(i,j,1), error(i,j,1)
        end do
        write(matunit1,*) ' '
    end do
120 format (5E26.16)

    close(matunit1)

end subroutine poisson_fort_output_ascii

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
!     path:      $Source: /storm/rc1/cvsroot/rc/rrtmg_sw/src/rrtmg_sw_setcoef.f90,v $
!     author:    $Author: mike $
!     revision:  $Revision: 1.2 $
!     created:   $Date: 2007/08/23 20:40:14 $

      module rrtmg_sw_setcoef

!  --------------------------------------------------------------------------
! |                                                                          |
! |  Copyright 2002-2007, Atmospheric & Environmental Research, Inc. (AER).  |
! |  This software may be used, copied, or redistributed as long as it is    |
! |  not sold and this copyright notice is reproduced on each copy made.     |
! |  This model is provided as is without any express or implied warranties. |
! |                       (http://www.rtweb.aer.com/)                        |
! |                                                                          |
!  --------------------------------------------------------------------------

! ------- Modules -------

      use shr_kind_mod, only: r8 => shr_kind_r8

!      use parkind, only : jpim, jprb
      use parrrsw, only : mxmol
      use rrsw_ref, only : pref, preflog, tref
      use rrsw_vsn, only : hvrset, hnamset

      implicit none

      contains

!----------------------------------------------------------------------------
      subroutine setcoef_sw(nlayers, pavel, tavel, pz, tz, tbound, coldry, wkl, &
                            laytrop, layswtch, laylow, jp, jt, jt1, &
                            co2mult, colch4, colco2, colh2o, colmol, coln2o, &
                            colo2, colo3, fac00, fac01, fac10, fac11, &
                            selffac, selffrac, indself, forfac, forfrac, indfor)
!----------------------------------------------------------------------------
!
! Purpose:  For a given atmosphere, calculate the indices and
! fractions related to the pressure and temperature interpolations.

! Modifications:
! Original: J. Delamere, AER, Inc. (version 2.5, 02/04/01)
! Revised: Rewritten and adapted to ECMWF F90, JJMorcrette 030224
! Revised: For uniform rrtmg formatting, MJIacono, Jul 2006

! ------ Declarations -------

! ----- Input -----
      integer, intent(in) :: nlayers         ! total number of layers

      real(kind=r8), intent(in) :: pavel(:)           ! layer pressures (mb)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(in) :: tavel(:)           ! layer temperatures (K)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(in) :: pz(0:)             ! level (interface) pressures (hPa, mb)
                                                        !    Dimensions: (0:nlayers)
      real(kind=r8), intent(in) :: tz(0:)             ! level (interface) temperatures (K)
                                                        !    Dimensions: (0:nlayers)
      real(kind=r8), intent(in) :: tbound             ! surface temperature (K)
      real(kind=r8), intent(in) :: coldry(:)          ! dry air column density (mol/cm2)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(in) :: wkl(:,:)           ! molecular amounts (mol/cm-2)
                                                        !    Dimensions: (mxmol,nlayers)

! ----- Output -----
      integer, intent(out) :: laytrop        ! tropopause layer index
      integer, intent(out) :: layswtch       !
      integer, intent(out) :: laylow         !

      integer, intent(out) :: jp(:)          !
                                                        !    Dimensions: (nlayers)
      integer, intent(out) :: jt(:)          !
                                                        !    Dimensions: (nlayers)
      integer, intent(out) :: jt1(:)         !
                                                        !    Dimensions: (nlayers)

      real(kind=r8), intent(out) :: colh2o(:)         ! column amount (h2o)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: colco2(:)         ! column amount (co2)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: colo3(:)          ! column amount (o3)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: coln2o(:)         ! column amount (n2o)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: colch4(:)         ! column amount (ch4)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: colo2(:)          ! column amount (o2)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: colmol(:)         !
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: co2mult(:)        !
                                                        !    Dimensions: (nlayers)

      integer, intent(out) :: indself(:)
                                                        !    Dimensions: (nlayers)
      integer, intent(out) :: indfor(:)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: selffac(:)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: selffrac(:)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: forfac(:)
                                                        !    Dimensions: (nlayers)
      real(kind=r8), intent(out) :: forfrac(:)
                                                        !    Dimensions: (nlayers)

      real(kind=r8), intent(out) :: &                 !
                         fac00(:), fac01(:), &          !    Dimensions: (nlayers)
                         fac10(:), fac11(:)

! ----- Local -----

      integer :: indbound
      integer :: indlev0
      integer :: lay
      integer :: jp1

      real(kind=r8) :: stpfac
      real(kind=r8) :: tbndfrac
      real(kind=r8) :: t0frac
      real(kind=r8) :: plog
      real(kind=r8) :: fp
      real(kind=r8) :: ft
      real(kind=r8) :: ft1
      real(kind=r8) :: water
      real(kind=r8) :: scalefac
      real(kind=r8) :: factor
      real(kind=r8) :: co2reg
      real(kind=r8) :: compfp

! Initializations
      stpfac = 296._r8/1013._r8

      indbound = tbound - 159._r8
      tbndfrac = tbound - int(tbound)
      indlev0  = tz(0) - 159._r8
      t0frac   = tz(0) - int(tz(0))

      laytrop  = 0
      layswtch = 0
      laylow   = 0

! Begin layer loop
      do lay = 1, nlayers
! Find the two reference pressures on either side of the
! layer pressure.  Store them in JP and JP1.  Store in FP the
! fraction of the difference (in ln(pressure)) between these
! two values that the layer pressure lies.

         plog = log(pavel(lay))
         jp(lay) = int(36._r8 - 5*(plog+0.04_r8))
         if (jp(lay) .lt. 1) then
            jp(lay) = 1
         elseif (jp(lay) .gt. 58) then
            jp(lay) = 58
         endif
         jp1 = jp(lay) + 1
         fp = 5._r8 * (preflog(jp(lay)) - plog)

! Determine, for each reference pressure (JP and JP1), which
! reference temperature (these are different for each
! reference pressure) is nearest the layer temperature but does
! not exceed it.  Store these indices in JT and JT1, resp.
! Store in FT (resp. FT1) the fraction of the way between JT
! (JT1) and the next highest reference temperature that the
! layer temperature falls.

         jt(lay) = int(3._r8 + (tavel(lay)-tref(jp(lay)))/15._r8)
         if (jt(lay) .lt. 1) then
            jt(lay) = 1
         elseif (jt(lay) .gt. 4) then
            jt(lay) = 4
         endif
         ft = ((tavel(lay)-tref(jp(lay)))/15._r8) - float(jt(lay)-3)
         jt1(lay) = int(3._r8 + (tavel(lay)-tref(jp1))/15._r8)
         if (jt1(lay) .lt. 1) then
            jt1(lay) = 1
         elseif (jt1(lay) .gt. 4) then
            jt1(lay) = 4
         endif
         ft1 = ((tavel(lay)-tref(jp1))/15._r8) - float(jt1(lay)-3)

         water = wkl(1,lay)/coldry(lay)
         scalefac = pavel(lay) * stpfac / tavel(lay)

! If the pressure is less than ~100mb, perform a different
! set of species interpolations.

         if (plog .le. 4.56_r8) go to 5300
         laytrop =  laytrop + 1
         if (plog .ge. 6.62_r8) laylow = laylow + 1

! Set up factors needed to separately include the water vapor
! foreign-continuum in the calculation of absorption coefficient.

         forfac(lay) = scalefac / (1.+water)
         factor = (332.0_r8-tavel(lay))/36.0_r8
         indfor(lay) = min(2, max(1, int(factor)))
         forfrac(lay) = factor - float(indfor(lay))

! Set up factors needed to separately include the water vapor
! self-continuum in the calculation of absorption coefficient.

         selffac(lay) = water * forfac(lay)
         factor = (tavel(lay)-188.0_r8)/7.2_r8
         indself(lay) = min(9, max(1, int(factor)-7))
         selffrac(lay) = factor - float(indself(lay) + 7)

! Calculate needed column amounts.

         colh2o(lay) = 1.e-20_r8 * wkl(1,lay)
         colco2(lay) = 1.e-20_r8 * wkl(2,lay)
         colo3(lay) = 1.e-20_r8 * wkl(3,lay)
!           colo3(lay) = 0._r8
!           colo3(lay) = colo3(lay)/1.16_r8
         coln2o(lay) = 1.e-20_r8 * wkl(4,lay)
         colch4(lay) = 1.e-20_r8 * wkl(6,lay)
         colo2(lay) = 1.e-20_r8 * wkl(7,lay)
         colmol(lay) = 1.e-20_r8 * coldry(lay) + colh2o(lay)
!           colco2(lay) = 0._r8
!           colo3(lay) = 0._r8
!           coln2o(lay) = 0._r8
!           colch4(lay) = 0._r8
!           colo2(lay) = 0._r8
!           colmol(lay) = 0._r8
         if (colco2(lay) .eq. 0._r8) colco2(lay) = 1.e-32_r8 * coldry(lay)
         if (coln2o(lay) .eq. 0._r8) coln2o(lay) = 1.e-32_r8 * coldry(lay)
         if (colch4(lay) .eq. 0._r8) colch4(lay) = 1.e-32_r8 * coldry(lay)
         if (colo2(lay) .eq. 0._r8) colo2(lay) = 1.e-32_r8 * coldry(lay)
! Using E = 1334.2 cm-1.
         co2reg = 3.55e-24_r8 * coldry(lay)
         co2mult(lay)= (colco2(lay) - co2reg) * &
               272.63_r8*exp(-1919.4_r8/tavel(lay))/(8.7604e-4_r8*tavel(lay))
         goto 5400

! Above laytrop.
 5300    continue

! Set up factors needed to separately include the water vapor
! foreign-continuum in the calculation of absorption coefficient.

         forfac(lay) = scalefac / (1.+water)
         factor = (tavel(lay)-188.0_r8)/36.0_r8
         indfor(lay) = 3
         forfrac(lay) = factor - 1.0_r8

! Calculate needed column amounts.

         colh2o(lay) = 1.e-20_r8 * wkl(1,lay)
         colco2(lay) = 1.e-20_r8 * wkl(2,lay)
         colo3(lay)  = 1.e-20_r8 * wkl(3,lay)
         coln2o(lay) = 1.e-20_r8 * wkl(4,lay)
         colch4(lay) = 1.e-20_r8 * wkl(6,lay)
         colo2(lay)  = 1.e-20_r8 * wkl(7,lay)
         colmol(lay) = 1.e-20_r8 * coldry(lay) + colh2o(lay)
         if (colco2(lay) .eq. 0._r8) colco2(lay) = 1.e-32_r8 * coldry(lay)
         if (coln2o(lay) .eq. 0._r8) coln2o(lay) = 1.e-32_r8 * coldry(lay)
         if (colch4(lay) .eq. 0._r8) colch4(lay) = 1.e-32_r8 * coldry(lay)
         if (colo2(lay)  .eq. 0._r8) colo2(lay)  = 1.e-32_r8 * coldry(lay)
         co2reg = 3.55e-24_r8 * coldry(lay)
         co2mult(lay)= (colco2(lay) - co2reg) * &
               272.63_r8*exp(-1919.4_r8/tavel(lay))/(8.7604e-4_r8*tavel(lay))

         selffac(lay) = 0._r8
         selffrac(lay)= 0._r8
         indself(lay) = 0

 5400    continue

! We have now isolated the layer ln pressure and temperature,
! between two reference pressures and two reference temperatures
! (for each reference pressure).  We multiply the pressure
! fraction FP with the appropriate temperature fractions to get
! the factors that will be needed for the interpolation that yields
! the optical depths (performed in routines TAUGBn for band n).

         compfp = 1._r8 - fp
         fac10(lay) = compfp * ft
         fac00(lay) = compfp * (1._r8 - ft)
         fac11(lay) = fp * ft1
         fac01(lay) = fp * (1._r8 - ft1)

! End layer loop
      enddo

      end subroutine setcoef_sw

!***************************************************************************
      subroutine swatmref
!***************************************************************************

      save

! These pressures are chosen such that the ln of the first pressure
! has only a few non-zero digits (i.e. ln(PREF(1)) = 6.96000) and
! each subsequent ln(pressure) differs from the previous one by 0.2.

      pref(:) = (/ &
          1.05363e+03_r8,8.62642e+02_r8,7.06272e+02_r8,5.78246e+02_r8,4.73428e+02_r8, &
          3.87610e+02_r8,3.17348e+02_r8,2.59823e+02_r8,2.12725e+02_r8,1.74164e+02_r8, &
          1.42594e+02_r8,1.16746e+02_r8,9.55835e+01_r8,7.82571e+01_r8,6.40715e+01_r8, &
          5.24573e+01_r8,4.29484e+01_r8,3.51632e+01_r8,2.87892e+01_r8,2.35706e+01_r8, &
          1.92980e+01_r8,1.57998e+01_r8,1.29358e+01_r8,1.05910e+01_r8,8.67114e+00_r8, &
          7.09933e+00_r8,5.81244e+00_r8,4.75882e+00_r8,3.89619e+00_r8,3.18993e+00_r8, &
          2.61170e+00_r8,2.13828e+00_r8,1.75067e+00_r8,1.43333e+00_r8,1.17351e+00_r8, &
          9.60789e-01_r8,7.86628e-01_r8,6.44036e-01_r8,5.27292e-01_r8,4.31710e-01_r8, &
          3.53455e-01_r8,2.89384e-01_r8,2.36928e-01_r8,1.93980e-01_r8,1.58817e-01_r8, &
          1.30029e-01_r8,1.06458e-01_r8,8.71608e-02_r8,7.13612e-02_r8,5.84256e-02_r8, &
          4.78349e-02_r8,3.91639e-02_r8,3.20647e-02_r8,2.62523e-02_r8,2.14936e-02_r8, &
          1.75975e-02_r8,1.44076e-02_r8,1.17959e-02_r8,9.65769e-03_r8 /)

      preflog(:) = (/ &
           6.9600e+00_r8, 6.7600e+00_r8, 6.5600e+00_r8, 6.3600e+00_r8, 6.1600e+00_r8, &
           5.9600e+00_r8, 5.7600e+00_r8, 5.5600e+00_r8, 5.3600e+00_r8, 5.1600e+00_r8, &
           4.9600e+00_r8, 4.7600e+00_r8, 4.5600e+00_r8, 4.3600e+00_r8, 4.1600e+00_r8, &
           3.9600e+00_r8, 3.7600e+00_r8, 3.5600e+00_r8, 3.3600e+00_r8, 3.1600e+00_r8, &
           2.9600e+00_r8, 2.7600e+00_r8, 2.5600e+00_r8, 2.3600e+00_r8, 2.1600e+00_r8, &
           1.9600e+00_r8, 1.7600e+00_r8, 1.5600e+00_r8, 1.3600e+00_r8, 1.1600e+00_r8, &
           9.6000e-01_r8, 7.6000e-01_r8, 5.6000e-01_r8, 3.6000e-01_r8, 1.6000e-01_r8, &
          -4.0000e-02_r8,-2.4000e-01_r8,-4.4000e-01_r8,-6.4000e-01_r8,-8.4000e-01_r8, &
          -1.0400e+00_r8,-1.2400e+00_r8,-1.4400e+00_r8,-1.6400e+00_r8,-1.8400e+00_r8, &
          -2.0400e+00_r8,-2.2400e+00_r8,-2.4400e+00_r8,-2.6400e+00_r8,-2.8400e+00_r8, &
          -3.0400e+00_r8,-3.2400e+00_r8,-3.4400e+00_r8,-3.6400e+00_r8,-3.8400e+00_r8, &
          -4.0400e+00_r8,-4.2400e+00_r8,-4.4400e+00_r8,-4.6400e+00_r8 /)

! These are the temperatures associated with the respective
! pressures for the MLS standard atmosphere.

      tref(:) = (/ &
           2.9420e+02_r8, 2.8799e+02_r8, 2.7894e+02_r8, 2.6925e+02_r8, 2.5983e+02_r8, &
           2.5017e+02_r8, 2.4077e+02_r8, 2.3179e+02_r8, 2.2306e+02_r8, 2.1578e+02_r8, &
           2.1570e+02_r8, 2.1570e+02_r8, 2.1570e+02_r8, 2.1706e+02_r8, 2.1858e+02_r8, &
           2.2018e+02_r8, 2.2174e+02_r8, 2.2328e+02_r8, 2.2479e+02_r8, 2.2655e+02_r8, &
           2.2834e+02_r8, 2.3113e+02_r8, 2.3401e+02_r8, 2.3703e+02_r8, 2.4022e+02_r8, &
           2.4371e+02_r8, 2.4726e+02_r8, 2.5085e+02_r8, 2.5457e+02_r8, 2.5832e+02_r8, &
           2.6216e+02_r8, 2.6606e+02_r8, 2.6999e+02_r8, 2.7340e+02_r8, 2.7536e+02_r8, &
           2.7568e+02_r8, 2.7372e+02_r8, 2.7163e+02_r8, 2.6955e+02_r8, 2.6593e+02_r8, &
           2.6211e+02_r8, 2.5828e+02_r8, 2.5360e+02_r8, 2.4854e+02_r8, 2.4348e+02_r8, &
           2.3809e+02_r8, 2.3206e+02_r8, 2.2603e+02_r8, 2.2000e+02_r8, 2.1435e+02_r8, &
           2.0887e+02_r8, 2.0340e+02_r8, 1.9792e+02_r8, 1.9290e+02_r8, 1.8809e+02_r8, &
           1.8329e+02_r8, 1.7849e+02_r8, 1.7394e+02_r8, 1.7212e+02_r8 /)

      end subroutine swatmref

      end module rrtmg_sw_setcoef

max_stars_repo_licenses: NCSA Unlicense


text:
C Compute the phase, incidence, and emission angles for a point on the ring.
C
      subroutine ring_light(vsc,vsun,r,lon,i,e,p)
      implicit none
      real*8 vsc(3)             !Vector from target center to S/C (RS vector)
      real*8 vsun(3)            !Vector from target center to sun
      real*8 r,lon		!Radius and longitude of ring point
      real*8 i,e,p		!Incidence, emission, and phase angles (deg)

      real*8 v(3),s(3),c(3),PI,rpd,dpr
      real*8 solrange,scrange,elon
      integer j

      PI = 3.141592653589793D0
      rpd = PI/180.d0		!Radians per degree
      dpr = 180.d0/PI           !Degrees per radian
      elon = 360.d0 - lon	!East longitude

C     ....Compute vector from planet center to ring point
      v(1) = r*dcos(elon*rpd)
      v(2) = r*dsin(elon*rpd)
      v(3) = 0.

      do j=1,3
         s(j) = vsun(j) - v(j)		!Vector from ring point to sun
         c(j) = vsc(j)  - v(j)		!Vector from ring point to S/C
      enddo

C     ....Convert to unit vectors
      solrange = dsqrt(s(1)**2+s(2)**2+s(3)**2)		!Sun range
      scrange  = dsqrt(c(1)**2+c(2)**2+c(3)**2)		!S/C range
      do j=1,3
         s(j) = s(j)/solrange
         c(j) = c(j)/scrange
      enddo

C     ....First assume the sun illuminates the northern half of the ring plane.
C     ....Unit normal to ring plane is n = (0,0,1)
C     ....Cosine of incidence angle is dot product of s and n
      i = dacos(s(3))*dpr
C     ....Cosine of emission angle is dot product of c and n
      e = dacos(c(3))*dpr
C     ....cosine of phase angle is dot product of s and c
      p = dacos(s(1)*c(1)+s(2)*c(2)+s(3)*c(3))*dpr
C     ....If sun illuminates southern half of ring plane, then
C     ....unit normal to ring plane is n = (0,0,-1)
      if (i .gt. 90.) then
         i = 180. - i
         e = 180. - e
      endif
      return
      end

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
module test_var
    use testdrive, only : new_unittest, unittest_type, error_type, check
    use stdlib_kinds, only: sp, dp, int32, int64
    use stdlib_stats, only: var
    use,intrinsic :: ieee_arithmetic, only : ieee_is_nan
    implicit none

    real(sp), parameter :: sptol = 1000 * epsilon(1._sp)
    real(dp), parameter :: dptol = 1000 * epsilon(1._dp)

    integer(int32), parameter :: i321(5) = [1, 2, 3, 4, 5]
    integer(int64), parameter :: i641(5) = [1, 2, 3, 4, 5]

    real(sp), parameter :: s1(5) = [1.0_sp, 2.0_sp, 3.0_sp, 4.0_sp, 5.0_sp]
    real(dp), parameter :: d1(5) = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp]

    real(dp), parameter :: d(4, 3) = reshape([1._dp, 3._dp, 5._dp, 7._dp,&
                                   2._dp, 4._dp, 6._dp, 8._dp,&
                                   9._dp, 10._dp, 11._dp, 12._dp], [4, 3])
    real(dp), parameter :: d3(4, 3, 3) = reshape([d, d*2, d*4], shape(d3))
    real(sp), parameter :: s(4, 3) = d
    real(sp), parameter :: s3(4, 3, 3) = reshape([s, s*2, s*4], shape(s3))

    integer(int32), parameter :: i32(4, 3) = d
    integer(int32), parameter :: i323(4, 3, 3) = d3
    integer(int64), parameter :: i64(4, 3) = d
    integer(int64), parameter :: i643(4, 3, 3) = d3

    complex(sp), parameter :: cs1(5) = [ cmplx(0.57706_sp, 0.00000_sp, sp),&
                            cmplx(0.00000_sp, 1.44065_sp, sp),&
                            cmplx(1.26401_sp, 0.00000_sp, sp),&
                            cmplx(0.00000_sp, 0.88833_sp, sp),&
                            cmplx(1.14352_sp, 0.00000_sp, sp)]
    complex(dp), parameter :: cd1(5) = [ cmplx(0.57706_dp, 0.00000_dp,kind=dp),&
                            cmplx(0.00000_dp, 1.44065_dp,kind=dp),&
                            cmplx(1.26401_dp, 0.00000_dp,kind=dp),&
                            cmplx(0.00000_dp, 0.88833_dp,kind=dp),&
                            cmplx(1.14352_dp, 0.00000_dp,kind=dp)]
    complex(sp), parameter :: cs(5,3) = reshape([cs1, cs1*3.0_sp, cs1*1.5_sp], shape(cs))
    complex(dp), parameter :: cd(5,3) = reshape([cd1, cd1*3.0_dp, cd1*1.5_dp], shape(cd))

contains

    !> Collect all exported unit tests
    subroutine collect_var(testsuite)
        !> Collection of tests
        type(unittest_type), allocatable, intent(out) :: testsuite(:)

        testsuite = [ &
            new_unittest("sp_1dim", test_sp_1dim), &
            new_unittest("sp_1dim_mask", test_sp_1dim_mask), &
            new_unittest("sp_2dim", test_sp_2dim), &
            new_unittest("sp_2dim_mask", test_sp_2dim_mask), &
            new_unittest("sp_2dim_mask_array", test_sp_2dim_mask_array), &
            new_unittest("sp_3dim", test_sp_3dim), &
            new_unittest("sp_3dim_mask", test_sp_3dim_mask), &
            new_unittest("sp_3dim_mask_array", test_sp_3dim_mask_array), &
            new_unittest("dp_1dim", test_dp_1dim), &
            new_unittest("dp_1dim_mask", test_dp_1dim_mask), &
            new_unittest("dp_1dim_mask_array", test_dp_1dim_mask_array), &
            new_unittest("dp_2dim", test_dp_2dim), &
            new_unittest("dp_2dim_mask", test_dp_2dim_mask), &
            new_unittest("dp_2dim_mask_array", test_dp_2dim_mask_array), &
            new_unittest("dp_3dim", test_dp_3dim), &
            new_unittest("dp_3dim_mask", test_dp_3dim_mask), &
            new_unittest("dp_3dim_mask_array", test_dp_3dim_mask_array), &
            new_unittest("int32_1dim", test_int32_1dim), &
            new_unittest("int32_1dim_mask", test_int32_1dim_mask), &
            new_unittest("int32_1dim_mask_array", test_int32_1dim_mask_array), &
            new_unittest("int32_2dim", test_int32_2dim), &
            new_unittest("int32_2dim_mask", test_int32_2dim_mask), &
            new_unittest("int32_2dim_mask_array", test_int32_2dim_mask_array), &
            new_unittest("int32_3dim", test_int32_3dim), &
            new_unittest("int32_3dim_mask", test_int32_3dim_mask), &
            new_unittest("int32_3dim_mask_array", test_int32_3dim_mask_array), &
            new_unittest("int64_1dim", test_int64_1dim), &
            new_unittest("int64_1dim_mask", test_int64_1dim_mask), &
            new_unittest("int641_1dim_mask_array", test_int641_1dim_mask_array), &
            new_unittest("int64_2dim", test_int64_2dim), &
            new_unittest("int64_2dim_mask", test_int64_2dim_mask), &
            new_unittest("int64_2dim_mask_array", test_int64_2dim_mask_array), &
            new_unittest("int64_3dim", test_int64_3dim), &
            new_unittest("int64_3dim_mask", test_int64_3dim_mask), &
            new_unittest("int64_3dim_mask_array", test_int64_3dim_mask_array), &
            new_unittest("csp_1dim", test_csp_1dim), &
            new_unittest("csp_1dim_mask", test_csp_1dim_mask), &
            new_unittest("csp_1dim_mask_array", test_csp_1dim_mask_array), &
            new_unittest("csp_2dim", test_csp_2dim), &
            new_unittest("csp_2dim_mask", test_csp_2dim_mask), &
            new_unittest("csp_2dim_mask_array", test_csp_2dim_mask_array), &
            new_unittest("cdp_1dim", test_cdp_1dim), &
            new_unittest("cdp_1dim_mask", test_cdp_1dim_mask), &
            new_unittest("cdp_1dim_mask_array", test_cdp_1dim_mask_array), &
            new_unittest("cdp_2dim", test_cdp_2dim), &
            new_unittest("cdp_2dim_mask", test_cdp_2dim_mask), &
            new_unittest("cdp_2dim_mask_array", test_cdp_2dim_mask_array) &
            ]
    end subroutine collect_var

    subroutine test_sp_1dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(s1) - 2.5) < sptol)
        call check(error, abs(var(s1, dim=1) - 2.5) < sptol)

    end subroutine test_sp_1dim

    subroutine test_sp_1dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(s1, .false.)))
        call check(error, ieee_is_nan(var(s1, 1, .false.)))

    end subroutine test_sp_1dim_mask

    subroutine test_sp_1dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(s1, s1 < 5) - 5./3.) < sptol)
        call check(error, ieee_is_nan((var(s1, s1 < 0.))))
        call check(error, ieee_is_nan((var(s1, s1 == 1.))))
        call check(error, abs(var(s1, 1, s1 < 5) - 5./3.) < sptol)

    end subroutine test_sp_1dim_mask_array

    subroutine test_sp_2dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(s) - 13) < sptol)
        call check(error, all( abs( var(s, 1) - [20. / 3., 20. / 3., 5. / 3.]) < sptol))
        call check(error, all( abs( var(s, 2) - [19.0, 43. / 3., 31. / 3. , 7.0]) < sptol))

    end subroutine test_sp_2dim

    subroutine test_sp_2dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(s, .false.)))
        call check(error, any(ieee_is_nan(var(s, 1, .false.))))
        call check(error, any(ieee_is_nan(var(s, 2, .false.))))

    end subroutine test_sp_2dim_mask

    subroutine test_sp_2dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(s, s < 11) - 27.5 / 3.) < sptol)
        call check(error, all( abs( var(s, 1, s < 11) - [20. / 3., 20. / 3., 0.5]) < sptol))
        call check(error, all( abs( var(s, 2, s < 11) - [19.0, 43. / 3., 0.5 , 0.5]) < sptol))

    end subroutine test_sp_2dim_mask_array

    subroutine test_sp_3dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(s3) - 153.4) < sptol)
        call check(error, all( abs( var(s3, 1) -&
                 reshape([20. / 3., 20. / 3., 5. / 3.,&
                          4* 20. / 3., 4* 20. / 3., 4* 5. / 3.,&
                          16* 20. / 3., 16* 20. / 3., 16* 5. / 3.],&
                          [size(s3,2), size(s3,3)]))&
                 < sptol))
        call check(error, all( abs( var(s3, 2) -&
                 reshape([19.0, 43. / 3., 31. / 3. , 7.0,&
                          4* 19.0, 4* 43. / 3., 4* 31. / 3. , 4* 7.0,&
                          16* 19.0, 16* 43. / 3., 16* 31. / 3. , 16* 7.0],&
                          [size(s3,1), size(s3,3)] ))&
                 < sptol))
        call check(error, all(abs( var(s3, 3) -&
                 reshape([ 7./3., 21., 175./3.,&
                           343./3., 28./3., 112./3.,&
                           84., 448./3., 189.,&
                           700./3., 847./3., 336.], [size(s3,1), size(s3,2)] ))&
                 < sptol))

    end subroutine test_sp_3dim

    subroutine test_sp_3dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(s3, .false.)))
        call check(error, any(ieee_is_nan(var(s3, 1, .false.))))
        call check(error, any(ieee_is_nan(var(s3, 2, .false.))))
        call check(error, any(ieee_is_nan(var(s3, 3, .false.))))

    end subroutine test_sp_3dim_mask

    subroutine test_sp_3dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(s3, s3 < 11) - 8.2205877_sp) < sptol)
        call check(error, all( abs( var(s3, 1, s3 < 45) -&
                  reshape([20./3., 20./3., 5./3.,  80./3., 80./3., 20./3.,&
                           320./3., 320./3., 16.],&
                  [size(s3, 2), size(s3, 3)])) < sptol ))
        call check(error, any( ieee_is_nan( var(s3, 2, s3 < 25))))
        call check(error, all( abs( var(s3, 3, s3 < 25) -&
                 reshape([ 7./3., 21., 175./3.,&
                           24.5, 28./3., 112./3.,&
                           84., 32., 40.5,&
                           50., 60.5, 72.], [size(s3,1), size(s3,2)] ))&
                 < sptol ))

    end subroutine test_sp_3dim_mask_array

    subroutine test_dp_1dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(d1) - 2.5) < dptol)
        call check(error, abs(var(d1, 1) - 2.5) < dptol)

    end subroutine test_dp_1dim

    subroutine test_dp_1dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(d1, .false.)))
        call check(error, ieee_is_nan(var(d1, 1, .false.)))

    end subroutine test_dp_1dim_mask

    subroutine test_dp_1dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(d1, d1 < 5) - 5._dp/3._dp) < dptol)
        call check(error, ieee_is_nan((var(d1, d1 < 0.))))
        call check(error, ieee_is_nan((var(d1, d1 == 1.))))
        call check(error, abs(var(d1, 1, d1 < 5) - 5._dp/3._dp) < dptol)

    end subroutine test_dp_1dim_mask_array

    subroutine test_dp_2dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(d) - 13) < dptol)
        call check(error, all( abs( var(d,1) -&
            [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol))
        call check(error, all( abs( var(d,2) -&
            [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol))

    end subroutine test_dp_2dim

    subroutine test_dp_2dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

    call check(error, ieee_is_nan(var(d, .false.)))
    call check(error, any(ieee_is_nan(var(d, 1, .false.))))
    call check(error, any(ieee_is_nan(var(d, 2, .false.))))

    end subroutine test_dp_2dim_mask

    subroutine test_dp_2dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

    call check(error, abs(var(d, d < 11) - 27.5_dp / 3._dp) < dptol)
    call check(error, all( abs( var(d, 1, d < 11) -&
                 [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol))
    call check(error, all( abs( var(d, 2, d < 11) -&
                 [19.0_dp, 43._dp / 3._dp, 0.5_dp, 0.5_dp]) < dptol))

    end subroutine test_dp_2dim_mask_array

    subroutine test_dp_3dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(d3) - 153.4_dp) < dptol)
        call check(error, all( abs( var(d3, 1) -&
                 reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,&
                          4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,&
                          16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],&
                          [size(d3,2), size(d3,3)]))&
                 < dptol))
        call check(error, all( abs( var(d3, 2) -&
                 reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,&
                          4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,&
                          16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,&
                          16* 7.0_dp],&
                          [size(d3,1), size(d3,3)] ))&
                 < dptol))
        call check(error, all(abs( var(d3, 3) -&
                 reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                           343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,&
                           84._dp, 448._dp/3._dp, 189._dp,&
                           700._dp/3._dp, 847._dp/3._dp, 336._dp],&
                           [size(d3,1), size(d3,2)] ))&
                 < dptol))

    end subroutine test_dp_3dim

    subroutine test_dp_3dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(d3, .false.)))
        call check(error, any(ieee_is_nan(var(d3, 1, .false.))))
        call check(error, any(ieee_is_nan(var(d3, 2, .false.))))
        call check(error, any(ieee_is_nan(var(d3, 3, .false.))))

    end subroutine test_dp_3dim_mask

    subroutine test_dp_3dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(d3, d3 < 25) - 46.041379310344823_dp) < dptol)
        call check(error, all( abs( var(d3, 1, d3 < 45) -&
                  reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,&
                           80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,&
                           320._dp/3._dp, 320._dp/3._dp, 16._dp],&
                           [size(d3, 2), size(d3, 3)]))&
                  < dptol ))
        call check(error, any( ieee_is_nan( var(d3, 2, d3 < 25))))
        call check(error, all( abs( var(d3, 3, d3 < 25) -&
                 reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                           24.5_dp, 28._dp/3._dp, 112._dp/3._dp,&
                           84._dp, 32._dp, 40.5_dp,&
                           50._dp, 60.5_dp, 72._dp],&
                           [size(d3,1), size(d3,2)] ))&
                 < dptol ))

    end subroutine test_dp_3dim_mask_array

    subroutine test_int32_1dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i321) - 2.5) < dptol)
        call check(error, abs(var(i321, 1) - 2.5) < dptol)

    end subroutine test_int32_1dim

    subroutine test_int32_1dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(i321, .false.)))
        call check(error, ieee_is_nan(var(i321, 1, .false.)))

    end subroutine test_int32_1dim_mask

    subroutine test_int32_1dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i321, i321 < 5) - 5._dp/3._dp) < dptol)
        call check(error, ieee_is_nan((var(i321, i321 < 0))))
        call check(error, ieee_is_nan((var(i321, i321 == 1))))
        call check(error, abs(var(i321, 1, i321 < 5) - 5._dp/3._dp) < dptol)

    end subroutine test_int32_1dim_mask_array

    subroutine test_int32_2dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i32) - 13) < dptol)
        call check(error, all( abs( var(i32,1) -&
            [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol))
        call check(error, all( abs( var(i32,2) -&
            [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol))

    end subroutine test_int32_2dim

    subroutine test_int32_2dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(i32, .false.)))
        call check(error, any(ieee_is_nan(var(i32, 1, .false.))))
        call check(error, any(ieee_is_nan(var(i32, 2, .false.))))

    end subroutine test_int32_2dim_mask

    subroutine test_int32_2dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i32, i32 < 11) - 27.5_dp / 3._dp) < dptol)
        call check(error, all( abs( var(i32, 1, i32 < 11) -&
            [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol))
        call check(error, all( abs( var(i32, 2, i32 < 11) -&
            [19.0_dp, 43._dp / 3._dp, 0.5_dp, 0.5_dp]) < dptol))

    end subroutine test_int32_2dim_mask_array

    subroutine test_int32_3dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i323) - 153.4_dp) < dptol)
        call check(error, all( abs( var(i323, 1) -&
                 reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,&
                          4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,&
                          16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],&
                          [size(i323,2), size(i323,3)]))&
                 < dptol))
        call check(error, all( abs( var(i323, 2) -&
                 reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,&
                          4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,&
                          16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,&
                          16* 7.0_dp],&
                          [size(i323,1), size(i323,3)] ))&
                 < dptol))
        call check(error, all(abs( var(i323, 3) -&
                 reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                           343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,&
                           84._dp, 448._dp/3._dp, 189._dp,&
                           700._dp/3._dp, 847._dp/3._dp, 336._dp],&
                           [size(i323,1), size(i323,2)] ))&
                 < dptol))

    end subroutine test_int32_3dim

    subroutine test_int32_3dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(i323, .false.)))
        call check(error, any(ieee_is_nan(var(i323, 1, .false.))))
        call check(error, any(ieee_is_nan(var(i323, 2, .false.))))
        call check(error, any(ieee_is_nan(var(i323, 3, .false.))))

    end subroutine test_int32_3dim_mask

    subroutine test_int32_3dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i323, i323 < 25) - 46.041379310344823_dp) < dptol)
        call check(error, all( abs( var(i323, 1, i323 < 45) -&
                  reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,&
                           80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,&
                           320._dp/3._dp, 320._dp/3._dp, 16._dp],&
                           [size(i323, 2), size(i323, 3)]))&
                  < dptol ))
        call check(error, any( ieee_is_nan( var(i323, 2, i323 < 25))))
        call check(error, all( abs( var(i323, 3, i323 < 25) -&
                 reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                           24.5_dp, 28._dp/3._dp, 112._dp/3._dp,&
                           84._dp, 32._dp, 40.5_dp,&
                           50._dp, 60.5_dp, 72._dp],&
                           [size(i323,1), size(i323,2)] ))&
                 < dptol ))

    end subroutine test_int32_3dim_mask_array

    subroutine test_int64_1dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i641) - 2.5) < dptol)
        call check(error, abs(var(i641, 1) - 2.5) < dptol)

    end subroutine test_int64_1dim

    subroutine test_int64_1dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(i641, .false.)))
        call check(error, ieee_is_nan(var(i641, 1, .false.)))

    end subroutine test_int64_1dim_mask

    subroutine test_int641_1dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i641, i641 < 5) - 5._dp/3._dp) < dptol)
        call check(error, ieee_is_nan((var(i641, i641 < 0))))
        call check(error, ieee_is_nan((var(i641, i641 == 1))))
        call check(error, abs(var(i641, 1, i641 < 5) - 5._dp/3._dp) < dptol)

    end subroutine test_int641_1dim_mask_array

    subroutine test_int64_2dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i64) - 13) < dptol)
        call check(error, all( abs( var(i64,1) -&
            [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol))
        call check(error, all( abs( var(i64,2) -&
            [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol))

    end subroutine test_int64_2dim

    subroutine test_int64_2dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(i64, .false.)))
        call check(error, any(ieee_is_nan(var(i64, 1, .false.))))
        call check(error, any(ieee_is_nan(var(i64, 2, .false.))))

    end subroutine test_int64_2dim_mask

    subroutine test_int64_2dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i64, i64 < 11) - 27.5_dp / 3._dp) < dptol)
        call check(error, all( abs( var(i64, 1, i64 < 11) -&
            [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol))
        call check(error, all( abs( var(i64, 2, i64 < 11) -&
            [19.0_dp, 43._dp / 3._dp, 0.5_dp, 0.5_dp]) < dptol))

    end subroutine test_int64_2dim_mask_array

    subroutine test_int64_3dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i643) - 153.4_dp) < dptol)
        call check(error, all( abs( var(i643, 1) -&
                 reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,&
                          4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,&
                          16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],&
                          [size(i643,2), size(i643,3)]))&
                 < dptol))
        call check(error, all( abs( var(i643, 2) -&
                 reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,&
                          4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,&
                          16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,&
                          16* 7.0_dp],&
                          [size(i643,1), size(i643,3)] ))&
                 < dptol))
        call check(error, all(abs( var(i643, 3) -&
                 reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                           343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,&
                           84._dp, 448._dp/3._dp, 189._dp,&
                           700._dp/3._dp, 847._dp/3._dp, 336._dp],&
                           [size(i643,1), size(i643,2)] ))&
                 < dptol))

    end subroutine test_int64_3dim

    subroutine test_int64_3dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(i643, .false.)))
        call check(error, any(ieee_is_nan(var(i643, 1, .false.))))
        call check(error, any(ieee_is_nan(var(i643, 2, .false.))))
        call check(error, any(ieee_is_nan(var(i643, 3, .false.))))

    end subroutine test_int64_3dim_mask

    subroutine test_int64_3dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(i643, i643 < 25) - 46.041379310344823_dp) < dptol)
        call check(error, all( abs( var(i643, 1, i643 < 45) -&
                  reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,&
                           80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,&
                           320._dp/3._dp, 320._dp/3._dp, 16._dp],&
                           [size(i643, 2), size(i643, 3)]))&
                  < dptol ))
        call check(error, any( ieee_is_nan( var(i643, 2, i643 < 25))))
        call check(error, all( abs( var(i643, 3, i643 < 25) -&
                 reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                           24.5_dp, 28._dp/3._dp, 112._dp/3._dp,&
                           84._dp, 32._dp, 40.5_dp,&
                           50._dp, 60.5_dp, 72._dp],&
                           [size(i643,1), size(i643,2)] ))&
                 < dptol ))

    end subroutine test_int64_3dim_mask_array

    subroutine test_csp_1dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(cs1) - (var(real(cs1)) + var(aimag(cs1)))) < sptol)
        call check(error, abs(var(cs1, dim=1) - (var(real(cs1),1) + var(aimag(cs1), 1)) ) < sptol)

    end subroutine test_csp_1dim

    subroutine test_csp_1dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(cs1, .false.)))
        call check(error, ieee_is_nan(var(cs1, 1, .false.)))

    end subroutine test_csp_1dim_mask

    subroutine test_csp_1dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(cs1, aimag(cs1) == 0) - var(real(cs1), aimag(cs1) == 0)) < sptol)
        call check(error, abs(var(cs1, 1, aimag(cs1) == 0) - var(real(cs1), 1, aimag(cs1) == 0)) < sptol)

    end subroutine test_csp_1dim_mask_array

    subroutine test_csp_2dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(cs) - (var(real(cs)) + var(aimag(cs)))) < sptol)
        call check(error, all( abs( var(cs, 1) - (var(real(cs), 1) + var(aimag(cs), 1))) < sptol))
        call check(error, all( abs( var(cs, 2) - (var(real(cs), 2) + var(aimag(cs), 2))) < sptol))

    end subroutine test_csp_2dim

    subroutine test_csp_2dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(cs, .false.)))
        call check(error, any(ieee_is_nan(var(cs, 1, .false.))))
        call check(error, any(ieee_is_nan(var(cs, 2, .false.))))

    end subroutine test_csp_2dim_mask

    subroutine test_csp_2dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(cs, aimag(cs) == 0) - var(real(cs), aimag(cs) == 0)) < sptol)
        call check(error, all( abs( var(cs, 1, aimag(cs) == 0) - var(real(cs), 1, aimag(cs) == 0)) < sptol))
        call check(error, any( ieee_is_nan( var(cs, 2, aimag(cs) == 0))))

    end subroutine test_csp_2dim_mask_array

    subroutine test_cdp_1dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(cd1) - (var(real(cd1)) + var(aimag(cd1)))) < dptol)
        call check(error, abs(var(cd1, dim=1) - (var(real(cd1),1) + var(aimag(cd1), 1)) ) < dptol)

    end subroutine test_cdp_1dim

    subroutine test_cdp_1dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(cd1, .false.)))
        call check(error, ieee_is_nan(var(cd1, 1, .false.)))

    end subroutine test_cdp_1dim_mask

    subroutine test_cdp_1dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(cd1, aimag(cd1) == 0) - var(real(cd1), aimag(cd1) == 0)) < dptol)
        call check(error, abs(var(cd1, 1, aimag(cd1) == 0) - var(real(cd1), 1, aimag(cd1) == 0)) < dptol)

    end subroutine test_cdp_1dim_mask_array

    subroutine test_cdp_2dim(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(cd) - (var(real(cd)) + var(aimag(cd)))) < dptol)
        call check(error, all( abs( var(cd, 1) - (var(real(cd), 1) + var(aimag(cd), 1))) < dptol))
        call check(error, all( abs( var(cd, 2) - (var(real(cd), 2) + var(aimag(cd), 2))) < dptol))

    end subroutine test_cdp_2dim

    subroutine test_cdp_2dim_mask(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, ieee_is_nan(var(cd, .false.)))
        call check(error, any(ieee_is_nan(var(cd, 1, .false.))))
        call check(error, any(ieee_is_nan(var(cd, 2, .false.))))

    end subroutine test_cdp_2dim_mask

    subroutine test_cdp_2dim_mask_array(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, abs(var(cd, aimag(cd) == 0) - var(real(cd), aimag(cd) == 0)) < dptol)
        call check(error, all( abs( var(cd, 1, aimag(cd) == 0) - var(real(cd), 1, aimag(cd) == 0)) < dptol))
        call check(error, any( ieee_is_nan( var(cd, 2, aimag(cd) == 0))))

    end subroutine test_cdp_2dim_mask_array

end module

program tester
    use, intrinsic :: iso_fortran_env, only : error_unit
    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
    use test_var, only : collect_var
    implicit none
    integer :: stat, is
    type(testsuite_type), allocatable :: testsuites(:)
    character(len=*), parameter :: fmt = '("#", *(1x, a))'

    stat = 0

    testsuites = [ &
        new_testsuite("var", collect_var) &
        ]

    do is = 1, size(testsuites)
        write(error_unit, fmt) "Testing:", testsuites(is)%name
        call run_testsuite(testsuites(is)%collect, error_unit, stat)
    end do

    if (stat > 0) then
        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
        error stop
    end if
end program

max_stars_repo_name: freevryheid/stdlib

max_issues_repo_name: seanwallawalla-forks/stdlib

max_forks_repo_name: seanwallawalla-forks/stdlib