text:
	SUBROUTINE IM_SIMG ( cproj, imgfil, iret )
C************************************************************************
C* IM_SIMG								*
C*									*
C* This subroutine reads an image file header and sets the navigation	*
C* for an image file.							*
C*									*
C* IM_SIMG ( CPROJ, IMGFIL, IRET )					*
C*									*
C* Input parameters:							*
C*	CPROJ		CHAR*		Projection name			*
C*	IMGFIL		CHAR*		Image file name			*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 +1 = No entry found for image	*
C*					 -1 = Image file not found	*
C*					 -3 = Invalid image file format *
C*					 -4 = Invalid image product	*
C*					 -5 = Invalid navigation	*
C*					 -7 = Could not open imtyp table*
C**									*
C* Log:									*
C* J. Cowie/COMET	 2/95	Modified from GG_SAT()			*
C* C. Lin/EAI		 6/95	Restructure				*
C* J. Cowie/COMET	10/95	Accomodate small image products		*
C* J. Cowie/COMET	11/95	Add GINI file check			*
C* S. Jacobs/NCEP	 7/96	Changed IM_ARHD -> IM_AR2GM and		*
C*				IM_GNHD -> IM_GI2GM; Modified checks	*
C*				for different image types; Added MV_BTOI*
C* S. Jacobs/NCEP	 7/96	Fixed call to IM_GI2GM			*
C* R. Lindsay/GSC	12/96	Added code for UNISYS NIDS composites	*
C* J. Cowie/COMET	 1/97	Changed IMGDEF common variable names	*
C* S. Jacobs/NCEP	 4/97	Added call to IM_RTBL to read image tbl	*
C* S. Chiswell/Unidata	12/97	Added byte swap check for nids for OSF	*
C* S. Jacobs/NCEP	 2/98	Removed check for valid natl/regl prods	*
C* S. Jacobs/NCEP	 4/98	Fixed array size for ipflip		*
C* S. Jacobs/NCEP	 5/99	Added byte swap check for Unisys format	*
C* S. Jacobs/NCEP	10/99	Added check for NetCDF format images	*
C* R. Curtis/EAI	 8/00   Updated for NetCDF files                *
C* S. Chiswell/Unidata	11/00	Added compressed NIDS format		*
C* S. Chiswell/Unidata	 2/02	Added imcmn.cmn for imcalbar init	*
C* J. Levit/NCEP        10/07   Increased amount of available radar     *
C*                              products from 149 to 999.               *
C* m.gamazaychikov/CWS	01/10	Add calls to IM_DBGINI2GM, IM_DBNIDS2GM,*
C*                              IM_DBMOSAIC2GM, IM_DBMCIDAS2GM providing*
C*				ability to display GINI, NIDS, MOSAIC   *
C                               MCIDAS images stored in AWIPS II	*
C*				database				*
C* M. James/Unidata      2/10   Modified imdoff def. for FOS/WMO hdrs   *
C* X. Guo/CWS           04/10   Added codes to support 94 product       *
C* X. Guo/CWS           04/10   Added codes to support 99,134 and 135   *
C*                              product. Also set imftyp to IFHINIDS    *
C* S. Jacobs/NCEP       4/11    Removed AWIPS II database access        *
C* M. James/Unidata     1/11    Add support for both 4 & 8-bit NIDS     *
C* M. James/Unidata     8/13    Removed IM_DBGINI2GM, IM_DBNIDS2GM,     *
C*                              IM_DBMOSAIC2GM, IM_DBMCIDAS2GM          *
C* S. Guan/NCEP         06/15   Added code for NETCDF4 file             *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'IMGDEF.CMN'
	INCLUDE		'imcmn.cmn'
C*
	CHARACTER*(*)	cproj, imgfil
C*
	CHARACTER	temfil*132, tmproj*4
	LOGICAL		exist, found, negflg
	REAL		rnvblk ( LLNNAV )
C*
	INTEGER		iarr4 ( 8 ), ignhdr ( 135 )
	INTEGER*2	iarr2 ( 16 ), ipflip ( 16 )
C*
	EQUIVALENCE	( iarr4, iarr2 )
C------------------------------------------------------------------------
	iret   = 0
	imftyp = IFINVD
	found  = .false.
C
C*	Check projection type.
C
	CALL ST_LCUC ( cproj, tmproj, ier )
	IF  ( tmproj .eq. 'SAT' ) THEN
	    imbank = 1
	  ELSE IF ( tmproj .eq. 'RAD' ) THEN
	    imbank = 2
	  ELSE
	    iret = -5
	    RETURN
	END IF
C
C*	Get the actual file name to use.
C
	CALL FL_INQR ( imgfil, exist, temfil, ier )
	IF  ( .not. exist ) THEN
	    iret = -1
	    RETURN
	END IF
C
C*	Open the file and read enough information to detect the format.
C
	CALL FL_DOPN ( temfil, 8, .false., lunif, iret )
	IF  ( iret .ne. 0 ) THEN
	    iret = -1
	    RETURN
	END IF
	CALL FL_READ ( lunif, 1, 8, iarr4, iret )
	CALL FL_CLOS ( lunif, ierr )
C
C*	Verify we were able to read a header OK.
C
	IF ( iret .ne. 0 ) THEN
	    iret = -1
	    RETURN
	END IF
C
C*	Initialize calibration flag for color bar
C
	imcalbar = 0
C
C*	Check for AREA file.	Rules: Word 1 = 0, Word 2 = 4 (might be
C*				byte-swapped).
C
	ier = MV_SWP4 ( 1, iarr4 (2), iswap2 )
	IF  ( (   iarr4 (1) .eq. 0 ) .and.
     +	      ( ( iarr4 (2) .eq. 4 ) .or. ( iswap2 .eq. 4 ) ) ) THEN
C
	    imftyp = IFAREA
	    found  = .true.
C
C*	    Read header and set navigation for AREA files.
C
	    CALL IM_AR2GM ( temfil, rnvblk, iret )
	END IF
C
	IF  ( .not. found ) THEN
C
C*	    Check for GINI images.	Rules: Look at the WMO header.
C*					Char 1 is T (ascii 84), and
C*					there should be a LF in byte
C*					21 (or 25).
C
	    nbytes = 1
	    negflg = .false.
	    istart = 1 - 1
	    CALL MV_BTOI ( iarr4, istart, nbytes, negflg, ihdr1, ier )
	    istart = 21 - 1
	    CALL MV_BTOI ( iarr4, istart, nbytes, negflg, ihdr21, ier )
	    istart = 25 - 1
	    CALL MV_BTOI ( iarr4, istart, nbytes, negflg, ihdr25, ier )
	    IF  ( (   ihdr1  .eq. 84 ) .and.
     +		  ( ( ihdr21 .eq. 10 ) .or. ( ihdr25 .eq. 10 ) ) ) THEN
C
		imftyp = IFGINI
		found  = .true.
C
C*		Read header and set navigation for GINI files.
C
		CALL IM_GI2GM ( temfil, ignhdr, rnvblk, iret )
C
C*		Use the satellite color bank for these images (even radar!)
C
		imbank = 1
C
	    END IF
	END IF
C
C*	S. Chiswell, Updated for new products, up to 299
C
        IF  ( ( .not. found ) .and.
     +	      ( ( iarr2(1)   .eq. iarr2(16) ) .and.
     +		( ( iarr2(1) .ge.        16 ) .and.
     +		  ( iarr2(1) .le.       149 ) ) .and.
     +		( iarr2(10)  .eq.        -1 ) .and.
     +		( iarr2(7)   .lt.     10000 ) ) ) THEN
C
C*	    Check for raw NIDS format.	Rules: halfword 1 = halfword 16,
C*	    in the valid NIDS product range, halfword 10 = -1, and
C*	    halfword 7 less than 10000.
C
	    imftyp = IFNIDS
	    found  = .true.
C
C*	    Read header and set navigation for NIDS files.
C
            IF ( (iarr2(1) .eq. 94) .or. (iarr2(1) .eq. 99)
     +           .or. (iarr2(1) .eq. 134) .or. (iarr2(1) .eq. 135)
     +           .or. (iarr2(1) .eq. 138)) THEN
                   imftyp = IFHINIDS
                   CALL ST_NULL ( temfil, temfil, lens, iret )
                   CALL IM_NEXBZ ( temfil, itype, ioff, iret )
               ELSE
                   CALL IM_NIHD ( temfil, iret )
            END IF
	END IF
C
	CALL MV_SWP2 ( 8, iarr4 (1), ipflip )
        IF  ( ( .not. found ) .and.
     +	      ( ( ipflip(1)   .eq. ipflip(16) ) .and.
     +		( ( ipflip(1) .ge.         16 ) .and.
     +		  ( ipflip(1) .le.        149 ) ) .and.
     +		( ipflip(10)  .eq.         -1 ) .and.
     +		( ipflip(7)   .lt.      10000 ) ) ) THEN
C
C*	    Check for raw NIDS format with the bytes flipped.
C*	    Rules: halfword 1 = halfword 16, in the valid NIDS product
C*	    range, halfword 10 = -1, and halfword 7 less than 10000.
C
	    imftyp = IFNIDS
	    found  = .true.
C
C*	    Read header and set navigation for NIDS files.
C
            IF ( (ipflip(1) .eq. 94) .or. (ipflip(1) .eq. 99)
     +         .or. (ipflip(1) .eq. 134) .or. (ipflip(1) .eq. 135)
     +           .or. (iarr2(1) .eq. 138)) THEN
                   imftyp = IFHINIDS
                   CALL ST_NULL ( temfil, temfil, lens, iret )
                   CALL IM_NEXBZ ( temfil, itype, ioff, iret )
               ELSE
	           CALL IM_NIHD ( temfil, iret )
            END IF
	END IF
C
	IF  ( ( .not. found ) .and.
     +	      ( ( iarr2(1)   .eq. iarr2(16) ) .and.
     +		( ( iarr2(1) .ge.         5 ) .and.
     +		  ( iarr2(1) .le.       999 ) ) .and.
     +		( iarr2(10)  .eq.        -1 ) .and.
     +		( iarr2(7)   .ge.     10000 ) ) ) THEN
C
C*	    Check for UNISYS Mosaics (Alaska, Hawaii, National, and Regional).
C*	    Rules:  halfword 1 = halfword 16, and in the valid NIDS
C*		    product range, and halfword 10 = -1 and halfword 7
C*		    greater than or equal 10000 and first 2 bytes are
C*		    a valid product ID (see DATA statement).
C
	    imftyp = IFNIDS
	    found  = .true.
C
C*	    Read header and set navigation for UNISYS NIDS products.
C
	    CALL IM_UNHD ( temfil, iret )
	END IF
C
	CALL MV_SWP2 ( 8, iarr4 (1), ipflip )
	IF  ( ( .not. found ) .and.
     +	      ( ( ipflip(1)   .eq. ipflip(16) ) .and.
     +		( ( ipflip(1) .ge.          5 ) .and.
     +		  ( ipflip(1) .le.        999 ) ) .and.
     +		( ipflip(10)  .eq.         -1 ) .and.
     +		( ipflip(7)   .ge.      10000 ) ) ) THEN
C
C*	    Check for UNISYS Mosaics (Alaska, Hawaii, National, and Regional).
C*	    Rules:  halfword 1 = halfword 16, and in the valid NIDS
C*		    product range, and halfword 10 = -1 and halfword 7
C*		    greater than or equal 10000 and first 2 bytes are
C*		    a valid product ID (see DATA statement).
C
	    imftyp = IFNIDS
	    found  = .true.
C
C*	    Read header and set navigation for UNISYS NIDS products.
C
	    CALL IM_UNHD ( temfil, iret )
	END IF
C
	IF  ( .not. found ) THEN
C
C*	    Check for NOWrad images.	Rules: First three bytes are
C*					0 (0x00), 240 (0xf0), and
C*					9 (0x09) to indicate the
C*					broadcast image header.
	    nbytes = 1
	    negflg = .false.
	    istart = 1 - 1
	    CALL MV_BTOI ( iarr4, istart, nbytes, negflg, ihdr1, ier )
	    istart = 2 - 1
	    CALL MV_BTOI ( iarr4, istart, nbytes, negflg, ihdr2, ier )
	    istart = 3 - 1
	    CALL MV_BTOI ( iarr4, istart, nbytes, negflg, ihdr3, ier )
	    IF  ( ( ihdr1 .eq. 0 ) .and. ( ihdr2 .eq. 240 ) .and.
     +		  ( ihdr3 .eq. 9 ) ) THEN
C
		imftyp = IFNOWR
		found  = .true.
C
C*		Read header and set navigation for WSI NOWrad products.
C
		CALL IM_NOHD ( temfil, iret )
	    END IF
	END IF
C
	IF  ( .not. found )  THEN
C
C*	    Check for NetCDF images.	Rules: First four bytes are
C*					'C', 'D', 'F', 0x01.
C
	    CALL MV_SWP2 ( 1, iarr4(1), ipflip )
	    IF  ( ( (  iarr2(1) .eq. 17220 ) .and.
     +		    (  iarr2(2) .eq. 17921 ) ) .or.
     +		  ( ( ipflip(1) .eq. 17220 ) .and.
     +		    ( ipflip(2) .eq. 17921 ) ) )  THEN
C
		imftyp = IFNCDF
		found  = .true.
C
C*		Read header and set navigation for NetCDF products.
C
		CALL IM_NCDF ( temfil, iret )
		IF  ( iret .ne. 0 )  THEN
		    imftyp = IFINVD
		    found  = .false.
		ELSE
