text:
MODULE ghex_defs
  integer, public, parameter :: GHEX_ANY_SOURCE = -1

  integer, public, parameter :: ghex_fp_kind = GHEX_FORTRAN_FP_KIND

  integer, public, parameter :: GhexDeviceUnknown = 0
  integer, public, parameter :: GhexDeviceCPU = 1
  integer, public, parameter :: GhexDeviceGPU = 2

  integer, public, parameter :: GhexLayoutFieldLast  = 1
  integer, public, parameter :: GhexLayoutFieldFirst = 2

  integer, public, parameter :: GhexBarrierGlobal = 1
  integer, public, parameter :: GhexBarrierThread = 2
  integer, public, parameter :: GhexBarrierRank   = 3

  integer, public, parameter :: GhexAllocatorHost    = 1
  integer, public, parameter :: GhexAllocatorDevice  = 2

END MODULE ghex_defs

max_stars_repo_name: tehrengruber/GHEX

max_issues_repo_name: tehrengruber/GHEX

max_forks_repo_name: tehrengruber/GHEX


text:
!--------------------------------------------------
!PHY1038 Assignment 2 - Series Computation
!URN 6309823 - Penguin Lab Group 1K
!March 10th 2015
!--------------------------------------------------

!A Taylor approximation will be produced which will look close to the actual function for the range of x values suggested
!However there will be some noticeable deviations and the ends of the range.
!This program, will output to a file called "taylor.dat".

!Compile this program using the command gfortran

PROGRAM Series_Computation

  IMPLICIT NONE
  REAL :: a,x,e,dx,pi,factorial
  INTEGER :: n,i,j,num

  OPEN(unit=20,file='taylor.dat')

  WRITE(6,*) "This program will calculate the Taylor Series expansion of the function f(x)=sin(x)*exp(-x/2)."
  WRITE(6,*) " "
  WRITE(6,*) "This describes an oscillation whose amplitude is exponentially decaying"
  WRITE(6,*) " "
  WRITE(6,*) " "
  WRITE(6,*) "Specify the order (n) of the Taylor approximation"
  !A reasonable value is 15
  READ(5,*) n
  WRITE(6,*) " "
  WRITE(6,*) "Specify the number of x points"
  !A reasonable value is 1000
  READ(5,*) num
  WRITE(6,*) " "

  pi=acos(-1.0)

  dx=4.0*pi/num

  !Where dx is the step value

  DO j=0,num
    e=0.0
    x=-2.0*pi+dx*j+1e-6
    a=sin(x)*exp(-x/2)
    factorial=1
    DO i=1,n

      !This DO loop calculates the taylor series values

      factorial=factorial*i
      e=((((-sqrt(5.0)/2.0)**i)*(sin(-i*atan(2.0))))/factorial)*x**i+e
    END DO
    WRITE (20,*) x,a,e

    !x is the values of the variable x
    !a is the value of produced by the original function
    !e is the value of each order of the taylor series

  END DO

  WRITE(6,*) "The data is writen to the file taylor.dat in the format:"
  WRITE(6,*) "x, f(x), Taylor Expansion"

  !Write in Terminal "gnuplot" then press ENTER
  !Write in Terminal "plot "taylor.dat" using 1:2 w l, "" using 1:3 w l" and press ENTER
  !Write in Terminal "set yrange [-12,12]" and press ENTER
  !Write in Terminal "replot" and press ENTER

END PROGRAM Series_Computation

max_stars_repo_name: WilliamHoltam/computational-physics

max_issues_repo_name: WilliamHoltam/computational-physics

max_forks_repo_name: WilliamHoltam/computational-physics


text:
C ****************************************************************************
C FILE: mpi_struct.f
C DESCRIPTION:
C   MPI tutorial example code: Structure Derived Datatype
C AUTHOR: Blaise Barney
C LAST REVISED: 04/02/05
C ****************************************************************************

      program struct
      include 'mpif.h'

      integer NELEM
      parameter(NELEM=25)
      integer numtasks, rank, source, dest, tag, i,  ierr
      integer stat(MPI_STATUS_SIZE)

      type Particle
      sequence
      real*4 x, y, z, velocity
      integer n, type
      end type Particle

      type (Particle) p(NELEM), particles(NELEM)
      integer particletype, oldtypes(0:1), blockcounts(0:1),
     &        offsets(0:1), extent

      call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, numtasks, ierr)

C     Setup description of the 4 MPI_REAL fields x, y, z, velocity
      offsets(0) = 0
      oldtypes(0) = MPI_REAL
      blockcounts(0) = 4

C     Setup description of the 2 MPI_INTEGER fields n, type
C     Need to first figure offset by getting size of MPI_REAL
      call MPI_TYPE_EXTENT(MPI_REAL, extent, ierr)
      offsets(1) = 4 * extent
      oldtypes(1) = MPI_INTEGER
      blockcounts(1) = 2

C     Now define structured type and commit it
      call MPI_TYPE_STRUCT(2, blockcounts, offsets, oldtypes,
     &                     particletype, ierr)
      call MPI_TYPE_COMMIT(particletype, ierr)

C     Initialize the particle array and then send it to each task
      tag = 1
      if (rank .eq. 0) then
         do 10 i=0, NELEM-1
         particles(i) = Particle ( 1.0*i, -1.0*i, 1.0*i,
     &                  0.25, i, mod(i,2) )
  10     continue

         do 20 i=0, numtasks-1
         call MPI_SEND(particles, NELEM, particletype, i, tag,
     &                 MPI_COMM_WORLD, ierr)
  20     continue
      endif

      source = 0
      call MPI_RECV(p, NELEM, particletype, source, tag,
     &              MPI_COMM_WORLD, stat, ierr)

      print *, 'rank= ',rank,' p(3)= ',p(3)
      call MPI_TYPE_FREE(particletype, ierr)
      call MPI_FINALIZE(ierr)
      end

ext: f

max_stars_repo_name: hixio-mh/HPC-Tutorials

max_issues_repo_name: hixio-mh/HPC-Tutorials

max_forks_repo_name: hixio-mh/HPC-Tutorials


text:
      MODULE nmtejj_I
      INTERFACE
!
      INTEGER FUNCTION NMTEJJ(I2Q,I2J,J,NK,ND)
      INTEGER, INTENT(IN) :: I2Q, I2J, J, NK, ND
      END FUNCTION
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
	SUBROUTINE GH_RDPB ( strnam, advno, wtype, timstr,
     +                       wname, wadnm, wwnd, gust, wdir, wsped,
     +                       wpres, rlat, rlon, iret )
C************************************************************************
C* GH_RDPB								*
C*									*
C* This subroutine reads the public advisory text bulletin for          *
C* tropical storms.							*
C*									*
C* GH_RDPB ( STRNAM, ADVNO, WTYPE, TIMSTR, WNAME, WADNM, 		*
C*           WWND, WDIR, WSPED, WPRES, RLAT, RLON,  IRET )              *
C*									*
C* Input parameters:							*
C*	STRNAM		CHAR*		Tropical storm name		*
C*	ADVNO		CHAR*           Advisory number			*
C*									*
C* Output parameters:							*
C* 	WTYPE 		CHAR*		Tropical storm type		*
C* 	TIMSTR 		CHAR*		Date/time time string 		*
C*	WNAME 		CHAR*		Tropical storm name		*
C*	WADNM 		CHAR*		Advisory number			*
C*	WWND 		CHAR*		Max sustained winds		*
C*	GUST 		CHAR*		Wind gust			*
C*	WDIR 		CHAR*		Movement direction		*
C*	WSPED 		CHAR*		Movement speed (kts)		*
C*	WPRES 		CHAR*		Central minimum pressure (mb)	*
C*	RLAT 		REAL		Current/forecasted Latitudes	*
C*	RLON 		REAL		Current/forecasted Longitudes	*
C*	IRET		INTEGER		Return code			*
C*                                        0 = normal return             *
C*                                       -1 = could not open public     *
C*                                            advisory                  *
C*                                       -2 = could not read public     *
C*                                            advisory header           *
C*                                       -3 = could not find public     *
C*                                            advisory                  *
C*									*
C**									*
C* Log:									*
C* S. Gilbert/NCEP	 1/06	Modified from GH_RDAD                   *
C* m.gamazaychikov/SAIC 04/06   Added idtmch flag to CTB_DTGET CS       *
C* S. Gilbert/NCEP	 6/06	Added call to FL_SCND to get file list  *
C* F. J. Yen/NCEP        4/08   Added bin mins & mstrct to CTB_DTGET CSC*
C* C. Lauer/NHC          3/09   Added call to GH_PSUM for new pub format*
C* X. Guo/CWS		02/10   Changed GH_PSUM to support the new	*
C*                              public format and called GH_PSUM09 to   *
C*                              read the old format                     *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	strnam, advno, wtype, wname, wpres,
     +                  wadnm,  wwnd, wdir, wsped, timstr, gust
        REAL            rlat, rlon
C*
	CHARACTER*(MXFLSZ)      filnam, files (MXNMFL)
	CHARACTER	path*25, templ*25, tfile*128,
     +                  adnum*4, record*10000, rec*128
C*
	INCLUDE		'ERMISS.FNC'
C-----------------------------------------------------------------------
	iret = 0
C
C*      Get path to public advisory files
C
        filnam = 'HCNPUB'
        path  = ' '
        templ = ' '
        CALL ST_NULL ( filnam, filnam, nf, ier )
	CALL CTB_DTGET ( filnam, path, templ, ic, is, if, ir, ii, ion,
     +          ihb, mnb, iha, mna, mstrct, idtmch, ier )
        CALL ST_RNUL ( path, path, lenp, ier )
        CALL ST_RNUL ( templ, templ, lent, ier )
C
C*      Get list of public advisory files
C
        CALL ST_NULL ( advno, advno, lena, ier )
        CALL GH_ADVN( advno, iadnm, iflag, ier )
        WRITE ( adnum, fmt='(I3.3)' ) iadnm
        IF ( iflag .eq. 1 ) THEN
            CALL ST_RPST ( templ, 'X', 'a', ipos, templ, ier )
        ELSE IF ( iflag .eq. 2 ) THEN
            CALL ST_RPST ( templ, 'X', 'b', ipos, templ, ier )
        END IF
        nexp   = MXNMFL
        iorder = 1
        CALL FL_SCND ( path, templ, iorder, nexp, files, nfiles, ier )
        tfile = ' '
        DO ifl = 1, nfiles
           IF ( files(ifl)(1:8) .eq. strnam .AND.
     +          files(ifl)(19:21) .eq. adnum ) THEN
              tfile = path(:lenp) // '/' // files(ifl)
           ENDIF
        ENDDO
        IF ( tfile .eq. ' ' ) THEN
            iret = -3
            RETURN
        ENDIF
C
C*	Open/read public advisory
C
        CALL FL_SOPN ( tfile, lunf, ier )
        IF ( ier .ne. 0 ) THEN
            iret = -1
            RETURN
        ENDIF
C
        lrec = 0
        iostat = 0
        DO WHILE  ( iostat .eq. 0 )
            READ  ( lunf, 1000, IOSTAT = iostat ) rec
 1000       FORMAT (A)
            CALL ST_LSTR( rec, lent, ier )
            IF ( lrec .eq. 0 ) THEN
                lrec = lent
                record = rec (:lent)
            ELSE
                CALL ST_LSTR( record, lrec, ier )
                record = record(:lrec) // ' ' // rec (:lent)
                CALL ST_LSTR( record, lrec, ier )
            END IF
        END DO
        CALL FL_CLOS ( lunf, ier )
C
	CALL ST_LSTR ( record, lrec1, ier )
	CALL ST_UNPR ( record, lrec1, record, lrec, ier )
	CALL ST_LCUC ( record, record, ier )
C
C*      Decode the advisory header lines.
C
        CALL GH_PUHD ( record, wtype, wname, wadnm, timstr, ier )
        IF ( ier .ne. 0 ) THEN
            iret = -2
        END IF
C
C*      Find and decode the Repeat Section. (GH_PRPT is old format)
C
 	CALL GH_PSUM ( record, lrec, rlat, rlon, wdir, wsped,
     +                 wpres, wwnd, gust, ier )
        IF ( ier .ne. 0 ) THEN
            CALL GH_PSUM09 ( record, lrec, rlat, rlon, wdir, wsped,
     +                 wpres, wwnd, gust, ier )
	    IF ( ier .ne. 0 ) THEN
	    	iret = -3
	    END IF

        END IF
C
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


!*==cbdsqr.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022