C
C*		    Use the satellite color bank for these images (even qpe radar!)
C
		    imbank = 1
		END IF
	    END IF
	END IF
C
        IF  ( .not. found )  THEN
C
C*          Check for NetCDF4 images. 211   H   D   F
C
            CALL MV_SWP2 ( 1, iarr4(1), ipflip )
            IF  ( ( (  iarr2(1) .eq. 18569 ) .and.
     +              (  iarr2(2) .eq. 17988 ) ) .or.
     +            ( ( ipflip(1) .eq. 18569 ) .and.
     +              ( ipflip(2) .eq. 17988 ) ) )  THEN
C
                imftyp = IFNCDF
                found  = .true.
C
C*              Read header and set navigation for NetCDF4 products.
C
                CALL IM_NCDF4 ( temfil, iret )
                IF  ( iret .ne. 0 )  THEN
                    imftyp = IFINVD
                    found  = .false.
                END IF
            END IF
        END IF
C
	IF  ( .not. found )  THEN
C
C*	    Check to see if we have a FOS header before the product
C*	    actually starts.
C*
C*	    First 4 bytes are SOH, cr, cr, nl
C*
C
	    IF  ( ( ( iarr2(1) .eq.  269 ) .and.
     +		    ( iarr2(2) .eq. 3338 ) ) .or.
     +		  ( ( iarr2(1) .eq. 3329 ) .and.
     +		    ( iarr2(2) .eq. 2573 ) ) )  THEN
C
		CALL ST_NULL ( temfil, temfil, lens, iret )
		CALL IM_NEXZ ( temfil, 4, itype, ioff, iret )
		IF  ( iret .eq. 0)  THEN
		    found = .true.
		    imftyp = itype
		    IF ( ( imftyp .eq. IFNIDS ) .or.
     +		    ( imftyp .eq. IFNEXZ ) ) THEN
		    	imdoff = imdoff + ioff
		    END IF
		END IF
C
	    END IF
	END IF
C
	IF  ( .not. found )  THEN
C
C*	    Check to see if we have a WMO header before the product
C*	    actually starts.
C*
C*	    First 4 bytes are 'S', 'D', 'U', 'S'
C*
C
	    IF  ( ( ( iarr2(1) .eq. 21316 ) .and.
     +		    ( iarr2(2) .eq. 21843 ) ) .or.
     +		  ( ( iarr2(1) .eq. 17491 ) .and.
     +		    ( iarr2(2) .eq. 21333 ) ) )  THEN
C
		CALL ST_NULL ( temfil, temfil, lens, iret )
		CALL IM_NEXZ ( temfil, 2, itype, ioff, iret )
		IF  ( iret .eq. 0)  THEN
		    found = .true.
		    imftyp = itype
		    IF (imftyp .eq. IFNIDS ) THEN
                        imdoff = imdoff + ioff
		    END IF
		END IF
C
	    END IF
	END IF
C
	IF  ( .not. found )  THEN
C
C*	    Check for NEXRAD Level II Archive format
C*
C*	    First 8 bytes are A   R   C   H   I   V   E   2
C*	    or
C*	    With Build5 ROC header are A   R   2   V  x   x   x   x
C
	    CALL MV_SWP2 ( 1, iarr4(1), ipflip )
            IF  ( ( (  iarr2(1) .eq. 16722 ) .and.
     +              (  iarr2(2) .eq. 17224 ) .and.
     +		    (  iarr2(3) .eq. 18774 ) .and.
     +		    (  iarr2(4) .eq. 17714 ) ) .or.
     +            ( ( ipflip(1) .eq. 16722 ) .and.
     +              ( ipflip(2) .eq. 17224 ) .and.
     +		    ( ipflip(3) .eq. 18774 ) .and.
     +              ( ipflip(4) .eq. 17714 ) ) .or.
     +		  ( (  iarr2(1) .eq. 16722 ) .and.
     +		    (  iarr2(2) .eq. 12886 ) ) .or.
     +            ( ( ipflip(1) .eq. 16722 ) .and.
     +              ( ipflip(2) .eq. 12886 ) ) )  THEN

		CALL IM_NEX2 ( temfil, iret )
		IF ( iret .eq. 0 ) THEN
		    found = .true.
		    imftyp = IFNEX2
C
C*		    Use the satellite color bank for these images
C
		    imbank = 1
		END IF
	    END IF
	END IF
C
C*	Else, don't know this format.
C
	IF  ( .not. found )  THEN
	    iret = -3
	  ELSE
C
C*	    Read the image table.
C
	    CALL IM_RTBL ( iret )
C
C*	    NIDS non graphic products OK if not in table
C
	    IF  ( ( iret .ne. 0 ) .and.
     +		( ( imradf .eq. 1 ) .and. ( imrdfl .eq. 2 ) ) ) THEN
		iret = 0
	    END IF
	END IF
C
C*	If successfully down to here, set CMFILE.
C
	cmfile = temfil
C
C*      Use sat color bank if imndlv exceeds 20 for radar
C
        IF ( ( ( imbank .eq. 2 ) .and. ( imndlv .gt. 20 ) ) .or.
     +	 ( imftyp .eq. IFHINIDS ) ) THEN
           imbank = 1
        END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program main
  use mpi
  implicit none
  integer, parameter :: n = 8, nb = 1
  integer :: provided, ierr, myrank, nprocs
  integer :: icontxt, npcol, nprow
  integer :: info
  double precision, allocatable :: d(:), e(:), w(:)
  integer :: lwork, liwork
  double precision, allocatable :: work(:)
  integer, allocatable :: iwork(:)
  integer :: il, iu, m, nsplit
  double precision :: vl, vu, abstol
  integer, allocatable :: isplit(:), iblock(:)

  call MPI_init_thread(MPI_THREAD_MULTIPLE, provided, ierr)
  call MPI_comm_rank(MPI_COMM_WORLD, myrank, ierr)
  call MPI_comm_size(MPI_COMM_WORLD, nprocs, ierr)
  nprow = int(sqrt(1. * nprocs + 0.5))
  npcol = int(1. * nprocs / nprow)
  if (myrank == 0) then
     print *, "n =", n
     print *, "nprocs =", nprocs
     print *, "nprow =", nprow
     print *, "npcol =", npcol
     if ((nprocs /= nprow * npcol) .or. (mod(n, nprow) /= 0) .or. &
          (mod(n, npcol) /= 0)) then
        print *, "incompatible matrix size and number of processes"
        call MPI_Abort(MPI_COMM_WORLD, 127, ierr)
     end if
  end if
  call BLACS_get(0, 0, icontxt)
  call BLACS_gridinit(icontxt, 'R', nprow, npcol)  ! not used explicitly, but needed for pdstebz. If you remove it, segmentation fault occurs.

  allocate(d(n), e(n))
  call generate_laplacian(n, d, e)

  allocate(w(n))

  abstol = 0d0
  lwork = -1
  liwork = -1
  allocate(work(1), iwork(1))
  allocate(isplit(n))
  allocate(iblock(n))

  call pdstebz(icontxt, 'A', 'E', n, &
       & vl, vu, il, iu, &
       & abstol, d, e, m, nsplit, &
       & w, iblock, isplit, &
       & work, lwork, iwork, liwork, &
       & info)
  lwork = int(work(1))
  liwork = iwork(1)
  deallocate(work, iwork)
  allocate(work(lwork), iwork(liwork))

  call pdstebz(icontxt, 'A', 'E', n, &
       & vl, vu, il, iu, &
       & abstol, d, e, m, nsplit, &
       & w, iblock, isplit, &
       & work, lwork, iwork, liwork, &
       & info)

  call mpi_barrier(mpi_comm_world, ierr)
  if (myrank == 0) then
     print *, "eigenvalues:", w(:)
  end if

  deallocate(isplit, iblock)
  deallocate(d, e, w, work, iwork)
  call BLACS_gridexit(icontxt)
  call MPI_finalize(ierr)

contains

  subroutine generate_laplacian(n, d, e)
    implicit none
    integer, intent(in) :: n
    double precision, intent(out) :: d(:), e(:)
    integer :: i

    d(1) = 1d0;
    do i = 2, n
       d(i) = 2d0
    enddo
    do i = 1, n-1
       e(i) = -1d0
    enddo
  end subroutine generate_laplacian

end program main

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
      recursive subroutine fpinst(iopt,t,n,c,k,x,l,tt,nn,cc,nest)
      implicit none
c  given the b-spline representation (knots t(j),j=1,2,...,n, b-spline
c  coefficients c(j),j=1,2,...,n-k-1) of a spline of degree k, fpinst
c  calculates the b-spline representation (knots tt(j),j=1,2,...,nn,
c  b-spline coefficients cc(j),j=1,2,...,nn-k-1) of the same spline if
c  an additional knot is inserted at the point x situated in the inter-
c  val t(l)<=x<t(l+1). iopt denotes whether (iopt.ne.0) or not (iopt=0)
c  the given spline is periodic. in case of a periodic spline at least
c  one of the following conditions must be fulfilled: l>2*k or l<n-2*k.
c
c  ..scalar arguments..
      integer k,n,l,nn,iopt,nest
      real*8 x
c  ..array arguments..
      real*8 t(nest),c(nest),tt(nest),cc(nest)
c  ..local scalars..
      real*8 fac,per,one
      integer i,i1,j,k1,m,mk,nk,nk1,nl,ll
c  ..
      one = 0.1e+01
      k1 = k+1
      nk1 = n-k1
c  the new knots
      ll = l+1
      i = n
      do 10 j=ll,n
         tt(i+1) = t(i)
         i = i-1
  10  continue
      tt(ll) = x
      do 20 j=1,l
         tt(j) = t(j)
  20  continue
c  the new b-spline coefficients
      i = nk1
      do 30 j=l,nk1
         cc(i+1) = c(i)
         i = i-1
  30  continue
      i = l
      do 40 j=1,k
         m = i+k1
         fac = (x-tt(i))/(tt(m)-tt(i))
         i1 = i-1
         cc(i) = fac*c(i)+(one-fac)*c(i1)
         i = i1
  40  continue
      do 50 j=1,i
         cc(j) = c(j)
  50  continue
      nn = n+1
      if(iopt.eq.0) return
c   incorporate the boundary conditions for a periodic spline.
      nk = nn-k
      nl = nk-k1
      per = tt(nk)-tt(k1)
      i = k1
      j = nk
      if(ll.le.nl) go to 70
      do 60 m=1,k
         mk = m+nl
         cc(m) = cc(mk)
         i = i-1
         j = j-1
         tt(i) = tt(j)-per
  60  continue
      return
  70  if(ll.gt.(k1+k)) return
      do 80 m=1,k
         mk = m+nl
         cc(mk) = cc(m)
         i = i+1
         j = j+1
         tt(j) = tt(i)+per
  80  continue
      return
      end

max_stars_repo_name: jcampbell/scipy

max_issues_repo_name: jcampbell/scipy

max_forks_repo_name: jcampbell/scipy


! { dg-do run }
!
! PR fortran/37336
!
! Check the scalarizer/array packing with strides
! in the finalization wrapper
!
module m
  implicit none

  type t1
    integer :: i = 1
  contains
    final :: fini_elem
  end type t1

  type, extends(t1) :: t1e
    integer :: j = 11
  contains
    final :: fini_elem2
  end type t1e

  type t2
    integer :: i = 2
  contains
    final :: fini_shape
  end type t2

  type, extends(t2) :: t2e
    integer :: j = 22
  contains
    final :: fini_shape2
  end type t2e

  type t3
    integer :: i = 3
  contains
    final :: fini_explicit
  end type t3

  type, extends(t3) :: t3e
    integer :: j = 33
  contains
    final :: fini_explicit2
  end type t3e

  integer :: cnt1, cnt1e, cnt2, cnt2e, cnt3, cnt3e

contains

  impure elemental subroutine fini_elem(x)
    type(t1), intent(inout) :: x
    integer :: i, j, i2, j2

    if (cnt1e /= 5*4) STOP 1
    j = mod (cnt1,5)+1
    i = cnt1/5 + 1
    i2 = (i-1)*3 + 1
    j2 = (j-1)*2 + 1
    if (x%i /= j2 + 100*i2) STOP 2
    x%i = x%i * (-13)
    cnt1 = cnt1 + 1
  end subroutine fini_elem

  impure elemental subroutine fini_elem2(x)
    type(t1e), intent(inout) :: x
    integer :: i, j, i2, j2

    j = mod (cnt1e,5)+1
    i = cnt1e/5 + 1
    i2 = (i-1)*3 + 1
    j2 = (j-1)*2 + 1
    if (x%i /= j2 + 100*i2) STOP 3
    if (x%j /= (j2 + 100*i2)*100) STOP 4
    x%j = x%j * (-13)
    cnt1e = cnt1e + 1
  end subroutine fini_elem2

  subroutine fini_shape(x)
    type(t2) :: x(:,:)
    if (cnt2e /= 1 .or. cnt2 /= 0) STOP 5
    call check_var_sec(x%i, 1)
    x%i = x%i * (-13)
    cnt2 = cnt2 + 1
  end subroutine fini_shape

  subroutine fini_shape2(x)
    type(t2e) :: x(:,:)
    call check_var_sec(x%i, 1)
    call check_var_sec(x%j, 100)
    x%j = x%j * (-13)
    cnt2e = cnt2e + 1
  end subroutine fini_shape2

  subroutine fini_explicit(x)
    type(t3) :: x(5,4)
    if (cnt3e /= 1 .or. cnt3 /= 0) STOP 6
    call check_var_sec(x%i, 1)
    x%i = x%i * (-13)
    cnt3 = cnt3 + 1
  end subroutine fini_explicit

  subroutine fini_explicit2(x)
    type(t3e) :: x(5,4)
    call check_var_sec(x%i, 1)
    call check_var_sec(x%j, 100)
    x%j = x%j * (-13)
    cnt3e = cnt3e + 1
  end subroutine fini_explicit2

  subroutine fin_test_1(x)
    class(t1), intent(out) :: x(5,4)
  end subroutine fin_test_1

  subroutine fin_test_2(x)
    class(t2), intent(out) :: x(:,:)
  end subroutine fin_test_2

  subroutine fin_test_3(x)
    class(t3), intent(out) :: x(:,:)
    if (any (shape(x) /= [5,4])) STOP 7
  end subroutine fin_test_3

  subroutine check_var_sec(x, factor)
    integer :: x(:,:)
    integer, value :: factor
    integer :: i, j, i2, j2

    do i = 1, 4
      i2 = (i-1)*3 + 1
      do j = 1, 5
        j2 = (j-1)*2 + 1
        if (x(j,i) /= (j2 + 100*i2)*factor) STOP 8
      end do
    end do
  end subroutine check_var_sec
end module m

program test
  use m
  implicit none

  class(t1), allocatable :: x(:,:)
  class(t2), allocatable :: y(:,:)
  class(t3), allocatable :: z(:,:)
  integer :: i, j

  cnt1 = 0; cnt1e = 0; cnt2 = 0; cnt2e = 0;  cnt3 = 0; cnt3e = 0

  allocate (t1e :: x(10,10))
  allocate (t2e :: y(10,10))
  allocate (t3e :: z(10,10))

  select type(x)
    type is (t1e)
      do i = 1, 10
        do j = 1, 10
          x(j,i)%i = j + 100*i
          x(j,i)%j = (j + 100*i)*100
        end do
      end do
  end select

  select type(y)
    type is (t2e)
      do i = 1, 10
        do j = 1, 10
          y(j,i)%i = j + 100*i
          y(j,i)%j = (j + 100*i)*100
        end do
      end do
  end select

  select type(z)
    type is (t3e)
      do i = 1, 10
        do j = 1, 10
          z(j,i)%i = j + 100*i
          z(j,i)%j = (j + 100*i)*100
        end do
      end do
  end select

  if (cnt1 + cnt1e + cnt2 + cnt2e + cnt3 + cnt3e /= 0) STOP 9

  call fin_test_1(x(::2,::3))
  if (cnt1 /= 5*4) STOP 10
  if (cnt1e /= 5*4) STOP 11
  cnt1 = 0; cnt1e = 0
  if (cnt2 + cnt2e + cnt3 + cnt3e /= 0) STOP 12

  call fin_test_2(y(::2,::3))
  if (cnt2 /= 1) STOP 13
  if (cnt2e /= 1) STOP 14
  cnt2 = 0; cnt2e = 0
  if (cnt1 + cnt1e + cnt3 + cnt3e /= 0) STOP 15

  call fin_test_3(z(::2,::3))
  if (cnt3 /= 1) STOP 16
  if (cnt3e /= 1) STOP 17
  cnt3 = 0; cnt3e = 0
  if (cnt1 + cnt1e + cnt2 + cnt2e /= 0) STOP 18

  select type(x)
    type is (t1e)
      call check_val(x%i, 1, 1)
      call check_val(x%j, 100, 11)
  end select

  select type(y)
    type is (t2e)
      call check_val(y%i, 1, 2)
      call check_val(y%j, 100, 22)
  end select

  select type(z)
    type is (t3e)
      call check_val(z%i, 1, 3)
      call check_val(z%j, 100, 33)
  end select

contains
  subroutine check_val(x, factor, val)
    integer :: x(:,:)
    integer, value :: factor, val
    integer :: i, j
    do i = 1, 10
      do j = 1, 10
        if (mod (j-1, 2) == 0 .and. mod (i-1, 3) == 0) then
          if (x(j,i) /= val) STOP 19
        else
          if (x(j,i) /= (j + 100*i)*factor) STOP 20
        end if
      end do
    end do
  end subroutine check_val
end program test


program omp_par_do
  implicit none

  integer, parameter :: n = 100
  real, dimension(n) :: dat, result
  integer :: i

  !$OMP PARALLEL DO
  do i = 1, n
     result(i) = my_function(dat(i))
  end do
  !$OMP END PARALLEL DO

contains

  function my_function(d) result(y)
    real, intent(in) :: d
    real :: y

    ! do something complex with data to calculate y
  end function my_function

end program omp_par_do


text:
C    @(#)fnd_ptin.f	20.1 7/18/96
C****************************************************************
C
C     File: fnd_ptin.f
C
C     Purpose: Routine to determin PTI number hash index given PTI numbe
C
c     Return code:  n > 0 : PTI number hash index
c                   n < 0 : PTI number hash index (duplicate entity)
c                   n = 0 : Error or termination flag
c
C     Author: Walt Powell  Date: 21 May 1996
C     Called by: load_pti.f
C
C****************************************************************
      integer function fnd_ptin (numpti)
      integer numpti

      include 'ipfinc/parametr.inc'

      include 'ipfinc/pti_data.inc'

      integer        num, h, p
      logical        found

      num = iabs (numpti)
      if (num .gt. 0) then
        h = mod (num, PTI_HASHSIZE) + 1
        p = htable_n(h)
        found = .false.
        do while (p .gt. 0 .and. .not. found)
          if (num .ne. pti_num(p)) then
            p = nextptr_n(p)
          else
            found = .true.
          endif
        enddo
        fnd_ptin = p
      else
        fnd_ptin = 0
      endif
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
! Copyright (C) 2017 A. Gil, G. Navas-Palencia, J. Segura and N. M. Temme

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

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

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

program test
    use gnstlib
    implicit none

    ! parameters

    integer :: err_id, i
    real(kind=dp), parameter :: x = 3.2_dp
    real(kind=dp), parameter :: xsmall = 3.2e-9_dp
    complex(kind=dp), parameter :: z = complex(3.0_dp, 4.0_dp)

    real(kind=dp) :: v(6)
    real(kind=dp), allocatable :: r(:)

    v(1:6) = (/3.2_dp, 6.3_dp, 0.4_dp, 10.5_dp, 100.6_dp, -2.7_dp/)

    err_id = 0

    write(*,*) "real functions"
    write(*,*) "x          =", x
    write(*,*) "exp(x)     =", gnstlib_exp(x)
    write(*,*) "exp2(x)    =", gnstlib_exp2(x)
    write(*,*) "exp10(x)   =", gnstlib_exp10(x)
    write(*,*) "expm1(x)   =", gnstlib_expm1(x)
    write(*,*) "expm1x(x)  =", gnstlib_expm1x(x)
    write(*,*) "expm1mx(x) =", gnstlib_expm1mx(x)
    write(*,*)
    write(*,*) "log(x)     =", gnstlib_log(x)
    write(*,*) "log2(x)    =", gnstlib_log2(x)
    write(*,*) "log10(x)   =", gnstlib_log10(x)
    write(*,*) "log1p(x)   =", gnstlib_log1p(xsmall)
    write(*,*) "log1pmx(x) =", gnstlib_log1pmx(x)
    write(*,*) "log1mexp(x)=", gnstlib_log1mexp(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "log1pexp(x)=", gnstlib_log1pexp(x)

    write(*,*)
    write(*,*) "complex functions"
    write(*,*) "z          =", z
    write(*,*) "exp(z)     =", gnstlib_cexp(z)
    write(*,*) "log(z)     =", gnstlib_clog(z)

    write(*,*)
    write(*,*) "power functions"
    err_id = 0
    write(*,*) "powm1      =", gnstlib_powm1(-3.0_dp, -1.2_dp, err_id)
    write(*,*) "err_id     =", err_id

    write(*,*) "gamma functions"
    err_id = 0
    write(*,*) "gamma(x)      =", gnstlib_cgamma(z, err_id)
    write(*,*) "err_id        =", err_id
    write(*,*) "stirling(x)   =", gnstlib_stirling(11.5_dp, err_id)
    write(*,*) "err_id        =", err_id
    write(*,*) "gammastar(x)  =", gnstlib_gammastar(3.5_dp, err_id)
    write(*,*) "err_id        =", err_id
    write(*,*) "qgamma(x, y)  =", gnstlib_qgamma(28.5_dp, 36.2_dp, err_id)
    write(*,*) "err_id        =", err_id

    write(*,*)
    write(*,*) "vectorized gamma"
    call gnstlib_gamma_vec(6, v, r, 0)

    do i = 1, 6
        write(*,*) "gamma(",v(i),")=", r(i)
    end do

    write(*,*)
    write(*,*) "ei(x) = ", gnstlib_ei(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "e1(x) = ", gnstlib_e1(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "li(x) = ", gnstlib_li(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "ci(x) = ", gnstlib_ci(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "si(x) = ", gnstlib_si(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "chi(x) = ", gnstlib_chi(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "shi(x) = ", gnstlib_shi(x, err_id)
    write(*,*) "err_id     =", err_id

    write(*,*) "invei(x) = ", gnstlib_invei(5.4_dp, err_id)
    write(*,*) "err_id     =", err_id

    write(*,*)
    write(*,*) "erf(x)     =", gnstlib_erf(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "erfc(x)    =", gnstlib_erfc(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "erfcx(x)   =", gnstlib_erfcx(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "erfi(x)    =", gnstlib_erfi(x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "voigt(x)   =", gnstlib_voigt_profile(x, x, x, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "inverf(x)  =", gnstlib_inverf(0.2_dp, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "inverfc(x) =", gnstlib_inverfc(0.2_dp, err_id)
    write(*,*) "err_id     =", err_id

    write(*,*)
    write(*,*) "erf(z)     =", gnstlib_cerf(z, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "erfc(z)    =", gnstlib_cerfc(z, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "erfcx(z)   =", gnstlib_cerfcx(z, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "erfi(z)    =", gnstlib_cerfi(z, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "faddeeva(z)=", gnstlib_faddeeva(z, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "dawson(z)  =", gnstlib_cdawson(z, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "fresnelc(z)=", gnstlib_cfresnelc(z, err_id)
    write(*,*) "err_id     =", err_id
    write(*,*) "fresnels(z)=", gnstlib_cfresnels(z, err_id)
    write(*,*) "err_id     =", err_id

    write(*,*)
    write(*,*) "gammainc_p(a, x) = ", gnstlib_gammainc_p(2.4_dp, 1.3_dp, err_id)
    write(*,*) "gammainc_q(a, x) = ", gnstlib_gammainc_q(2.4_dp, 1.3_dp, err_id)
    write(*,*) "expint(v, x)     = ", gnstlib_expint(2.4_dp, 1.3_dp, err_id)
    write(*,*) "expint_acc(v, x) = ", gnstlib_expint_acc(2, 1e-14_dp, 1e-10_dp, err_id)

end program

max_stars_repo_name: guillermo-navas-palencia/gnstlib

max_issues_repo_name: guillermo-navas-palencia/gnstlib

max_forks_repo_name: guillermo-navas-palencia/gnstlib


text:
OPEN(ACCESS='SEQUENTIAL', DECIMAL='POINT', BLANK="ZERO", UNIT=1)
OPEN(STATUS='new', UNIT=2, FILE='hunkydory.txt')
END PROGRAM

max_stars_repo_name: ouankou/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


submodule (io) milestone

use timeutils, only : date_filename,dateinc
use h5fortran, only : h5exist
use pathlib, only : get_suffix

implicit none (type,external)   !! external procedures must be explicitly denoted thusly

contains

module procedure find_milestone

!> search path having output rate cadence (s) and find the last file that is a milestone.
integer, dimension(3) :: ymd
real(wp) :: UTsec
character(:), allocatable :: fn, suffix
logical :: exists
real(wp) :: tsim

tsim = 0
tmile = 0

ymd = cfg%ymd0
UTsec = cfg%UTsec0
ymdmile = cfg%ymd0
UTsecmile = cfg%UTsec0

suffix = get_suffix(cfg%indatsize)
filemile = date_filename(cfg%outdir, ymd, UTsec) // suffix
!! This presumes the first file output is a milestone.
!! We don't test the situation wheere a first output was not produced.
!! User should not be restarting in that case.

if (cfg%mcadence <= 0 .and. cfg%flagoutput/=1) then      !okay for milestone if full output specified
!! milestone was not in config.nml
  inquire(file=filemile, exist=exists)
 ! error stop filemile
  if (exists) error stop 'a fresh simulation should not have data in output directory: ' // filemile
  return
endif

if (suffix /= '.h5') return

milesearch : do
  !! new filename, add the 1 if it is the first
  fn = date_filename(cfg%outdir, ymd, UTsec) // suffix

  inquire(file=fn, exist=exists)
  if ( .not. exists ) exit milesearch
  !! last output file

  if (h5exist(fn, '/nsall')) then
    !! this file is milestone
    ymdmile=ymd
    UTsecmile=UTsec
    filemile=fn
    tmile=tsim
  end if

  !! next time
  call dateinc(cfg%dtout, ymd,UTsec)
  tsim = tsim + cfg%dtout
end do milesearch

end procedure find_milestone

end submodule milestone


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20090924 22:38:17 nbsp  Welcome to the Wiki. You can establish identity here in a variety of ways. Doing so will give more weight to your edit on the Café Méditerranée page. Youve already chosen your Importance of using your RealName username, so the next step would be to introduce yourself and make some positive contributions like adding business reviews, Photo Requests uploading photos, or adding to Seed Pages List pages that need more information. Is there a particular reason that you removed Users/Aaron.Curtin Aarons comment? Users/JasonAller

20100510 19:23:33 nbsp This is possibly an employee of Cafe Mediterranee, since she removed two pieces of information from that page that reflected poorly on the establishment. Users/DavisKing

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! Variables that all files may need access too
! For instance the outfile unit for writing to the main file
!

module global_parameters_M

    implicit none

!---Define floats
    integer, parameter :: dp=kind(0.d0) !---double precision
    integer, parameter :: sp=kind(1.0)  !--- single precision

!---Nonlinear variables
    integer :: max_iter = 1 ! max num of nonlinear iterations to do
    integer :: max_nl_iter  ! numer of nonllinear iterations to do
    real(dp) :: tolerance = 0.001 ! prescribed tolerance
!---
    integer :: NumTimeSteps

!---Math constants
    real(dp) :: pi
    parameter (pi = 3.1415926535897932_dp)
!---THIS ASSUMES 1 CM MESH SIZING
    real(dp), dimension(3) :: elem_vol_int_fe
    data elem_vol_int_fe / 0.16666667198838597,0.66666665602322805, &
                          0.16666667198838597 /

    integer :: beta_counter

end module global_parameters_M


text:
MODULE mod_blas

  IMPLICIT NONE

  INTERFACE

     ! BLAS Level 1 function
     ! Double precision Euclidean norm
     DOUBLE PRECISION FUNCTION DNRM2( n, X, incx )
       INTEGER n, incx
       DOUBLE PRECISION X(*)
     END FUNCTION DNRM2

     ! BLAS Level 2 subroutine
     ! Double precision matrix-vector multiply C = alpha * A x + beta * y
     SUBROUTINE DGEMV( trans, m, n, &
                       alpha, A, lda, &
                              x, incx, &
                       beta,  y, incy )
       CHARACTER*1 trans
       INTEGER m, n, lda, incx, incy
       DOUBLE PRECISION alpha, beta
       DOUBLE PRECISION A(lda,*), X(*), Y(*)
     END SUBROUTINE DGEMV

  END INTERFACE

END MODULE mod_blas

max_stars_repo_name: wyphan/testcodes

max_issues_repo_name: wyphan/testcodes

max_forks_repo_name: wyphan/testcodes


! { dg-do compile { target i?86-*-* x86_64-*-* } }
! { dg-options "-O3 -mavx -mvzeroupper -mtune=generic -dp" }

  SUBROUTINE func(kts, kte, qrz, qiz, rho)
  IMPLICIT NONE
  INTEGER, INTENT(IN)               :: kts, kte
  REAL,    DIMENSION(kts:kte), INTENT(INOUT) :: qrz, qiz, rho
  INTEGER                              :: k
  REAL, DIMENSION(kts:kte)    ::  praci, vtiold
  REAL                          :: fluxout
  INTEGER                       :: min_q, max_q, var
  do k=kts,kte
    praci(k)=1.0
  enddo
  min_q=kte
  max_q=kts-1
  DO var=1,20
    do k=max_q,min_q,-1
       fluxout=rho(k)*qrz(k)
    enddo
    qrz(min_q-1)=qrz(min_q-1)+fluxout
  ENDDO
  DO var=1,20
    do k=kts,kte-1
      vtiold(k)= (rho(k))**0.16
    enddo
  ENDDO
  STOP
  END SUBROUTINE func

! { dg-final { scan-assembler "avx_vzeroupper" } }


text:
      SUBROUTINE DT_SIGGA(Nti,Xi,Q2i,Ecmi,Xnui,Stot,Etot,Sin,Ein,Stot0)

C***********************************************************************
C Total/inelastic photon-nucleus cross sections.                       *
C     !!!! Overwrites SHMAKI-initialization. Do not use it during      *
C          production runs !!!!                                        *
C This version dated 27.03.96 is written by S. Roesler                 *
C***********************************************************************

      IMPLICIT NONE
      DOUBLE PRECISION AMPROT , DLARGE , ecm , Ecmi , Ein , Etot ,
     &                 OHALF , ONE , Q2i , Sin , Stot , Stot0 , TINY10 ,
     &                 TINY2 , x , Xi , xnu , Xnui , ZERO
      INTEGER nt , Nti
      SAVE

      INCLUDE 'inc/dtflka'

      PARAMETER (TINY10=1.0D-10,TINY2=1.0D-2,ZERO=0.0D0,DLARGE=1.0D10,
     &           OHALF=0.5D0,ONE=1.0D0)
      PARAMETER (AMPROT=0.938D0)

C emulsion treatment
      INCLUDE 'inc/dtcomp'

C Glauber formalism: cross sections
      INCLUDE 'inc/dtglxs'

      nt = Nti
      x = Xi
      Q2 = Q2i
      ecm = Ecmi
      xnu = Xnui
      IF ( (Ecmi.LE.ZERO) .AND. (Xnui.GT.ZERO) )
     &     ecm = SQRT(AMPROT**2-Q2+2.0D0*Xnui*AMPROT)
      CALL DT_XSGLAU(1,nt,7,x,Q2,ecm,1,1,-1)
      Stot = XSTot(1,1,1)
      Etot = XETot(1,1,1)
      Sin = XSPro(1,1,1)
      Ein = XEPro(1,1,1)

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
       LOGICAL FUNCTION IW3MAT(L1, L2, N)
C$$$   SUBPROGRAM  DOCUMENTATION  BLOCK
C
C FUNCTION: IW3MAT           TEST FOR MATCH TWO INTEGER ARRAYS
C   AUTHOR: STACKPOLE. J.D.  ORG: W342       DATE: 86-01-13
C
C ABSTACT: TEST N WORDS STARTING AT L1, L2 FOR EQUALITY, RETURN .TRUE.
C     IF ALL EQUAL; OTHERWISE .FALSE.
C
C PROGRAM HISTORY LOG:
C   86-01-13  J.D.STACKPOLE
C   90-03-15  R.E.JONES       CONVERT TO CRAY CFT77 FORTRAN
C
C USAGE:  II = IW3MAT(L1,L2,N)
C
C   INPUT VARIABLES:
C     NAMES  INTERFACE DESCRIPTION OF VARIABLES AND TYPES
C     ------ --------- -----------------------------------------------
C     L1     ARG LIST  INTEGER ARRAY TO MATCH WITH L2
C     L2     ARG LIST  INTEGER ARRAY TO MATCH WITH L1
C     N      ARG LIST  NUMBER OF INTEGER WORDS TO TEST FOR MATCH
C
C   OUTPUT VARIABLES:
C     NAMES  INTERFACE DESCRIPTION OF VARIABLES AND TYPES
C     ------ --------- -----------------------------------------------
C     IW3MAT FUNCTION  LOGICAL .TRUE. IF L1 AND L2 MATCH ON ALL WORDS,
C                      LOGICAL .FALSE. IF NOT MATCH ON ANY WORD
C
C ATTRIBUTES:
C   LANGUAGE: CRAY CFT77 FORTRAN
C   MACHINE:  CRAY Y-MP8/864, CRAY Y-MP EL@/256
C
C$$$
C
       INTEGER L1(*)
       INTEGER L2(*)
C
          IW3MAT = .TRUE.
          DO 10 I = 1,N
            IF (L1(I).NE.L2(I))  GO TO 20
   10     CONTINUE
          RETURN
C
   20     CONTINUE
          IW3MAT = .FALSE.
          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:
C
C
C +-----------------------------------------------------------------+
C |                                                                 |
C |                Copyright (C) 1986 by UCAR                       |
C |        University Corporation for Atmospheric Research          |
C |                    All Rights Reserved                          |
C |                                                                 |
C |                 NCARGRAPHICS  Version 1.00                      |
C |                                                                 |
C +-----------------------------------------------------------------+
C
C
C ---------------------------------------------------------------------
C
      SUBROUTINE AGSETI (TPID,IUSR)
C
      CHARACTER*(*) TPID
      DIMENSION FURA(1)
C
C The routine AGSETI may be used to set the integer-equivalent value of
C any single AUTOGRAPH control parameter.
C
      FURA(1)=FLOAT(IUSR)
      CALL AGSETP (TPID,FURA,1)
      RETURN
C
      END

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
C	THIS FUNCTION IS TO CALCULATE THE NUCLEAR PROFILE FUNCTION
C       OF THE  COLLIDERING SYSTEM (IN UNITS OF 1/mb)
	FUNCTION FLAP(X)
        REAL HIPR1, HINT1
        COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
        REAL BB, B1, PHI, Z1
        COMMON/PACT_HIJING/BB,B1,PHI,Z1
        EXTERNAL FGP1
        B1=X
        FLAP=GAUSS2(FGP1,0.0,2.0*HIPR1(40),0.01)
	RETURN
	END
C
        FUNCTION FLAP1(X)
        REAL BB, B1, PHI, Z1
        COMMON/PACT_HIJING/BB,B1,PHI,Z1
        R=SQRT(BB**2+X**2)
        FLAP1=WDSAX1(R)
        RETURN
        END
C
        FUNCTION FLAP2(X)
        REAL BB, B1, PHI, Z1
        COMMON/PACT_HIJING/BB,B1,PHI,Z1
        R=SQRT(BB**2+X**2)
        FLAP2=WDSAX2(R)
        RETURN
        END
C
	FUNCTION FGP1(X)
        REAL HIPR1, HINT1
        COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
        REAL BB, B1, PHI, Z1
        COMMON/PACT_HIJING/BB,B1,PHI,Z1
        EXTERNAL FGP2
        PHI=X
        FGP1=2.0*GAUSS3(FGP2,0.0,HIPR1(34),0.01)
	RETURN
	END
C
        FUNCTION FGP2(X)
        REAL HIPR1, HINT1
        COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
        REAL BB, B1, PHI, Z1
        COMMON/PACT_HIJING/BB,B1,PHI,Z1
        EXTERNAL FGP3
        Z1=X
        FGP2=2.0*GAUSS4(FGP3,0.0,HIPR1(35),0.01)
	RETURN
	END
C
        FUNCTION FGP3(X)
        REAL BB, B1, PHI, Z1
        COMMON/PACT_HIJING/BB,B1,PHI,Z1
        R1=SQRT(B1**2+Z1**2)
        R2=SQRT(BB**2+B1**2-2.0*B1*BB*COS(PHI)+X**2)
        FGP3=B1*WDSAX1(R1)*WDSAX2(R2)
        RETURN
        END
C
	FUNCTION WDSAX1(X)
        REAL HIPR1, HINT1
        COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
      	WDSAX1=HINT1(75)*(1.+HINT1(74)*(X/HINT1(72))**2)/
     &       (1+EXP((X-HINT1(72))/HINT1(73)))
       	IF (HINT1(74).LT.0.) THEN
       		IF (X.GE.HINT1(72)/SQRT(ABS(HINT1(74)))) WDSAX1=0.
       	ENDIF
      	RETURN
      	END
C
	FUNCTION WDSAX2(X)
        REAL HIPR1, HINT1
        COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
      	WDSAX2=HINT1(79)*(1.+HINT1(78)*(X/HINT1(76))**2)/
     &       (1+EXP((X-HINT1(76))/HINT1(77)))
       	IF (HINT1(78).LT.0.) THEN
       		IF (X.GE.HINT1(76)/SQRT(ABS(HINT1(78)))) WDSAX2=0.
       	ENDIF
      	RETURN
      	END
C
	FUNCTION  PROFILE(XB)
        REAL HIPR1, HINT1
        COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
        REAL BB, B1, PHI, Z1
        COMMON/PACT_HIJING/BB,B1,PHI,Z1
	SAVE
	EXTERNAL FLAP, FLAP1, FLAP2
C
        BB=XB
        PROFILE=1.0
        IF(IHNT2(1).GT.1 .AND. IHNT2(3).GT.1) THEN
           PROFILE=float(IHNT2(1))*float(IHNT2(3))*0.1*
     &          GAUSS1(FLAP,0.0,HIPR1(34),0.01)
        ELSE IF(IHNT2(1).EQ.1 .AND. IHNT2(3).GT.1) THEN
           PROFILE=0.2*float(IHNT2(3))*
     &          GAUSS1(FLAP2,0.0,HIPR1(35),0.001)
        ELSE IF(IHNT2(1).GT.1 .AND. IHNT2(3).EQ.1) THEN
           PROFILE=0.2*float(IHNT2(1))*
     &          GAUSS1(FLAP1,0.0,HIPR1(34),0.001)
        ENDIF
	RETURN
	END

max_stars_repo_name: AllaMaevskaya/AliRoot

max_issues_repo_name: AllaMaevskaya/AliRoot

max_forks_repo_name: AllaMaevskaya/AliRoot


text:
! { dg-do run }
IF (T1(1.0,1.0) .NE. (1.0,1.0) ) CALL ABORT()
IF (T1(1.0) .NE. (1.0,0.0)) CALL ABORT()
IF (M1(1,2,3) .NE. 3) CALL ABORT()
IF (M1(1,2,A4=4) .NE. 4) CALL ABORT()
CONTAINS

COMPLEX FUNCTION T1(X,Y)
  REAL :: X
  REAL, OPTIONAL :: Y
  T1=CMPLX(X,Y)
END FUNCTION T1

INTEGER FUNCTION M1(A1,A2,A3,A4)
  INTEGER :: A1,A2
  INTEGER, OPTIONAL :: A3,A4
  M1=MAX(A1,A2,A3,A4)
END FUNCTION M1

END

max_stars_repo_name: maurizioabba/rose


text:
      program fixpoint05

C     Suppose J is known to be greater than 5 (in fact, equal to 5
C     because of the previous loop) and the loop certainly entered: the
C     information is lost!

C     The problem seems to be due to the precondition equation for
C     entered loops: p' = t(t*(p)). The certainly executed iteration is
C     applied too late, when the information brought by p is lost thru
C     t*. However, t*(t(p)) would still not be satisfactory. t* should
C     be computed for the subset defined by p.

      if(j.ge.5) then

         do i = 1, 5
            if(j.ge.5) then
               j = 5
            endif
         enddo

         print *, j

      endif

      end


text:
!
! Copyright (c) 2018, 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 first_submod35
character (len = 10) :: word
integer :: arr(2:7, 3:18)
integer :: arr2(4:17)
interface
    module subroutine printWord2
    end subroutine
    module subroutine printWord3
    end subroutine
    module subroutine printWord4
    end subroutine
    module subroutine printWord5
    end subroutine
end interface
contains
    subroutine printWord1
        print *, word;
    end subroutine
end module first_submod35

submodule (first_submod35) second_submod35
character (len = 10) word2
contains
   module procedure printWord2
        call set
        print *, word, word2
        print *, "shape of arr is ", shape(arr)
    end procedure
    subroutine set
        word2 = "second"
    end subroutine
end submodule second_submod35

submodule (first_submod35:second_submod35) third_submod35
character (len = 10) word3
contains
   module procedure printWord3
        call set
        print *, word, word2, word3
        print *, "lbounds of arr are: ", lbound(arr)
    end procedure
    subroutine set
        word3 = "third"
    end subroutine

end submodule third_submod35

submodule (first_submod35:third_submod35) fourth_submod35
character (len = 10) word4
contains
   module procedure printWord4
        call set
        print *, word, word2, word3,  word4
        print *, "ubounds of arr are: ", ubound(arr)
        ! check the upper bounds of arr in the fourth submodule
        if ( .NOT. ((ubound(arr, dim=1) .EQ. 7) .AND. (ubound(arr, dim=2) .EQ. 18) )) then
            print *, "TEST FAILED upper bounds incorrect"
        end if
    end procedure
    subroutine set
        word4 = "fourth"
    end subroutine

end submodule fourth_submod35

submodule (first_submod35:fourth_submod35) fifth_submod35
character (len = 10) word5
contains
   module procedure printWord5
        call set
        print *, word, word2, word3, word4, word5
        print *, "shape of result of shape of arr is: ", shape(shape(arr))
        if ( (word .EQ. "first") .AND. (word2 .EQ. "second") .AND. (word3 .EQ. "third") .AND. (word4 .EQ. "fourth") .AND. (word5 .EQ. "fifth") ) then
            print *, " PASS"
        else
            print *, "TEST FAILED submodule host association not working"
        end if
        ! check the upper bounds of arr in the fifth submodule
        if ( .NOT. ((ubound(arr, dim=1) .EQ. 7) .AND. (ubound(arr, dim=2) .EQ. 18) )) then
            print *, "TEST FAILED upper bounds incorrect"
        end if

    end procedure
    subroutine set
        word5 = "fifth"
    end subroutine
end submodule fifth_submod35

program dothis
use first_submod35
implicit none
word = 'first'

call printWord1
call printWord2
call printWord3
call printWord4
call printWord5
        ! check the upper bounds of arr in the main program
        if ( .NOT. ((ubound(arr, dim=1) .EQ. 7) .AND. (ubound(arr, dim=2) .EQ. 18) )) then
            print *, "TEST FAILED upper bounds incorrect"
        end if

end program dothis

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
c
c
        subroutine model(bgmodel,grid_lat,grid_lon,lat_s,lon_s,elev_s,
     &                   i4time_cur,laps_cycle_time,n_obs,ni,nj,nk,
     &                   mxsta,dt,dtd,du,dv,dpmsl,dalt,badflag)
c
c
c*********************************************************************
c
c     Original: John McGinley, NOAA/FSL  Spring 1998
c     Changes:
c       24 Aug 1998  Peter Stamus, NOAA/FSL
c          Make code dynamic, housekeeping changes, for use in LAPS.
c       14 Dec 1999  John McGinley and Peter Stamus, NOAA/FSL
c          Completely new version.
c       06 Jan 2000  Peter Stamus, NOAA/FSL
c          Yet another new version.  This one uses actual model data.
c
c     Notes: Currently (Jan 2000), if any of the LGA fields is missing
c            we set all of the trends to zero.  Plans are to check for
c            each variable separately, and to use a simple model (like
c            diurnal cycle) to estimate the trend better.
c
c*********************************************************************
c
c.....  Arrays for the LAPS domain lat/lon
c
        real grid_lat(ni,nj), grid_lon(ni,nj)
c
c.....  Arrays for 3d fields.
c
        real heights_3d(ni,nj,nk), p_3d(ni,nj,nk)
        real u_3d(ni,nj,nk), v_3d(ni,nj,nk)
        real t_3d(ni,nj,nk), q_3d(ni,nj,nk)
        real p_1d(nk)
        character var_lga*3
c
c.....  Arrays for current run interpolation from 3d to station.
c
        real t_cur(mxsta), td_cur(mxsta)
        real u_cur(mxsta), v_cur(mxsta)
        real psfc_cur(mxsta), pmsl_cur(mxsta), alt_cur(mxsta)
c
c.....  Arrays for previous run interpolation from 3d to station.
c
        real t_1(mxsta), td_1(mxsta)
        real u_1(mxsta), v_1(mxsta)
        real psfc_1(mxsta), pmsl_1(mxsta), alt_1(mxsta)
c
c.....  Arrays for trends of each variable over the past cycle.
c.....  This is the output.
c
        real dt(mxsta), dtd(mxsta), du(mxsta), dv(mxsta)
        real dpmsl(mxsta), dalt(mxsta)
c
c.....  Arrays for station location info (lat/lon/elevation and i,j)
c
        real lat_s(mxsta), lon_s(mxsta), elev_s(mxsta)
        integer ii(mxsta), jj(mxsta)
c
c.....  Other stuff
c
        real    lapse_t, lapse_td
        integer bgmodel
c
c
c.....  Start.  Fill the trend arrays with something.
c
        do i=1,mxsta
           dt(i)    = badflag
           dtd(i)   = badflag
           du(i)    = badflag
           dv(i)    = badflag
           dpmsl(i) = badflag
           dalt(i)  = badflag
        enddo !i
c
        t_ref = -47.           !for specific humidity to Td calc
        lapse_t = -0.01167     !std atm lapse rate
        lapse_td = -0.007      !       "
c
c.....  Find and store the i,j location of each station in the LAPS
c.....  grid.
c
        do ista=1,n_obs
          call latlon_to_rlapsgrid(lat_s(ista),lon_s(ista),grid_lat,
     &       grid_lon,ni,nj,rii,rjj,istatus)
          ii(ista) = rii
          jj(ista) = rjj
        enddo !ista
c
c.....  Start with the current hour.
c.....  Get the 3d fields for the current time.
c
        print *,' Getting LAPS pressure levels'
        call get_pres_3d(i4time_cur,ni,nj,nk,p_3d,istatus)
        do k=1,nk  !for now...grab a column of pressures
           p_1d(k) = p_3d(1,1,k) * 0.01  !conv from Pa to mb too...
        enddo !k
c
        print *,' Getting LGA 3-d Heights'
        var_lga = 'HT '
        call get_modelfg_3d(i4time_cur,var_lga,ni,nj,nk,heights_3d,
     &                      istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA HT not available. '
           go to 800
        else

        endif
c
        print *,' Getting LGA 3-d Temperatures'
        var_lga = 'T3 '
        call get_modelfg_3d(i4time_cur,var_lga,ni,nj,nk,t_3d,istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA T not available. '
           go to 800
        else

        endif
c
        print *,' Getting LGA 3-d Specific Humidity'
        var_lga = 'SH '  ! specific humidity
        call get_modelfg_3d(i4time_cur,var_lga,ni,nj,nk,q_3d,istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA SH not available. '
           go to 800
        else

        endif
c
        print *,' Getting LGA 3-d U-wind component'
        var_lga = 'U3 '
        call get_modelfg_3d(i4time_cur,var_lga,ni,nj,nk,u_3d,istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA U not available. '
           go to 800
        else

        endif
c
        print *,' Getting LGA 3-d V-wind component'
        var_lga = 'V3 '
        call get_modelfg_3d(i4time_cur,var_lga,ni,nj,nk,v_3d,istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA V not available. '
           go to 800
        else

        endif
c
c.....  Now, for the current time, interpolate the 3-d variables to the
c.....  proper location of each station.
c
        do ista=1,n_obs
c
c.....  Is the station within the LAPS grid?  If so, continue.  Otherwise,
c.....  set the current variables to badflag.
c
           if(ii(ista).ge.1 .and. ii(ista).le.ni .and.
     &        jj(ista).ge.1 .and. jj(ista).le.nj) then
c
c.....  Find the vertical index for the station elevation.
c
              iii = ii(ista)
              jjj = jj(ista)
              zlow = height_to_zcoord2(elev_s(ista),heights_3d,ni,nj,nk,
     &                              iii,jjj,istatus)
              if(istatus .ne. 1)then
                 write(6,*)
     &           ' Error in height_to_zcoord2 in interp_to_sfc', istatus
                 write(6,*) iii,jjj,zlow,elev_s(ista),
     &                                 (heights_3d(iii,jjj,k),k=1,nk)
                 return
              endif
c
              k = max(zlow, 1.)
c
              t_in = t_3d(iii,jjj,k)   !K          -- temp
              q_in = q_3d(iii,jjj,k)   !kg/kg      -- specific hum in
              p_in = 1.                !dummy arg  -- nothing needed here
c
              call compute_sfc_bgfields(bgmodel,ni,nj,nk,iii,jjj,
     &               k,elev_s(ista),heights_3d,t_3d,p_1d,
     &               q_3d,t_ref,p_in,t_in,q_in)
c
              t_cur(ista)    = t_in !K   -- temp
              td_cur(ista)   = q_in !K   -- Td returned here
              psfc_cur(ista) = p_in !Pa  -- calc sfc p
c
              call interp_to_elevation(elev_s(ista),iii,jjj,
     &                u_3d,heights_3d,ni,nj,nk,badflag,u_cur(ista))
              call interp_to_elevation(elev_s(ista),iii,jjj,
     &                v_3d,heights_3d,ni,nj,nk,badflag,v_cur(ista))
c
c.....  Convert units to what LAPS is using.
c
              t_cur(ista) = ((t_cur(ista) - 273.17) * 1.8) + 32.   !K to F
              td_cur(ista) = ((td_cur(ista) - 273.17) * 1.8) + 32. !K to F
              u_cur(ista) = u_cur(ista) * 1.94254                  !m/s to kt
              v_cur(ista) = v_cur(ista) * 1.94254                  !m/s to kt
              psfc_cur(ista) = psfc_cur(ista) * 0.01               !Pa to mb
c
c.....  Use the interpolated values to calculate a MSL pressure and
c.....  altimeter.  Note that the routine used is the same as that
c.....  in the LAPS surface analysis, which may or may not be giving
c.....  good values in all areas (e.g., high terrain, very cold areas).
c
              call reduce_p(t_cur(ista),td_cur(ista),psfc_cur(ista),
     &                   elev_s(ista),lapse_t,lapse_td,pmsl_cur(ista),
     &                   0.,badflag)
c
c.....  Altimeter is defined as being 10 feet (3.048 meters) above the
c.....  station elevation.  So, we'll add that height and use the
c.....  same routine.  There are other differences in the way that
c.....  MSL and altimeter are calculated, but we'll ignore them here.
c.....  The bias introduced should be small.
c
              elev_alt = elev_s(ista) + 3.048 !meters
              call reduce_p(t_cur(ista),td_cur(ista),psfc_cur(ista),
     &                   elev_alt,lapse_t,lapse_td,alt_cur(ista),
     &                   0.,badflag)
c
           else  !station outside of grid
c
              t_cur(ista)    = badflag
              td_cur(ista)   = badflag
              u_cur(ista)    = badflag
              v_cur(ista)    = badflag
              psfc_cur(ista) = badflag
              pmsl_cur(ista) = badflag
              alt_cur(ista)  = badflag
c
           endif
c
        enddo !ista
c
c.....  Now do all that again, but for the previous cycle.
c
        i4time_pre = i4time_cur - laps_cycle_time
c
c.....  Get the 3d fields for the previous cycle.
c
        print *,' Getting LAPS pressure levels'
        call get_pres_3d(i4time_pre,ni,nj,nk,p_3d,istatus)
        do k=1,nk  !for now...grab a column of pressures
           p_1d(k) = p_3d(1,1,k) * 0.01  !conv from Pa to mb too...
        enddo !k
c
        print *,' Getting LGA 3-d Heights'
        var_lga = 'HT '
        call get_modelfg_3d(i4time_pre,var_lga,ni,nj,nk,heights_3d,
     &                      istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA HT not available. '
           go to 800
        else

        endif
c
        print *,' Getting LGA 3-d Temperatures'
        var_lga = 'T3 '
        call get_modelfg_3d(i4time_pre,var_lga,ni,nj,nk,t_3d,istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA T not available. '
           go to 800
        else

        endif
c
        print *,' Getting LGA 3-d Specific Humidity'
        var_lga = 'SH '  ! specific humidity
        call get_modelfg_3d(i4time_pre,var_lga,ni,nj,nk,q_3d,istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA SH not available. '
           go to 800
        else

        endif
c
        print *,' Getting LGA 3-d U-wind component'
        var_lga = 'U3 '
        call get_modelfg_3d(i4time_pre,var_lga,ni,nj,nk,u_3d,istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA U not available. '
           go to 800
        else

        endif
c
        print *,' Getting LGA 3-d V-wind component'
        var_lga = 'V3 '
        call get_modelfg_3d(i4time_pre,var_lga,ni,nj,nk,v_3d,istatus)
c
        if(istatus .ne. 1)  then
           print *,' LGA V not available. '
           go to 800
        else

        endif
c
c.....  Now, for the previous cycle, interpolate the 3-d variables to the
c.....  proper location of each station.
c
        do ista=1,n_obs
c
c.....  Is the station within the LAPS grid?  If so, continue.  Otherwise,
c.....  set the previous variables to badflag.
c
           if(ii(ista).ge.1 .and. ii(ista).le.ni .and.
     &        jj(ista).ge.1 .and. jj(ista).le.nj) then
c
c.....  Find the vertical index for the station elevation.
c
              iii = ii(ista)
              jjj = jj(ista)
c
              zlow = height_to_zcoord2(elev_s(ista),heights_3d,ni,nj,nk,
     &                              iii,jjj,istatus)
              if(istatus .ne. 1)then
                 write(6,*)
     &           ' Error in height_to_zcoord2 in interp_to_sfc', istatus
                 write(6,*) iii,jjj,zlow,elev_s(ista),
     &                               (heights_3d(iii,jjj,k),k=1,nk)
                 return
              endif
c
              k = max(zlow, 1.)
c
              t_in = t_3d(iii,jjj,k)    !K          -- temp
              q_in = q_3d(iii,jjj,k)    !kg/kg      -- specific hum in
              p_in = 1.                 !dummy arg  -- nothing needed here
c
              call compute_sfc_bgfields(bgmodel,ni,nj,nk,iii,jjj,
     &               k,elev_s(ista),heights_3d,t_3d,p_1d,
     &               q_3d,t_ref,p_in,t_in,q_in)
c
              t_1(ista)    = t_in      !K   -- temp
              td_1(ista)   = q_in      !K   -- Td returned here
              psfc_1(ista) = p_in      !Pa  -- calc sfc p
c
              call interp_to_elevation(elev_s(ista),iii,jjj,
     &                   u_3d,heights_3d,ni,nj,nk,badflag,u_1(ista))
              call interp_to_elevation(elev_s(ista),iii,jjj,
     &                   v_3d,heights_3d,ni,nj,nk,badflag,v_1(ista))
c
c.....  Convert units to what LAPS is using.
c
              t_1(ista) = ((t_1(ista) - 273.17) * 1.8) + 32.   !K to F
              td_1(ista) = ((td_1(ista) - 273.17) * 1.8) + 32. !K to F
              u_1(ista) = u_1(ista) * 1.94254                  !m/s to kt
              v_1(ista) = v_1(ista) * 1.94254                  !m/s to kt
              psfc_1(ista) = psfc_1(ista) * 0.01               !Pa to mb
c
c.....  Use the interpolated values to calculate a MSL pressure and
c.....  altimeter.  Note that the routine used is the same as that
c.....  in the LAPS surface analysis, which may or may not be giving
c.....  good values in all areas (e.g., high terrain, very cold areas).
c
              call reduce_p(t_1(ista),td_1(ista),psfc_1(ista),
     &                   elev_s(ista),lapse_t,lapse_td,pmsl_1(ista),
     &                   0.,badflag)
c
c.....  Altimeter is defined as being 10 feet (3.048 meters) above the
c.....  station elevation.  So, we'll add that height and use the
c.....  same routine.  There are other differences in the way that
c.....  MSL and altimeter are calculated, but we'll ignore them here.
c.....  The bias introduced should be small.
c
              elev_alt = elev_s(ista) + 3.048 !meters
              call reduce_p(t_1(ista),td_1(ista),psfc_1(ista),
     &                   elev_alt,lapse_t,lapse_td,alt_1(ista),
     &                   0.,badflag)
c
           else  !station outside of grid
c
              t_1(ista)    = badflag
              td_1(ista)   = badflag
              u_1(ista)    = badflag
              v_1(ista)    = badflag
              psfc_1(ista) = badflag
              pmsl_1(ista) = badflag
              alt_1(ista)  = badflag
c
           endif
c

         enddo !ista

c
c.....  Now calculate the trends for each variable at each station.
c.....  Trend is defined as var(current) minus var(previous).
c.....  For now, if either estimate is bad, set the trend to zero
c.....  (which is the same as persistance).
c
         do i=1,n_obs
            if(t_cur(i).le.badflag .or. t_1(i).le.badflag) then
               dt(i) = 0.
            else
               dt(i) = t_cur(i) - t_1(i)                   !deg F
            endif
c
            if(td_cur(i).le.badflag .or. td_1(i).le.badflag) then
               dtd(i) = 0.
            else
               dtd(i) = td_cur(i) - td_1(i)                !deg F
            endif
c
            if(u_cur(i).le.badflag .or. u_1(i).le.badflag) then
               du(i) = 0.
            else
               du(i) = u_cur(i) - u_1(i)                   !kt
            endif
c
            if(v_cur(i).le.badflag .or. v_1(i).le.badflag) then
               dv(i) = 0.
            else
               dv(i) = v_cur(i) - v_1(i)                   !kt
            endif
c
            if(pmsl_cur(i).le.badflag .or. pmsl_1(i).le.badflag) then
               dpmsl(i) = 0.
            else
               dpmsl(i) = pmsl_cur(i) - pmsl_1(i)          !mb
            endif
c
            if(alt_cur(i).le.badflag .or. alt_1(i).le.badflag) then
               dalt(i) = 0.
            else
               dalt(i) = alt_cur(i) - alt_1(i)             !mb
            endif
         enddo !i
c
c..... That's all.
c
        return
c
c..... If missing the LGA from one or the other times, do something
c..... else to get the trends.  For now, just set them equal to
c..... zero (persistance forecast).
c
 800    continue
        print *,
     &    ' **WARNING. Problem getting LGA. Setting trends to zero.'
        do i=1,n_obs
           dt(i)    = 0.
           dtd(i)   = 0.
           du(i)    = 0.
           dv(i)    = 0.
           dpmsl(i) = 0.
           dalt(i)  = 0.
        enddo !i
        return
c
        end
c
c
        subroutine interp_to_elevation(pt_elev,pt_i,pt_j,field_3d,
     &                         heights_3d,ni,nj,nk,badflag,interp)
c
c=================================================================
c
c     Routine to interpolate a 3-d field to a single point in 3-d
c     space.  The point can be the actual "surface" (the
c     topography), or any other point within the grid.
c
c     Based on code in the LAPS wind analysis, Steve Albers, FSL,
c     and the "interp_to_sfc" routine.
c
c     Original: 01-06-00   Peter A. Stamus, NOAA/FSL
c     Changes:
c
c=================================================================
c
        real pt_elev                 !elevation of desired point
        real interp                  !output interpolated value
        real field_3d(ni,nj,nk),     !3d array of desired variable
     &       heights_3d(ni,nj,nk)    !3d height field
        integer pt_i, pt_j           !i,j of desired point
c
cc      i_sfc_bad = 0
c
c..... Interpolate from the 3-d grid to the elevation at the point.
c
        zlow = height_to_zcoord2(pt_elev,heights_3d,ni,nj,nk,
     &                           pt_i,pt_j,istatus)
        if(istatus .ne. 1)then
           write(6,*)
     &       ' Error in height_to_zcoord2 in interp_to_sfc', istatus
           write(6,*) pt_i,pt_j,zlow,pt_elev,(heights_3d(i,j,k),k=1,nk)
           return
        endif
c
        klow = max(zlow, 1.)
        khigh = klow + 1
        fraclow = float(khigh) - zlow
        frachigh = 1.0 - fraclow
c
        if( field_3d(pt_i,pt_j,klow)  .eq. badflag .or.
     &      field_3d(pt_i,pt_j,khigh) .eq. badflag) then
           write(6,3333)
 3333      format(' Warning: cannot interpolate to the point')
c          i_sfc_bad = 1
           interp = badflag
        else
           interp = field_3d(pt_i,pt_j,klow ) * fraclow  +
     &              field_3d(pt_i,pt_j,khigh) * frachigh
        endif
c
c..... That's all.
c
        return
        end
c
c
	Subroutine perturb(ta,tb,dta,imax,m,offset,onoff)
c
c*********************************************************************
c
c     Original: John McGinley, NOAA/FSL  Spring 1998
c     Changes:
c       24 Aug 1998  Peter Stamus, NOAA/FSL
c          Make code dynamic, housekeeping changes, for use in LAPS.
c
c     Notes:
c
c*********************************************************************
c
	parameter(badflag=-99.9)
	real ta(m),tb(m),dta(m)
	integer onoff
c
	if (onoff.eq.1) then
	   do i=1,imax
	      if(ta(i).ne.badflag.and.tb(i).ne.badflag) then
		 dta(i)=ta(i)-tb(i)+offset
	      else
		 dta(i)=badflag
	      endif
	   enddo !i
	else
	   do i=1,imax
	      ta(i)=tb(i)+(dta(i)-offset)
	   enddo !i
	endif
c
	return
	end
c
c
	Subroutine weights(mwt,dwt,ua,va,theta,lat,lon,icycle,imax,m,
     &                      hcutoff,vcutoff)
c
c*********************************************************************
c
c     Original: John McGinley, NOAA/FSL  Spring 1998
c     Changes:
c       24 Aug 1998  Peter Stamus, NOAA/FSL
c          Make code dynamic, housekeeping changes, for use in LAPS.
c       14 Dec 1999  Peter Stamus, NOAA/FSL
c          New version (minor changes from old).
c
c     Notes:
c
c*********************************************************************
c
	common tab(10000),pi,re,rdpdg,reorpd
	real mwt(m,m),dwt(m,m),ua(m),va(m),theta(m),lat(m),lon(m)
        integer icycle
c
        cycle = float( icycle )
	expsclr=1000.
c       hcutoff= specify  !dist in m for wt e**-1
c       vcutoff=specify	 !deg C for vertical wt to go to e**-1
	C3=2.		 !wt on wind factor
	c1=sqrt(expsclr)/hcutoff**2
	c1m=sqrt(expsclr)/hcutoff**2
	c2=hcutoff**2/vcutoff**2
	c2m=hcutoff**2/vcutoff**2
c
c corrects for wind component along station/station vector
c
	do i=1,imax
	do j=1,imax
	   ang=(lat(i)+lat(j))*.5*rdpdg
	   dy=(lat(j)-lat(i))*reorpd
	   dx=(lon(j)-lon(i))*reorpd*cos(ang)
	   rr=(dx*dx+dy*dy)+1.e-20
	   r=sqrt(rr)
	   uav=dx*(ua(i)+ua(j))*.5/r
	   vav=dy*(va(i)+va(i))*.5/r
	   dxx=dx+uav*cycle/c3
	   dyy=dy+vav*cycle/c3
	   dz=theta(j)-theta(i)
	   iii=int(c1*((dx*dx+dy*dy) +c2*dz*dz))+1
	   iiii=int(c1m*((dxx*dxx+dyy*dyy) +c2m*dz*dz))+1
	   if(iii.gt.10000) iii=10000
	   if(iiii.gt.10000) iiii=10000
	   dwt(i,j)=tab(iii)
	   mwt(i,j)=tab(iiii)
	enddo !j
	enddo !i
c
c normalize dwt,mwt
c
	do i=1,imax
	   sum=0.
	   summ=0.
	   do j=1,imax
	      sum=dwt(i,j)+sum
	      summ=mwt(i,j)+summ
	   enddo !on j
	   do j=1,imax
	      dwt(i,j)=dwt(i,j)/sum
	      mwt(i,j)=mwt(i,j)/summ
	   enddo !on j
	enddo !on i
c
	return
	end
c
c
        subroutine project(ord,y,monster,nv,nvar,maxsta,m,ncycles,
     &              nn,atime,it,icyc,oberr,badthr,bmonster,ihr)
c
c*********************************************************************
c
c     This routine uses Taylor series of order 'ord' to forecast
c     the variable 'nv' one cycle.
c
c     Original: John McGinley, NOAA/FSL  Spring 1998
c     Changes:
c       24 Aug 1998  Peter Stamus, NOAA/FSL
c          Make code dynamic, housekeeping changes, for use in LAPS.
c       14 Dec 1999  John McGinley and Peter Stamus, NOAA/FSL
c          Completely new version.
c       Bias correction (bmonster) added to self trend. Includes
c      a cycle dependent bias correction based on z time (ihr)
c      and model performance
c
c     Notes:
c          For the time being (14 Dec 1999) 'ord' is limited to 2.
c
c*********************************************************************
c
        parameter(badflag=-99.9)
        common tab(10000),pi,re,rdpdg,reorpd
        integer ord
        real y(m), monster(m,ncycles,nvar),bmonster(m,24,nvar)
        character*24 atime
c
c.....  Apply a filter to each derivative
c
        em1=exp(-1.)
        em1s=em1*em1
        do i=1,maxsta
           a=1.
           b=1.
           if(monster(i,2,nv).eq.badflag) a=0.
           if(monster(i,3,nv).eq.badflag) b=0.
           sum0 = monster(i,1,nv)*(1.+a*em1+b*em1s*.5) +
     &            a*monster(i,2,nv)*(-a*em1-b*em1s) +
     &            monster(i,3,nv)*(b*em1s*.5)
           y(i)=sum0-bmonster(i,ihr,nv)
        enddo !on i
c
        return
        end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
!##############################################################################
!# ****************************************************************************
!# <name> cc2dminim2boundary </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains the definition of the analytic boundary conditions
!# as well as discretisation routines for the boundary.
!#
!# The following files can be found here:
!#
!# 1.) c2d2_initAnalyticBC
!#     -> Initialise analytic boundary conditions
!#
!# 2.) c2d2_initDiscreteBC
!#     -> Discretise analytic boundary conditions, create discrete
!#        boundary conditions
!#
!# 3.) c2d2_updateDiscreteBC
!#     -> Reassemble discrete boundary conditions
!#
!# 4.) c2d2_doneBC
!#     -> Release discrete and analytic boundary conditions
!#
!# 5.) c2d2_implementBC
!#     -> Implement discrete boundary conditions into solution/RHS vectors
!#        and matrix on finest level
!#
!# </purpose>
!##############################################################################

module cc2dmediumm2boundary

  use fsystem
  use storage
  use linearsolver
  use boundary
  use bilinearformevaluation
  use linearformevaluation
  use cubature
  use matrixfilters
  use vectorfilters
  use bcassembly
  use triangulation
  use spatialdiscretisation
  use coarsegridcorrection
  use spdiscprojection
  use nonlinearsolver
  use paramlist

  use collection
  use convection

  use cc2dmediumm2basic
  use cc2dmediumm2boundarydef
  use cc2dmedium_callback

  implicit none

contains

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

!<subroutine>

  subroutine c2d2_initAnalyticBC (rproblem)

!<description>
  ! This initialises the analytic boundary conditions of the problem
  ! and saves them to the problem structure.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem
!</inputoutput>

!</subroutine>

    ! local variables

    ! A pointer to the discretisation structure with the data.
    type(t_blockDiscretisation), pointer :: p_rdiscretisation

    integer :: i

    ! Initialise the boundary condition by parsing the parameter files.
    ! This initialises rproblem%p_rboundaryConditions by parsing DAT-file
    ! parameters in the parameter list rproblem%rparamList
    call c2d2_parseBDconditions (rproblem)

    ! Initialise the boundary conditions of fictitious boundary components
    call c2d2_parseFBDconditions (rproblem)

  end subroutine

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

!<subroutine>

  subroutine c2d2_initDiscreteBC (rproblem,rvector,rrhs)

!<description>
  ! This calculates the discrete version of the boundary conditions and
  ! assigns it to the system matrix and RHS vector.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem

  ! A vector structure for the solution vector. The discrete BC structures are
  ! attached to that.
  type(t_vectorBlock), intent(INOUT) :: rvector

  ! A vector structure for the RHS vector. The discrete BC structures are
  ! attached to that.
  type(t_vectorBlock), intent(INOUT) :: rrhs
!</inputoutput>

!</subroutine>

  ! local variables
  integer :: i

  ! A pointer to the system matrix and the RHS vector as well as
  ! the discretisation
  type(t_blockDiscretisation), pointer :: p_rdiscretisation

  ! Pointer to structure for saving discrete BC's:
  type(t_discreteBC), pointer :: p_rdiscreteBC
  type(t_discreteFBC), pointer :: p_rdiscreteFBC

    ! Initialise the collection for the assembly process with callback routines.
    ! Basically, this stores the simulation time in the collection if the
    ! simulation is nonstationary.
    call c2d2_initCollectForAssembly (rproblem,rproblem%rcollection)

    do i=rproblem%NLMIN,rproblem%NLMAX

      ! From the matrix or the RHS we have access to the discretisation and the
      ! analytic boundary conditions.
      p_rdiscretisation => rproblem%RlevelInfo(i)%p_rdiscretisation

      ! For the discrete problem, we need a discrete version of the above
      ! boundary conditions. So we have to discretise them.
      ! The following routine gives back p_rdiscreteBC, a pointer to a
      ! discrete version of the boundary conditions. Remark that
      ! the pointer has to be nullified before calling the routine,
      ! otherwise, the routine tries to update the boundary conditions
      ! in p_rdiscreteBC!
      ! getBoundaryValues is a callback routine that specifies the
      ! values on the boundary. We pass our collection structure as well
      ! to this routine, so the callback routine has access to everything what is
      ! in the collection.
      !
      ! On maximum level, discretrise everything. On lower level, discretise
      ! only for the implementation into the matrices and defect vector.
      ! That's enough, as the lower levels are only used for preconditioning
      ! of defect vectors.

      nullify(rproblem%RlevelInfo(i)%p_rdiscreteBC)
      if (i .eq. rproblem%NLMAX) then
        call bcasm_discretiseBC (p_rdiscretisation, &
                                 rproblem%RlevelInfo(i)%p_rdiscreteBC, &
                                .false.,getBoundaryValues, &
                                rproblem%rcollection)
      else
        call bcasm_discretiseBC (p_rdiscretisation, &
                                 rproblem%RlevelInfo(i)%p_rdiscreteBC, &
                                .false.,getBoundaryValues, &
                                rproblem%rcollection,BCASM_DISCFORDEFMAT)
      end if

      ! The same way, discretise the fictitious boundary conditions and hang
      ! them in into the matrix/vectors.
      nullify(rproblem%RlevelInfo(i)%p_rdiscreteFBC)
      if (i .eq. rproblem%NLMAX) then
        call bcasm_discretiseFBC (p_rdiscretisation,&
                                  rproblem%RlevelInfo(i)%p_rdiscreteFBC,.false., &
                                  getBoundaryValuesFBC,rproblem%rcollection)
      else
        call bcasm_discretiseFBC (p_rdiscretisation,&
                                  rproblem%RlevelInfo(i)%p_rdiscreteFBC,.false., &
                                  getBoundaryValuesFBC,rproblem%rcollection,&
                                  BCASM_DISCFORDEFMAT)
      end if

      ! Hang the pointer into the the matrix. That way, these
      ! boundary conditions are always connected to that matrix and that
      ! vector.
      p_rdiscreteBC => rproblem%RlevelInfo(i)%p_rdiscreteBC

      rproblem%RlevelInfo(i)%rpreallocatedSystemMatrix%p_rdiscreteBC => p_rdiscreteBC

      ! Also hang in the boundary conditions into the temporary vector that is
      ! used for the creation of solutions on lower levels.
      ! This allows us to filter this vector when we create it.
      rproblem%RlevelInfo(i)%rtempVector%p_rdiscreteBC => p_rdiscreteBC

      ! The same for the fictitious boundary boundary conditions.
      p_rdiscreteFBC => rproblem%RlevelInfo(i)%p_rdiscreteFBC
      rproblem%RlevelInfo(i)%rpreallocatedSystemMatrix%p_rdiscreteBCfict => &
          p_rdiscreteFBC
      rproblem%RlevelInfo(i)%rtempVector%p_rdiscreteBCfict => p_rdiscreteFBC

    end do

    ! On the finest level, attach the discrete BC also
    ! to the solution and RHS vector. They need it to be compatible
    ! to the matrix on the finest level.
    p_rdiscreteBC => rproblem%RlevelInfo(rproblem%NLMAX)%p_rdiscreteBC

    rrhs%p_rdiscreteBC => p_rdiscreteBC
    rvector%p_rdiscreteBC => p_rdiscreteBC

    ! The same with the fictitious boundary BC's
    p_rdiscreteFBC => rproblem%RlevelInfo(rproblem%NLMAX)%p_rdiscreteFBC
    rrhs%p_rdiscreteBCfict => p_rdiscreteFBC
    rvector%p_rdiscreteBCfict => p_rdiscreteFBC

    ! Clean up the collection (as we are done with the assembly, that's it.
    call c2d2_doneCollectForAssembly (rproblem,rproblem%rcollection)

  end subroutine

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

!<subroutine>

  subroutine c2d2_updateDiscreteBC (rproblem, bforce)

!<description>
  ! This updates the discrete version of the boundary conditions. The BC's
  ! are reassembled according to the current situation of the simulation.
!</description>

!<input>
  ! Normally, only those BC's are update that are marked as 'nonstatic'
  ! (e.g. for time-dependent simulation). When bforce=TRUE, all discrete
  ! BC's are rebuild, which might be necessary in optimisation environments
  ! when testing multiple stationary configurations.
  logical, intent(IN) :: bforce
!</input>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem
!</inputoutput>

!</subroutine>

    ! local variables
    integer :: i

    ! A pointer to the system matrix and the RHS vector as well as
    ! the discretisation
    type(t_blockDiscretisation), pointer :: p_rdiscretisation

    ! Initialise the collection for the assembly process with callback routines.
    ! Basically, this stores the simulation time in the collection if the
    ! simulation is nonstationary.
    call c2d2_initCollectForAssembly (rproblem,rproblem%rcollection)

    do i=rproblem%NLMIN,rproblem%NLMAX

      ! From the matrix or the RHS we have access to the discretisation and the
      ! analytic boundary conditions.
      p_rdiscretisation => rproblem%RlevelInfo(i)%p_rdiscretisation

      ! For the discrete problem, we need a discrete version of the above
      ! boundary conditions. So we have to discretise them.
      ! The following routine gives back p_rdiscreteBC, a pointer to a
      ! discrete version of the boundary conditions. Remark that
      ! the pointer has to be nullified before calling the routine,
      ! otherwise, the routine tries to update the boundary conditions
      ! in p_rdiscreteBC!
      ! getBoundaryValues is a callback routine that specifies the
      ! values on the boundary. We pass our collection structure as well
      ! to this routine, so the callback routine has access to everything what is
      ! in the collection.
      !
      ! On maximum level, discretrise everything. On lower level, discretise
      ! only for the implementation into the matrices and defect vector.
      ! That's enough, as the lower levels are only used for preconditioning
      ! of defect vectors.
      !
      ! Note that we do not NULLIFY the pointer to the discrete BC's here.
      ! That's the only difference to the initDiscreteBC routine above,
      ! as this instructs bcasm_discretiseBC to update nonstatic BC's.

      if (i .eq. rproblem%NLMAX) then
        call bcasm_discretiseBC (p_rdiscretisation, &
                                 rproblem%RlevelInfo(i)%p_rdiscreteBC, &
                                 bforce,getBoundaryValues, &
                                 rproblem%rcollection)
      else
        call bcasm_discretiseBC (p_rdiscretisation, &
                                 rproblem%RlevelInfo(i)%p_rdiscreteBC, &
                                 bforce,getBoundaryValues, &
                                 rproblem%rcollection,BCASM_DISCFORDEFMAT)
      end if

      ! The same way, discretise the fictitious boundary conditions and hang
      ! them in into the matrix/vectors.
      if (i .eq. rproblem%NLMAX) then
        call bcasm_discretiseFBC (p_rdiscretisation,&
                                  rproblem%RlevelInfo(i)%p_rdiscreteFBC,bforce, &
                                  getBoundaryValuesFBC,rproblem%rcollection)
      else
        call bcasm_discretiseFBC (p_rdiscretisation,&
                                  rproblem%RlevelInfo(i)%p_rdiscreteFBC,bforce, &
                                  getBoundaryValuesFBC,rproblem%rcollection,&
                                  BCASM_DISCFORDEFMAT)
      end if

    end do

    ! Clean up the collection (as we are done with the assembly, that's it.
    call c2d2_doneCollectForAssembly (rproblem,rproblem%rcollection)

  end subroutine

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

!<subroutine>

  subroutine c2d2_implementBC (rproblem,rvector,rrhs,&
      bsolvector,brhsvector)

!<description>
  ! Implements boundary conditions into the RHS and into a given solution vector.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem

  ! A vector structure for the solution vector. The discrete BC's are implemented
  ! into that.
  type(t_vectorBlock), intent(INOUT) :: rvector

  ! A vector structure for the RHS vector. The discrete BC's are implamented into that.
  type(t_vectorBlock), intent(INOUT) :: rrhs

  ! Whether to implement the BC's into the solution vector or not
  logical, intent(IN) :: bsolvector

  ! Whether to implement the BC's into the solution vector or not
  logical, intent(IN) :: brhsvector
!</inputoutput>

!</subroutine>

  ! local variables
  integer :: i,ilvmax

  ! A pointer to the system matrix and the RHS vector as well as
  ! the discretisation
  type(t_matrixBlock), pointer :: p_rmatrix

    ! Get our the right hand side and solution from the problem structure
    ! on the finest level
    ilvmax = rproblem%NLMAX

    if (bsolvector) then

      ! Implement discrete boundary conditions into solution vector by
      ! filtering the vector.
      call vecfil_discreteBCsol (rvector)

      ! Implement discrete boundary conditions of fictitioous boundary comnponents
      ! into solution vector by filtering the vector.
      call vecfil_discreteFBCsol (rvector)

    end if

    if (brhsvector) then

      ! Implement pressure drop boundary conditions into RHS vector
      ! if there are any.
      call vecfil_discreteNLPDropBCrhs (rrhs)

      ! Implement discrete boundary conditions into RHS vector by
      ! filtering the vector.
      call vecfil_discreteBCrhs (rrhs)

      ! Implement discrete boundary conditions of fictitious boundary components
      ! into RHS vector by filtering the vector.
      call vecfil_discreteFBCrhs (rrhs)

    end if

    ! Implementation of boundary conditions into matrices deactivated.
    ! It's simply not used as the global matrices are actually not used outside
    ! of the nonlinear iteration!
    !
    !  IF (bmatrices) THEN
    !
    !    ! Implement discrete boundary conditions into the matrices on all
    !    ! levels, too.
    !    ! In fact, this modifies the B-matrices. The A-matrices are overwritten
    !    ! later and must then be modified again!
    !    DO i=rproblem%NLMIN ,rproblem%NLMAX
    !      p_rmatrix => rproblem%RlevelInfo(i)%rpreallocatedSystemMatrix
    !      CALL matfil_discreteBC (p_rmatrix)  ! standard boundary conditions
    !      CALL matfil_discreteFBC (p_rmatrix)  ! fictitious boundary boundary conditions
    !    END DO
    !
    !  END IF

  end subroutine

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

!<subroutine>

  subroutine c2d2_doneBC (rproblem)

!<description>
  ! Releases discrete and analytic boundary conditions from the heap.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(INOUT), target :: rproblem
!</inputoutput>

!</subroutine>

  ! local variables
  integer :: i
  type(t_fparser), pointer :: p_rparser

    do i=rproblem%NLMAX,rproblem%NLMIN,-1
      ! Release our discrete version of the boundary conditions
      call bcasm_releaseDiscreteBC (rproblem%RlevelInfo(i)%p_rdiscreteBC)

      ! as well as the discrete version of the BC's for fictitious boundaries
      call bcasm_releaseDiscreteFBC (rproblem%RlevelInfo(i)%p_rdiscreteFBC)

      ! ...and also the corresponding analytic description.
      call bcond_doneBC (rproblem%p_rboundaryConditions)
    end do

    ! Release the parser object with all the expressions to be evaluated
    ! on the boundary.
    p_rparser => collct_getvalue_pars (rproblem%rcollection, BDC_BDPARSER, &
                                       0, SEC_SBDEXPRESSIONS)
    call fparser_release (p_rparser)
    deallocate (p_rparser)

    ! Remove the boundary value parser from the collection
    call collct_deletevalue (rproblem%rcollection, BDC_BDPARSER)

    ! Remove the boundary-expression section we added earlier,
    ! with all their content.
    call collct_deletesection(rproblem%rcollection,SEC_SBDEXPRESSIONS)

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
subroutine SHReturnTapersM(theta0, lmax, m, tapers, eigenvalues, shannon, &
                           degrees, ntapers, exitstatus)
!------------------------------------------------------------------------------
!
!   This subroutine will return all the eigenvalues and eigenfunctions for the
!   space-concentration problem of a spherical cap of angular radius theta0.
!   The returned eigenfunctions correspond to "geodesy" normalized spherical
!   harmonic coefficients, and the eigenfunctions are further normalized such
!   that they have unit power (i.e., the integral of the function squared over
!   the sphere divided by 4 pi is 1, and the sum of the squares of their
!   coefficients is 1). If the optional vector DEGREES is specified, then the
!   eigenfunctions will be computed using only those degrees where DEGREES(l+1)
!   is not zero.
!
!   When possible, the eigenfunctions are calculated using the kernel of
!   Grunbaum et al. 1982 and the eigenvalues are then calculated by integration
!   using the definition of the space-concentration problem. Use of the
!   Grunbaum et al. kernel is prefered over the space-concentration kernel as
!   the eigenfunctions of the later are unreliable when there are several
!   eigenvalues identical (within machine precision) to either 1 or zero. If,
!   the optional parameter DEGREES is specified, and at least one element is
!   zero for degrees greater or equal to abs(m), then the eigenfunctions and
!   eigenvalues will instead be computed directly using the space-concentration
!   kernel.
!
!   Calling Parameters
!
!       IN
!           theta0          Angular radius of spherical cap in RADIANS.
!           lmax            Maximum spherical harmonic degree
!                           for the concentration problem.
!           m               Angular order of the concentration
!                           problem (m=0 corresponds to isotropic case).
!
!       OUT
!           tapers          An (lmax+1) by (lmax+1-abs(m)) array containing
!                           all the eigenfunctions of the space-
!                           concentration kernel. Eigenfunctions
!                           are listed by columns in decreasing order
!                           corresponding to value of their eigenvalue.
!                           Only the first ntapers columns are non-zero.
!           eigenvalues     A vector of length lmax+1-abs(m) containing the
!                           eigenvalued corresponding to the individual
!                           eigenfunctions.
!
!       OPTIONAL (IN)
!           degrees     Specify those degrees of the coupling matrix to
!                       compute. If degrees(l+1) is zero, degree l will not
!                       be computed.
!
!       OPTIONAL (OUT)
!           shannon     Shannon number as calculated from the trace of the
!                       kernel.
!           ntapers     The number of non-zero tapers.
!           exitstatus  If present, instead of executing a STOP when an error
!                       is encountered, the variable exitstatus will be
!                       returned describing the error.
!                       0 = No errors;
!                       1 = Improper dimensions of input array;
!                       2 = Improper bounds for input variable;
!                       3 = Error allocating memory;
!                       4 = File IO error.
!
!   Copyright (c) 2005-2019, SHTOOLS
!   All rights reserved.
!
!------------------------------------------------------------------------------
    use SHTOOLS, only: ComputeDG82, ComputeDm, EigValVecSymTri, EigValVecSym, &
                       PreGLQ, PlmBar, PlmIndex
    use ftypes

    implicit none

    real(dp), intent(in) :: theta0
    integer, intent(in) :: lmax, m
    real(dp), intent(out) :: tapers(:,:), eigenvalues(:)
    real(dp), intent(out), optional :: shannon
    integer, intent(in), optional :: degrees(:)
    integer, intent(out), optional :: ntapers
    integer, intent(out), optional :: exitstatus
    integer :: l, n, n_int, j, i, astat(3), ind(lmax+1), use_dg82
    real(dp) :: eval(lmax+1), pi, upper, lower, zero(lmax+1), w(lmax+1), h
    real(dp), allocatable :: evec(:,:), dllmtri(:,:), p(:), dllm(:,:)

    use_dg82 = 1

    if (present(exitstatus)) exitstatus = 0

    if (size(tapers(:,1)) < (lmax+1) .or. &
            size(tapers(1,:)) < (lmax+1-abs(m)) ) then
        print*, "Error --- SHReturnTapersM"
        print*, "TAPERS must be dimensioned as (LMAX+1, LMAX+1-ABS(M)) " // &
                "where LMAX is ", lmax
        print*, "Input array is dimensioned as ", size(tapers(:,1)), &
                size(tapers(1,:))
        if (present(exitstatus)) then
            exitstatus = 1
            return
        else
            stop
        end if

    else if(size(eigenvalues) < (lmax+1-abs(m)) ) then
        print*, "Error --- SHReturnTapersM"
        print*, "EIGENVALUES must be dimensioned as (LMAX+1-ABS(m)) " // &
                "where LMAX is ", lmax
        print*, "Input array is dimensioned as ", size(eigenvalues)
        if (present(exitstatus)) then
            exitstatus = 1
            return
        else
            stop
        end if

    else if (m > lmax) then
        print*, "Error --- SHReturnTapersM"
        print*, "M must be less than or equal to LMAX."
        print*, "M = ", m
        print*, "LMAX = ", lmax
        if (present(exitstatus)) then
            exitstatus = 2
            return
        else
            stop
        end if

    end if

    if (present(degrees)) then
        if (size(degrees) < lmax+1) then
            print*, "Error --- SHReturnTapersM"
            print*, "DEGREES must have dimension LMAX+1, where LMAX is ", lmax
            print*, "Input array is dimensioned as ", size(degrees)
            if (present(exitstatus)) then
                exitstatus = 1
                return
            else
                stop
            end if
        end if
    end if

    if (present(degrees)) then
        if (sum(degrees(1+abs(m):lmax+1)) /= lmax + 1 - abs(m)) use_dg82 = 0
    end if

    allocate (evec(lmax+1, lmax+1-abs(m)), stat = astat(1))
    allocate (dllmtri(lmax+1, lmax+1), stat = astat(2))
    allocate (p((lmax+1)*(lmax+2)/2), stat = astat(3))

    if (astat(1) /= 0 .or. astat(2) /= 0 .or. astat(3) /= 0) then
        print*, "Error --- SHReturnTapersM"
        print*, "Problem allocating arrays EVEC, DLLMTRI, and P", &
            astat(1), astat(2), astat(3)
        if (present(exitstatus)) then
            exitstatus = 3
            return
        else
            stop
        end if

    end if

    pi = acos(-1.0_dp)

    tapers = 0.0_dp
    eigenvalues = 0.0_dp
    eval = 0.0_dp
    evec = 0.0_dp

    if (use_dg82 == 1) then
        !----------------------------------------------------------------------
        !
        !   Calculate space-concentration kernel and the corresponding
        !   eigenfunctions of the Grunbaum et al. kernel.
        !
        !----------------------------------------------------------------------
        n = lmax + 1 - abs(m)

        if (present(exitstatus)) then
            call ComputeDG82(dllmtri(1:n,1:n), lmax, m, theta0, &
                             exitstatus = exitstatus)
        if (exitstatus /= 0) return
            call EigValVecSymTri(dllmtri(1:n,1:n), n, eval(1:n), &
                                 tapers(1+abs(m):lmax+1,1:n), &
                                 exitstatus = exitstatus)
            if (exitstatus /= 0) return
        else
            call ComputeDG82(dllmtri(1:n,1:n), lmax, m, theta0)
            call EigValVecSymTri(dllmtri(1:n,1:n), n, eval(1:n), &
                                 tapers(1+abs(m):lmax+1,1:n))
        end if

        !----------------------------------------------------------------------
        !
        !   Calculate true eigenvalues
        !
        !----------------------------------------------------------------------
        upper = 1.0_dp
        lower = cos(theta0)
        n_int = lmax + 1

        if (present(exitstatus)) then
            call PreGLQ(lower, upper, n_int, zero, w, exitstatus = exitstatus)
            if (exitstatus /= 0) return
        else
            call PreGLQ(lower, upper, n_int, zero, w)
        end if

        do i=1, n_int
            if (present(exitstatus)) then
                call PlmBar(p, lmax, zero(i), exitstatus = exitstatus)
                if (exitstatus /= 0) return
            else
                call PlmBar(p, lmax, zero(i))
            end if

            do j = 1, n
                h = 0.0_dp

                do l = abs(m), lmax
                    h = h + p(PlmIndex(l, abs(m))) * tapers(l+1, j)

                end do

                eigenvalues(j) = eigenvalues(j) + w(i) * h**2

            end do
        end do

        if (m == 0) then
            eigenvalues(1:n) = eigenvalues(1:n) / 2.0_dp

        else
            eigenvalues(1:n) = eigenvalues(1:n) / 4.0_dp

        end if

    else
        !----------------------------------------------------------------------
        !
        !   Calculate space-concentration kernel and the corresponding
        !   eigenfunctions.
        !
        !----------------------------------------------------------------------
        ind = 0
        i = 0

        do l=abs(m), lmax, 1
            if (degrees(l+1) /= 0) then
                i = i + 1
                ind(i) = l + 1
            end if
        end do

        n = i

        if (n /= 0) then

            allocate (dllm(n, n), stat = astat(1))
            if (astat(1) /= 0) then
                print*, "Error --- SHReturnTapersM"
                print*, "Problem allocating array DLLM ", astat(1)
                if (present(exitstatus)) then
                    exitstatus = 3
                    return
                else
                    stop
                end if
            end if

            dllm = 0.0_dp

            if (present(exitstatus)) then
                call ComputeDm(dllmtri, lmax, m, theta0, degrees=degrees, &
                               exitstatus=exitstatus)
                if (exitstatus /= 0) return
            else
                call ComputeDm(dllmtri, lmax, m, theta0, degrees=degrees)
            end if

            do i=1, n
                do j=1, n
                    dllm(i, j) = dllmtri(ind(i), ind(j))
                end do
            end do

            if (present(exitstatus)) then
                call EigValVecSym(dllm(1:n, 1:n), n, &
                                  eval(1:n), evec(1:n,1:n), &
                                exitstatus = exitstatus)
                if (exitstatus /= 0) return
            else
                call EigValVecSym(dllm, n, eval(1:n),&
                                  evec(1:n,1:n))
            end if

            do i=1, n
                tapers(ind(i), 1:n) = evec(i, 1:n)
            end do

            eigenvalues(1:n) = eval(1:n)

        end if

    end if

    if (present(shannon)) then
        shannon = sum(eigenvalues(1:lmax+1-abs(m)))
    end if

    if (present(ntapers)) then
        ntapers = n
    end if

    ! deallocate memory
    call PlmBar (p, -1, zero(1))
    deallocate (evec)
    deallocate (dllmtri)
    deallocate (p)
    if (use_dg82 == 0 .and. n /= 0) deallocate (dllm)

end subroutine SHReturnTapersM

max_stars_repo_name: dilkins/SHTOOLS

max_issues_repo_name: dilkins/SHTOOLS

max_forks_repo_name: dilkins/SHTOOLS