!> \brief \b CBDSQR
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download CBDSQR + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
!                          LDU, C, LDC, RWORK, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          UPLO
!       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
!       ..
!       .. Array Arguments ..
!       REAL               D( * ), E( * ), RWORK( * )
!       COMPLEX            C( LDC, * ), U( LDU, * ), VT( LDVT, * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> CBDSQR computes the singular values and, optionally, the right and/or
!> left singular vectors from the singular value decomposition (SVD) of
!> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
!> zero-shift QR algorithm.  The SVD of B has the form
!>
!>    B = Q * S * P**H
!>
!> where S is the diagonal matrix of singular values, Q is an orthogonal
!> matrix of left singular vectors, and P is an orthogonal matrix of
!> right singular vectors.  If left singular vectors are requested, this
!> subroutine actually returns U*Q instead of Q, and, if right singular
!> vectors are requested, this subroutine returns P**H*VT instead of
!> P**H, for given complex input matrices U and VT.  When U and VT are
!> the unitary matrices that reduce a general matrix A to bidiagonal
!> form: A = U*B*VT, as computed by CGEBRD, then
!>
!>    A = (U*Q) * S * (P**H*VT)
!>
!> is the SVD of A.  Optionally, the subroutine may also compute Q**H*C
!> for a given complex input matrix C.
!>
!> See "Computing  Small Singular Values of Bidiagonal Matrices With
!> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
!> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
!> no. 5, pp. 873-912, Sept 1990) and
!> "Accurate singular values and differential qd algorithms," by
!> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
!> Department, University of California at Berkeley, July 1992
!> for a detailed description of the algorithm.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] UPLO
!> \verbatim
!>          UPLO is CHARACTER*1
!>          = 'U':  B is upper bidiagonal;
!>          = 'L':  B is lower bidiagonal.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix B.  N >= 0.
!> \endverbatim
!>
!> \param[in] NCVT
!> \verbatim
!>          NCVT is INTEGER
!>          The number of columns of the matrix VT. NCVT >= 0.
!> \endverbatim
!>
!> \param[in] NRU
!> \verbatim
!>          NRU is INTEGER
!>          The number of rows of the matrix U. NRU >= 0.
!> \endverbatim
!>
!> \param[in] NCC
!> \verbatim
!>          NCC is INTEGER
!>          The number of columns of the matrix C. NCC >= 0.
!> \endverbatim
!>
!> \param[in,out] D
!> \verbatim
!>          D is REAL array, dimension (N)
!>          On entry, the n diagonal elements of the bidiagonal matrix B.
!>          On exit, if INFO=0, the singular values of B in decreasing
!>          order.
!> \endverbatim
!>
!> \param[in,out] E
!> \verbatim
!>          E is REAL array, dimension (N-1)
!>          On entry, the N-1 offdiagonal elements of the bidiagonal
!>          matrix B.
!>          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
!>          will contain the diagonal and superdiagonal elements of a
!>          bidiagonal matrix orthogonally equivalent to the one given
!>          as input.
!> \endverbatim
!>
!> \param[in,out] VT
!> \verbatim
!>          VT is COMPLEX array, dimension (LDVT, NCVT)
!>          On entry, an N-by-NCVT matrix VT.
!>          On exit, VT is overwritten by P**H * VT.
!>          Not referenced if NCVT = 0.
!> \endverbatim
!>
!> \param[in] LDVT
!> \verbatim
!>          LDVT is INTEGER
!>          The leading dimension of the array VT.
!>          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
!> \endverbatim
!>
!> \param[in,out] U
!> \verbatim
!>          U is COMPLEX array, dimension (LDU, N)
!>          On entry, an NRU-by-N matrix U.
!>          On exit, U is overwritten by U * Q.
!>          Not referenced if NRU = 0.
!> \endverbatim
!>
!> \param[in] LDU
!> \verbatim
!>          LDU is INTEGER
!>          The leading dimension of the array U.  LDU >= max(1,NRU).
!> \endverbatim
!>
!> \param[in,out] C
!> \verbatim
!>          C is COMPLEX array, dimension (LDC, NCC)
!>          On entry, an N-by-NCC matrix C.
!>          On exit, C is overwritten by Q**H * C.
!>          Not referenced if NCC = 0.
!> \endverbatim
!>
!> \param[in] LDC
!> \verbatim
!>          LDC is INTEGER
!>          The leading dimension of the array C.
!>          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
!> \endverbatim
!>
!> \param[out] RWORK
!> \verbatim
!>          RWORK is REAL array, dimension (4*N)
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  If INFO = -i, the i-th argument had an illegal value
!>          > 0:  the algorithm did not converge; D and E contain the
!>                elements of a bidiagonal matrix which is orthogonally
!>                similar to the input matrix B;  if INFO = i, i
!>                elements of E have not converged to zero.
!> \endverbatim
!
!> \par Internal Parameters:
!  =========================
!>
!> \verbatim
!>  TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8)))
!>          TOLMUL controls the convergence criterion of the QR loop.
!>          If it is positive, TOLMUL*EPS is the desired relative
!>             precision in the computed singular values.
!>          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
!>             desired absolute accuracy in the computed singular
!>             values (corresponds to relative accuracy
!>             abs(TOLMUL*EPS) in the largest singular value.
!>          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
!>             between 10 (for fast convergence) and .1/EPS
!>             (for there to be some accuracy in the results).
!>          Default is to lose at either one eighth or 2 of the
!>             available decimal digits in each computed singular value
!>             (whichever is smaller).
!>
!>  MAXITR  INTEGER, default = 6
!>          MAXITR controls the maximum number of passes of the
!>          algorithm through its inner loop. The algorithms stops
!>          (and so fails to converge) if the number of passes
!>          through the inner loop exceeds MAXITR*N**2.
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup complexOTHERcomputational
!
!  =====================================================================
      SUBROUTINE CBDSQR(Uplo,N,Ncvt,Nru,Ncc,D,E,Vt,Ldvt,U,Ldu,C,Ldc,    &
     &                  Rwork,Info)
      IMPLICIT NONE
!*--CBDSQR227
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Uplo
      INTEGER Info , Ldc , Ldu , Ldvt , N , Ncc , Ncvt , Nru
!     ..
!     .. Array Arguments ..
      REAL D(*) , E(*) , Rwork(*)
      COMPLEX C(Ldc,*) , U(Ldu,*) , Vt(Ldvt,*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      REAL ZERO
      PARAMETER (ZERO=0.0E0)
      REAL ONE
      PARAMETER (ONE=1.0E0)
      REAL NEGONE
      PARAMETER (NEGONE=-1.0E0)
      REAL HNDRTH
      PARAMETER (HNDRTH=0.01E0)
      REAL TEN
      PARAMETER (TEN=10.0E0)
      REAL HNDRD
      PARAMETER (HNDRD=100.0E0)
      REAL MEIGTH
      PARAMETER (MEIGTH=-0.125E0)
      INTEGER MAXITR
      PARAMETER (MAXITR=6)
!     ..
!     .. Local Scalars ..
      LOGICAL lower , rotate
      INTEGER i , idir , isub , iter , j , ll , lll , m , maxit , nm1 , &
     &        nm12 , nm13 , oldll , oldm
      REAL abse , abss , cosl , cosr , cs , eps , f , g , h , mu ,      &
     &     oldcs , oldsn , r , shift , sigmn , sigmx , sinl , sinr ,    &
     &     sll , smax , smin , sminl , sminoa , sn , thresh , tol ,     &
     &     tolmul , unfl
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      REAL SLAMCH
      EXTERNAL LSAME , SLAMCH
!     ..
!     .. External Subroutines ..
      EXTERNAL CLASR , CSROT , CSSCAL , CSWAP , SLARTG , SLAS2 ,        &
     &         SLASQ1 , SLASV2 , XERBLA
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , MAX , MIN , REAL , SIGN , SQRT
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters.
!
      Info = 0
      lower = LSAME(Uplo,'L')
      IF ( .NOT.LSAME(Uplo,'U') .AND. .NOT.lower ) THEN
         Info = -1
      ELSEIF ( N<0 ) THEN
         Info = -2
      ELSEIF ( Ncvt<0 ) THEN
         Info = -3
      ELSEIF ( Nru<0 ) THEN
         Info = -4
      ELSEIF ( Ncc<0 ) THEN
         Info = -5
      ELSEIF ( (Ncvt==0 .AND. Ldvt<1) .OR. (Ncvt>0 .AND. Ldvt<MAX(1,N)) &
     &         ) THEN
         Info = -9
      ELSEIF ( Ldu<MAX(1,Nru) ) THEN
         Info = -11
      ELSEIF ( (Ncc==0 .AND. Ldc<1) .OR. (Ncc>0 .AND. Ldc<MAX(1,N)) )   &
     &         THEN
         Info = -13
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('CBDSQR',-Info)
         RETURN
      ENDIF
      IF ( N==0 ) RETURN
      IF ( N==1 ) GOTO 300
!
!     ROTATE is true if any singular vectors desired, false otherwise
!
      rotate = (Ncvt>0) .OR. (Nru>0) .OR. (Ncc>0)
!
!     If no singular vectors desired, use qd algorithm
!
      IF ( .NOT.rotate ) THEN
         CALL SLASQ1(N,D,E,Rwork,Info)
!
!     If INFO equals 2, dqds didn't finish, try to finish
!
         IF ( Info/=2 ) RETURN
         Info = 0
      ENDIF
!
      nm1 = N - 1
      nm12 = nm1 + nm1
      nm13 = nm12 + nm1
      idir = 0
!
!     Get machine constants
!
      eps = SLAMCH('Epsilon')
      unfl = SLAMCH('Safe minimum')
!
!     If matrix lower bidiagonal, rotate to be upper bidiagonal
!     by applying Givens rotations on the left
!
      IF ( lower ) THEN
         DO i = 1 , N - 1
            CALL SLARTG(D(i),E(i),cs,sn,r)
            D(i) = r
            E(i) = sn*D(i+1)
            D(i+1) = cs*D(i+1)
            Rwork(i) = cs
            Rwork(nm1+i) = sn
         ENDDO
!
!        Update singular vectors if desired
!
         IF ( Nru>0 ) CALL CLASR('R','V','F',Nru,N,Rwork(1),Rwork(N),U, &
     &                           Ldu)
         IF ( Ncc>0 ) CALL CLASR('L','V','F',N,Ncc,Rwork(1),Rwork(N),C, &
     &                           Ldc)
      ENDIF
!
!     Compute singular values to relative accuracy TOL
!     (By setting TOL to be negative, algorithm will compute
!     singular values to absolute accuracy ABS(TOL)*norm(input matrix))
!
      tolmul = MAX(TEN,MIN(HNDRD,eps**MEIGTH))
      tol = tolmul*eps
!
!     Compute approximate maximum, minimum singular values
!
      smax = ZERO
      DO i = 1 , N
         smax = MAX(smax,ABS(D(i)))
      ENDDO
      DO i = 1 , N - 1
         smax = MAX(smax,ABS(E(i)))
      ENDDO
      sminl = ZERO
      IF ( tol>=ZERO ) THEN
!
!        Relative accuracy desired
!
         sminoa = ABS(D(1))
         IF ( sminoa/=ZERO ) THEN
            mu = sminoa
            DO i = 2 , N
               mu = ABS(D(i))*(mu/(mu+ABS(E(i-1))))
               sminoa = MIN(sminoa,mu)
               IF ( sminoa==ZERO ) EXIT
            ENDDO
         ENDIF
         sminoa = sminoa/SQRT(REAL(N))
         thresh = MAX(tol*sminoa,MAXITR*N*N*unfl)
      ELSE
!
!        Absolute accuracy desired
!
         thresh = MAX(ABS(tol)*smax,MAXITR*N*N*unfl)
      ENDIF
!
!     Prepare for main iteration loop for the singular values
!     (MAXIT is the maximum number of passes through the inner
!     loop permitted before nonconvergence signalled.)
!
      maxit = MAXITR*N*N
      iter = 0
      oldll = -1
      oldm = -1
!
!     M points to last element of unconverged part of matrix
!
      m = N
!
!     Begin main iteration loop
!
!
!     Check for convergence or exceeding iteration count
!
 100  DO WHILE ( m>1 )
         IF ( iter>maxit ) THEN
!
!     Maximum number of iterations exceeded, failure to converge
!
            Info = 0
            DO i = 1 , N - 1
               IF ( E(i)/=ZERO ) Info = Info + 1
            ENDDO
            GOTO 99999
         ELSE
!
!     Find diagonal block of matrix to work on
!
            IF ( tol<ZERO .AND. ABS(D(m))<=thresh ) D(m) = ZERO
            smax = ABS(D(m))
            smin = smax
            DO lll = 1 , m - 1
               ll = m - lll
               abss = ABS(D(ll))
               abse = ABS(E(ll))
               IF ( tol<ZERO .AND. abss<=thresh ) D(ll) = ZERO
               IF ( abse<=thresh ) GOTO 150
               smin = MIN(smin,abss)
               smax = MAX(smax,abss,abse)
            ENDDO
            ll = 0
            GOTO 200
         ENDIF
 150     E(ll) = ZERO
!
!     Matrix splits since E(LL) = 0
!
         IF ( ll==m-1 ) THEN
!
!        Convergence of bottom singular value, return to top of loop
!
            m = m - 1
            CYCLE
         ENDIF
         GOTO 200
      ENDDO
      GOTO 300
 200  ll = ll + 1
!
!     E(LL) through E(M-1) are nonzero, E(LL-1) is zero
!
      IF ( ll==m-1 ) THEN
!
!        2 by 2 block, handle separately
!
         CALL SLASV2(D(m-1),E(m-1),D(m),sigmn,sigmx,sinr,cosr,sinl,cosl)
         D(m-1) = sigmx
         E(m-1) = ZERO
         D(m) = sigmn
!
!        Compute singular vectors, if desired
!
         IF ( Ncvt>0 ) CALL CSROT(Ncvt,Vt(m-1,1),Ldvt,Vt(m,1),Ldvt,cosr,&
     &                            sinr)
         IF ( Nru>0 ) CALL CSROT(Nru,U(1,m-1),1,U(1,m),1,cosl,sinl)
         IF ( Ncc>0 ) CALL CSROT(Ncc,C(m-1,1),Ldc,C(m,1),Ldc,cosl,sinl)
         m = m - 2
         GOTO 100
      ENDIF
!
!     If working on new submatrix, choose shift direction
!     (from larger end diagonal element towards smaller)
!
      IF ( ll>oldm .OR. m<oldll ) THEN
         IF ( ABS(D(ll))>=ABS(D(m)) ) THEN
!
!           Chase bulge from top (big end) to bottom (small end)
!
            idir = 1
         ELSE
!
!           Chase bulge from bottom (big end) to top (small end)
!
            idir = 2
         ENDIF
      ENDIF
!
!     Apply convergence tests
!
      IF ( idir==1 ) THEN
!
!        Run convergence test in forward direction
!        First apply standard test to bottom of matrix
!
         IF ( ABS(E(m-1))<=ABS(tol)*ABS(D(m)) .OR.                      &
     &        (tol<ZERO .AND. ABS(E(m-1))<=thresh) ) THEN
            E(m-1) = ZERO
            GOTO 100
         ENDIF
!
         IF ( tol>=ZERO ) THEN
!
!           If relative accuracy desired,
!           apply convergence criterion forward
!
            mu = ABS(D(ll))
            sminl = mu
            DO lll = ll , m - 1
               IF ( ABS(E(lll))<=tol*mu ) THEN
                  E(lll) = ZERO
                  GOTO 100
               ENDIF
               mu = ABS(D(lll+1))*(mu/(mu+ABS(E(lll))))
               sminl = MIN(sminl,mu)
            ENDDO
         ENDIF
!
      ELSE
!
!        Run convergence test in backward direction
!        First apply standard test to top of matrix
!
         IF ( ABS(E(ll))<=ABS(tol)*ABS(D(ll)) .OR.                      &
     &        (tol<ZERO .AND. ABS(E(ll))<=thresh) ) THEN
            E(ll) = ZERO
            GOTO 100
         ENDIF
!
         IF ( tol>=ZERO ) THEN
!
!           If relative accuracy desired,
!           apply convergence criterion backward
!
            mu = ABS(D(m))
            sminl = mu
            DO lll = m - 1 , ll , -1
               IF ( ABS(E(lll))<=tol*mu ) THEN
                  E(lll) = ZERO
                  GOTO 100
               ENDIF
               mu = ABS(D(lll))*(mu/(mu+ABS(E(lll))))
               sminl = MIN(sminl,mu)
            ENDDO
         ENDIF
      ENDIF
      oldll = ll
      oldm = m
!
!     Compute shift.  First, test if shifting would ruin relative
!     accuracy, and if so set the shift to zero.
!
      IF ( tol>=ZERO .AND. N*tol*(sminl/smax)<=MAX(eps,HNDRTH*tol) )    &
     &     THEN
!
!        Use a zero shift to avoid loss of relative accuracy
!
         shift = ZERO
      ELSE
!
!        Compute the shift from 2-by-2 block at end of matrix
!
         IF ( idir==1 ) THEN
            sll = ABS(D(ll))
            CALL SLAS2(D(m-1),E(m-1),D(m),shift,r)
         ELSE
            sll = ABS(D(m))
            CALL SLAS2(D(ll),E(ll),D(ll+1),shift,r)
         ENDIF
!
!        Test if shift negligible, and if so set to zero
!
         IF ( sll>ZERO ) THEN
            IF ( (shift/sll)**2<eps ) shift = ZERO
         ENDIF
      ENDIF
!
!     Increment iteration count
!
      iter = iter + m - ll
!
!     If SHIFT = 0, do simplified QR iteration
!
      IF ( shift==ZERO ) THEN
         IF ( idir==1 ) THEN
!
!           Chase bulge from top to bottom
!           Save cosines and sines for later singular vector updates
!
            cs = ONE
            oldcs = ONE
            DO i = ll , m - 1
               CALL SLARTG(D(i)*cs,E(i),cs,sn,r)
               IF ( i>ll ) E(i-1) = oldsn*r
               CALL SLARTG(oldcs*r,D(i+1)*sn,oldcs,oldsn,D(i))
               Rwork(i-ll+1) = cs
               Rwork(i-ll+1+nm1) = sn
               Rwork(i-ll+1+nm12) = oldcs
               Rwork(i-ll+1+nm13) = oldsn
            ENDDO
            h = D(m)*cs
            D(m) = h*oldcs
            E(m-1) = h*oldsn
!
!           Update singular vectors
!
            IF ( Ncvt>0 ) CALL CLASR('L','V','F',m-ll+1,Ncvt,Rwork(1),  &
     &                               Rwork(N),Vt(ll,1),Ldvt)
            IF ( Nru>0 ) CALL CLASR('R','V','F',Nru,m-ll+1,Rwork(nm12+1)&
     &                              ,Rwork(nm13+1),U(1,ll),Ldu)
            IF ( Ncc>0 ) CALL CLASR('L','V','F',m-ll+1,Ncc,Rwork(nm12+1)&
     &                              ,Rwork(nm13+1),C(ll,1),Ldc)
!
!           Test convergence
!
            IF ( ABS(E(m-1))<=thresh ) E(m-1) = ZERO
!
         ELSE
!
!           Chase bulge from bottom to top
!           Save cosines and sines for later singular vector updates
!
            cs = ONE
            oldcs = ONE
            DO i = m , ll + 1 , -1
               CALL SLARTG(D(i)*cs,E(i-1),cs,sn,r)
               IF ( i<m ) E(i) = oldsn*r
               CALL SLARTG(oldcs*r,D(i-1)*sn,oldcs,oldsn,D(i))
               Rwork(i-ll) = cs
               Rwork(i-ll+nm1) = -sn
               Rwork(i-ll+nm12) = oldcs
               Rwork(i-ll+nm13) = -oldsn
            ENDDO
            h = D(ll)*cs
            D(ll) = h*oldcs
            E(ll) = h*oldsn
!
!           Update singular vectors
!
            IF ( Ncvt>0 ) CALL CLASR('L','V','B',m-ll+1,Ncvt,           &
     &                               Rwork(nm12+1),Rwork(nm13+1),       &
     &                               Vt(ll,1),Ldvt)
            IF ( Nru>0 ) CALL CLASR('R','V','B',Nru,m-ll+1,Rwork(1),    &
     &                              Rwork(N),U(1,ll),Ldu)
            IF ( Ncc>0 ) CALL CLASR('L','V','B',m-ll+1,Ncc,Rwork(1),    &
     &                              Rwork(N),C(ll,1),Ldc)
!
!           Test convergence
!
            IF ( ABS(E(ll))<=thresh ) E(ll) = ZERO
         ENDIF
!
!        Use nonzero shift
!
      ELSEIF ( idir==1 ) THEN
!
!           Chase bulge from top to bottom
!           Save cosines and sines for later singular vector updates
!
         f = (ABS(D(ll))-shift)*(SIGN(ONE,D(ll))+shift/D(ll))
         g = E(ll)
         DO i = ll , m - 1
            CALL SLARTG(f,g,cosr,sinr,r)
            IF ( i>ll ) E(i-1) = r
            f = cosr*D(i) + sinr*E(i)
            E(i) = cosr*E(i) - sinr*D(i)
            g = sinr*D(i+1)
            D(i+1) = cosr*D(i+1)
            CALL SLARTG(f,g,cosl,sinl,r)
            D(i) = r
            f = cosl*E(i) + sinl*D(i+1)
            D(i+1) = cosl*D(i+1) - sinl*E(i)
            IF ( i<m-1 ) THEN
               g = sinl*E(i+1)
               E(i+1) = cosl*E(i+1)
            ENDIF
            Rwork(i-ll+1) = cosr
            Rwork(i-ll+1+nm1) = sinr
            Rwork(i-ll+1+nm12) = cosl
            Rwork(i-ll+1+nm13) = sinl
         ENDDO
         E(m-1) = f
!
!           Update singular vectors
!
         IF ( Ncvt>0 ) CALL CLASR('L','V','F',m-ll+1,Ncvt,Rwork(1),     &
     &                            Rwork(N),Vt(ll,1),Ldvt)
         IF ( Nru>0 ) CALL CLASR('R','V','F',Nru,m-ll+1,Rwork(nm12+1),  &
     &                           Rwork(nm13+1),U(1,ll),Ldu)
         IF ( Ncc>0 ) CALL CLASR('L','V','F',m-ll+1,Ncc,Rwork(nm12+1),  &
     &                           Rwork(nm13+1),C(ll,1),Ldc)
!
!           Test convergence
!
         IF ( ABS(E(m-1))<=thresh ) E(m-1) = ZERO
!
      ELSE
!
!           Chase bulge from bottom to top
!           Save cosines and sines for later singular vector updates
!
         f = (ABS(D(m))-shift)*(SIGN(ONE,D(m))+shift/D(m))
         g = E(m-1)
         DO i = m , ll + 1 , -1
            CALL SLARTG(f,g,cosr,sinr,r)
            IF ( i<m ) E(i) = r
            f = cosr*D(i) + sinr*E(i-1)
            E(i-1) = cosr*E(i-1) - sinr*D(i)
            g = sinr*D(i-1)
            D(i-1) = cosr*D(i-1)
            CALL SLARTG(f,g,cosl,sinl,r)
            D(i) = r
            f = cosl*E(i-1) + sinl*D(i-1)
            D(i-1) = cosl*D(i-1) - sinl*E(i-1)
            IF ( i>ll+1 ) THEN
               g = sinl*E(i-2)
               E(i-2) = cosl*E(i-2)
            ENDIF
            Rwork(i-ll) = cosr
            Rwork(i-ll+nm1) = -sinr
            Rwork(i-ll+nm12) = cosl
            Rwork(i-ll+nm13) = -sinl
         ENDDO
         E(ll) = f
!
!           Test convergence
!
         IF ( ABS(E(ll))<=thresh ) E(ll) = ZERO
!
!           Update singular vectors if desired
!
         IF ( Ncvt>0 ) CALL CLASR('L','V','B',m-ll+1,Ncvt,Rwork(nm12+1),&
     &                            Rwork(nm13+1),Vt(ll,1),Ldvt)
         IF ( Nru>0 ) CALL CLASR('R','V','B',Nru,m-ll+1,Rwork(1),       &
     &                           Rwork(N),U(1,ll),Ldu)
         IF ( Ncc>0 ) CALL CLASR('L','V','B',m-ll+1,Ncc,Rwork(1),       &
     &                           Rwork(N),C(ll,1),Ldc)
      ENDIF
!
!     QR iteration finished, go back and check convergence
!
      GOTO 100
!
!     All singular values converged, so make them positive
!
 300  DO i = 1 , N
         IF ( D(i)<ZERO ) THEN
            D(i) = -D(i)
!
!           Change sign of singular vectors, if desired
!
            IF ( Ncvt>0 ) CALL CSSCAL(Ncvt,NEGONE,Vt(i,1),Ldvt)
         ENDIF
      ENDDO
!
!     Sort the singular values into decreasing order (insertion sort on
!     singular values, but only one transposition per singular vector)
!
      DO i = 1 , N - 1
!
!        Scan for smallest D(I)
!
         isub = 1
         smin = D(1)
         DO j = 2 , N + 1 - i
            IF ( D(j)<=smin ) THEN
               isub = j
               smin = D(j)
            ENDIF
         ENDDO
         IF ( isub/=N+1-i ) THEN
!
!           Swap singular values and vectors
!
            D(isub) = D(N+1-i)
            D(N+1-i) = smin
            IF ( Ncvt>0 ) CALL CSWAP(Ncvt,Vt(isub,1),Ldvt,Vt(N+1-i,1),  &
     &                               Ldvt)
            IF ( Nru>0 ) CALL CSWAP(Nru,U(1,isub),1,U(1,N+1-i),1)
            IF ( Ncc>0 ) CALL CSWAP(Ncc,C(isub,1),Ldc,C(N+1-i,1),Ldc)
         ENDIF
      ENDDO
!
!     End of CBDSQR
!
99999 END SUBROUTINE CBDSQR


text:
      program tiling01

C     Check tiling: use integer loop bounds, mutiple of the tile size

      real x(0:15, 0:10)

      do 100 i = 1, 15
         do 200 j = 1, 10
            x(i,j) = float(i+j)
 200     continue
 100  continue

      print *, x

      end


text:
!##############################################################################
!# ****************************************************************************
!# <name> structuresoptcontrol </name>
!# ****************************************************************************
!#
!# <purpose>
!# Underlying structures of the space-time optimal control solver.
!#
!# Routines in this module:
!#
!# 1.) soptc_initParOptControl
!#     -> Reads the parameters for the optimal control from the parameter list
!#        and stores them in the structure.
!#
!# 2.) soptc_doneParOptControl
!#     -> Releases data allocated in soptc_initParOptControl.
!#
!# 3.) soptc_initStepPathFollower
!#     -> Adapts the parameters in the optimal control structure
!#        according to path following information.
!# </purpose>
!##############################################################################

module structuresoptcontrol

  use fsystem
  use storage
  use genoutput
  use boundary
  use triangulation
  use paramlist
  use discretebc
  use discretefbc
  use fparser

  use collection
  use paramlist
  use linearsystemblock

  use analyticsolution

  use timediscretisation

  use structuresdiscretisation

  implicit none

  private

!<types>

!<typeblock>

!  ! Type block encapsuling constraints, only in space.
!  type t_optcconstraintsSpace
!
!    ! Type of constraints to apply to the control u.
!    ! =0: No constraints.
!    ! =1: Constant constraints on u active:
!    !     rconstraintsDistCtrl%dmin1 <= u_1 <= rconstraintsDistCtrl%dmax1, rconstraintsDistCtrl%dmin2 <= u_2 <= rconstraintsDistCtrl%dmax2
!    !     Implementation by DOF.
!    ! =2: Constant constraints on u active:
!    !     rconstraintsDistCtrl%dmin1 <= u_1 <= rconstraintsDistCtrl%dmax1, rconstraintsDistCtrl%dmin2 <= u_2 <= rconstraintsDistCtrl%dmax2
!    !     Implemented by cubature point.
!    ! =3: Constant constraints on u active:
!    !     rconstraintsDistCtrl%dmin1 <= u_1 <= rconstraintsDistCtrl%dmax1, rconstraintsDistCtrl%dmin2 <= u_2 <= rconstraintsDistCtrl%dmax2.
!    !     Implemented by DOF. Newton matrix is approximative.
!    ! =4: Constant constraints on u active:
!    !     rconstraintsDistCtrl%dmin1 <= u_1 <= rconstraintsDistCtrl%dmax1, rconstraintsDistCtrl%dmin2 <= u_2 <= rconstraintsDistCtrl%dmax2.
!    !     Implemented by cubature point (more exact). Adaptive integration.
!    integer :: ccontrolConstraints = 0
!
!    ! Type of definition of the constraints if ccontrolConstraints <> 0.
!    ! =0: constants specified in rconstraintsDistCtrl%dmin1/2, rconstraintsDistCtrl%dmax1/2.
!    ! =1: analytical functions defined in p_rumin1/2, p_rumax1/2.
!    integer :: ccontrolConstraintsType = 0
!
!    ! Constraints on u_1
!    real(DP) :: rconstraintsDistCtrl%dmin1 = -1.0E10
!    real(DP) :: rconstraintsDistCtrl%dmax1 = 1.0E10
!
!    ! Constraints in u_2
!    real(DP) :: rconstraintsDistCtrl%dmin2 = -1.0E10
!    real(DP) :: rconstraintsDistCtrl%dmax2 = 1.0E10
!
!    ! Analytical constraints for u_1
!    type(t_anSolution), pointer :: p_rumin1 => null()
!    type(t_anSolution), pointer :: p_rumax1 => null()
!
!    ! Analytical constraints for u_2
!    type(t_anSolution), pointer :: p_rumin2 => null()
!    type(t_anSolution), pointer :: p_rumax2 => null()
!
!    ! Discrete constraints for u_1, u_2
!    type(t_vectorBlock), pointer :: p_rvectorumin => null()
!    type(t_vectorBlock), pointer :: p_rvectorumax => null()
!
!
!    ! Type of constraints to apply to the state y.
!    ! =0: No constraints.
!    ! =1: Constant constraints on u active:
!    !     rconstraintsState%dmin1 <= y_1 <= rconstraintsState%dmax1, rconstraintsState%dmin2 <= y_2 <= rconstraintsState%dmax2
!    !     Implementation by DOF.
!    ! =2: Constant constraints on y active:
!    !     rconstraintsState%dmin1 <= y_1 <= rconstraintsState%dmax1, rconstraintsState%dmin2 <= y_2 <= rconstraintsState%dmax2
!    !     Implemented by cubature point.
!    integer :: rconstraintsState%cconstrType%cconstraints = 0
!
!    ! Type of definition of the constraints if ccontrolConstraints <> 0.
!    ! =0: constants specified in rconstraintsDistCtrl%dmin1/2, rconstraintsDistCtrl%dmax1/2.
!    ! =1: analytical functions defined in p_rumin1/2, p_rumax1/2.
!    integer :: rconstraintsState%cconstrType%cconstraintsType = 0
!
!    ! Regularisation parameter for the Moreau-Yosida regularisation
!    real(DP) :: rconstraintsState%dregularisation = 1.0_DP
!
!    ! Constraints on u_1
!    real(DP) :: rconstraintsState%dmin1 = -1.0E10
!    real(DP) :: rconstraintsState%dmax1 = 1.0E10
!
!    ! Constraints in u_2
!    real(DP) :: rconstraintsState%dmin2 = -1.0E10
!    real(DP) :: rconstraintsState%dmax2 = 1.0E10
!
!    ! Analytical constraints for u_1
!    type(t_anSolution), pointer :: p_rymin1 => null()
!    type(t_anSolution), pointer :: p_rymax1 => null()
!
!    ! Analytical constraints for u_2
!    type(t_anSolution), pointer :: p_rymin2 => null()
!    type(t_anSolution), pointer :: p_rymax2 => null()
!
!    ! Discrete constraints for u_1, u_2
!    type(t_vectorBlock), pointer :: p_rvectorymin => null()
!    type(t_vectorBlock), pointer :: p_rvectorymax => null()
!
!
!    ! Time used for the evaluation of the analytical control constraints
!    ! p_rumin1/2, p_rumax1/2.
!    real(DP) :: dstateConstrTime = 0.0_DP
!
!    ! Time used for the evaluation of the analytical state constraints
!    ! p_rymin1/2, p_rymax1/2.
!    real(DP) :: dcontrolConstrTime = 0.0_DP
!
!  end type
!
!!</types>
!
!  public :: t_optcconstraintsSpace

!<types>

!<typeblock>

  ! Defines a set of contraints to any function
  type t_optcConstraintsDef

    ! Type of constraints to apply.
    ! =0: No constraints.
    ! =1: Constant constraints active:
    !     dmin1 <= u_1 <= dmax1, dmin2 <= u_2 <= dmax2
    !     Implementation by DOF.
    ! =2: Constant constraints active:
    !     dmin1 <= u_1 <= dmax1, dmin2 <= u_2 <= dmax2
    !     Implemented by cubature point.
    integer :: cconstraints = 0

    ! Type of definition of the constraints if ccontrolConstraints <> 0.
    ! =0: constants specified in dmin1/2, dmax1/2.
    ! =1: analytical functions defined in p_rmin1/2, p_rmax1/2.
    integer :: cconstrType = 0

    ! Constraints on u_1
    real(DP) :: dmin1 = -1.0E10_DP
    real(DP) :: dmax1 = 1.0E10_DP

    ! Section defining the minimum/maximum for the 1st component
    character(LEN=SYS_STRLEN) :: ssectionMin1 = ""
    character(LEN=SYS_STRLEN) :: ssectionMax1 = ""

    ! Analytical constraints for u_1.
    type(t_anSolution), pointer :: p_rmin1 => null()
    type(t_anSolution), pointer :: p_rmax1 => null()

    ! Constraints in u_2
    real(DP) :: dmin2 = -1.0E10_DP
    real(DP) :: dmax2 = 1.0E10_DP

    ! Section defining the minimum/maximum for the 1st component
    character(LEN=SYS_STRLEN) :: ssectionMin2 = ""
    character(LEN=SYS_STRLEN) :: ssectionMax2 = ""

    ! Analytical constraints for u_2
    type(t_anSolution), pointer :: p_rmin2 => null()
    type(t_anSolution), pointer :: p_rmax2 => null()

    ! Regularisation parameter. Not used for all types of constraints
    real(DP) :: dregularisation = 0.0_DP

  end type

!</typeblock>

  public :: t_optcConstraintsDef

!<typeblock>

  ! Type block encapsuling constraints, in space/time.
  type t_optcconstraintsSpaceTime

    ! Constraints in the distributed control
    type(t_optcConstraintsDef) :: rconstraintsDistCtrl

    ! Constraints in the L2 boundary control
    type(t_optcConstraintsDef) :: rconstraintsL2BdC

    ! Constraints in the $H^{1/2}$ boundary control
    type(t_optcConstraintsDef) :: rconstraintsH12BdC

    ! Constraints in the state
    type(t_optcConstraintsDef) :: rconstraintsState

  end type

!</typeblock>

  public :: t_optcconstraintsSpaceTime

!<typeblock>

  ! This structure collects a set of path following data.
  ! During the nonlinear iteration, this data is automatically
  ! adapted from nonlinear to nonlinear step.
  type t_optcpathfollowingdata

    ! List of $\alpha$ parameter of the optimal control functional parameters.
    ! Distributed control.
    ! If this is not =>NULL, a path following algorithm will be applied.
    ! Every nonlinear step, the next value is taken.
    real(DP), dimension(:), pointer :: p_dalphaDistC => null()

    ! List of $\alpha$ parameter of the optimal control functional parameters.
    ! $L_2$ boundary control.
    ! If this is not =>NULL, a path following algorithm will be applied.
    ! Every nonlinear step, the next value is taken.
    real(DP), dimension(:), pointer :: p_dalphaL2BdC => null()

    ! List of $\alpha$ parameter of the optimal control functional parameters.
    ! $H^{1/2}$ boundary control.
    ! If this is not =>NULL, a path following algorithm will be applied.
    ! Every nonlinear step, the next value is taken.
    real(DP), dimension(:), pointer :: p_dalphaH12BdC => null()

    ! List of $\gamma$ parameter of the optimal control functional parameters.
    ! If this is not =>NULL, a path following algorithm will be applied.
    ! Every nonlinear step, the next value is taken.
    real(DP), dimension(:), pointer :: p_dalphaEndTimeC => null()

  end type

!</typeblock>

  public :: t_optcpathfollowingdata

!<typeblock>

  ! This block saves parameters of the optimal control problem
  type t_settings_optcontrol

    ! $\alpha$ parameter of the optimal control functional
    real(DP) :: dalphaDistC = 1.0_DP

    ! $\alpha$ parameter of the optimal control functional via L2 boundary control
    real(DP) :: dalphaL2BdC = 1.0_DP

    ! $\alpha$ parameter of the optimal control functional via $H^{1/2}$ boundary control
    real(DP) :: dalphaH12BdC = 1.0_DP

    ! Penalty parameter for the dirichlet boundary control
    real(DP) :: ddirichletBCPenalty = 100.0_DP

    ! Penalty parameter for the end time condition in the dual equation.
    real(DP) :: dendTimeCondDualPenalty = 0.0_DP

    ! $\gamma$ parameter of the nonstationary optimal control functional
    real(DP) :: dalphaEndTimeC = 0.0_DP

    ! $\delta$ parameter of the terminal condition
    real(DP) :: ddeltaC = 0.0_DP

    ! $\sigma$-parameter that controls the blending in
    ! $F(u) = u - P( u-\sigma(\alpha u + \lambda) )$. Defaults to 1.0.
    ! A value of sigma=1/alpha would give
    ! $F(u) = u - P( -1/alpha \lambda) )$ which seems to be more instable.
    real(DP) :: dsigmaC = 1.0_DP

    ! Formulation of the Space-time problem.
    ! =0: usual formulation as specified in the DFG applicance
    ! =1: Formulation for the generation of reference results from papers
    ! The two formulations differ in a "-"-sign in front of the dual velocity.
    integer :: ispaceTimeFormulation = 0

    ! Whether to treat the convection explicitly or implicitly.
    ! =0: Treat the convection implicitely.
    ! =1: Treat the convection explicitely. (This is the formulation of the paper
    !     of Baerwolff and Hinze!)
    integer :: iconvectionExplicit = 0

    !<!-- --------------- -->
    !<!-- TARGET FUNCTION -->
    !<!-- --------------- -->

    ! Analytic solution defining the target function
    type(t_anSolution) :: rtargetFunction

    !<!-- ----------- -->
    !<!-- CONSTRAINTS -->
    !<!-- ----------- -->

    ! Possible constraints in the problem.
    type(t_optcconstraintsSpaceTime) :: rconstraints

    !<!-- ---------------- -->
    !<!-- OBSERVATION AREA -->
    !<!-- ---------------- -->

    ! Observation area. If this points to NULL(), the whole
    ! domain is observed. Otherwise, this is an array with four entries
    ! in the form
    !   p_DobservationArea = (x1 y1 x2 y2)
    ! specifying the area to be observed.
    real(DP), dimension(:), pointer :: p_DobservationArea => null()

    !<!-- ---------------- -->
    !<!-- OBSERVATION AREA -->
    !<!-- ---------------- -->

    ! Contains data for path following strategies
    type(t_optcpathfollowingdata) :: rpathfollowingdata
  end type

!</typeblock>

  public :: t_settings_optcontrol

!</types>

  public :: soptc_initParOptControl
  public :: soptc_doneParOptControl
  public :: soptc_initStepPathFollower

contains

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

!<subroutine>

  subroutine sys_getParameter (sstring,stoken,istart)

!<description>
  ! Gets a parameter out of a list of parameters.
!</description>

!<input>
  ! A string with parameters, e.g. "0.1 0.2 0.3",...
  character(len=*), intent(in) :: sstring

  ! String that receives the parameter.
  character(len=*), intent(out) :: stoken
!</input>

!<inputoutput>
  ! On input: Must be set to =1 for the first call.
  !   Position from where to search for the next parameter
  ! On output:
  !   =0, if this was the last parameter
  !   Otherwise: position in sstring where the next parameter starts.
  integer, intent(inout) :: istart
!</inputoutput>

!</subroutine>

    integer :: i,slen

    ! If istart=0, do not do anything. There is nothing in the string
    stoken = ""
    if (istart .le. 0) return

    ! If the string length is =0, finish.
    slen = len_trim(sstring)
    if (slen .eq. 0) return

    ! Skip whitespaces.
    do while (istart .le. slen)
      if (sstring(istart:istart) .eq. " ") then
        istart = istart + 1
      else
        exit
      end if
    end do

    ! End of the string?
    i = istart
    do while (i .le. slen)
      if (sstring(i:i) .ne. " ") then
        i = i + 1
      else
        exit
      end if
    end do

    ! Copy
    if (istart .le. slen) then
      stoken = sstring(istart:i)
    end if

    ! Put istart behind the parameter or set it to 0.
    if (i .ge. slen) then
      ! No more parameters
      istart = 0
    else
      ! Put istart behind the parameter -- on the whitespace which
      ! is skipped in the next loop.
      istart = i+1
    end if

  end subroutine

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

!<subroutine>

  subroutine sys_countParameters (sstring,ntokens)

!<description>
  ! Returns the number of parameters in sstring
!</description>

!<input>
  ! A string with parameters, e.g. "0.1 0.2 0.3",...
  character(len=*), intent(in) :: sstring
!</input>

!<output>
  ! Number of parameters in sstring
  integer, intent(out) :: ntokens
!</output>

!</subroutine>

    integer :: slen, istart

    ! If the string length is =0, finish.
    ntokens = 0

    slen = len_trim(sstring)
    if (slen .eq. 0) return

    ! Find all substrings
    istart = 1

    do
      ! Cancel if nothing left.
      if (istart .gt. slen) exit

      ! Skip whitespaces.
      do while (istart .le. slen)
        if (sstring(istart:istart) .eq. " ") then
          istart = istart + 1
        else
          exit
        end if
      end do

      ! Cancel if we reached the string end
      if (istart .gt. slen) exit

      ! End of the string?
      do while (istart .le. slen)
        if (sstring(istart:istart) .ne. " ") then
          istart = istart + 1
        else
          exit
        end if
      end do

      ! One more
      ntokens = ntokens + 1

    end do

  end subroutine

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

!<subroutine>

  subroutine soptc_initParOptControl (rparlist,ssectionOptC,roptcontrol)

!<description>
  ! Reads the parameters for the optimal control from the parameter list
  ! and stores them in the structure.
  !
  ! NOTE: Analytical constraints are not read from disc!
  ! The caller has to do that!
!</description>

!<input>
  ! Parameter list
  type(t_parlist), intent(in) :: rparlist

  ! Section where the parameters of the optimal control can be found.
  character(len=*), intent(in) :: ssectionOptC
!</input>

!<inputoutput>
  ! Structure defining the parameters for the optimal control.
  type(t_settings_optcontrol), intent(inout) :: roptcontrol
!</inputoutput>

!</subroutine>

    character(len=SYS_STRLEN) :: sstring,stoken
    integer :: ntokens,istart,itoken

    ! Sigma-parameter that controls the formulation.
    call parlst_getvalue_double (rparlist,ssectionOptC,&
        "dsigmaC",roptcontrol%dsigmaC,roptcontrol%dsigmaC)

    ! Alpha/Gamma parameters.
    !
    ! These parameters may have two types: either it is a single
    ! parameter or a list of parameters. In the second case,
    ! this defines a path following strategy, i.e., the parameters
    ! change according to the list in every nonlinear iteration.
    call parlst_getvalue_string (rparlist,ssectionOptC,&
        "dalphaDistC",sstring,"-1.0",bdequote=.true.)

    ! Is this a list?
    call sys_countParameters (sstring,ntokens)

    if (ntokens .eq. 1) then

      read(sstring,*) roptcontrol%dalphaDistC

    else if (ntokens .gt. 1) then
      ! Read the list
      roptcontrol%dalphaDistC = 0.0_DP
      allocate(roptcontrol%rpathfollowingdata%p_dalphaDistC(ntokens))

      ! Loop through all tokens and save them
      istart = 1
      do itoken = 1,ntokens
        call sys_getParameter (sstring,stoken,istart)
        read(stoken,*) roptcontrol%rpathfollowingdata%p_dalphaDistC(itoken)
      end do
    end if

    call parlst_getvalue_string (rparlist,ssectionOptC,&
        "dalphaL2BdC",sstring,"-1.0",bdequote=.true.)

    ! Is this a list?
    call sys_countParameters (sstring,ntokens)

    if (ntokens .eq. 1) then

      read(sstring,*) roptcontrol%dalphaL2BdC

    else if (ntokens .gt. 1) then
      ! Read the list
      roptcontrol%dalphaL2BdC = 0.0_DP
      allocate(roptcontrol%rpathfollowingdata%p_dalphaL2BdC(ntokens))

      ! Loop through all tokens and save them
      istart = 1
      do itoken = 1,ntokens
        call sys_getParameter (sstring,stoken,istart)
        read(stoken,*) roptcontrol%rpathfollowingdata%p_dalphaL2BdC(itoken)
      end do
    end if

    call parlst_getvalue_string (rparlist,ssectionOptC,&
        "dalphaH12BdC",sstring,"-1.0",bdequote=.true.)

    ! Is this a list?
    call sys_countParameters (sstring,ntokens)

    if (ntokens .eq. 1) then

      read(sstring,*) roptcontrol%dalphaH12BdC

    else if (ntokens .gt. 1) then
      ! Read the list
      roptcontrol%dalphaH12BdC = 0.0_DP
      allocate(roptcontrol%rpathfollowingdata%p_dalphaH12BdC(ntokens))

      ! Loop through all tokens and save them
      istart = 1
      do itoken = 1,ntokens
        call sys_getParameter (sstring,stoken,istart)
        read(stoken,*) roptcontrol%rpathfollowingdata%p_dalphaH12BdC(itoken)
      end do
    end if

    call parlst_getvalue_string (rparlist,ssectionOptC,&
        "dalphaEndTimeC",sstring,"0.0",bdequote=.true.)

    ! Is this a list?
    call sys_countParameters (sstring,ntokens)

    if (ntokens .eq. 1) then

      read(sstring,*) roptcontrol%dalphaEndTimeC

    else if (ntokens .gt. 1) then
      ! Read the list
      roptcontrol%dalphaEndTimeC = 0.0_DP
      allocate(roptcontrol%rpathfollowingdata%p_dalphaEndTimeC(ntokens))

      ! Loop through all tokens and save them
      istart = 1
      do itoken = 1,ntokens
        call sys_getParameter (sstring,stoken,istart)
        read(stoken,*) roptcontrol%rpathfollowingdata%p_dalphaEndTimeC(itoken)
      end do
    end if

    ! Penalty parameter for the dirichlet boundary control
    call parlst_getvalue_double (rparlist,ssectionOptC,&
        "ddirichletBCPenalty",roptcontrol%ddirichletBCPenalty,100.0_DP)

    # Penalty parameter for the end time condition in the dual equation.
    call parlst_getvalue_double (rparlist,ssectionOptC,&
        "dendTimeCondDualPenalty",roptcontrol%dendTimeCondDualPenalty,0.0_DP)

    ! Type of the formulation
    call parlst_getvalue_int (rparlist,ssectionOptC,&
        "ispaceTimeFormulation",roptcontrol%ispaceTimeFormulation,0)

    call parlst_getvalue_int (rparlist,ssectionOptC,&
        "iconvectionExplicit",roptcontrol%iconvectionExplicit,0)

    ! Observation area
    call parlst_getvalue_string (rparlist,ssectionOptC,&
        "DobservationArea",sstring,"")

    if (sstring .ne. "") then
      ! Read the observation area. This is a box, the parameters
      ! have the format "x1, y1, x2, y2".
      allocate (roptcontrol%p_DobservationArea(4))
      read (sstring,*) &
          roptcontrol%p_DobservationArea(1), &
          roptcontrol%p_DobservationArea(2), &
          roptcontrol%p_DobservationArea(3), &
          roptcontrol%p_DobservationArea(4)
    end if

    ! -------------------------------------------
    ! Read the constraints
    ! -------------------------------------------

    ! Distributed control
    call parlst_getvalue_string (rparlist,ssectionOptC,&
        "ssectionConstrDistCtrl",sstring,"",bdequote=.true.)

    call soptc_getConstraints (rparlist,sstring,&
        roptcontrol%rconstraints%rconstraintsDistCtrl)

    ! L2 boundary control
    call parlst_getvalue_string (rparlist,ssectionOptC,&
        "ssectionConstrL2BdC",sstring,"",bdequote=.true.)

    call soptc_getConstraints (rparlist,sstring,&
        roptcontrol%rconstraints%rconstraintsL2BdC)

    ! H^^1/2 boundary control
    call parlst_getvalue_string (rparlist,ssectionOptC,&
        "ssectionConstrH12BdC",sstring,"",bdequote=.true.)

    call soptc_getConstraints (rparlist,sstring,&
        roptcontrol%rconstraints%rconstraintsH12BdC)

    ! State
    call parlst_getvalue_string (rparlist,ssectionOptC,&
        "ssectionConstrState",sstring,"",bdequote=.true.)

    call soptc_getConstraints (rparlist,sstring,&
        roptcontrol%rconstraints%rconstraintsState)

  end subroutine

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

!<subroutine>

  subroutine soptc_getConstraints (rparlist,ssection,rconstraints)

!<description>
  ! Reads the parameters that defines a set of constraints.
  !
  ! NOTE: Analytical constraints are not read from disc!
!</description>

!<input>
  ! Parameter list
  type(t_parlist), intent(in) :: rparlist

  ! Section where the parameters of the constraints can be found.
  character(len=*), intent(in) :: ssection
!</input>

!<inputoutput>
  ! Constraints-structure, to be set up.
  type(t_optcConstraintsDef), intent(out) :: rconstraints
!</inputoutput>

!</subroutine>

    call parlst_getvalue_int (rparlist,ssection,&
        "cconstraints",rconstraints%cconstraints,0)

    call parlst_getvalue_int (rparlist,ssection,&
        "cconstrType",rconstraints%cconstrType,0)

    call parlst_getvalue_double (rparlist,ssection,&
        "dmin1",rconstraints%dmin1,-1.0E10_DP)

    call parlst_getvalue_double (rparlist,ssection,&
        "dmax1",rconstraints%dmax1,1.0E10_DP)

    call parlst_getvalue_double (rparlist,ssection,&
        "dmin2",rconstraints%dmin2,-1.0E10_DP)

    call parlst_getvalue_double (rparlist,ssection,&
        "dmax2",rconstraints%dmax2,1.0E10_DP)

    call parlst_getvalue_double (rparlist,ssection,&
        "dregularisation",rconstraints%dregularisation,1.0_DP)

    call parlst_getvalue_string (rparlist,ssection,&
        "ssectionMin1",rconstraints%ssectionMin1,"",bdequote=.true.)

    call parlst_getvalue_string (rparlist,ssection,&
        "ssectionMax1",rconstraints%ssectionMax1,"",bdequote=.true.)

    call parlst_getvalue_string (rparlist,ssection,&
        "ssectionMin2",rconstraints%ssectionMin2,"",bdequote=.true.)

    call parlst_getvalue_string (rparlist,ssection,&
        "ssectionMax2",rconstraints%ssectionMax2,"",bdequote=.true.)

  end subroutine

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

!<subroutine>

  subroutine soptc_doneParOptControl (roptcontrol)

!<description>
  ! Cleans up information in the optimal control structure.
!</description>

!<inputoutput>
  ! Structure defining the parameters for the optimal control.
  type(t_settings_optcontrol), intent(inout) :: roptcontrol
!</inputoutput>

!</subroutine>

    ! Release the observation area
    if (associated(roptcontrol%p_DobservationArea)) then
      deallocate(roptcontrol%p_DobservationArea)
    end if

    ! Release path following data
    if (associated(roptcontrol%rpathfollowingdata%p_dalphaDistC)) then
      deallocate(roptcontrol%rpathfollowingdata%p_dalphaDistC)
    end if

    if (associated(roptcontrol%rpathfollowingdata%p_dalphaH12BdC)) then
      deallocate(roptcontrol%rpathfollowingdata%p_dalphaH12BdC)
    end if

    if (associated(roptcontrol%rpathfollowingdata%p_dalphaL2BdC)) then
      deallocate(roptcontrol%rpathfollowingdata%p_dalphaL2BdC)
    end if

    if (associated(roptcontrol%rpathfollowingdata%p_dalphaEndTimeC)) then
      deallocate(roptcontrol%rpathfollowingdata%p_dalphaEndTimeC)
    end if

  end subroutine

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

!<subroutine>

  subroutine soptc_initStepPathFollower (roptcontrol,istep,bprint)

!<description>
  ! Adapts the parameters in the optimal control structure
  ! according to path following information.
  ! The step istep defines the current setting of the parameters
  ! which are adapted via path following methods.
!</description>

!<input>
  ! Current step in the path following strategy
  integer, intent(in) :: istep

  ! Whether or not to print information about the path following parameters.
  logical, intent(in) :: bprint
!</input>

!<inputoutput>
  ! Structure defining the parameters for the optimal control.
  ! The parameters in this structure which are influenced
  ! via the path following data are modified according to
  ! step istep.
  type(t_settings_optcontrol), intent(inout) :: roptcontrol
!</inputoutput>

!</subroutine>

    ! Change the data in roptcontrol.
    ! istep defines the array position of the element to take.
    if (associated(roptcontrol%rpathfollowingdata%p_dalphaDistC)) then
      roptControl%dalphaDistC = roptControl%rpathfollowingdata%p_dalphaDistC(&
              min(size(roptcontrol%rpathfollowingdata%p_dalphaDistC),istep))
      if (bprint) then
        call output_line ("Path-Follower: Parameter ALPHA_DIST = "//&
            trim(sys_sdEL(roptcontrol%dalphaDistC,10)))
      end if
    end if

    if (associated(roptcontrol%rpathfollowingdata%p_dalphaH12BdC)) then
      roptcontrol%dalphaL2BdC = roptControl%rpathfollowingdata%p_dalphaH12BdC(&
              min(size(roptcontrol%rpathfollowingdata%p_dalphaH12BdC),istep))
      if (bprint) then
        call output_line ("Path-Follower: Parameter ALPHA_L2BDC  = "//&
            trim(sys_sdEL(roptcontrol%dalphaL2BdC,10)))
      end if
    end if

    if (associated(roptcontrol%rpathfollowingdata%p_dalphaH12BdC)) then
      roptcontrol%dalphaH12BdC = roptControl%rpathfollowingdata%p_dalphaH12BdC(&
              min(size(roptcontrol%rpathfollowingdata%p_dalphaH12BdC),istep))
      if (bprint) then
        call output_line ("Path-Follower: Parameter ALPHA_H12BDC  = "//&
            trim(sys_sdEL(roptcontrol%dalphaH12BdC,10)))
      end if
    end if

    if (associated(roptcontrol%rpathfollowingdata%p_dalphaEndTimeC)) then
      roptcontrol%dalphaEndTimeC = roptControl%rpathfollowingdata%p_dalphaEndTimeC(&
              min(size(roptcontrol%rpathfollowingdata%p_dalphaEndTimeC),istep))
      if (bprint) then
        call output_line ("Path-Follower: Parameter GAMMA = "//&
            trim(sys_sdEL(roptcontrol%dalphaEndTimeC,10)))
      end if
    end if

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
c> @file
c> Fortran subroutines for metric terms

c ---------------------------------------------------------
c> @brief @copybrief ::fclaw2d_fort_compute_mesh_t
c>
c> Default implimentation
c>
c> @details @copydetails ::fclaw2d_fort_compute_mesh_t
c ---------------------------------------------------------
      subroutine fclaw2d_fort_compute_mesh(mx,my,mbc,
     &      xlower,ylower, dx,dy,blockno,xp,yp,zp,xd,yd,zd)
      implicit none

      integer mx,my, mbc, blockno
      double precision dx,dy,xlower,ylower

      double precision xp(-mbc:mx+mbc+1,-mbc:my+mbc+1)
      double precision yp(-mbc:mx+mbc+1,-mbc:my+mbc+1)
      double precision zp(-mbc:mx+mbc+1,-mbc:my+mbc+1)

      double precision xd(-mbc:mx+mbc+2,-mbc:my+mbc+2)
      double precision yd(-mbc:mx+mbc+2,-mbc:my+mbc+2)
      double precision zd(-mbc:mx+mbc+2,-mbc:my+mbc+2)

      integer i,j
      double precision dxf,dyf, xc,yc,xd1, yd1,zd1

      integer*8 map_context_ptr, get_context

      map_context_ptr = get_context()

c     # We need both cell centered and node locations to
c     # compute the normals at cell edges.
      dxf = dx/2.d0
      dyf = dy/2.d0
      do i = -2*mbc-1,2*(mx+mbc+1)+1
         do j = -2*mbc-1,2*(my+mbc+1)+1
            if (abs(mod(i,2)) .ne. abs(mod(j,2))) then
               cycle
            endif
            xc = xlower + (i-1)*dxf
            yc = ylower + (j-1)*dyf

c            call mapc2m(xc,yc ,xd1,yd1,zd1)

            call fclaw2d_map_c2m(map_context_ptr,
     &            blockno,xc,yc,xd1,yd1,zd1)

            if (abs(mod(i,2)) .eq. 1) then
c              # Physical location of cell vertices
               xd((i-1)/2 + 1, (j-1)/2 + 1) = xd1
               yd((i-1)/2 + 1, (j-1)/2 + 1) = yd1
               zd((i-1)/2 + 1, (j-1)/2 + 1) = zd1
            else
c              # Physical locations of cell centers
               xp(i/2,j/2) = xd1
               yp(i/2,j/2) = yd1
               zp(i/2,j/2) = zd1
            endif
         enddo
      enddo
      end

c ----------------------------------------------------------
c> @brief @copybrief ::fclaw2d_fort_compute_area_t
c>
c> Default implimentation
c>
c> @details @copydetails ::fclaw2d_fort_compute_area_t
c ---------------------------------------------------------
      subroutine fclaw2d_fort_compute_area(mx,my,mbc,dx,dy,
     &      xlower, ylower, blockno,area, quadsize, quadstore,
     &      ghost_only)
      implicit none

      integer mx,my,mbc,blockno, quadsize
      double precision quadstore(0:quadsize,0:quadsize,3)
      double precision dx,dy, xlower, ylower
      integer ghost_only
      double precision area(-mbc:mx+mbc+1,-mbc:my+mbc+1)

      integer*8 cont, get_context
      logical isaffine

      cont = get_context()
      if (isaffine()) then
c        # We don't need to compute areas all the way to the
c        # finest level.
         call fclaw2d_fort_compute_area_affine(mx,my,mbc,dx,dy,
     &         xlower, ylower, blockno,area, ghost_only)
      else
         call fclaw2d_fort_compute_area_general(mx,my,mbc,dx,dy,
     &         xlower, ylower, blockno,area,quadsize,quadstore,
     &         ghost_only)
      endif

      end

c ---------------------------------------------------------
c> @brief Compute the area for each cell for general mappings
c>
c> @param[in] mx, my the number of cells in the x and y directions
c> @param[in] mbc the number of ghost cells
c> @param[in] dx, dy the spacings in the x and y direcitons
c> @param[in] xlower, ylower the lower left coordinate of the patch
c> @param[in] blockno the block number
c> @param[out] area the area of each cell
c> @param[in] quadsize the length of the quad
c> @param[in] quadstore stores a group of cell values
c> @param[in] ghost_only
c ---------------------------------------------------------
      subroutine fclaw2d_fort_compute_area_general(mx,my,mbc,
     &      dx,dy, xlower, ylower, blockno,area,
     &      quadsize, quadstore,ghost_only)
      implicit none

      integer mx,my,mbc,blockno, quadsize
      double precision dx,dy, xlower, ylower
      double precision quadstore(0:quadsize,0:quadsize,3)
      integer ghost_only

      double precision area(-mbc:mx+mbc+1,-mbc:my+mbc+1)

      integer i,j,ii,jj, rfactor, icell, jcell
      double precision sum_area
      double precision xp1,yp1,zp1
      double precision quad(0:1,0:1,3)
      double precision get_area_approx

      double precision dxf, dyf
      double precision xef, yef, xe,ye
      logical is_area_interior

      integer*8 cont, get_context

      cont = get_context()

      rfactor = quadsize
      dxf = dx/rfactor
      dyf = dy/rfactor

c     # Primary cells.  Note that we don't do anything special
c     # in the diagonal cells - so the areas there are less accurate
c     # than in the rest of the mesh.
      do j = -mbc,my+mbc+1
         do i = -mbc,mx+mbc+1
            if (is_area_interior(mx,my,i,j) .and.
     &            ghost_only .eq. 1) then
               cycle
            endif
            xe = xlower + (i-1)*dx
            ye = ylower + (j-1)*dy

            do ii = 0,rfactor
               do jj = 0,rfactor
                  xef = xe + ii*dxf
                  yef = ye + jj*dyf

                  call fclaw2d_map_c2m(cont,
     &                  blockno,xef,yef,xp1,yp1,zp1)

                  quadstore(ii,jj,1) = xp1
                  quadstore(ii,jj,2) = yp1
                  quadstore(ii,jj,3) = zp1
               enddo
            enddo

            sum_area = 0.d0
            do ii = 0,rfactor-1
               do jj = 0,rfactor-1
                  do icell = 0,1
                     do jcell = 0,1
                        quad(icell,jcell,1) =
     &                        quadstore(ii+icell,jj+jcell,1)
                        quad(icell,jcell,2) =
     &                        quadstore(ii+icell,jj+jcell,2)
                        quad(icell,jcell,3) =
     &                        quadstore(ii+icell,jj+jcell,3)
                     enddo
                  enddo
                  sum_area = sum_area + get_area_approx(quad)
               enddo
            enddo
            area(i,j) = sum_area
         enddo
      enddo

      end

c ---------------------------------------------------------
c> @brief Compute the area for each cell for affine mappings.
c>
c> If the mapping is affine, (e.g. Ax + b) then we don't need to sum
c> the finer level areas.
c>
c> @param[in] mx, my the number of cells in the x and y directions
c> @param[in] mbc the number of ghost cells
c> @param[in] dx, dy the spacings in the x and y direcitons
c> @param[in] xlower, ylower the lower left coordinate of the patch
c> @param[in] blockno the block number
c> @param[out] area the area of each cell
c> @param[in] quadsize the length of the quad
c> @param[in] quadstore stores a group of cell values
c> @param[in] ghost_only
c ---------------------------------------------------------

      subroutine fclaw2d_fort_compute_area_affine(mx,my,mbc,dx,dy,
     &      xlower, ylower, blockno,area,ghost_only)
      implicit none

      integer mx,my,mbc,blockno
      double precision dx,dy, xlower, ylower
      integer ghost_only

      double precision area(-mbc:mx+mbc+1,-mbc:my+mbc+1)

      integer i,j, icell, jcell
      double precision xcorner, ycorner
      double precision xe,ye, xp1,yp1,zp1
      double precision quad(0:1,0:1,3)
      double precision get_area_approx
      logical is_area_interior

      integer*8 map_context_ptr, get_context

      map_context_ptr = get_context()

      do j = -mbc,my+mbc+1
         do i = -mbc,mx+mbc+1
            if (is_area_interior(mx,my,i,j) .and.
     &            ghost_only .eq. 1) then
               cycle
            endif
            xe = xlower + (i-1)*dx
            ye = ylower + (j-1)*dy
            do icell = 0,1
               do jcell = 0,1
                  xcorner = xe + icell*dx
                  ycorner = ye + jcell*dy
                  call fclaw2d_map_c2m(map_context_ptr,
     &                  blockno,xcorner,ycorner,xp1,yp1,zp1)
                  quad(icell,jcell,1) = xp1
                  quad(icell,jcell,2) = yp1
                  quad(icell,jcell,3) = zp1
               enddo
            enddo
            area(i,j) = get_area_approx(quad)
         enddo
      enddo

      end

c ---------------------------------------------------------
c> @brief Check if the index is an interior index
c>
c> @param[in] mx, my the number of cell sin the x and y directions
c> @param[in] i, j the index
c> @return true if the index is interior
c ---------------------------------------------------------
      logical function is_area_interior(mx,my,i,j)
      implicit none
      integer mx,my,i,j

      is_area_interior = i .ge. 0 .and. i .le. mx+1 .and.
     &      j .ge. 0 .and. j .le. my+1

      end

c ---------------------------------------------------------
c> @brief Approximate the area based on the corners of the mesh cell
c>
c> This is the area element based on a bilinear approximation to
c> the surface defined by the corners of the mesh cell.
c>
c> @param[in] quad the coordinates of the cell nodes
c> @return the approximate area
c ---------------------------------------------------------
      double precision function get_area_approx(quad)
      implicit none

      double precision quad(0:1,0:1,3)

      double precision a00(3), a01(3),a10(3), a11(3)
      double precision Txi(3), Teta(3)
      double precision norm_cross, xi,eta

      integer m

c     # Coefficients for bilinear approximation to surface
      do m = 1,3
         a00(m) = quad(0,0,m)
         a01(m) = quad(1,0,m) - quad(0,0,m)
         a10(m) = quad(0,1,m) - quad(0,0,m)
         a11(m) = quad(1,1,m) - quad(1,0,m) -
     &         quad(0,1,m) + quad(0,0,m)
      enddo

      eta = 0.5d0
      xi  = 0.5d0

c     # Mesh square is approximated by
c     #       T(xi,eta) = a00 + a01*xi + a10*eta + a11*xi*eta
c     #
c     # Area is the length of the cross product of T_xi and T_eta.
c     # computed at the center of the mesh cell.
      do m = 1,3
         txi(m)  = (a01(m) + a11(m)*eta)
         teta(m) = (a10(m) + a11(m)*xi)
      enddo

      get_area_approx = norm_cross(txi,teta)

      end

c ---------------------------------------------------------
c> @brief Get the l2norm of the vector cross product
c>
c> @param[in] u, v the vectors
c> @return the norm of the cross product
c ---------------------------------------------------------
      double precision function norm_cross(u,v)
      implicit none

      double precision u(3), v(3), w(3)

      w(1) =  u(2)*v(3) - u(3)*v(2)
      w(2) = -u(1)*v(3) + u(3)*v(1)
      w(3) =  u(1)*v(2) - u(2)*v(1)

      norm_cross = sqrt(w(1)*w(1) + w(2)*w(2) + w(3)*w(3))

      end

c ---------------------------------------------------------
c> @brief @copybrief ::fclaw2d_fort_compute_normals_t
c>
c> Default implimentation
c>
c> @details @copydetails ::fclaw2d_fort_compute_normals_t
c ---------------------------------------------------------
      subroutine fclaw2d_fort_compute_normals(mx,my,mbc,
     &      xp,yp,zp,xd,yd,zd,xnormals,ynormals)
      IMPLICIT NONE

      INTEGER mx,my,mbc

      double precision xp(-mbc:mx+mbc+1,-mbc:my+mbc+1)
      double precision yp(-mbc:mx+mbc+1,-mbc:my+mbc+1)
      double precision zp(-mbc:mx+mbc+1,-mbc:my+mbc+1)

      double precision xd(-mbc:mx+mbc+2,-mbc:my+mbc+2)
      double precision yd(-mbc:mx+mbc+2,-mbc:my+mbc+2)
      double precision zd(-mbc:mx+mbc+2,-mbc:my+mbc+2)

      double precision xnormals(-mbc:mx+mbc+2,-mbc:my+mbc+2,3)
      double precision ynormals(-mbc:mx+mbc+2,-mbc:my+mbc+2,3)

      INTEGER i,j,m
      DOUBLE PRECISION taud(3),taup(3),nv(3), sp

c     # Compute normals at all interior edges.

c     # Get x-face normals
      DO j =  -mbc,my+mbc+1
         DO i = 1-mbc,mx+mbc+1

            taud(1) = xp(i,j) - xp(i-1,j)
            taud(2) = yp(i,j) - yp(i-1,j)
            taud(3) = zp(i,j) - zp(i-1,j)

            taup(1) = xd(i,j+1) - xd(i,j)
            taup(2) = yd(i,j+1) - yd(i,j)
            taup(3) = zd(i,j+1) - zd(i,j)

            CALL get_normal(taup,taud,nv,sp)

            DO m = 1,3
               xnormals(i,j,m) = nv(m)
            ENDDO
         ENDDO
      ENDDO

      DO j = 1-mbc,my+mbc+1
         DO i =  -mbc,mx+mbc+1
c           # Now do y-faces
            taud(1) = xp(i,j) - xp(i,j-1)
            taud(2) = yp(i,j) - yp(i,j-1)
            taud(3) = zp(i,j) - zp(i,j-1)

            taup(1) = xd(i+1,j) - xd(i,j)
            taup(2) = yd(i+1,j) - yd(i,j)
            taup(3) = zd(i+1,j) - zd(i,j)

            call get_normal(taup,taud,nv,sp)

c           # nv has unit length
            do m = 1,3
               ynormals(i,j,m) = nv(m)
            enddo
         enddo
      enddo

      end subroutine

c ---------------------------------------------------------
c> @brief @copybrief ::fclaw2d_fort_compute_tangents_t
c>
c> Default implimentation
c>
c> @details @copydetails ::fclaw2d_fort_compute_tangents_t
c ---------------------------------------------------------
      subroutine fclaw2d_fort_compute_tangents(mx,my,mbc,
     &      xd,yd,zd, xtangents,ytangents,edge_lengths)
      IMPLICIT NONE

      INTEGER mx,my,mbc

      double precision xd(-mbc:mx+mbc+2,-mbc:my+mbc+2)
      double precision yd(-mbc:mx+mbc+2,-mbc:my+mbc+2)
      double precision zd(-mbc:mx+mbc+2,-mbc:my+mbc+2)

      double precision xtangents(-mbc:mx+mbc+2,-mbc:my+mbc+2,3)
      double precision ytangents(-mbc:mx+mbc+2,-mbc:my+mbc+2,3)
      double precision edge_lengths(-mbc:mx+mbc+2,-mbc:my+mbc+2,2)

      INTEGER i,j,m
      DOUBLE PRECISION taup(3),tlen

c     # Compute normals at all interior edges.

c     # Get x-face normals
      DO j = -mbc,my+mbc+1
         DO i = -mbc,mx+mbc+2

            taup(1) = xd(i,j+1) - xd(i,j)
            taup(2) = yd(i,j+1) - yd(i,j)
            taup(3) = zd(i,j+1) - zd(i,j)
            tlen = sqrt(taup(1)**2 + taup(2)**2 + taup(3)**2)

            DO m = 1,3
c               xtangents(i,j,m) = taup(m)/tlen
               xtangents(i,j,m) = taup(m)
            ENDDO
            edge_lengths(i,j,1) = tlen
         ENDDO
      ENDDO

      DO j = -mbc,my+mbc+2
         DO i = -mbc,mx+mbc+1
c           # Now do y-faces
            taup(1) = xd(i+1,j) - xd(i,j)
            taup(2) = yd(i+1,j) - yd(i,j)
            taup(3) = zd(i+1,j) - zd(i,j)
            tlen = sqrt(taup(1)**2 + taup(2)**2 + taup(3)**2)

            do m = 1,3
c               ytangents(i,j,m) = taup(m)/tlen
               ytangents(i,j,m) = taup(m)
            enddo
            edge_lengths(i,j,2) = tlen
         enddo
      enddo

      end subroutine

c ---------------------------------------------------------
c> Compute an approximate unit normal to cell edge
c>
c> @param[in] taud     Vector between adjoining cell centers ("dual" cell edge)
c> @param[in] taup     Vector between edge nodes ("primal" cell edge)
c> @param[out] nv       Unit vector vector normal to taup, and tangent to the surface
c> @param[out] sp       Length of taup
c>
c> Idea is to construct normal to edge vector using
c>
c>     t^i = a^{i1} t_1 + a^{i2} t_2
c>
c> where t_1, t_2 are coordinate basis vectors T_xi, T_eta
c> and t^i dot t_j = 1 (i == j), 0 (i /= j)
c>
c> At left edge, we have
c>       taud ~ T_xi = t_1  and   taup ~ T_eta = t_2
c>
c> At bottom edge, we have
c>       taup ~ T_xi = t_1  and   taud ~ T_eta = t_2
c>
c> Metric a_{ij} = t_i dot t_j
c> Metric a^{ij} is the inverse of a_{ij}
c>
c> Both versions use this idea.  That is, both compute
c> coefficients c1, c2 so that the unit normal at the edge
c> is expressed as
c>
c>        n = c1*taup + c2*taud
c>
c> In version 1, trig. identities are used to express these coefficients
c> in terms of the angle between taup and taud.  In version 2, entries of
c> the conjugate tensor are used directly.   Both compute the same vector, but
c> version 2 is about 20%-30% faster.
c>
c> Formulas work for both left and bottom edges.
c>
c> Neither version depends on any knowledge of the surface normals.
c ---------------------------------------------------------
      subroutine get_normal(taup,taud,nv,sp)
      implicit none

      double precision taup(3),taud(3),nv(3), sp
      integer version

      data version /2/

      if (version .eq .1) then
         call get_normal_ver1(taup,taud,nv)
      elseif (version .eq. 2) then
         call get_normal_ver2(taup,taud,nv,sp)
      endif

      end

c ---------------------------------------------------------
c> @brief see get_normal()
c>
c> @param[in] taud     Vector between adjoining cell centers ("dual" cell edge)
c> @param[in] taup     Vector between edge nodes ("primal" cell edge)
c> @param[out] nv       Unit vector vector normal to taup, and tangent to the surface
c ---------------------------------------------------------
      subroutine get_normal_ver1(taup,taud,nv)
      implicit none

      double precision taup(3),taud(3),nv(3)
      double precision sp,sd,ct,st,c1,c2
      double precision tp(3), td(3), nlen

      integer m

      sp = 0
      sd = 0
      do m = 1,3
         sp = sp + taup(m)*taup(m)
         sd = sd + taud(m)*taud(m)
      enddo
      sp = sqrt(sp)
      sd = sqrt(sd)

c     # st = sin(theta)
c     # ct = cos(theta)
c     # theta = angle between taup and taud
      ct = 0
      do m = 1,3
         tp(m) = taup(m)/sp
         td(m) = taud(m)/sd
         ct = ct + tp(m)*td(m)
      enddo
      st = sqrt(1.d0-ct*ct)

      c1 = 1.d0/st
      c2 = -ct/st
      do m = 1,3
         nv(m) = c1*td(m) + c2*tp(m)
         nlen = nlen + nv(m)*nv(m)
      enddo
c      if (nlen .eq. 0) then
c         write(6,*) 'nlen == 0'
c      endif

      end

c ---------------------------------------------------------
c> @brief see get_normal()
c>
c> @param[in] taud     Vector between adjoining cell centers ("dual" cell edge)
c> @param[in] taup     Vector between edge nodes ("primal" cell edge)
c> @param[out] nv       Unit vector vector normal to taup, and tangent to the surface
c> @param[out] sp       Length of taup
c ---------------------------------------------------------
      subroutine get_normal_ver2(taup,taud,nv,sp)
      implicit none

      double precision taup(3),taud(3),nv(3), sp

      integer m
      double precision nlen,aii,ai2

c     # we leave out any scaling by the determinant of the metric,
c     # since we just normalize the vector anyway.
      aii = 0
      ai2 = 0
      do m = 1,3
         aii = aii + taup(m)*taup(m)
         ai2 = ai2 - taud(m)*taup(m)
      enddo

      nlen = 0
      do m = 1,3
         nv(m) = aii*taud(m) + ai2*taup(m)
         nlen = nlen + nv(m)*nv(m)
      enddo
      nlen = sqrt(nlen)

c      c1 = aii/nlen
c      c2 = ai2/nlen

      do m = 1,3
c         nv(m) = c1*taud(m) + c2*taup(m)
         nv(m) = nv(m)/nlen
      enddo
      sp = sqrt(aii)

      end

c ---------------------------------------------------------
c> @brief @copybrief ::fclaw2d_fort_compute_surf_normals_t
c>
c> Default implimentation
c>
c> @details @copydetails ::fclaw2d_fort_compute_surf_normals_t
c ---------------------------------------------------------
      subroutine fclaw2d_fort_compute_surf_normals(mx,my,mbc,
     &      xnormals,ynormals,edge_lengths,curvature,
     &      surfnormals,area)
      IMPLICIT NONE

      INTEGER mx,my,mbc

      double precision     xnormals(-mbc:mx+mbc+2,-mbc:my+mbc+2,3)
      double precision     ynormals(-mbc:mx+mbc+2,-mbc:my+mbc+2,3)
      double precision edge_lengths(-mbc:mx+mbc+2,-mbc:my+mbc+2,2)
      double precision  surfnormals(-mbc:mx+mbc+1,-mbc:my+mbc+1,3)
      double precision    curvature(-mbc:mx+mbc+1,-mbc:my+mbc+1)
      double precision         area(-mbc:mx+mbc+1,-mbc:my+mbc+1)

      integer i,j,m
      double precision nv(3), nlen, av(3), bv(3), kappa

      logical isflat

      do i = 1-mbc,mx+mbc
         do j = 1-mbc,my+mbc
            if (isflat()) then
c              # The surface normal can be computed from the cross product of
c              # an xnormal and a ynormal.  The curvature is zero in this
c              # case
               do m = 1,3
c                 # It shouldn't matter which ones we pick.
                  av(m) = xnormals(i,j,m)
                  bv(m) = ynormals(i,j,m)
               enddo
c              # construct cross product
               nv(1) =   av(2)*bv(3) - av(3)*bv(2)
               nv(2) = -(av(1)*bv(3) - av(3)*bv(1))
               nv(3) =   av(1)*bv(2) - av(2)*bv(1)
               nlen = sqrt(nv(1)*nv(1) + nv(2)*nv(2) + nv(3)*nv(3))
               kappa = 0
            else
               nlen = 0
               do m = 1,3
                  nv(m) = 0
               enddo
               do m = 1,3
                  nv(m) =
     &                  edge_lengths(i+1,j,1)*xnormals(i+1,j,m) -
     &                  edge_lengths(i,j,1)*xnormals(i,j,m) +
     &                  edge_lengths(i,j+1,2)*ynormals(i,j+1,m) -
     &                  edge_lengths(i,j,2)*ynormals(i,j,m)
                  nlen = nlen + nv(m)*nv(m)
               enddo
               nlen = sqrt(nlen)
               kappa = 0.5d0*nlen/area(i,j)
            endif
            do m = 1,3
               surfnormals(i,j,m) = nv(m)/nlen
            enddo
            curvature(i,j) = kappa
         enddo
      enddo

      end

c ---------------------------------------------------------
c> @brief Averages the area from a fine grid to a coarse grid
c>
c> @param[in] mx, my the number of cells in the x and y directions
c> @param[in] mbc the number of ghost cells
c> @param[in] areafine the cell areas of the fine patch
c> @param[out] areacoarse the cell areas of the coarse patch
c> @param[in] igrid the index of the fine grid in the child array
c ---------------------------------------------------------
      subroutine fclaw2d_fort_average_area(mx,my,mbc,
     &      areacoarse, areafine, igrid)
      implicit none

      integer mx,my,mbc,p4est_refineFactor, refratio, igrid

c     # these will be empty if we are not on a manifold.
      double precision areacoarse(-mbc:mx+mbc+1,-mbc:my+mbc+1)
      double precision   areafine(-mbc:mx+mbc+1,-mbc:my+mbc+1)

      integer i,j, ig, jg, ic_add, jc_add, ii, jj
      double precision sum

c     # This should be refratio*refratio.
      integer i1,j1, r2, m
      integer rr2
      parameter(rr2 = 4)
      integer i2(0:rr2-1),j2(0:rr2-1)
      double precision kf

      p4est_refineFactor = 2
      refratio = 2

c     # 'iface' is relative to the coarse grid

      r2 = refratio*refratio
      if (r2 .ne. rr2) then
         write(6,*) 'average_face_ghost (claw2d_utils.f) ',
     &         '  Refratio**2 is not equal to rr2'
         stop
      endif

c     # Get (ig,jg) for grid from linear (igrid) coordinates
      ig = mod(igrid,refratio)
      jg = (igrid-ig)/refratio

c     # Get rectangle in coarse grid for fine grid.
      ic_add = ig*mx/p4est_refineFactor
      jc_add = jg*mx/p4est_refineFactor

      r2 = refratio*refratio
      do j = 0,my/p4est_refineFactor+1
         do i = 0,mx/p4est_refineFactor +1
            i1 = i+ic_add
            j1 = j+jc_add
            m = 0
            do jj = 1,refratio
               do ii = 1,refratio
                  i2(m) = (i-1)*refratio + ii
                  j2(m) = (j-1)*refratio + jj
                  m = m + 1
               enddo
            enddo
            sum = 0
            do m = 0,r2-1
               kf = areafine(i2(m),j2(m))
               sum = sum + kf
            enddo
            areacoarse(i1,j1) = sum
         enddo
      enddo

c     # Compute area in the ghost cells

      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


module fiducial
  use settings,  only: mcp        ! real precision
  use constants, only: const_pi  ! pi
  use numrec
  implicit none
  !values of dm and H used to take into account the fiducial cosmology.
  !Assumes that the fiducial cosmology corresponds to a simple flat LCDM model
  !as the one used for all BOSS measurements.
  real(mcp)       :: save_om_m !used to pass Omega_m to the integrand of D_M(z)
  contains

    subroutine compute_scaling(om_m_fid, h0_fid, z, d_m, hub)!Ariel: added this subroutine
      use camb
      implicit none
      real(kind=mcp), intent(in)  :: om_m_fid, h0_fid, z
      real(kind=mcp), intent(out) :: d_m, hub
      hub  = e(om_m_fid, z)!*h0_fid  !H(z)/100.
      d_m  = dm(om_m_fid, z)!/h0_fid !D_M(z)
    end subroutine compute_scaling

    function e(om_m, z)    !E(z)=H(z)/H_0
      implicit none
      real(kind=mcp)   :: om_m, z
      real(kind=mcp)   :: e2, e, zplus, a
      zplus  = 1._mcp + z
      a      = 1._mcp / zplus
      e2 = om_m/a**3 + (1._mcp- om_m)
      e = sqrt(e2)
    end function e

    function einv(z)
      implicit none
      real(kind=mcp) :: einv, z
      einv = 1._mcp/e(save_om_m, z)
    end function einv

    function dm(om_m, z) !for a flat cosmology dm = r(z)
      implicit none
      real(kind=mcp)   :: dm, om_m, z, ans
      save_om_m = om_m
      call qromb(einv, 0._mcp, z, ans)
      dm = ans*2997.92458_mcp
    end function dm

end module fiducial


text:
! Test the MVBITS intrinsic subroutine
INTEGER*4 :: from, to, result
integer*8 :: from8, to8

DATA from / z'0003FFFC' /
DATA to / z'77760000' /
DATA result / z'7777FFFE' /

CALL mvbits(from, 2, 16, to, 1)
if (to /= result) CALL abort()

to8 = 0_8
from8 = b'1011'*2_8**32
call mvbits (from8, 33, 3, to8, 2)
if (to8 /= b'10100') call abort
end

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


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

!----------------------------------------------------------------------
!> Utilities that can be used by any of the obs_def_xxx_mods.
!> In a separate module to obs_def_mod to avoid a circular dependency.
module obs_def_utilities_mod

use types_mod, only : r8, MISSING_R8

implicit none

private

public :: track_status, set_debug_fwd_op

logical :: debug = .false. ! Flag for debugging

contains

!----------------------------------------------------------------------
!> This routine can be used to set the debugging level of track_status.
!> It is called from filter_mod.
subroutine set_debug_fwd_op(setme)

logical, intent(in) :: setme

debug = setme

end subroutine set_debug_fwd_op

!----------------------------------------------------------------------
!> track_status can be used to keep track of the status of
!> each ensemble member during multiple calls to model_interpolate
!> for a given obs_def.
!> It assumes that you are starting with istatus(:) = 0
!> If debugging, return_now is only set to true if all istatuses are non-zero
!> If not debugging, return_now is set to true if any istatues are non-zero
!>  and any remaining zero istatues are set to 1.
subroutine track_status(ens_size, val_istatus, val_data, istatus, return_now)

integer,  intent(in)    :: ens_size
integer,  intent(in)    :: val_istatus(ens_size)
real(r8), intent(inout) :: val_data(ens_size) !! expected_obs for obs_def
integer,  intent(inout) :: istatus(ens_size) !! istatus for obs_def
logical,  intent(out)   :: return_now

where (istatus == 0) istatus = val_istatus
where (istatus /= 0) val_data = MISSING_R8

return_now = .false.
if (debug) then
   if( all(istatus /= 0))then
      return_now = .true.
      val_data(:) = missing_r8
   endif
else
   if( any(istatus /= 0)) then
      return_now = .true.
      val_data(:) = missing_r8
      where (istatus == 0) istatus = 1
   endif
endif

end subroutine track_status
!----------------------------------------------------------------------

end module obs_def_utilities_mod


! chk_quadcomplex.f90 --
!     Check: does the compiler support complex numbers with larger precision than double precision?
!
!     Note: if not, then the compilation will fail (kind < 0)
!
program chk_quadcomplex
    implicit none

    integer, parameter :: quadruple = selected_real_kind(precision(1.0d0)+1)
    complex(quadruple) :: quadruple_value = (1.0_quadruple,1.0_quadruple)

    write(*,'(a,i0)') 'Precision for "quadruple" complex numbers: ', precision(quadruple_value)
    write(*,'(a,i0)') 'Range for "quadruple" complex numbers:     ', range(quadruple_value)
end program chk_quadcomplex


text:
c# ioncap.f
      subroutine voacap(igraph)
c******************************************************************
c            this routine executes VOACAP
c******************************************************************
      use voacapl_defs
      use crun_directory
      INCLUDE 'ficepac.hdr'
      common /cbotlines/ nbotlines,linesbot(14)
      common /ctoplines/ ntoplines,linestop( 7)
c jw      common /crun_directory/ run_directory
c jw         character run_directory*50
c******************************************************************
      character tns*1,tew*1,rns*1,rew*1,alf*80,antfile_name*21
      character path(2)*1,coefflist(2)*4
      data path/'S','L'/
      data coefflist/'CCIR','URSI'/

      nch_run=lcount(run_directory,50)
      if(igraph.eq.0) then
         open(31,file=run_directory(1:nch_run)//PATH_SEPARATOR//'voacapx.dat')
      else
         open(31,file=run_directory(1:nch_run)//PATH_SEPARATOR//'voacapg.dat')
      end if
      rewind(31)
c***************************************************************
c          Copy defaults from VOACAP.DEF
      open(32,file=run_directory(1:nch_run-3)//'database'//PATH_SEPARATOR//'voacap.def',
     +     status='old',err=100)
90    read(32,'(a)',end=99) alf
      nch=lcount(alf,80)
      write(31,'(a)') alf(1:nch)
      go to 90
99    close(32)
c***************************************************************
100   write(31,101) coefflist(icoeffs)(1:4)
101   format('COEFFS    ',a4)
c***************************************************************
      if(igraph.eq.0) then               !  circuit
         if(method.eq.11 .or. method.eq.28) then
            write(31,3) 1,24,1,1         !  MUF-LUF plots
         else
            iflag=-1                     !  LMT
            if(itimecode.eq.1) iflag=1   !  UT
            write(31,3) ihr1,ihr2,ihrinc,iflag
3           format('TIME      ',4i5)
         end if
      else                               !  graph, force 24 hours UT
         write(31,3) 1,24,1,1
      end if

      do 120 n=1,10
      if(montha(n).eq.0.) go to 125
120   continue
      n=11
125   n=n-1
ccc      write(31,4) year,(int(montha(i)),i=1,n)
ccc4     format('MONTH     ',11i5)
      write(31,4) year,(montha(i),i=1,n)
4     format('MONTH     ',i5,10f5.2)
      write(31,2) (float(ssna(i)),i=1,n)
2     format('SUNSPOT   ',12f5.0)

      call getTR(tlat,tlon,rlat,rlon)
      call benchfix(rname,rlat,rlon)      !  fix coordinates for benchmark
      write(31,5) tname,rname
5     format('LABEL     ',2a20)
      tlatd=tlat
      tns='N'
      if(tlatd.lt.0.) tns='S'
      tlond=tlon
      if(tlond.gt.180.) tlond=tlond-360.
      tew='E'
      if(tlond.lt.0.) tew='W'
      rlatd=rlat
      rns='N'
      if(rlatd.lt.0.) rns='S'
      rlond=rlon
      if(rlond.gt.180.) rlond=rlond-360.
      rew='E'
      if(rlond.lt.0.) rew='W'
      write(31,6) abs(tlatd),tns,abs(tlond),tew,
     +            abs(rlatd),rns,abs(rlond),rew,path(ipath),ipath-1
6     format('CIRCUIT   ',f5.2,a1,3(f9.2,a1),2x,a1,1x,i5)
      xnoise=noise
      if(noise.eq.1) xnoise=-1.       !  -140.4
      if(noise.eq.2) xnoise=-2.       !  -144.7
      if(noise.eq.3) xnoise=-3.       !  -150.0
      if(noise.eq.4) xnoise=-4.       !  -164.1
      write(31,111) 1.0,xnoise,amind,float(xlufp),float(rsn),pmp,dmpx
111   format('SYSTEM    ',f5.0,f5.0,f5.2,f5.0,3f5.2)
      write(31,13) fprob
13    format('FPROB     ',4f5.2)
      do 10 i=1,numants
      call antfile(xmtr_dir(i),xmtr_file(i),antfile_name)
10    write(31,11) 1,i,minfreq(i),maxfreq(i),design_freq(i),
     +             antfile_name,beam_main(i),TxPower(i)
11    format('ANTENNA   ',4i5,f10.3,1h[,a21,1h],f5.1,f10.4)
      call antfile(rec_dir,rec_file,antfile_name)
      write(31,11) 2,numants+1,2,30,0.,antfile_name,rec_bear,rec_gain

      if(igraph.eq.0) then
         write(31,7) (float(freq(i))/1000.,i=1,11)
7        format('FREQUENCY ',11f5.2)
      else
            write(31,17) (freq(i),i=1,11)
17          format('COMMENT   FREQ ',11i5)
ccc         write(31,8) 26,0                 !  method 26 is LUF-MUF table
ccc         write(31,'(7HEXECUTE)')
         write(31,7) -1.                  !  do all frequencies for graph
      end if

      npage0=0
      if(igraph.eq.0) then
         write(31,8) method,npage0
8        format('METHOD    ',2i5)
      else
         met=method
         if(met.lt.16 .or. met.gt.23) met=16
         if(method.eq.30) met=30      !  method 30 is allowed for graphs
         write(31,8) met,npage0
      end if
      if(method.eq.23) then
         if(nbotlines.eq.0) then
            nbotlines=7
            linesbot(1)=1         !  MODE
            linesbot(2)=2         !  ANGLE
            linesbot(3)=19        !  Tgain
            linesbot(4)=7         !  DBU
            linesbot(5)=10        !  SNR
            linesbot(6)=21        !  SNR90
            linesbot(7)=12        !  REL
         end if
         write(31,23) (linesbot(i),i=1,nbotlines)
23       format('BOTLINES  ',14i5)
         if(ntoplines.eq.0) then
            ntoplines=7
            do 123 i=1,7
123         linestop(i)=i        !  turn on all
         end if
         write(31,24) (linestop(i),i=1,ntoplines)
24       format('TOPLINES  ',7i5)

      end if

      write(31,'(7HEXECUTE)')

      write(31,'(4HQUIT)')
      close(31)
ccc      if(igraph.eq.0) then
ccc         if(method.eq.11 .or.method.eq.28) stop 25   !  graphics plots
ccc         stop 5
ccc      else
ccc         stop 15
ccc      end if
      end
c------------------------------------------------------------------

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl