text:
C    @(#)calcintflo.f	20.4 7/18/96
C****************************************************************
C
C   File: calcintflo.f
C
C   Purpose: Routine to calculate selected area intertie "I" flows.
c
C   Author: Walt Powell
C   Date: 5 March 1996
C   Called by: gtoutput.f
C
C****************************************************************
C
        integer function calcintflo (ka1, ka2, sched_export,
     &      actual_export, int_switch)

        include 'ipfinc/parametr.inc'
        include 'ipfinc/blank.inc'
        include 'ipfinc/bus.inc'
        include 'ipfinc/area.inc'
        include 'ipfinc/alpha.inc'
        include 'ipfinc/arcntl.inc'
        include 'ipfinc/dc2t.inc'
        include 'ipfinc/dcmt.inc'
        include 'ipfinc/intchg.inc'
        include 'ipfinc/slnopt.inc'
        include 'ipfinc/prt.inc'
        include 'ipfinc/lfiles.inc'
        include 'ipfinc/ordsta.inc'

        common /is_batch / is_batch

        common /area_data/ gentot, lodtot, lostot, arcact(4*MAXCAR),
     &                     arcflg(4*MAXCAR)
        real gentot, lodtot, lostot, arcact
        integer arcflg

        logical found

        sched_export = 0.0
        actual_export = 0.0
        int_switch = 0
        calcintflo = 0

        jt = 1
        do while (jt .le. ntotic .and.
     &            arcint(1,jt) .ne. arcnam(ka1))
           jt = jt + 1
        enddo
        if (jt .gt. ntotic .or.
     &      arcint(1,jt) .ne. arcnam(ka1)) go to 100
        do while (jt .le. ntotic .and.
     &            arcint(1,jt) .eq. arcnam(ka1) .and.
     &            arcint(2,jt) .ne. arcnam(ka2))
           jt = jt + 1
        enddo
        if (jt .le. ntotic .and.
     &      arcint(1,jt) .eq. arcnam(ka1) .and.
     &      arcint(2,jt) .eq. arcnam(ka2)) then
          sched_export = arcinp(jt)
          actual_export = arcact(jt)
          int_switch = 1
        endif
  100   return
        end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
c=======================================================================
c
c     subroutine MULTIVOLRFR
c
c     Risk Budgeting allocation strategy with a risk-free asset
c
c     w(i), i=1,...,n the riskly assets (portfolio weights)
c     w(0), the risk-free asset
c
c     Max[ w'*rho - k*w'*Q*w]
c      s.t.
c     sqrt[w'*cov_i*w] <= sigma(i)  (i-th risk budgeting constraint)
c     C*w <= b                      (general linear constraints)
c     Cinf(i) <= w(i) <= Csup(i)    (riskly assets lower/upper bounds)
c     Cinf(0) <= w(0) <= Csup(0)    (risk-free asset lower/upper bounds)
c
c        Q   : covariance matrix
c        rho : expected returns
c        rfr : risk-free rate
c
c-----------------------------------------------------------------------
      SUBROUTINE multivolrfr ( n, cov, Q, kappa, rho, rfr,
     &                         neq, nin, ccst, bcst,
     &                         cinf, csup, cinfrfr, csuprfr,
     &                         nbc, class, sigma,
     &                         iwork, dwork, wopt, wrf, lambda, info )
c-----------------------------------------------------------------------
c
c     INPUT
c       n      : portfolio size                                  integer
c       cov    : covariance matrix (n*n)                          double
c       Q      : risk aversion matrix (n*n)                       double
c       kappa  : risk aversion coefficient                        double
c       rho    : expected returns vector (n)                      double
c       rfr    : expected risk-free rate                          double
c       neq    : number equality constraints                     integer
c       nin    : number inequality constraints                   integer
c       ccst   : matrix of constraints (nasset*(neq+nin))         double
c       bcst   : vector initial of constraints (neq+nin)          double
c       cinf   : lower bound (n)                                  double
c       csup   : upper bound (n)                                  double
c       cinfrfr: risk-free lower bound                            double
c       csuprfr: risk-free upper bound                            double
c       nbc    : number of risk budgeting constraints (>=1)      integer
c       class  : block definition (n)                            integer
c       sigma  : volatilities budgets constraints (nbc)           double
c
c     WORKSPACE
c       iwork  : 2*nbc + 5*n + neq + 2*tnin + 7                  integer
c              = 2*nbc + 5*n + neq + 2*nin + 11
c       dwork  : nbc*(nbc+21)/2
c              + n*(8*n+2*neq+2*tnin+22)+5*neq+7*tnin+16          double
c              = nbc*(nbc+21)/2
c              + n*(8*n+2*neq+2*nin+26)+ 5*neq+7*nin +30
c
c     OUTPUT
c       wopt   : optimal portfolio (n)                            double
c       wrf    : optimal risk-free asset                          double
c       lambda : optimal lambda (nbc)                             double
c       info   : diagnostic argument                             integer
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, info, neq, nin, nbc
      INTEGER class(*)
      DOUBLE PRECISION wrf, rfr, kappa, cinfrfr, csuprfr
      DOUBLE PRECISION cov(*), Q(*), rho(*), cinf(*), csup(*),
     &                 ccst(n,*), bcst(*), wopt(*), sigma(*), lambda(*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER i, j, npk, totitr, totsim, ncst, ntot, tnin,
     &        infoBFGS, infot, IZERO, test
      INTEGER piwbf, pisim, pinfo, pin, pineq, pinin, piclass, piqp,
     &        pydual, pbinf, pbsup, pgrad, pwbf, pdsim,
     &        pdcov, pdrho, pdcinf, pdcsup, pdsigma, pdc, pdb,
     &        pdcov1, pdQ, pdkappa, pdcinf0, pdcsup0, pdqp, pdlagr
      DOUBLE PRECISION funct, epsbfg, sum, scal, tmp
      DOUBLE PRECISION INFINI, MINFINI, ZERO, EPS, ONE, COEF, TOL
      PARAMETER (IZERO=0, ZERO=0.E0,EPS=1.E-15,ONE=1.E0,COEF=-0.5E0)
      PARAMETER (TOL=1.E-15, INFINI=1.E+15, MINFINI = -1.E+10)
c
c     external subroutines
      EXTERNAL simmvolrfr, gesterr, testsdp, OVTMCV, IVX, YV,
     &         bfgsmvol, checklinbox, initfeas, IVI, SEV
c
c     intrinsic functions
c
c-----------------------------------------------------------------------
c
c
c     initializations
      info = 0          ! diagnostic argument
      ntot = neq + nin  ! number total of input constraint
      tnin = nin + 2    ! add risk-free constraints (lower/upper)
      ncst = neq + tnin ! number total of constrains
c
c     pointers for integer work space : iwork
c     ---------------------------------------
      piwbf = 1
c     piwbf : pointer for BFGSBOX who needs (2*nbc + 1)
      pisim  = piwbf + (2*nbc + 1)
c     pisim : pointer for simulator SIMMVOLRFR (5*n + neq + 2*tnin + 6)
      pinfo = pisim
c     pinfo : pointer for info (errors of simul), so 1 more
      pin   = pinfo + ( 1 )
c     pin   : pointer for value of n, so 1 more
      pineq = pin + ( 1 )
c     pineq : pointer for value of neq, so 1 more
      pinin = pineq + ( 1 )
c     pinin : pointer for value of nin, so 1 more
      piclass = pinin + ( 1 )
c     piclass: pointer for value of class, so n more
      piqp = piclass + ( n )
c     piqp : pointer for QP (SIMMVOLRFR)
c
c     Total size of dwork array = (2*nbc + 1) + (5*n + neq + 2*tnin + 6)
c                               = 2*nbc + 5*n + neq + 2*nin + 10
c
c     pointers for double precision work space  : dwork
c     -------------------------------------------------
      pydual = 1
c     pydual: pointer for initial value of ydual (nbc)
      pbinf = pydual + ( nbc )
c     pbinf : pointer for inferior bounds binf (nbc)
      pbsup = pbinf + ( nbc )
c     pbsup : pointer for superior bounds bsup (nbc)
      pgrad = pbsup + ( nbc )
c     pgrad : pointer for gradient at solution (nbc)
      pwbf = pgrad + ( nbc )
c     pwbf  : pointer for BFGSMVOL ( nbc*(nbc+11)/2 )
      pdsim = pwbf + ( nbc*(nbc+11)/2 )
c     pdsim  : pointer for the simulator SIMMVOLRFR
c              n*(7*n+2*neq+2*tnin+21)+nbc+4*neq+6*tnin+16
      pdcov   = pdsim
c     pdcov : pointer for cov. matrix (n*n)
      pdQ     = pdcov + ( n*n )
c     pdQ   : pointer for Q matrix (n*n)
      pdkappa = pdQ + ( n*n )
c     pdcov : pointer for kappa (1)
      pdrho   = pdkappa + ( 1 )
c     pdrho : pointer for expected returns (n)
      pdcinf  = pdrho + ( n )
c     pdcinf : pointer for lower bounds (n)
      pdcsup  = pdcinf + ( n )
c     pdcsup : pointer for upper bounds (n)
      pdcinf0 = pdcsup + ( n )
c     pdcinf0 : pointer for risk-free lower bounds (1)
      pdcsup0 = pdcinf0 + ( 1 )
c     pdcsup0 : pointer for risk-free upper bounds (1)
      pdsigma = pdcsup0 + ( 1 )
c     pdsigma : pointer for risk budget constraint (nbc)
      pdc     = pdsigma + ( nbc )
c     pdc : pointer for matrix constraint (n*(neq + tnin))
      pdb     = pdc + ( n*(neq + tnin) )
c     pdb : pointer for vector constraint (neq + tnin)
      pdqp    = pdb + ( neq + tnin )
c     pdqp : pointer for QP (SIMMVOLRFR), n*(neq+tnin+2*n+15)+2*neq +4*tnin+13
      pdlagr  = pdqp + ( n*(5*n+neq+tnin+18)+3*neq+5*tnin+13 )
c     pdlagr : pointer for QP Lagrange parameter (n + neq + tnin)
      pdcov1  = pdlagr + ( n + neq + tnin )
c     + n*n
c
c     Total size of dwork array = 4*nbc + ( nbc*(nbc+11)/2 )
c                               + n*(7*n+2*neq+2*tnin+21)+nbc+4*neq+6*tnin+16 (SIMMVOLRFR)
c                               + n*n
c                               + n+neq+tnin
c              = nbc*(nbc+21)/2 + n*(8*n+2*neq+2*tnin+22)+5*neq+7*tnin+16
c-----------------------------------------------------------------------
c
c      TODO - test if classes are not all null
c
c-----------------------------------------------------
c     satisfiability conditions (well-defined problem)
c-----------------------------------------------------
c
c     test if the covariance matrix is SDP
      CALL testsdp (n, cov, iwork(pisim), dwork(pdsim), info)
      IF (info .NE. 0) THEN
         info = -108
         RETURN
      ENDIF
      IF (kappa .LT. EPS) THEN
        CALL YV ( n*n, cov, dwork(pdQ) )
      ELSE
        CALL YV ( n*n, Q, dwork(pdQ) )
        CALL testsdp (n, dwork(pdQ), iwork(pisim), dwork(pdsim), info)
        IF (info .NE. 0) THEN
            info = -109
            RETURN
        ENDIF
      ENDIF
c
c     case sigma(1)=0
      IF (nbc .EQ. 1) THEN
        IF (ABS(sigma(1)) .LT. EPS) THEN
            CALL IVX ( n, wopt, ZERO )
            IF (csuprfr .GT. (1.0-10*EPS)) THEN
                wrf = 1.0
            ELSE
                wrf  = cinfrfr
                info = 110
            ENDIF
            RETURN
        ENDIF
      ENDIF
c
c     problem initialization and feasibility check
      CALL initfeas ( n, cinf, csup, neq, nin, ccst, bcst,
     &               iwork(piqp), dwork(pdqp), wopt, info )
      IF (info .EQ. 1111) THEN
        info = -100
        RETURN
      ENDIF
c
c-----------------------------------------------------------------
c     construction of work vectors (communication with simmvolrfr)
c-----------------------------------------------------------------
c
c     initialization
      CALL IVX ( nbc, dwork(pydual), ONE )   ! dual
      CALL IVX ( nbc, dwork(pbinf), ZERO )   ! lower bounds
      CALL IVX ( nbc, dwork(pbsup), INFINI ) ! upper bounds
c
c     construct of iwork (communication with simmvol)
      iwork(pinfo) = 0       ! diagnostic argument
      iwork(pin)   = n       ! number of asset(s)
      iwork(pineq) = neq     ! number of eq. constraint(s)
      iwork(pinin) = nin     ! number of ineq. constraint(s)
      CALL YVI ( n, class, iwork(piclass) ) ! cov. matrix partition
c
c     construct dwork (communication with simmvol)
      CALL YV ( n*n, cov, dwork(pdcov) )     ! cov. matrix
      CALL YV ( n*n, Q, dwork(pdQ) )         ! Q matrix (risk aversion)
      CALL YV ( n, rho, dwork(pdrho) )       ! rho
      DO i = 1,n
        dwork(pdrho + i - 1) = dwork(pdrho + i - 1) - rfr
      ENDDO
      dwork(pdkappa) = kappa                 ! kappa coef.
      CALL YV ( n, cinf, dwork(pdcinf) )     ! lower bounds
      CALL YV ( n, csup, dwork(pdcsup) )     ! upper bounds
      CALL YV ( 1, cinfrfr, dwork(pdcinf0) ) ! risk-free lower bound
      CALL YV ( 1, csuprfr, dwork(pdcsup0) ) ! risk-free upper bound
      CALL YV ( nbc, sigma, dwork(pdsigma) ) ! risk budget(s)
c
c     case rho(1) = rho(2) = ... = rho(n)
      tmp = rho(1)
      test = 0
      DO i = 2,n
       IF (ABS(tmp - rho(i)) .GT. EPS) THEN
        test = 1
       ENDIF
      ENDDO
      IF (test .EQ. 0) THEN
        dwork(pdrho + 1) = dwork(pdrho + 1)
     &                   + 1.E-4*(MAX(dwork(pdrho + 1),1.E-6))
      ENDIF
c
c     case ncst <> 0 (neq=0, nin=0)
      IF (ntot .NE. 0) THEN
        CALL YM ( n, ntot, ccst, dwork(pdc) )  ! constraint matrix
        CALL YV ( ntot, bcst, dwork(pdb) )    ! constraint vector
      ENDIF
c
c     linear constraints (C*w <= b)
      DO i = 1,n
        dwork(pdc + ntot*n + i - 1)     =  1.0
        dwork(pdc + (ntot+1)*n + i - 1) = -1.0
      ENDDO
      dwork(pdb + ntot)     = 1.0 - cinfrfr
      dwork(pdb + ntot + 1) = csuprfr - 1.0
c
c-------------------------------------------------------------------
c     Non Linear optimization (BFGS)
c-------------------------------------------------------------------
      epsbfg = 1.E-10   ! precision stop test
      CALL bfgsmvol ( simmvolrfr, gesterr, nbc, dwork(pydual), epsbfg,
     &                dwork(pbinf), dwork(pbsup),
     &                iwork(piwbf), dwork(pwbf),
     &                funct, dwork(pgrad), totitr, totsim, info)
      CALL YV ( nbc, dwork(pydual), lambda ) ! ydual -> lambda
      infoBFGS = info ! BFGS diagnostic argument
c
c-------------------------------------------------
c     Optimize with optimal dual solution (lambda)
c-------------------------------------------------
c
      CALL IMX ( n, n, dwork(pdcov), ZERO ) ! initialize pdcov
c
c     construct cov. matrix with optimal dual solution
      DO i = 1,nbc
         CALL IVI(n, iwork(piwbf), IZERO) ! initialize piwbf
c
c        dimension of block i
         npk = 0
         DO j = 1,n
            IF (class(j) .EQ. i) THEN
               npk = npk + 1
               iwork(piwbf + npk - 1) = j
            ENDIF
         ENDDO
         IF (npk .GT. 1) THEN
c
c           initialize temporary cov. matrix
            CALL IMX ( npk, npk, dwork(pdcov1), ZERO )
c
c           extract block i
            CALL YMCPI ( n, cov, npk, iwork(piwbf), dwork(pdcov1), info)
c
c           product lambda(i)*cov(i)
            CALL PMX (npk, npk, dwork(pdcov1), lambda(i), dwork(pdcov1))
c
c           sub-block -> cov. matrix
            CALL YMCPIR (npk, dwork(pdcov1), n, iwork(piwbf),
     &                   dwork(pdcov), info)
         ENDIF
      ENDDO
c
c     product kappa*Q -> pdcov1
      CALL PMX (n, n, Q, kappa, dwork(pdcov1))
c
c     kappa*Q + Sum(Cov_i)
      CALL SM ( n, n, dwork(pdcov), dwork(pdcov1), dwork(pdcov))
c
c     linear part: l(i) = -0.5*rho(i)
      CALL PVX ( n, rho, COEF, dwork(pdrho) )
c
c     problem initialization and feasibility check
      CALL initfeas ( n, cinf, csup, neq, tnin, dwork(pdc), dwork(pdb),
     &                iwork(piqp), dwork(pdqp), wopt, info)
      IF (info .EQ. 1111) THEN
        info = -100
        RETURN
      ENDIF
c
c     quadratic solver
      CALL qp ( n,dwork(pdcov), dwork(pdrho), neq, tnin,
     &          dwork(pdc), dwork(pdb), cinf, csup,
     &          iwork(piqp), dwork(pdqp),
     &          dwork(pdlagr), wopt, info )
c
c     risk-free asset
      CALL SEV ( n, wopt, sum )
      wrf = 1.0 - sum
c
c-------------------------------------------------------------
c     test if the ith-volatility budget constraint is too high/small
c     (budget constraint not attainable) -> info = 111/110
c-------------------------------------------------------------
c
      IF (info .NE. 1001) THEN
        DO i = 1,nbc
            IF (ABS(lambda(i)) .LT. 1.E-15) THEN
                info = 111
            ENDIF
            IF (ABS(lambda(i)) .GT. 1.E+14) THEN
                info = 110
            ENDIF
        ENDDO
      ENDIF
c
c     checking non-emptiness of linear constraints
      CALL checklinbox ( n, wopt, neq, tnin,
     &                   dwork(pdc), dwork(pdb),
     &                   cinf, csup,
     &                   dwork(pwbf), infot )
      IF (infot .LT. 0) THEN
        info = infot
        RETURN
      ENDIF
c
c     construct cov. matrix with optimal dual solution
      DO i = 1,nbc
         DO j = 1,n
            iwork(piwbf + i - 1) = 0
         ENDDO
c
c        dimension of the i-th block
         npk = 0
         DO j = 1,n
            IF (class(j) .EQ. i) THEN
               npk = npk + 1
               iwork(piwbf + npk - 1) = j
            ENDIF
         ENDDO
         IF (npk .GT. 1) THEN
c
c           initialize temporary cov. matrix
            CALL IMX ( npk, npk, dwork(pdcov1), ZERO )
            CALL IMX ( n, n, dwork(pdcov), ZERO )
c
c           extracting the i-th block
            CALL YMCPI(n, cov, npk, iwork(piwbf), dwork(pdcov1), infot)
c
c           sub-block -> cov. matrix
            CALL YMCPIR (npk, dwork(pdcov1), n, iwork(piwbf),
     &                   dwork(pdcov), infot)
c
c           variance of the i-th block, w'*Cov*w
            CALL OVTMCV(n, dwork(pdcov), wopt, scal)
            scal = SQRT(scal)
            IF (scal .GT. (sigma(i)+1.E-4)) THEN
                info = -120
                RETURN
            ENDIF
         ENDIF
      ENDDO
c
c     BFGS diagnostic (very complex case)
      IF ((infoBFGS .NE. 0) .AND. (info .EQ. 0)) THEN
        info = infoBFGS
      ENDIF
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
program P10799
    implicit none
    character(len = 100000) :: parens
    integer :: stack, sum
    integer :: i

    read *, parens
    sum = 0
    stack = 0
    do i = 1, len_trim(parens)
        if (parens(i:i) == '(') then
            stack = stack + 1
        else if (parens(i:i) == ')') then
            stack = stack - 1
            if (parens((i - 1):(i - 1)) == '(') then
                sum = sum + stack
            else
                sum = sum + 1
            end if
        end if
    end do
    write (*, '(I0)') sum
end program P10799

max_stars_repo_name: daily-boj/kiwiyou

max_issues_repo_name: daily-boj/kiwiyou

max_forks_repo_name: daily-boj/kiwiyou


text:
C
C $Id: agsrch.f,v 1.8 2008-07-27 00:14:35 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE AGSRCH (TPID,IPID,IKWL,UKWL,LKWL)
C
      CHARACTER*(*) TPID,UKWL,LKWL
C
C The routine AGSRCH is used by AGSCAN to search a parameter identifier
C for the next keyword and return the index of that keyword in a list of
C keywords.  It has the following arguments.
C
C -- TPID is the parameter identifier, a character string.
C
C -- IPID is the index of the last character examined in TPID.  It is
C    updated by AGSRCH to point to the first slash or period following
C    the next keyword.
C
C -- IKWL is returned containing the index (in the keyword list) of the
C    next keyword in the parameter identifier (list length, plus one,
C    if the keyword is not found in the list).
C
C -- UKWL is the keyword list - 4*NOKW characters in all.
C
C -- LKWL is identical to UKWL, but in lower case.
C
C ICHR is used to hold up to four characters of a keyword.
C
      CHARACTER*4 ICHR
C
C LPID is the assumed maximum length of a parameter identifier.
C
      DATA LPID / 100 /
C
C Compute the number of 4-character keywords in the keyword list.
C
      NOKW=LEN(UKWL)/4
C
C Find the next non-blank in the parameter identifier.
C
  101 IPID=IPID+1
      IF (IPID.GT.LPID) GO TO 107
      IF (TPID(IPID:IPID).EQ.' ') GO TO 101
C
C Pick up at most four characters of the keyword, stopping on the first
C blank, slash, or period encountered.
C
      NCHR=0
C
  102 IF (TPID(IPID:IPID).EQ.' '.OR.
     +    TPID(IPID:IPID).EQ.'/'.OR.
     +    TPID(IPID:IPID).EQ.'.') GO TO 103
C
      NCHR=NCHR+1
      ICHR(NCHR:NCHR)=TPID(IPID:IPID)
C
      IPID=IPID+1
C
      IF (NCHR.LT.4) GO TO 102
C
C If the keyword found has zero length, error.
C
  103 IF (NCHR.EQ.0) GO TO 107
C
C Scan ahead for the next slash or period.
C
  104 IF (TPID(IPID:IPID).EQ.'/'.OR.TPID(IPID:IPID).EQ.'.') GO TO 105
C
      IPID=IPID+1
      IF (IPID.GT.LPID) GO TO 107
      GO TO 104
C
C Search the keyword list for the keyword found.
C
  105 DO 106 I=1,NOKW
        IKWL=I
        ISTR=(I-1)*4+1
        IEND=(I-1)*4+NCHR
        IF (ICHR(1:NCHR).EQ.UKWL(ISTR:IEND).OR.
     +      ICHR(1:NCHR).EQ.LKWL(ISTR:IEND)) RETURN
  106 CONTINUE
C
C Keyword not found - set IKWL to impossible value and return.
C
  107 IKWL=NOKW+1
      RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
!-----------------------------------------------------------------------
! Skeleton 2D Electrostatic OpenMP/Vector PIC code
! written by Viktor K. Decyk, UCLA and Ricardo Fonseca, ISCTE
      program vmpic2
      use iso_c_binding
      use sselib2_c
      use ssempush2_c
      use vmpush2_h
      use omplib_h
      implicit none
! indx/indy = exponent which determines grid points in x/y direction:
! nx = 2**indx, ny = 2**indy.
      integer, parameter :: indx =   9, indy =   9
! npx/npy = number of electrons distributed in x/y direction.
      integer, parameter :: npx =  3072, npy =   3072
! ndim = number of velocity coordinates = 2
      integer, parameter :: ndim = 2
! tend = time at end of simulation, in units of plasma frequency.
! dt = time interval between successive calculations.
! qme = charge on electron, in units of e.
      real, parameter :: tend = 10.0, dt = 0.1, qme = -1.0
! vtx/vty = thermal velocity of electrons in x/y direction
! vx0/vy0 = drift velocity of electrons in x/y direction.
      real, parameter :: vtx = 1.0, vty = 1.0, vx0 = 0.0, vy0 = 0.0
! ax/ay = smoothed particle size in x/y direction
      real :: ax = .912871, ay = .912871
! idimp = number of particle coordinates = 4
! ipbc = particle boundary condition: 1 = periodic
      integer :: idimp = 4, ipbc = 1
! wke/we/wt = particle kinetic/electric field/total energy
      real, target :: wke = 0.0, we = 0.0
      real :: wt = 0.0
! mx/my = number of grids in x/y in sorting tiles
      integer :: mx = 16, my = 16
! xtras = fraction of extra particles needed for particle management
      real :: xtras = 0.2
! kvec = (1,2) = run (autovector,SSE2) version
      integer :: kvec = 1

! declare scalars for standard code
      integer :: np, nx, ny, nxh, nyh, nxe, nye, nxeh, nxyh, nxhy
      integer :: mx1, my1, mxy1, ntime, nloop, isign, lvect
      integer :: irc = 0
      real :: qbme, affp
!
! declare scalars for OpenMP code
      integer :: nppmx, nppmx0, ntmax, npbmx
      integer :: nvp
!
! declare arrays for standard code:
! part = original particle array
      real, dimension(:,:), pointer :: part => null()
! qe = electron charge density with guard cells
      real, dimension(:,:), pointer :: qe => null()
! fxye = smoothed electric field with guard cells
      real, dimension(:,:,:), pointer :: fxye => null()
! ffc = form factor array for poisson solver
      complex, dimension(:,:), pointer :: ffc => null()
! mixup = bit reverse table for FFT
      integer, dimension(:), pointer :: mixup => null()
! sct = sine/cosine table for FFT
      complex, dimension(:), pointer :: sct => null()
!
! declare arrays for OpenMP (tiled) code:
! ppartt = tiled particle array
! ppbuff = buffer array for reordering tiled particle array
      real, dimension(:,:,:), pointer :: ppartt => null()
      real, dimension(:,:,:), pointer :: ppbuff => null()
! kpic = number of particles in each tile
      integer, dimension(:), pointer :: kpic => null()
! ncl = number of particles departing tile in each direction
      integer, dimension(:,:), pointer :: ncl => null()
! ihole = location/destination of each particle departing tile
      integer, dimension(:,:,:), pointer :: ihole => null()
! kp = original location of reordered particle
      integer, dimension(:,:), pointer :: kp => null()
!
! declare and initialize timing data
      real :: time
      type, bind(C) :: timeval
         integer(c_long) :: tv_sec
         integer(c_long) :: tv_usec
      end type
      type (timeval) :: itime
      double precision :: dtime
      real :: tdpost = 0.0, tguard = 0.0, tfft = 0.0, tfield = 0.0
      real :: tpush = 0.0, tsort = 0.0
!
      interface
         subroutine dtimer(time,itime,icntrl) bind(C,name='dtimer')
         use iso_c_binding
         import :: timeval
         implicit none
         real(c_double) :: time
         type (timeval) :: itime
         integer(c_int), value :: icntrl
         end subroutine
      end interface
!
      irc = 0
! nvp = number of shared memory nodes (0=default)
      nvp = 0
!     write (*,*) 'enter number of nodes:'
!     read (5,*) nvp
! initialize for shared memory parallel processing
      call INIT_OMP(nvp)
!
! initialize scalars for standard code
! np = total number of particles in simulation
! nx/ny = number of grid points in x/y direction
      np = npx*npy; nx = 2**indx; ny = 2**indy; nxh = nx/2; nyh = ny/2
      nxe = nx + 2; nye = ny + 1; nxeh = nxe/2
      nxyh = max(nx,ny)/2; nxhy = max(nxh,ny)
! mx1/my1 = number of tiles in x/y direction
      mx1 = (nx - 1)/mx + 1; my1 = (ny - 1)/my + 1; mxy1 = mx1*my1
! nloop = number of time steps in simulation
! ntime = current time step
      nloop = tend/dt + .0001; ntime = 0
      qbme = qme
      affp = real(nx*ny)/real(np)
!
! allocate data for standard code
      allocate(part(idimp,np))
      allocate(mixup(nxhy),sct(nxyh))
      allocate(kpic(mxy1))
!
      lvect = 4
! allocate vector field data
      nxe = lvect*((nxe - 1)/lvect + 1)
      nxeh = nxe/2
      call sse_f2allocate(qe,nxe,nye,irc)
      call sse_f3allocate(fxye,ndim,nxe,nye,irc)
      call sse_c2allocate(ffc,nxh,nyh,irc)
      if (irc /= 0) then
         write (*,*) 'aligned field allocation error: irc = ', irc
      endif
!
! prepare fft tables
      call WFFT2RINIT(mixup,sct,indx,indy,nxhy,nxyh)
! calculate form factors
      isign = 0
      call VMPOIS22(qe,fxye,isign,ffc,ax,ay,affp,we,nx,ny,nxeh,nye,nxh, &
     &nyh)
! initialize electrons
      call DISTR2(part,vtx,vty,vx0,vy0,npx,npy,idimp,np,nx,ny,ipbc)
!
! find number of particles in each of mx, my tiles: updates kpic, nppmx
      call DBLKP2L(part,kpic,nppmx,idimp,np,mx,my,mx1,mxy1,irc)
      if (irc /= 0) then
         write (*,*) 'DBLKP2L error, irc=', irc
         stop
      endif
! allocate vector particle data
      nppmx0 = (1.0 + xtras)*nppmx
      ntmax = xtras*nppmx
      npbmx = xtras*nppmx
! align data for Vector Processor
      nppmx0 = lvect*((nppmx0 - 1)/lvect + 1)
      ntmax = lvect*(ntmax/lvect + 1)
      npbmx = lvect*((npbmx - 1)/lvect + 1)
      call sse_f3allocate(ppartt,nppmx0,idimp,mxy1,irc)
      call sse_f3allocate(ppbuff,npbmx,idimp,mxy1,irc)
      allocate(ncl(8,mxy1))
      allocate(ihole(2,ntmax+1,mxy1))
      allocate(kp(nppmx0,mxy1))
      if (irc /= 0) then
         write (*,*) 'aligned particle allocation error: irc = ', irc
      endif
!
! copy ordered particle data for OpenMP: updates ppartt, kpic, and kp
      call PPMOVIN2LTP(part,ppartt,kpic,kp,nppmx0,idimp,np,mx,my,mx1,   &
     &mxy1,irc)
      if (irc /= 0) then
         write (*,*) 'PPMOVIN2LTP overflow error, irc=', irc
         stop
      endif
! sanity check
      call PPCHECK2LT(ppartt,kpic,idimp,nppmx0,nx,ny,mx,my,mx1,my1,irc)
      if (irc /= 0) then
         write (*,*) 'PPCHECK2LT error: irc=', irc
         stop
      endif
!
! * * * start main iteration loop * * *
!
  500 if (nloop <= ntime) go to 2000
!     write (*,*) 'ntime = ', ntime
!
! deposit charge with OpenMP: updates qe
      call dtimer(dtime,itime,-1)
      qe = 0.0
      if (kvec==1) then
         call VGPPOST2LT(ppartt,qe,kpic,qme,nppmx0,idimp,mx,my,nxe,nye, &
     &mx1,mxy1)
! SSE2 function
      else if (kvec==2) then
         call csse2gppost2lt(c_loc(ppartt(1,1,1)),c_loc(qe(1,1)),       &
     &c_loc(kpic(1)),qme,nppmx0,idimp,mx,my,nxe,nye,mx1,mxy1)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tdpost = tdpost + time
!
! add guard cells with OpenMP: updates qe
      call dtimer(dtime,itime,-1)
      if (kvec==1) then
         call AGUARD2L(qe,nx,ny,nxe,nye)
! SSE2 function
      else if (kvec==2) then
         call csse2aguard2l(c_loc(qe(1,1)),nx,ny,nxe,nye)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tguard = tguard + time
!
! transform charge to fourier space with OpenMP: updates qe
      call dtimer(dtime,itime,-1)
      isign = -1
      if (kvec==1) then
         call WFFT2RVMX(qe,isign,mixup,sct,indx,indy,nxeh,nye,nxhy,nxyh)
! SSE2 function
      else if (kvec==2) then
         call csse2wfft2rmx(c_loc(qe(1,1)),isign,c_loc(mixup(1)),       &
     &c_loc(sct(1)),indx,indy,nxeh,nye,nxhy,nxyh)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfft = tfft + time
!
! calculate force/charge in fourier space with OpenMP: updates fxye, we
      call dtimer(dtime,itime,-1)
      isign = -1
      if (kvec==1) then
         call VMPOIS22(qe,fxye,isign,ffc,ax,ay,affp,we,nx,ny,nxeh,nye,  &
     &nxh,nyh)
! SSE2 function
      else if (kvec==2) then
         call csse2mpois22(c_loc(qe(1,1)),c_loc(fxye(1,1,1)),isign,     &
     &c_loc(ffc(1,1)),ax,ay,affp,c_loc(we),nx,ny,nxeh,nye,nxh,nyh)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfield = tfield + time
!
! transform force to real space with OpenMP: updates fxye
      call dtimer(dtime,itime,-1)
      isign = 1
      if (kvec==1) then
         call WFFT2RVM2(fxye,isign,mixup,sct,indx,indy,nxeh,nye,nxhy,   &
     &nxyh)
! SSE2 function
      else if (kvec==2) then
         call csse2wfft2rm2(c_loc(fxye(1,1,1)),isign,c_loc(mixup(1))    &
     &,c_loc(sct(1)),indx,indy,nxeh,nye,nxhy,nxyh)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfft = tfft + time
!
! copy guard cells with OpenMP: updates fxye
      call dtimer(dtime,itime,-1)
      if (kvec==1) then
         call CGUARD2L(fxye,nx,ny,nxe,nye)
! SSE2 function
      else if (kvec==2) then
         call csse2cguard2l(c_loc(fxye(1,1,1)),nx,ny,nxe,nye)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tguard = tguard + time
!
! push particles with OpenMP:
      wke = 0.0
      call dtimer(dtime,itime,-1)
! updates ppartt, wke
!     if (kvec==1) then
!        call VGPPUSH2LT(ppartt,fxye,kpic,qbme,dt,wke,idimp,nppmx0,nx,ny&
!    &,mx,my,nxe,nye,mx1,mxy1,ipbc)
! SSE2 function
!     else if (kvec==2) then
!        call csse2gppush2lt(c_loc(ppartt(1,1,1)),c_loc(fxye(1,1,1)),   &
!    &c_loc(kpic(1)),qbme,dt,c_loc(wke),idimp,nppmx0,nx,ny,mx,my,nxe,nye&
!    &,mx1,mxy1,ipbc)
!     endif
! updates ppartt, ncl, ihole, wke, irc
      if (kvec==1) then
         call VGPPUSHF2LT(ppartt,fxye,kpic,ncl,ihole,qbme,dt,wke,idimp, &
     &nppmx0,nx,ny,mx,my,nxe,nye,mx1,mxy1,ntmax,irc)
! SSE2 function
      else if (kvec==2) then
         call csse2gppushf2lt(c_loc(ppartt(1,1,1)),c_loc(fxye(1,1,1)),  &
     &c_loc(kpic(1)),c_loc(ncl(1,1)),c_loc(ihole(1,1,1)),qbme,dt,       &
     &c_loc(wke),idimp,nppmx0,nx,ny,mx,my,nxe,nye,mx1,mxy1,ntmax,irc)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tpush = tpush + time
      if (irc /= 0) then
         write (*,*) 'VGPPUSHF2LT error: irc=', irc
         stop
      endif
!
! reorder particles by tile with OpenMP:
      call dtimer(dtime,itime,-1)
! updates ppartt, ppbuff, kpic, ncl, ihole, and irc
!     if (kvec==1) then
!        call VPPORDER2LT(ppartt,ppbuff,kpic,ncl,ihole,idimp,nppmx0,nx, &
!    &ny,mx,my,mx1,my1,npbmx,ntmax,irc)
! SSE2 function
!     else if (kvec==2) then
!        call csse2pporder2lt(c_loc(ppartt(1,1,1)),c_loc(ppbuff(1,1,1)),&
!    &c_loc(kpic(1)),c_loc(ncl(1,1)),c_loc(ihole(1,1,1)),idimp,nppmx0,nx&
!    &,ny,mx,my,mx1,my1,npbmx,ntmax,irc)
!     endif
! updates ppartt, ppbuff, kpic, ncl, and irc
      if (kvec==1) then
         call VPPORDERF2LT(ppartt,ppbuff,kpic,ncl,ihole,idimp,nppmx0,mx1&
     &,my1,npbmx,ntmax,irc)
! SSE2 function
      else if (kvec==2) then
         call csse2pporderf2lt(c_loc(ppartt(1,1,1)),c_loc(ppbuff(1,1,1))&
     &,c_loc(kpic(1)),c_loc(ncl(1,1)),c_loc(ihole(1,1,1)),idimp,nppmx0, &
     &mx1,my1,npbmx,ntmax,irc)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tsort = tsort + time
      if (irc /= 0) then
         write (*,*) 'VPPORDERF2LT error: ntmax, irc=', ntmax, irc
         stop
      endif
!
      if (ntime==0) then
         write (*,*) 'Initial Field, Kinetic and Total Energies:'
         write (*,'(3e14.7)') we, wke, wke + we
      endif
      ntime = ntime + 1
      go to 500
 2000 continue
!
! * * * end main iteration loop * * *
!
      write (*,*) 'ntime = ', ntime, 'kvec = ', kvec
      write (*,*) 'Final Field, Kinetic and Total Energies:'
      write (*,'(3e14.7)') we, wke, wke + we
!
      write (*,*)
      write (*,*) 'deposit time = ', tdpost
      write (*,*) 'guard time = ', tguard
      write (*,*) 'solver time = ', tfield
      write (*,*) 'fft time = ', tfft
      write (*,*) 'push time = ', tpush
      write (*,*) 'sort time = ', tsort
      tfield = tfield + tguard + tfft
      write (*,*) 'total solver time = ', tfield
      time = tdpost + tpush + tsort
      write (*,*) 'total particle time = ', time
      wt = time + tfield
      write (*,*) 'total time = ', wt
      write (*,*)
!
      wt = 1.0e+09/(real(nloop)*real(np))
      write (*,*) 'Push Time (nsec) = ', tpush*wt
      write (*,*) 'Deposit Time (nsec) = ', tdpost*wt
      write (*,*) 'Sort Time (nsec) = ', tsort*wt
      write (*,*) 'Total Particle Time (nsec) = ', time*wt
      write (*,*)
!
      call sse_deallocate(c_loc(ppartt(1,1,1)),irc); nullify(ppartt)
      call sse_deallocate(c_loc(ppbuff(1,1,1)),irc); nullify(ppbuff)
      call sse_deallocate(c_loc(ffc(1,1)),irc); nullify(ffc)
      call sse_deallocate(c_loc(fxye(1,1,1)),irc); nullify(fxye)
      call sse_deallocate(c_loc(qe(1,1)),irc); nullify(qe)
!
      stop
      end program

max_stars_repo_name: gcasabona/cuda

max_issues_repo_name: gcasabona/cuda

max_forks_repo_name: gcasabona/cuda


PROGRAM MAIN

  IMPLICIT NONE

  ! Loop counter and maximum value
  INTEGER :: i, max_val
  ! Variables for current value
  ! I choose to use a real. THIS WILL GET A SLIGHTLY WRONG ANSWER BUT
  ! it wont overflow nearly as soon. I use a REAL64 here
  INTEGER, PARAMETER :: mykind = SELECTED_REAL_KIND(6, 37)
  !INTEGER, PARAMETER :: mykind = SELECTED_REAL_KIND(15, 307)
  REAL(KIND=mykind) :: factorial

  ! Do the first 50 terms
  max_val = 50

  ! Optional last but to add- sanity check on max_val
  IF(max_val < 1 .OR. max_val > 50) THEN
    PRINT*, "Sequence is too short or long. max_val should be between 1 and 50"
  END IF

  factorial = 1.0

  ! Loop from first value to max_value, calculating and printing the term
  DO i = 1, max_val

    IF(factorial .GT. HUGE(factorial)/i) THEN
      PRINT*, "Value Overflowing!"
      STOP
    END IF

    factorial = factorial * REAL(i)
    ! Print the calculated value
    PRINT*, "n= ", i, " n!=", factorial

  END DO

END PROGRAM


text:
C
C**************************  MOISTABPIC ********************************
C
C  calculate characteristic curves for unsaturated zone
C  by interpolating the moisture tables (Picard case)
C
C***********************************************************************
C
      SUBROUTINE MOISTABPIC(NLKP,NT,NTRI,TETRA,PTNEW,PNODI,SNODI,
     1                      SWE,CKRWE,ETAE)
C
      IMPLICIT NONE
      INTEGER  IEL,ISTR,IR,MTYPE,I,INOD
      INTEGER  NLKP,NT,NTRI
      INTEGER  TETRA(5,*)
      REAL*8   PAVG,SSAVG,FIAVG,KRLOC,ETALOC,SWLOC
      REAL*8   PTNEW(*),PNODI(*),SNODI(*)
      REAL*8   SWE(*),CKRWE(*),ETAE(*)
      real*8   zero,quarter
      parameter (zero=0.0d0,quarter=0.25d0)
C
      DO IEL=1,NT
         ISTR=1+IEL/(NTRI*3)
         IR=MOD(IEL,NTRI*3)
         IF (IR .EQ. 0) ISTR=ISTR-1
         MTYPE=TETRA(5,IEL)
         PAVG = ZERO
         SSAVG = ZERO
         FIAVG = ZERO
         DO I=1,4
            INOD = TETRA(I,IEL)
            PAVG = PAVG + PTNEW(INOD)
            SSAVG = SSAVG + SNODI(INOD)
            FIAVG = FIAVG + PNODI(INOD)
         END DO
         PAVG = PAVG*QUARTER
         SSAVG = SSAVG*QUARTER
         FIAVG = FIAVG*QUARTER
         CALL INTERPPIC(NLKP,ISTR,MTYPE,
     1                  PAVG,SSAVG,FIAVG,KRLOC,ETALOC,SWLOC)
         CKRWE(IEL)=KRLOC
         ETAE(IEL)=ETALOC
         SWE(IEL)=SWLOC
      END DO
C
      RETURN
      END


PROGRAM bad
IMPLICIT NONE
INTEGER :: i
INTEGER, DIMENSION(3) :: subs = (/ 1, 8, 11 /)
INTEGER, DIMENSION(16), TARGET :: info = (/ (i, i=1,16) /)
INTEGER, DIMENSION(:), POINTER :: ptr1
ptr1 => info(subs)
WRITE (*,'(A,16I3)') ' ptr1 = ', ptr1
END PROGRAM bad


text:
      SUBROUTINE BESSEL (Z,J0,J0P)
C
C     BESSEL EVALUATES THE ZERO-ORDER BESSEL FUNCTION AND ITS DERIVATIVE
C     FOR COMPLEX ARGUMENT Z.
C
      COMPLEX J0,J0P,P0Z,P1Z,Q0Z,Q1Z,Z,ZI,ZI2,ZK,FJ,CZ,SZ,J0X,J0PX
      DIMENSION M(101), A1(25), A2(25), FJX(2)
      EQUIVALENCE (FJ,FJX)
      DATA PI,C3,P10,P20,Q10,Q20/3.141592654,.7978845608,.0703125,.11215
     120996,.125,.0732421875/
      DATA P11,P21,Q11,Q21/.1171875,.1441955566,.375,.1025390625/
      DATA POF,INIT/.7853981635,0/,FJX/0.,1./
      IF (INIT.EQ.0) GO TO 5
1     ZMS=Z*CONJG(Z)
      IF (ZMS.GT.1.E-12) GO TO 2
      J0=(1.,0.)
      J0P=-.5*Z
      RETURN
2     IB=0
      IF (ZMS.GT.37.21) GO TO 4
      IF (ZMS.GT.36.) IB=1
C     SERIES EXPANSION
      IZ=1.+ZMS
      MIZ=M(IZ)
      J0=(1.,0.)
      J0P=J0
      ZK=J0
      ZI=Z*Z
      DO 3 K=1,MIZ
      ZK=ZK*A1(K)*ZI
      J0=J0+ZK
3     J0P=J0P+A2(K)*ZK
      J0P=-.5*Z*J0P
      IF (IB.EQ.0) RETURN
      J0X=J0
      J0PX=J0P
C     ASYMPTOTIC EXPANSION
4     ZI=1./Z
      ZI2=ZI*ZI
      P0Z=1.+(P20*ZI2-P10)*ZI2
      P1Z=1.+(P11-P21*ZI2)*ZI2
      Q0Z=(Q20*ZI2-Q10)*ZI
      Q1Z=(Q11-Q21*ZI2)*ZI
      ZK=CEXP(FJ*(Z-POF))
      ZI2=1./ZK
      CZ=.5*(ZK+ZI2)
      SZ=FJ*.5*(ZI2-ZK)
      ZK=C3*CSQRT(ZI)
      J0=ZK*(P0Z*CZ-Q0Z*SZ)
      J0P=-ZK*(P1Z*SZ+Q1Z*CZ)
      IF (IB.EQ.0) RETURN
      ZMS=COS((SQRT(ZMS)-6.)*31.41592654)
      J0=.5*(J0X*(1.+ZMS)+J0*(1.-ZMS))
      J0P=.5*(J0PX*(1.+ZMS)+J0P*(1.-ZMS))
      RETURN
C     INITIALIZATION OF CONSTANTS
5     DO 6 K=1,25
      A1(K)=-.25/(K*K)
6     A2(K)=1./(K+1.)
      DO 8 I=1,101
      TEST=1.
      DO 7 K=1,24
      INIT=K
      TEST=-TEST*I*A1(K)
      IF (TEST.LT.1.E-6) GO TO 8
7     CONTINUE
8     M(I)=INIT
      GO TO 1
      END

max_stars_repo_name: ViperScience/ViperNEC

max_issues_repo_name: ViperScience/ViperNEC

max_forks_repo_name: ViperScience/ViperNEC


text:
module probdata_module

  ! make the model name enter through the probin file
  use amrex_fort_module, only : rt => amrex_real
  character (len=80), save  :: model_name

  ! arrange storage for read_in model-- not worrying about efficiency,
  ! since this will only be called once
  real(rt)        , allocatable, save ::  hse_r(:), hse_rho(:)
  real(rt)        , allocatable, save ::  hse_t(:), hse_p(:)
  real(rt)        , allocatable, save ::  hse_s(:,:)

  ! hold the state at the top of the initial model for the boundary
  ! conditions
  real(rt)        , save              :: hse_rho_top, hse_T_top
  real(rt)        , save              :: hse_p_top, hse_eint_top
  real(rt)        , allocatable, save :: hse_X_top(:)

  real(rt)        , save :: xmin, xmax, ymin, ymax, zmin, zmax

  real(rt)        , save :: heating_time, heating_rad, &
                            heating_peak, heating_sigma

  ! the prob_type matches MAESTRO test_basestate
  integer         , save :: prob_type

end module probdata_module

max_stars_repo_name: taehoryu/Castro

max_issues_repo_name: taehoryu/Castro

max_forks_repo_name: taehoryu/Castro


text:
! Module for timing passages of code.
!
! Example:
!
!  use timing
!  implicit none
!  timer_state::my_timer
!
!  call timer_init(my_timer)
!  call timer_start(my_timer)
!   ... do stuff ...
!  call timer_stop(my_timer)
!  call timer_output(<file>, my_timer, "my_timer")

module timing
  implicit none

  type timer_state
     integer::count
     real(kind=8)::start
     real(kind=8)::sum
  end type timer_state

end module timing

subroutine timer_init(timer)
  use timing
  implicit none
  type(timer_state)::timer

  timer%count = 0
  timer%start = 0.0
  timer%sum = 0.0
end subroutine

subroutine timer_start(timer)
  use timing
  implicit none
  include 'mpif.h'
  type(timer_state)::timer
  integer::info

  timer%start = MPI_WTIME(info)
end subroutine

subroutine timer_stop(timer)
  use timing
  implicit none
  include 'mpif.h'
  type(timer_state)::timer
  real(kind=8)::end
  integer::info

  if (timer%start.le.0.0) return

  end = MPI_WTIME(info)
  timer%sum = timer%sum + end - timer%start
  timer%count = timer%count + 1
end subroutine

subroutine timer_inc_count(timer, delta)
  use timing
  implicit none
  include 'mpif.h'
  type(timer_state)::timer
  integer::delta

  if (timer%start.le.0.0) return
  timer%count = timer%count + delta
end subroutine

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse


SUBROUTINE printPsi

 !---------------------------------------------------------------
 !
 ! Prints the stored nodal and edge spatial angular flux moments
 !
 !---------------------------------------------------------------
USE invar
USE solvar
IMPLICIT NONE
INTEGER :: i, j, k, l, n, m, g
INTEGER :: ydir,xdir,ys,ye,incy,xs,xe,incx
REAL*8,ALLOCATABLE :: psi(:,:,:,:,:,:)
REAL*8,ALLOCATABLE :: LR (:,:,:,:,:)
REAL*8,ALLOCATABLE :: BT (:,:,:,:,:)
!
! Angular Flux File
!
OPEN ( UNIT = 31 , FILE = "scratch8" , STATUS = "OLD" , ACTION = "READ" )
OPEN ( UNIT = 22 , FILE = "scratch5" , STATUS = "NEW" , ACTION = "WRITE")
!
WRITE(22,121) "Angular Flux Solution:",title
! k: k-th x-moment
! l: l-th y-moment
! i: x-dimension
! j: y-dimension
! n: # of Angle
! m: 1 - Both positve  2 - mu negative  3 - eta negative  4 - both negative
! g: group

ALLOCATE( psi(nx,ny,0:lambda,0:lambda,apo,4) )
!
DO g = 1,ng
   !
   ! Reading from scratch 8
   !
   DO n = 1, apo
      DO ydir = 1, 2
         IF (ydir == 1) THEN
            ys = ny
            ye = 1
            incy = -1
         ELSE IF (ydir == 2) THEN
            ys = 1
            ye = ny
            incy = 1
         END IF
         DO j = ys, ye, incy
            DO xdir = 1, 2
               IF (xdir == 1) THEN
                  xs = nx
                  xe = 1
                  incx = -1
               ELSE IF (xdir == 2) THEN
                  xs = 1
                  xe = nx
                  incx = 1
               END IF
               DO i = xs, xe, incx
                  DO k = 0, lambda
                     DO l = 0, lambda
                        IF      (incy == 1 .and. incx ==  1) THEN
                           READ(31,122) psi(i,j,k,l,n,1)
                        ELSE IF (incy == 1 .and. incx == -1) THEN
                           READ(31,122) psi(i,j,k,l,n,2)
                        ELSE IF (incy == -1 .and. incx ==  1) THEN
                           READ(31,122) psi(i,j,k,l,n,3)
                        ELSE IF (incy == -1 .and. incx == -1) THEN
                           READ(31,122) psi(i,j,k,l,n,4)
                        END IF
                     END DO
                  END DO
               END DO
            END DO
         END DO
      END DO
   END DO
   !
   ! Writing to scratch 5
   !
   DO m = 1, 4
      DO n = 1, apo
         DO j = 1, ny
            DO i = 1, nx
               DO l = 0, lambda
                  DO k = 0, lambda
                     WRITE(22,122) psi(i,j,k,l,n,m)
                  END DO
               END DO
            END DO
         END DO
      END DO
   END DO
   !
END DO
!
DEALLOCATE(psi)
CLOSE ( UNIT = 22 )
CLOSE ( UNIT = 31)
!
! LR
!
OPEN ( UNIT = 23 , FILE = "scratch6" , STATUS = "NEW" , ACTION = "WRITE")
OPEN ( UNIT = 32 , FILE = "scratch9" , STATUS = "OLD" , ACTION = "READ" )
WRITE(23,121) "Angular Flux RL Edge Solution:",title
ALLOCATE( LR(nx,ny,0:lambda,n,m) )
DO g = 1, ng
   !
   ! READ
   !
   DO n = 1, apo
      DO ydir = 1, 2
         IF (ydir == 1) THEN
            ys = ny
            ye = 1
            incy = -1
         ELSE IF (ydir == 2) THEN
            ys = 1
            ye = ny
            incy = 1
         END IF
         DO j = ys, ye, incy
            DO xdir = 1, 2
               IF (xdir == 1) THEN
                  xs = nx
                  xe = 1
                  incx = -1
               ELSE IF (xdir == 2) THEN
                  xs = 1
                  xe = nx
                  incx = 1
               END IF
               DO i = xs, xe, incx
                  DO l = 0, lambda
                        IF      (incy == 1 .and. incx ==  1) THEN
                           READ(32,122) LR(i,j,l,n,1)
                        ELSE IF (incy == 1 .and. incx == -1) THEN
                           READ(32,122) LR(i,j,l,n,2)
                        ELSE IF (incy == -1 .and. incx ==  1) THEN
                           READ(32,122) LR(i,j,l,n,3)
                        ELSE IF (incy == -1 .and. incx == -1) THEN
                           READ(32,122) LR(i,j,l,n,4)
                        END IF
                  END DO
               END DO
            END DO
         END DO
      END DO
   END DO
   !
   ! Write
   !
   DO m = 1, 4
      DO n = 1, apo
         DO j = 1, ny
            DO i = 1, nx
               DO l = 0, lambda
                  WRITE(23,122) LR(i,j,l,n,m)
               END DO
            END DO
         END DO
      END DO
   END DO
   !
END DO
DEALLOCATE(LR)
CLOSE( UNIT = 23)
!
! Bottom Top
!
OPEN ( UNIT = 24 , FILE = "scratch7" , STATUS = "NEW" , ACTION = "WRITE")
OPEN ( UNIT = 33 , FILE = "scratch10" , STATUS = "OLD" , ACTION = "READ" )
WRITE(24,121) "Angular Flux RL Edge Solution:",title
ALLOCATE( BT(nx,ny,0:lambda,n,m) )
DO g = 1, ng
   !
   ! READ
   !
   DO n = 1, apo
      DO ydir = 1, 2
         IF (ydir == 1) THEN
            ys = ny
            ye = 1
            incy = -1
         ELSE IF (ydir == 2) THEN
            ys = 1
            ye = ny
            incy = 1
         END IF
         DO j = ys, ye, incy
            DO xdir = 1, 2
               IF (xdir == 1) THEN
                  xs = nx
                  xe = 1
                  incx = -1
               ELSE IF (xdir == 2) THEN
                  xs = 1
                  xe = nx
                  incx = 1
               END IF
               DO i = xs, xe, incx
                  DO l = 0, lambda
                        IF      (incy == 1 .and. incx ==  1) THEN
                           READ(33,122) BT(i,j,l,n,1)
                        ELSE IF (incy == 1 .and. incx == -1) THEN
                           READ(33,122) BT(i,j,l,n,2)
                        ELSE IF (incy == -1 .and. incx ==  1) THEN
                           READ(33,122) BT(i,j,l,n,3)
                        ELSE IF (incy == -1 .and. incx == -1) THEN
                           READ(33,122) BT(i,j,l,n,4)
                        END IF
                  END DO
               END DO
            END DO
         END DO
      END DO
   END DO
   !
   ! Write
   !
   DO m = 1, 4
      DO n = 1, apo
         DO j = 1, ny
            DO i = 1, nx
               DO l = 0, lambda
                  WRITE(24,122) BT(i,j,l,n,m)
               END DO
            END DO
         END DO
      END DO
   END DO
   !
END DO
DEALLOCATE(BT)
CLOSE( UNIT = 24)
CLOSE( UNIT = 33)

121 FORMAT(1X,A,1X,A)
122 FORMAT(1ES24.15)

RETURN
END SUBROUTINE printPsi


text:
subroutine nsc_EndBouope(a)
!-----------------------------------------------------------------------
!
! DESCRIPTION
!    The routine calculates forces on bodies
!
!-----------------------------------------------------------------------
   use typre
   use Mod_Mesh
   use Mod_Element
   use Mod_Memor
   use Mod_iofile
   use Mod_int2str
   use Mod_NSCompressible
   use Mod_nsc_nscbc_ex
   use MPI

   implicit none
   class(NSCompressibleProblem) :: a
   class(FiniteElement), pointer :: e => NULL()

   real(rp), allocatable :: tract(:),forceR(:,:),momenR(:,:)
   real(rp), allocatable :: gpcod(:)
   real(rp), allocatable :: boden(:), bomom(:,:), boene(:)
   real(rp), allocatable :: elden(:), elmom(:,:), elene(:)
   real(rp), allocatable :: gpmom(:)
   real(rp), allocatable :: grden(:), grmom(:,:), grene(:)
   real(rp), allocatable :: qflx(:), mgrm(:)
   real(rp), allocatable :: incre(:)
   real(rp), allocatable :: bomom0(:,:), gpmom0(:)
   real(rp), allocatable :: boden0(:),boene0(:)
   real(rp), allocatable :: borhs(:,:)

   character(150) :: fil_force
   real(rp)    :: weightfactor
   real(rp)    :: dsurf,dsurf0
   real(rp)    :: acvis,actco,accph,accvh
   real(rp)    :: aux,invcvh,invgpd
   real(rp)    :: gpden(1),gpene(1),gppre,divm,mgrd
   real(rp)    :: gpden0(1),gpene0(1)
   real(rp)    :: tempn0(a%ndofn)
   integer(ip) :: idime,iboun,ndime,igaub,jdime,nbody,nboun,nelem,ibody
   integer(ip), pointer :: lbody => NULL()
   integer(ip) :: icount,npoinLocal,inodb,ipoin,idofn,ierr,npoin

   logical            :: isALE

   !Initializations
   call a%Mesh%GetNelem(nelem)
   call a%Mesh%GetNdime(ndime)
   call a%Mesh%GetNboun(nboun)
   call a%Mesh%GetNbody(nbody)
   call a%Mesh%GetNpoin(npoin)
   call a%Mesh%GetNpoinLocal(npoinLocal)
   call a%Mesh%GetALE(isALE)

   if(a%kfl_openforce) then
      if (a%MPIrank == a%MPIroot) then
         do ibody=1,nbody
            fil_force = trim(a%OutputFolder)//'/'//adjustl(trim(a%namda))//adjustl(trim(int2str(a%MPIrank)))//'.'//adjustl(trim(a%exmod))//adjustl(trim(int2str(ibody)))//'.frc'
            call iofile(zero,a%lun_force(ibody),fil_force,adjustl(trim(a%exmod))//'FORCES & MOMENTS')
            if (ndime==2) write (a%lun_force(ibody),12)
            if (ndime==3) write (a%lun_force(ibody),13)
         end do
         a%kfl_openforce = .false.
      endif
   end if

   call a%Mesh%ElementAlloc(e,a%Memor,'DefaultRule','nsc_EndBouope')

   !Allocation

   call a%Memor%alloc(ndime+1,tract,'tract','nsc_EndBouope')
   call a%Memor%alloc(ndime,gpcod,'gpcod','nsc_EndBouope')
   call a%Memor%alloc(e%mnodb,boden,'boden','nsc_EndBouope')
   call a%Memor%alloc(e%ndime,e%mnodb,bomom,'bomom','nsc_EndBouope')
   call a%Memor%alloc(e%mnodb,boene,'boene','nsc_EndBouope')
   call a%Memor%alloc(e%mnode,elden,'elden','nsc_EndBouope')
   call a%Memor%alloc(e%ndime,e%mnode,elmom,'elmom','nsc_EndBouope')
   call a%Memor%alloc(e%mnode,elene,'elene','nsc_EndBouope')
   call a%Memor%alloc(e%ndime,gpmom,'gpmom','nsc_EndBouope')
   call a%Memor%alloc(e%ndime,grden,'grden','nsc_EndBouope')
   call a%Memor%alloc(e%ndime,e%ndime,grmom,'grmom','nsc_EndBouope')
   call a%Memor%alloc(e%ndime,grene,'grene','nsc_EndBouope')
   call a%Memor%alloc(e%ndime,qflx,'qflx','nsc_EndBouope')
   call a%Memor%alloc(e%ndime,mgrm,'mgrm','nsc_EndBouope')
   call a%Memor%alloc(a%ndofn,incre,'incre','nsc_bouope_im')
   call a%Memor%alloc(e%ndime,e%mnodb,bomom0,'bomom0','nsc_EndBouope')
   call a%Memor%alloc(e%ndime,gpmom0,'gpmom0','nsc_EndBouope')
   call a%Memor%alloc(e%mnodb,boden0,'boden0','nsc_EndBouope')
   call a%Memor%alloc(e%mnodb,boene0,'boene0','nsc_EndBouope')
   !nscnbc(boundary problem)
   call a%Memor%alloc(a%ndofn,npoin,borhs,'borhs','nsc_EndBouope')
   !working variables (reuced variables)
   call a%Memor%alloc(ndime+1,nbody,forceR,'forceR','nsc_EndBouope')
   call a%Memor%alloc(3,nbody,momenR,'momenR','nsc_EndBouope')

   !Inializations
   a%btraction=0.0_rp
   !Force and moments fields
   if(a%npp_stepi(11) /= 0) then
      a%forfl=0.0_rp
      a%momfl=0.0_rp
      a%hfxfl=0.0_rp
   end if
   if(a%kfl_outfm==1) then
   !Force and moments
      a%force=0.0_rp
      a%formo=0.0_rp
   end if
   !Reduced values
   forceR=0.0_rp
   momenR=0.0_rp
   !Working variables
   ibody=0
   !nscnbc
   borhs=0.0_rp

   ! Loop over boundaries
   boundaries: do iboun=1,nboun
      !Load Element
      call a%Mesh%BoundaryLoad(iboun,e)

      call a%Mesh%GetLbody(lbody,iboun)

      ibody=lbody

      !to delete repeated entities
      icount = 0
      do inodb = 1,e%pnodb
         if (e%lnodb(inodb) <= npoinLocal) then
            icount = icount +1
         endif
      enddo
      weightfactor = real(icount)/real(e%pnodb)

      if((a%kfl_fixbo(iboun)==10) .or. (a%kfl_fixbo(iboun)==12) .or. ((ibody>0).and.(a%kfl_outfm==1)))then
         !Physical Parameters
         call a%GetPhysicalParameters(acvis,actco,accph,accvh)

         call e%elmdel

         !Gathers in boundary elements
         call e%gatherb(1_ip,boden,a%densf(:,1))
         call e%gatherb(ndime,bomom,a%momen(:,:,1))
         call e%gatherb(1_ip,boene,a%energ(:,1))
         call e%gather(1_ip,elden,a%densf(:,1))
         call e%gather(ndime,elmom,a%momen(:,:,1))
         call e%gather(1_ip,elene,a%energ(:,1))
         call e%gatherb(1_ip,boden0,a%densf(:,3))
         call e%gatherb(ndime,bomom0,a%momen(:,:,3))
         call e%gatherb(1_ip,boene0,a%energ(:,3))

         dsurf0 = 0.0_rp

         !Gauss-Point Loop
         do igaub=1,e%pgaub
            e%igaub = igaub

            !Traction at gauss point
            tract=0.0_rp
            divm=0.0_rp

            !Calculate exterior Normal
            call e%bounor

            dsurf=e%weigb(e%igaub)*e%eucta
            dsurf0 = dsurf0 + dsurf

            call e%interpb(ndime,e%bocod,gpcod)
            call e%interpb(1_ip,boden,gpden)
            call e%interpb(ndime,bomom,gpmom)
            call e%interpb(1_ip,boene,gpene)

            !Derivatives at the boundary
            call e%elmderb
            call e%gradientb(1_ip,elden,grden)
            call e%gradientb(ndime,elmom,grmom)
            call e%gradientb(1_ip,elene,grene)

            !Momentum divergence
            do idime =1,e%ndime
              divm = divm + grmom(idime,idime)
            end do

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

            !Auxiliary variables
            invcvh = 1.0_rp / accvh
            invgpd = 1.0_rp / gpden(1)
            aux = dot_product(gpmom,gpmom)

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

            !Traction value
            !  t=n·S  S=-p*I+tau
            if (a%lawde == 1) then
               !  Ideal state law pressure
               !  p=(gamma-1)*(s- m·m/2rho)

               gppre = (accph*invcvh - 1.0_rp)*(gpene(1) - (aux*invgpd*0.5_rp))

            else if (a%lawde /= 1) then
               call runend('Nsc_EndBouope: Non-ideal state law not ready')
            endif

            !  tau=(2mu/rho)*gradS(m)-(2mu/3rho)*div(m)*I
            !      -(2mu/rho^2)*gradS(rho)+(2mu/3rho^2)*m·grad(rho)*I

            mgrd = dot_product(gpmom,grden)

            do idime=1,ndime
                tract(idime)= tract(idime) - e%baloc(idime,ndime)*gppre&
                -(2/3)*acvis*invgpd*divm*e%baloc(idime,ndime)&
                +(2/3)*acvis*invgpd*invgpd*mgrd*e%baloc(idime,ndime)
               do jdime=1,ndime
                  tract(idime) = tract(idime) &
                     + acvis*invgpd*e%baloc(jdime,ndime)*(grmom(idime,jdime)+grmom(jdime,idime))&
                     - acvis*invgpd*invgpd*e%baloc(jdime,ndime)*(gpmom(idime)*grden(jdime)&
                           + gpmom(jdime)*grden(idime))
               end do
            end do

            !Now we fill the traction vector
            do inodb = 1,e%pnodb
                ipoin= e%lnodB(inodb)
                a%btraction(1:ndime,ipoin) = a%btraction(1:ndime,ipoin) + tract(1:ndime)*e%shapb(inodb,e%igaub)*dsurf*weightfactor
            end do

            !Heat flux:
            !   Q=q·n
            if (a%lawde == 1) then
               !   Ideal state law heat flux
               !   q=-actco*(1/rhoCv)*[grad(s)-(1/rho)m·grad(m)^T
               !             +(m·m/rho^2)*grad(rho)-(s/rho)*grad(rho)]

               mgrm = matmul(gpmom,grmom)
               qflx = -actco*invgpd*invcvh*(grene - invgpd*mgrm&
                    + (aux*invgpd*invgpd - gpene(1)*invgpd)*grden)

            else if (a%lawde /= 1) then
               call runend('Nsc_EndBouope: Non-ideal state law not ready')
            endif

            tract(ndime+1) = dot_product(qflx,e%baloc(:,e%ndime))

            !Forces and moments to global fields
            if(a%npp_stepi(11) /= 0) then
               do inodb=1,e%pnodb
                  ipoin = e%lnodb(inodb)
                  do idime=1,ndime
                   a%forfl(idime,ipoin)=a%forfl(idime,ipoin)+dsurf*tract(idime)*e%shapb(inodb,e%igaub)*weightfactor
                  end do
                  if(ndime==2)then
                     a%momfl(1,ipoin) = a%momfl(1,ipoin) &
                        + tract(2)*(gpcod(1)-a%origm(1))*dsurf*e%shapb(inodb,e%igaub)*weightfactor &
                        - tract(1)*(gpcod(2)-a%origm(2))*dsurf*e%shapb(inodb,e%igaub)*weightfactor
                  elseif(ndime==3)then
                     a%momfl(1,ipoin) = a%momfl(1,ipoin) &
                        + tract(3)*(gpcod(2)-a%origm(2))*dsurf*e%shapb(inodb,e%igaub)*weightfactor &
                        - tract(2)*(gpcod(3)-a%origm(3))*dsurf*e%shapb(inodb,e%igaub)*weightfactor
                     a%momfl(2,ipoin) = a%momfl(2,ipoin) &
                        + tract(1)*(gpcod(3)-a%origm(3))*dsurf*e%shapb(inodb,e%igaub)*weightfactor &
                        - tract(3)*(gpcod(1)-a%origm(1))*dsurf*e%shapb(inodb,e%igaub)*weightfactor
                     a%momfl(3,ipoin) = a%momfl(3,ipoin) &
                        + tract(2)*(gpcod(1)-a%origm(1))*dsurf*e%shapb(inodb,e%igaub)*weightfactor &
                        - tract(1)*(gpcod(2)-a%origm(2))*dsurf*e%shapb(inodb,e%igaub)*weightfactor
                  end if
                  a%hfxfl(ipoin)=a%hfxfl(ipoin)+dsurf*tract(ndime+1)*e%shapb(inodb,e%igaub)*weightfactor
               end do
            end if

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

            tract = tract*dsurf*weightfactor

            !Forces
            if(ibody>0)then

               do idime=1,ndime+1
                  a%force(idime,ibody)=a%force(idime,ibody) + tract(idime)
               end do

               !Moments
               if(ndime==2)then
                  a%formo(3,ibody) = a%formo(3,ibody) &
                     + tract(2)*(gpcod(1)-a%origm(1)) &
                     - tract(1)*(gpcod(2)-a%origm(2))
               elseif(ndime==3)then
                  a%formo(1,ibody) = a%formo(1,ibody) &
                     + tract(3)*(gpcod(2)-a%origm(2)) &
                     - tract(2)*(gpcod(3)-a%origm(3))
                  a%formo(2,ibody) = a%formo(2,ibody) &
                     + tract(1)*(gpcod(3)-a%origm(3)) &
                     - tract(3)*(gpcod(1)-a%origm(1))
                  a%formo(3,ibody) = a%formo(3,ibody) &
                     + tract(2)*(gpcod(1)-a%origm(1)) &
                     - tract(1)*(gpcod(2)-a%origm(2))
               end if

            end if

            if(a%kfl_fixbo(iboun)==10) then
               incre = 0.0_rp
               call nsc_nscbc_ex(e,accph,accvh,gpden(1),gpmom,gpene(1),grden,grmom,grene,a%bvnat(iboun)%a,incre)

               do inodb=1,e%pnodb
                  ipoin = e%lnodb(inodb)
                  borhs(1:a%ndofn,ipoin) = borhs(1:a%ndofn,ipoin) + incre(1:a%ndofn)*e%shapb(inodb,e%igaub)*weightfactor*dsurf
               end do

            elseif(a%kfl_fixbo(iboun)==12) then
               call e%interpb(1_ip,boden0,gpden0)
               call e%interpb(ndime,bomom0,gpmom0)
               call e%interpb(1_ip,boene0,gpene0)
               incre = 0.0_rp
               call nsc_nscbc_density_ex(e,a%dtinv,accph,accvh,gpden(1),gpmom,gpene(1),gpden0(1),gpmom0,gpene0(1),grden,grmom,grene,incre)

               do inodb=1,e%pnodb
                  ipoin = e%lnodb(inodb)
                  borhs(1,ipoin) = borhs(1,ipoin) + incre(1)*e%shapb(inodb,e%igaub)*weightfactor*dsurf
               end do

            end if

         end do

      end if

   end do boundaries

  if(nbody>0 .and. a%kfl_outfm==1) then
      !Collect Values
      call  MPI_REDUCE(a%force, forceR, (ndime+1)*nbody, MPI_REAL8, MPI_SUM, a%MPIroot,a%MPIcomm, ierr )
      a%force = forceR
      call  MPI_REDUCE(a%formo, momenR, 3*nbody, MPI_REAL8, MPI_SUM, a%MPIroot,a%MPIcomm, ierr )
      a%formo = momenR
  end if

   ! Output forces
   ! Applied to body => Negative sign.
   if(a%kfl_outfm==1) then
      if (a%MPIrank == a%MPIroot) then

         if(ndime==2) then
            do ibody = 1,nbody
               write(a%lun_force(ibody),20) a%ctime,      &
                  -a%force(1,ibody)/a%adimf(1),    &
                  -a%force(2,ibody)/a%adimf(2),    &
                  -a%formo(3,ibody)/a%adimm(3),    &
                  -a%force(3,ibody)/a%adimf(4)
            end do
         else if(ndime==3) then
            do ibody = 1,nbody
               write(a%lun_force(ibody),30) a%ctime,      &
                  -a%force(1,ibody)/a%adimf(1),    &
                  -a%force(2,ibody)/a%adimf(2),    &
                  -a%force(3,ibody)/a%adimf(3),    &
                  -a%formo(1,ibody)/a%adimm(1),    &
                  -a%formo(2,ibody)/a%adimm(2),    &
                  -a%formo(3,ibody)/a%adimm(3),    &
                  -a%force(4,ibody)/a%adimf(4)
            end do
         end if

      end if
   end if

   if (a%kfl_flush == 1) then
      do ibody = 1,nbody
         call flush(a%lun_force(ibody))
      end do
   endif

   if( a%kfl_nscnbc==1) then
      ! Loop over points
      do ipoin = 1,npoin
         tempn0(1:a%ndofn)=a%bvess(1:a%ndofn,ipoin,2)
         a%bvess(1:a%ndofn,ipoin,2)=(4.0_rp*a%bvess(1:a%ndofn,ipoin,1)-tempn0(1:a%ndofn))/3.0_rp
         a%bvess(1:a%ndofn,ipoin,1)=a%bvess(1:a%ndofn,ipoin,2) - a%bvmass(ipoin)*borhs(1:a%ndofn,ipoin)/a%dtinv
      enddo
      !Communicate between subdomains
      call a%Mesh%ArrayCommunicator%GhostCommunicate(a%ndofn,a%bvess(:,:,1))
   end if

   !Dellocation

   call a%Memor%dealloc(ndime+1,tract,'tract','nsc_EndBouope')
   call a%Memor%dealloc(ndime,gpcod,'gpcod','nsc_EndBouope')
   call a%Memor%dealloc(e%mnodb,boden,'boden','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,e%mnodb,bomom,'bomom','nsc_EndBouope')
   call a%Memor%dealloc(e%mnodb,boene,'boene','nsc_EndBouope')
   call a%Memor%dealloc(e%mnode,elden,'elden','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,e%mnode,elmom,'elmom','nsc_EndBouope')
   call a%Memor%dealloc(e%mnode,elene,'elene','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,gpmom,'gpmom','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,grden,'grden','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,e%ndime,grmom,'grmom','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,grene,'grene','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,qflx,'qflx','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,mgrm,'mgrm','nsc_EndBouope')
   call a%Memor%dealloc(a%ndofn,incre,'incre','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,e%mnodb,bomom0,'bomom0','nsc_EndBouope')
   call a%Memor%dealloc(e%ndime,gpmom0,'gpmom0','nsc_EndBouope')
   call a%Memor%dealloc(e%mnodb,boden0,'boden0','nsc_EndBouope')
   call a%Memor%dealloc(e%mnodb,boene0,'boene0','nsc_EndBouope')
   !nscnbc(boundary problem)
   call a%Memor%dealloc(a%ndofn,npoin,borhs,'borhs','nsc_EndBouope')
   !working variables (reduced variables)
   call a%Memor%dealloc(ndime+1,nbody,forceR,'forceR','nsc_EndBouope')
   call a%Memor%dealloc(3,nbody,momenR,'momenR','nsc_EndBouope')

   !Element Deallocation
   call a%Mesh%ElementDealloc(e,a%Memor,'DefaultRule','nsc_EndBouope')

   !Formats

   12   format( &
          '$        Forces and Moments over bodies          ',/,          &
          '$    Time',15x,'Fx',15x,'Fy',15x,'Mz',15x'Q')
   13   format( &
          '$          Forces and Moments over bodies            ',/,      &
          '$ Time  Fx  Fy  Fz  Mx  My  Mz  Q')

   20   format(1x,e12.6,2x,4(2x,e15.8))
   30   format(1x,e12.6,2x,7(2x,e15.8))

   50   format(a1,i1)
   51   format(a1,i2)
   52   format(a1,i3)
   53   format(5x,'WARNING: CANNOT WRITE MORE THAN 999 bodies')

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
      program while04

C     Check the controlizer

      real x(100)

      i = 100
      do while(i.gt.1)
         if(x(i).gt.10.) go to 200
         x(i) = 0.
         i = i - 1
      enddo

 200  continue

      end


text:
module gravity_1D_module

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

  public

contains

  subroutine ca_test_residual(lo, hi, &
       rhs, rhl1, rhh1,  &
       ecx, ecxl1, ecxh1, &
       dx, problo, coord_type) bind(C, name="ca_test_residual")

    use amrex_constants_module

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer , intent(in   ) :: lo(1),hi(1)
    integer , intent(in   ) :: rhl1, rhh1
    integer , intent(in   ) :: ecxl1, ecxh1
    integer , intent(in   ) :: coord_type
    real(rt), intent(inout) :: rhs(rhl1:rhh1)
    real(rt), intent(in   ) :: ecx(ecxl1:ecxh1)
    real(rt), intent(in   ) :: dx(1),problo(1)

    real(rt)         :: lapphi
    real(rt)         :: rlo,rhi,rcen
    integer          :: i

    ! Cartesian
    if (coord_type .eq. 0) then

       do i=lo(1),hi(1)
          lapphi = (ecx(i+1)-ecx(i)) / dx(1)
          rhs(i) = rhs(i) - lapphi
       enddo

       ! r-z
    else if (coord_type .eq. 1) then

       do i=lo(1),hi(1)
          rlo  = problo(1) + dble(i)*dx(1)
          rhi  = rlo + dx(1)
          rcen = HALF * (rlo+rhi)
          lapphi = (rhi*ecx(i+1)-rlo*ecx(i)) / (rcen*dx(1))
          rhs(i) = rhs(i) - lapphi
       enddo

       ! spherical
    else if (coord_type .eq. 2) then

       do i=lo(1),hi(1)
          rlo  = problo(1) + dble(i)*dx(1)
          rhi  = rlo + dx(1)
          rcen = HALF * (rlo+rhi)
          lapphi = (rhi**2 * ecx(i+1)-rlo**2 * ecx(i)) / (rcen**2 * dx(1))
          rhs(i) = rhs(i) - lapphi
       enddo

    else
       print *,'Bogus coord_type in test_residual ' ,coord_type
       call amrex_error("Error:: Gravity_1d.f90 :: ca_test_residual")
    end if

  end subroutine ca_test_residual

  subroutine ca_compute_radial_mass (lo,hi,dx,dr,&
                                     state,r_l1,r_h1, &
                                     radial_mass,radial_vol,problo, &
                                     n1d,drdxfac,level) bind(C, name="ca_compute_radial_mass")

    use amrex_constants_module, only: HALF, FOUR3RD, M_PI
    use prob_params_module, only: center, Symmetry, physbc_lo, coord_type
    use meth_params_module, only: NVAR, URHO

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer , intent(in   ) :: lo(1), hi(1)
    real(rt), intent(in   ) :: dx(1), dr
    real(rt), intent(in   ) :: problo(1)

    integer , intent(in   ) :: n1d, drdxfac, level
    real(rt), intent(inout) :: radial_mass(0:n1d-1)
    real(rt), intent(inout) :: radial_vol (0:n1d-1)

    integer , intent(in   ) :: r_l1, r_h1
    real(rt), intent(in   ) :: state(r_l1:r_h1,NVAR)

    integer          :: i, index
    integer          :: ii
    real(rt)         :: r
    real(rt)         :: dx_frac, fac, vol
    real(rt)         :: lo_i, rlo, rhi

    if (physbc_lo(1) .ne. Symmetry) then
       call amrex_error("Error: Gravity_1d.f90 :: 1D gravity assumes symmetric lower boundary.")
    endif

    if (coord_type .ne. 2) then
       call amrex_error("Error: Gravity_1d.f90 :: 1D gravity assumes spherical coordinates.")
    endif

    fac = dble(drdxfac)
    dx_frac = dx(1) / fac

    do i = lo(1), hi(1)

       r = abs(problo(1) + (dble(i) + HALF) * dx(1) - center(1))

       index = int(r / dr)

       if (index .gt. n1d-1) then

          if (level .eq. 0) then
             print *,'   '
             print *,'>>> Error: Gravity_1d::ca_compute_radial_mass ', i
             print *,'>>> ... index too big: ', index,' > ',n1d-1
             print *,'>>> ... at i     : ', i
             print *,'    '
             call amrex_error("Error:: Gravity_1d.f90 :: ca_compute_radial_mass")
          end if

       else

          ! Note that we assume we are in spherical coordinates in 1D or we wouldn't be
          ! doing monopole gravity.

          lo_i = problo(1) + dble(i) * dx(1) - center(1)

          do ii = 0, drdxfac-1

             r   = abs(lo_i + (dble(ii  ) + HALF) * dx_frac)
             rlo = abs(lo_i +  dble(ii  )         * dx_frac)
             rhi = abs(lo_i +  dble(ii+1)         * dx_frac)

             vol = FOUR3RD * M_PI * (rhi**3 - rlo**3)

             index = int(r / dr)

             if (index .le. n1d-1) then
                radial_mass(index) = radial_mass(index) + vol * state(i,URHO)
                radial_vol (index) = radial_vol (index) + vol
             end if

          enddo

       endif

    enddo

  end subroutine ca_compute_radial_mass

  subroutine ca_put_radial_grav(lo,hi,dx,dr,&
                                grav,g_l1,g_h1, &
                                radial_grav,problo,n1d,level) bind(C, name="ca_put_radial_grav")

    use prob_params_module, only: center
    use amrex_constants_module, only: HALF, TWO

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer , intent(in   ) :: lo(1), hi(1)
    real(rt), intent(in   ) :: dx(1), dr
    real(rt), intent(in   ) :: problo(1)

    integer , intent(in   ) :: n1d, level
    real(rt), intent(in   ) :: radial_grav(0:n1d-1)

    integer , intent(in   ) :: g_l1, g_h1
    real(rt), intent(inout) :: grav(g_l1:g_h1)

    integer          :: i, index
    real(rt)         :: r, mag_grav
    real(rt)         :: cen, xi, slope, glo, gmd, ghi, minvar, maxvar

    ! Note that we are interpolating onto the entire range of grav,
    ! including the ghost cells. Taking the absolute value of r ensures
    ! that we will get the correct behavior even for the ghost zones with
    ! negative indices, which have a reflecting boundary condition.

    do i = lo(1), hi(1)

       r = abs(problo(1) + (dble(i) + HALF) * dx(1) - center(1))

       index = int(r / dr)

       cen = (dble(index) + HALF) * dr
       xi = r - cen

       if (index == 0) then

          ! Linear interpolation or extrapolation
          slope = ( radial_grav(index+1) - radial_grav(index) ) / dr
          mag_grav = radial_grav(index) + slope * xi

       else if (index == n1d-1) then

          ! Linear interpolation or extrapolation
          slope = ( radial_grav(index) - radial_grav(index-1) ) / dr
          mag_grav = radial_grav(index) + slope * xi

       else if (index .gt. n1d-1) then

          if (level .eq. 0) then
             print *,'PUT_RADIAL_GRAV: INDEX TOO BIG ', index, ' > ', n1d-1
             print *,'AT i ', i
             print *,'R / DR IS ', r, dr
             call amrex_error("Error:: Gravity_1d.f90 :: ca_put_radial_grav")
          else
             ! NOTE: we don't do anything to this point if it's outside the
             !       radial grid and level > 0
          end if

       else

          ! Quadratic interpolation
          ghi = radial_grav(index+1)
          gmd = radial_grav(index  )
          glo = radial_grav(index-1)
          mag_grav = ( ghi -  TWO*gmd + glo)*xi**2/(TWO*dr**2) + &
                     ( ghi       - glo     )*xi   /(TWO*dr   ) + &
                     (-ghi + 26.e0_rt*gmd - glo)/24.e0_rt

          minvar = min(gmd, min(glo,ghi))
          maxvar = max(gmd, max(glo,ghi))
          mag_grav = max(mag_grav,minvar)
          mag_grav = min(mag_grav,maxvar)

       end if

       if (index .le. n1d-1) then

          grav(i) = mag_grav

       end if

    enddo

  end subroutine ca_put_radial_grav

  subroutine ca_put_radial_phi (lo,hi,domlo,domhi,dx,dr,&
                                phi,p_l1,p_h1, &
                                radial_phi,problo, &
                                numpts_1d,fill_interior) bind(C, name="ca_put_radial_phi")

    use prob_params_module, only: center
    use amrex_constants_module, only: HALF, TWO

    use amrex_fort_module, only : rt => amrex_real
    implicit none
    integer , intent(in   ) :: lo(1), hi(1)
    integer , intent(in   ) :: domlo(1), domhi(1)
    real(rt), intent(in   ) :: dx(1), dr
    real(rt), intent(in   ) :: problo(1)

    integer , intent(in   ) :: numpts_1d
    real(rt), intent(in   ) :: radial_phi(0:numpts_1d-1)
    integer , intent(in   ) :: fill_interior

    integer , intent(in   ) :: p_l1, p_h1
    real(rt), intent(inout) :: phi(p_l1:p_h1)

    integer          :: i, index
    real(rt)         :: r
    real(rt)         :: cen, xi, slope, p_lo, p_md, p_hi, minvar, maxvar

    ! Note that when we interpolate into the ghost cells we use the
    ! location of the edge, not the cell center.

    do i = lo(1), hi(1)
       if (i .gt. domhi(1)) then
          r = problo(1) + (dble(i  )     ) * dx(1) - center(1)
       else if (i .lt. domlo(1)) then
          r = problo(1) + (dble(i+1)     ) * dx(1) - center(1)
       else
          r = problo(1) + (dble(i  )+HALF) * dx(1) - center(1)
       end if

       index = int(r / dr)

       if (index .gt. numpts_1d-1) then
          print *,'PUT_RADIAL_PHI: INDEX TOO BIG ', index, ' > ', numpts_1d-1
          print *,'AT (i) ',i
          print *,'R / DR IS ', r, dr
          call amrex_error("Error:: Gravity_1d.f90 :: ca_put_radial_phi")
       end if

       if ( (fill_interior .eq. 1) .or. (i .lt. domlo(1) .or. i.gt.domhi(1)) ) then

          cen = (dble(index) + HALF) * dr
          xi = r - cen

          if (index == 0) then

             ! Linear interpolation or extrapolation
             slope = ( radial_phi(index+1) - radial_phi(index) ) / dr
             phi(i) = radial_phi(index) + slope * xi

          else if (index == numpts_1d-1) then

             ! Linear interpolation or extrapolation
             slope = ( radial_phi(index) - radial_phi(index-1) ) / dr
             phi(i) = radial_phi(index) + slope * xi

          else

             ! Quadratic interpolation
             p_hi = radial_phi(index+1)
             p_md = radial_phi(index  )
             p_lo = radial_phi(index-1)
             phi(i) = ( p_hi -  TWO*p_md + p_lo)*xi**2/(TWO*dr**2) + &
                      ( p_hi       - p_lo      )*xi   /(TWO*dr   ) + &
                      (-p_hi + 26.e0_rt*p_md - p_lo)/24.e0_rt
             minvar = min(p_md, min(p_lo,p_hi))
             maxvar = max(p_md, max(p_lo,p_hi))
             phi(i) = max(phi(i),minvar)
             phi(i) = min(phi(i),maxvar)

          end if

       end if

    enddo

  end subroutine ca_put_radial_phi

end module gravity_1D_module

max_stars_repo_name: harpolea/Castro

max_issues_repo_name: harpolea/Castro

max_forks_repo_name: harpolea/Castro


text:
c     Subroutine to convert variables
c     AJ_Kettle, 26Feb2019
c     04Jul2019: modified to include all 6 variables
c     14Mar2020: modified for USAF update

      SUBROUTINE convert_variables2(f_ndflag,
     +  l_rgh_datalines,l_data,
     +  s_vec_latitude,s_vec_longitude,s_vec_platformheight,
     +  s_vec_winddirection_deg,s_vec_windspeed_ms,
     +  s_vec_airtemperature_c,s_vec_dewpointtemperature_c,
     +  s_vec_sealevelpressure_hpa,s_vec_stationpressure_hpa,

     +  f_vec_latitude,f_vec_longitude,f_vec_platformheight,
     +  f_vec_winddirection_deg,f_vec_windspeed_ms,
     +  f_vec_airtemperature_c,f_vec_dewpointtemperature_c,
     +  f_vec_sealevelpressure_hpa,f_vec_stationpressure_hpa)

c     +  s_vec_latitude,s_vec_longitude,s_vec_platformheight,
c     +  s_vec_windspeed_ms,
c     +  s_vec_sealevelpressure_hpa,s_vec_stationpressure_hpa,

c     +  f_vec_latitude,f_vec_longitude,f_vec_platformheight,
c     +  f_vec_windspeed_ms,
c     +  f_vec_sealevelpressure_hpa,f_vec_stationpressure_hpa)

      IMPLICIT NONE
c************************************************************************
c     Declare variables passed into program

      INTEGER             :: l_rgh_datalines
      INTEGER             :: l_data

      REAL                :: f_ndflag

      CHARACTER(LEN=*)    :: s_vec_latitude(l_rgh_datalines)
      CHARACTER(LEN=*)    :: s_vec_longitude(l_rgh_datalines)
      CHARACTER(LEN=*)    :: s_vec_platformheight(l_rgh_datalines)
      CHARACTER(LEN=*)    :: s_vec_winddirection_deg(l_rgh_datalines)
      CHARACTER(LEN=*)    :: s_vec_windspeed_ms(l_rgh_datalines)
      CHARACTER(LEN=*)    :: s_vec_airtemperature_c(l_rgh_datalines)
      CHARACTER(LEN=*)   :: s_vec_dewpointtemperature_c(l_rgh_datalines)
      CHARACTER(LEN=*)    :: s_vec_sealevelpressure_hpa(l_rgh_datalines)
      CHARACTER(LEN=*)    :: s_vec_stationpressure_hpa(l_rgh_datalines)

      REAL                :: f_vec_latitude(l_rgh_datalines)
      REAL                :: f_vec_longitude(l_rgh_datalines)
      REAL                :: f_vec_platformheight(l_rgh_datalines)
      REAL                :: f_vec_winddirection_deg(l_rgh_datalines)
      REAL                :: f_vec_windspeed_ms(l_rgh_datalines)
      REAL                :: f_vec_airtemperature_c(l_rgh_datalines)
      REAL               :: f_vec_dewpointtemperature_c(l_rgh_datalines)
      REAL                :: f_vec_sealevelpressure_hpa(l_rgh_datalines)
      REAL                :: f_vec_stationpressure_hpa(l_rgh_datalines)
c*****
c     Variables used within program

      INTEGER             :: i,j,k,ii,jj,kk

      INTEGER             :: i_len

      CHARACTER(LEN=100)  :: s_input
      REAL                :: f_output
      INTEGER             :: i_badflag
c************************************************************************
c      print*,'just inside convert_variables'

c     Do conversions
      DO i=1,l_data
c*****
c      Initialize variables with ndflags
       f_vec_latitude(i)            =f_ndflag
       f_vec_longitude(i)           =f_ndflag
       f_vec_platformheight(i)      =f_ndflag
       f_vec_winddirection_deg(i)   =f_ndflag
       f_vec_windspeed_ms(i)        =f_ndflag
       f_vec_airtemperature_c(i)    =f_ndflag
       f_vec_dewpointtemperature_c(i)=f_ndflag
       f_vec_sealevelpressure_hpa(i)=f_ndflag
       f_vec_stationpressure_hpa(i) =f_ndflag
c*****
c*****
c      Latitude
       s_input=s_vec_latitude(i)
       i_len=LEN_TRIM(s_input)
       IF (i_len.GT.0) THEN
        CALL convert_single_string_to_float(s_input,f_ndflag,
     +    f_output,i_badflag)
         f_vec_latitude(i)=f_output

        IF (i_badflag.EQ.1) THEN
         print*,'bad latitude',s_input
         STOP 'convert_variables'
        ENDIF
       ENDIF
c*****
c      Longitude
       s_input=s_vec_longitude(i)
       i_len=LEN_TRIM(s_input)
       IF (i_len.GT.0) THEN
        CALL convert_single_string_to_float(s_input,f_ndflag,
     +    f_output,i_badflag)
         f_vec_longitude(i)=f_output

        IF (i_badflag.EQ.1) THEN
         print*,'bad longitude',s_input
         STOP 'convert_variables'
        ENDIF
       ENDIF
c*****
c      Platformheight
       s_input=s_vec_platformheight(i)
       i_len=LEN_TRIM(s_input)
       IF (i_len.GT.0) THEN
        CALL convert_single_string_to_float(s_input,f_ndflag,
     +    f_output,i_badflag)
         f_vec_platformheight(i)=f_output

        IF (i_badflag.EQ.1) THEN
         print*,'bad platformheight',s_input
         STOP 'convert_variables'
        ENDIF
       ENDIF
c*****
c      1. Winddirection
       s_input=s_vec_winddirection_deg(i)
       i_len=LEN_TRIM(s_input)
       IF (i_len.GT.0) THEN
        CALL convert_single_string_to_float(s_input,f_ndflag,
     +   f_output,i_badflag)
         f_vec_winddirection_deg(i)=f_output

        IF (i_badflag.EQ.1) THEN
         print*,'bad winddirection',s_input
         STOP 'convert_variables'
        ENDIF
       ENDIF
c*****
c      2. Windspeed
       s_input=s_vec_windspeed_ms(i)
       i_len=LEN_TRIM(s_input)
       IF (i_len.GT.0) THEN
        CALL convert_single_string_to_float(s_input,f_ndflag,
     +    f_output,i_badflag)
         f_vec_windspeed_ms(i)=f_output

        IF (i_badflag.EQ.1) THEN
         print*,'bad windspeed',s_input
         STOP 'convert_variables'
        ENDIF
       ENDIF
c*****
c      3. Airtemperature
       s_input=s_vec_airtemperature_c(i)
       i_len=LEN_TRIM(s_input)
       IF (i_len.GT.0) THEN
        CALL convert_single_string_to_float(s_input,f_ndflag,
     +    f_output,i_badflag)
         f_vec_airtemperature_c(i)=f_output

        IF (i_badflag.EQ.1) THEN
         print*,'bad airtemperature',s_input
         STOP 'convert_variables'
        ENDIF
       ENDIF
c*****
c      4. Dewpointtemperature
       s_input=s_vec_dewpointtemperature_c(i)
       i_len=LEN_TRIM(s_input)
       IF (i_len.GT.0) THEN
        CALL convert_single_string_to_float(s_input,f_ndflag,
     +    f_output,i_badflag)
         f_vec_dewpointtemperature_c(i)=f_output

        IF (i_badflag.EQ.1) THEN
         print*,'bad dewpointtemperature',s_input
         STOP 'convert_variables'
        ENDIF
       ENDIF
c*****
c      5. Sealevelpressure
       s_input=s_vec_sealevelpressure_hpa(i)
       i_len=LEN_TRIM(s_input)
       IF (i_len.GT.0) THEN
        CALL convert_single_string_to_float(s_input,f_ndflag,
     +    f_output,i_badflag)
        f_vec_sealevelpressure_hpa(i)=f_output

        IF (i_badflag.EQ.1) THEN
         print*,'bad sealevelpressure',s_input
         STOP 'convert_variables'
        ENDIF
       ENDIF
c*****
c      6. Station pressure
       s_input=s_vec_stationpressure_hpa(i)
       i_len=LEN_TRIM(s_input)
       IF (i_len.GT.0) THEN
        CALL convert_single_string_to_float(s_input,f_ndflag,
     +    f_output,i_badflag)
        f_vec_stationpressure_hpa(i)=f_output

        IF (i_badflag.EQ.1) THEN
         print*,'bad stationpressure',s_input
         STOP 'convert_variables'
        ENDIF
       ENDIF
c*****

      ENDDO

c      print*,'just leaving convert_variables'

c      STOP 'convert_variables'

      RETURN
      END

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


text: #redirect The Creatives Dept.  Advertising, Marketing & Design

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module d_subspace_solv_sl2_module

  implicit none

  private
  public :: d_subspace_solv_sl2

contains

  subroutine d_subspace_solv_sl2( sl, eig )
    use sl_variables, only: slinfo2, sl2, Dsub
    use sl_tools_module, only: slinfo, gather_matrix, distribute_matrix
    use rsdft_allreduce_module, only: rsdft_allreduce
    use rsdft_bcast_module, only: d_rsdft_bcast
    implicit none
    type(slinfo2),intent(inout) :: sl
    real(8),intent(inout) :: eig(:)
    integer :: ierr,nb,i
    integer :: LRWORK,LIWORK,TRILWMIN,NP,NQ
    integer :: itmp(1)
    integer,allocatable :: iwork(:)
    real(8) :: rtmp(1)
    real(8),allocatable :: rwork(:)
    real(8),allocatable :: Vsub(:,:)
    real(8),parameter :: z0=0.0d0
    character(100) :: msg
    character(1) :: ul
    logical :: disp_on
#ifndef _LAPACK_
    integer,external :: NUMROC
    type(slinfo) :: sl_

    call write_border( 1, " d_subspace_solv_sl2(start)" )
    call check_disp_switch( disp_on, 0 )

    ierr = 0
    nb = sl%nband
    ul = sl%uplo

    select case( sl%idiag )
    case( 'PDSYEVD' )

      allocate( Vsub(size(Dsub,1),size(Dsub,2)) ); Vsub=z0

      if ( sl%LWORK == 0 ) then
        call PDSYEVD('V',ul,nb,Dsub,1,1,sl%desca,eig,Vsub,1,1,sl%descz,rtmp,-1,itmp,-1,ierr)
        LRWORK = nint( rtmp(1) )
        LIWORK = itmp(1)
        NP = NUMROC(nb,sl%mbsize,sl%myrow,0,sl%nprow)
        NQ = NUMROC(nb,sl%nbsize,sl%mycol,0,sl%npcol)
        TRILWMIN = 3*nb + max( sl%mbsize*(NP+1), 3*sl%mbsize )
        if ( disp_on ) then
          write(*,'("(PDSYEVD) Work-array sizes (by query)   : LRWORK,LIWORK=",3i8)') LRWORK,LIWORK
          write(*,'("(PDSYEVD) Work-array sizes (in document): LRWORK,LIWORK=",3i8)') &
          max( 1+6*nb+2*NP*NQ, TRILWMIN ) + 2*nb, 7*nb + 8*sl%npcol + 2
        end if
        LRWORK = max( LRWORK, max(1+6*nb+2*NP*NQ,TRILWMIN)+2*nb )
        LIWORK = max( LIWORK, 7*nb+8*sl%npcol+2 )
        if ( disp_on ) then
          write(*,'("(PDSYEVD) Work-array sizes (allocate)   : LRWORK,LIWORK=",3i8)') LRWORK,LIWORK
        end if
        sl%lwork = LRWORK
        sl%lrwork = LRWORK
        sl%liwork = LIWORK
      else
        LRWORK = sl%lrwork
        LIWORK = sl%liwork
      end if

      allocate( rwork(LRWORK) ); rwork=z0
      allocate( iwork(LIWORK) ); iwork=0

      call PDSYEVD('V',ul,nb,Dsub,1,1,sl%desca,eig,Vsub,1,1,sl%descz,rwork,LRWORK,iwork,LIWORK,ierr)

      deallocate( iwork )
      deallocate( rwork )

      Dsub = Vsub

      deallocate( Vsub )

    end select

    if ( ierr /= 0 ) then
      write(msg,'("ierr=",i4,"(stop@d_subspace_solv_sl2)")') ierr
      call stop_program( msg )
    end if

    ! ---
    !
    call d_rsdft_bcast( eig, size(eig), 0, comm_label='grid' )
    call d_rsdft_bcast( eig, size(eig), 0, comm_label='band' )
    !
    ! ---
    !
    if ( allocated(Dsub) ) then

      sl_%nprow  = sl%nprow
      sl_%npcol  = sl%npcol
      sl_%myrow  = sl%myrow
      sl_%mycol  = sl%mycol
      sl_%mbsize = sl%mbsize
      sl_%nbsize = sl%nbsize

      allocate( Vsub(nb,nb) ); Vsub=z0
      call gather_matrix( sl_, Dsub, Vsub )

      deallocate( Dsub )

    end if
    !
    ! ---

    if ( .not.allocated(Vsub) ) then
      allocate( Vsub(nb,nb) ); Vsub=z0
    end if

    call rsdft_allreduce( Vsub, 'g' )
    call rsdft_allreduce( Vsub, 'b' )

    ! ---
    !
    allocate( Dsub(sl2%ldr,sl2%ldc) ); Dsub=z0

    sl_%nprow  = sl2%nprow
    sl_%npcol  = sl2%npcol
    sl_%myrow  = sl2%myrow
    sl_%mycol  = sl2%mycol
    sl_%mbsize = sl2%mbsize
    sl_%nbsize = sl2%nbsize

    call distribute_matrix( sl_, Vsub, Dsub )
    !
    ! ---

    call write_border( 1, " d_subspace_solv_sl2(end)" )

    return

#endif

  end subroutine d_subspace_solv_sl2

end module d_subspace_solv_sl2_module


text:
/*****************************************************************************
 * uvm_stim_factory_test_base.f
 *****************************************************************************/

+UVM_TESTNAME=uvm_stim_factory_test_base


text:
      program prog
      character s1*2, s2*3
      equivalence(s1(1:1), s2(2:2))
      s2 = 'pQr'
      write(6, '(a)') s1
      s1 = 'Ab'
      write(6, '(a)') s2
      end


text:
      subroutine fmulb1(n,h,x,hx,tabaux,nmisaj,prosca,izs,rzs,dzs)
      implicit double precision (a-h,o-z)
      external prosca
c
c
c parametres
      double precision   un     , deux
      parameter (        un=1.0d+0, deux=2.0d+0 )
c declarations
      double precision   h(*), x(n), hx(n), tabaux(n), dzs(*)
      real rzs(*)
      integer  izs(*)
      double precision   uscalx, sscalx, nu, eta, gamma, mu, sigma
      integer      n, nmisaj, memsup, ptnu, compt, iu, is, k
c
c calcul de la longueur d'un bloc
      memsup=2*n+2
c calcul de  h0*x=x=x
      do 1000 k=1,n
      hx(k)=x(k)
1000  continue
c
      if (nmisaj.eq.0) then
      return
      else
      ptnu=1
      compt=1
      endif
c
2000  iu=ptnu+1
      is=iu+n
      do 3000 k=1,n
      tabaux(k)=h(iu+k)
3000  continue
      call prosca(n,tabaux,x,uscalx,izs,rzs,dzs)
      do 4000 k=1,n
      tabaux(k)=h(is+k)
4000  continue
      call prosca(n,tabaux,x,sscalx,izs,rzs,dzs)
      nu=h(ptnu)
      eta=h(ptnu+1)
c calcul du nouveau terme et addition dans hx
      if (compt.eq.1) then
      gamma=eta / nu
      do 5000 k=1,n
      hx(k)=gamma * hx(k)
5000  continue
      mu=sscalx / nu
      sigma=-(deux * sscalx / eta)+(uscalx / nu)
      else
      mu=sscalx / eta
      sigma=-(un + nu / eta)* mu + uscalx / eta
      endif
c
      do 6000 k=1,n
      hx(k)=hx(k) - mu * h(iu+k) - sigma * h(is+k)
6000  continue
c
      compt=compt+1
      if (compt.le.nmisaj) then
      ptnu=ptnu+memsup
      goto 2000
      else
      return
      endif
      end


text:
program fortran_example
    use caliper_mod

    implicit none

    type(Annotation)      :: g_ann
    type(ScopeAnnotation) :: f_ann, w_ann
    type(ConfigManager)   :: mgr

    integer               :: i, count, argc, prop

    logical               :: ret
    character(len=:), allocatable :: errmsg
    character(len=256)    :: arg

    call cali_init()

    mgr = ConfigManager_new()

    argc = command_argument_count()
    if (argc .ge. 1) then
        call get_command_argument(1, arg)
        call mgr%add(arg)
        ret = mgr%error()
        if (ret) then
            errmsg = mgr%error_msg()
            write(*,*) 'ConfigManager: ', errmsg
        endif
    endif

    call mgr%start

    f_ann = ScopeAnnotation_begin('main')
    w_ann = ScopeAnnotation_begin('work')

    prop  = IOR(CALI_ATTR_GLOBAL, CALI_ATTR_ASVALUE)
    g_ann = Annotation_new_with_properties('global.int', prop)
    call g_ann%set_int(42)
    call Annotation_delete(g_ann)

    call ScopeAnnotation_end(w_ann)

    call cali_begin_region('foo')
    call cali_end_region('foo')

    call ScopeAnnotation_end(f_ann)

    call mgr%flush
    call ConfigManager_delete(mgr)
end program fortran_example

max_stars_repo_name: zfrye-llnl/Caliper

max_issues_repo_name: zfrye-llnl/Caliper

max_forks_repo_name: zfrye-llnl/Caliper


! A fortran95 program for G95
! By WQY
program RainDataClean
  implicit none
  integer y(40000), m(40000), d(40000), n, y0, i, j, k1, k2, a, b, c, ndays
  real p(40000)
  character file_in*30, file_out*30

        WRITE(*,*)'******************************************'
        WRITE(*,*)'               RainDataClean              '
        WRITE(*,*)'                                          '
        WRITE(*,*)'        Select event day only using       '
        WRITE(*,*)'      output data from RainData.f90       '
        WRITE(*,*)'                                          '
        WRITE(*,*)'        Technical University Berlin       '
        WRITE(*,*)'           Institute of Ecology           '
        WRITE(*,*)'                                          '
        WRITE(*,*)'            Pedro Alencar, 2019           '
        WRITE(*,*)'                                          '
        WRITE(*,*)'******************************************'

write(*,'(a)')'Insert the name of the input file: '
read(*,'(a30)')file_in
write(*,'(a)')'Insert the name of the output file: '
read(*,'(a30)')file_out
write(*,'(a)')'How many years long is the series?'
read(*,*)n
write(*,'(a)')'What is the startig year'
read(*,*)y0

j = mod(y0,4)
k1 = mod(n,4)
k2 = (n-k1)/4

if (j+k1.gt.4) then
    b = k2+1
    else
        b = k2
end if

ndays = b*366 + (n-b)*365
c = ndays+1

i=1
open(60, file=file_in)
open(70, file=file_out)

do while (i.le.c)
        read(60,*) y(i), m(i), d(i), p(i)
        if  (p(i) .gt. 0.) then
            write(70,*) y(i), m(i), d(i), p(i)
        end if
        i = i+1
end do
close(60)
close(70)

Write(*,*)' '
Write(*,*)' '
Write(*,*)'                                 ***************************'
Write(*,*)'                                 *  Check the output file! *'
Write(*,*)'                                 ***************************'
Write(*,*)' '
Write(*,*)' '

pause
end


text:
subroutine b4step1(mbc,mx,meqn,q,xlower,dx,t,dt,maux,aux)

    ! Called before each call to step1.
    ! Use to set time-dependent aux arrays or perform other tasks.
    !
    ! This default version does nothing.

    implicit none
    integer, intent(in) :: mbc,mx,meqn,maux
    real(kind=8), intent(in) :: xlower,dx,t,dt
    real(kind=8), intent(inout) :: q(meqn,1-mbc:mx+mbc)
    real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc)

end subroutine b4step1

max_stars_repo_name: ClimateImpactLab/classic

max_issues_repo_name: ClimateImpactLab/classic

max_forks_repo_name: cjvogl/classic


text:
Comfort Suites is a hotels hotel located in South Davis.

Hotel Features & Amenities:

  Complimentary Breakfast
  Free HighSpeed Internet Access (Wired & Wireless)
  Outdoor Heated Pool & Hot Tub
  Complimentary Business Center w/ Printer
  100% Smoke Free Hotel
  Meeting Space
  32 LCD Televisions, Cable, Free HBO
  Refrigerator, Microwave, Coffee Maker  all rooms
  Hair Dryer, Iron & Board  all rooms
  Room Types: Single King, Double Queen, Accessible Single Queen, Family Suites
  Close to UC Davis & Downtown Davis
  Recently Renovated  all rooms (2010)

  Looking for a place to stay (or for a guest to stay)?  Davis has many hotels hotel options!

  For even more Davis hotels information (and much more), check out the Yolo County Visitors Bureau!

20120202 23:31:32   (I borrowed a friends laptop, but, I am the one making this post.) I had a very negative experience with the General Manager while staying at the Comfort Suites UC (CA642)from January 13th to Feb 02, 2012. The General Manager Adam had the attitude, and even stated  the room is my room, not yours, I can enter it anytime I want. I had told the front desk clerk, Fouzy, that I thought one of the electrical outlets under the desk was not working. But, also said, I dont like anyone to enter the room when I am not in, when I have the DO NOT DISTURB sign on the door. I left to go to the IHOP restaurant with a neiighbor from the LA QUINTA HOTEL. When I returned the General Manager Adam told me, “you forgot to take your “Do not Disturb” sign off. We checked all the outlets,and they have 120 Vs. You should read the instruction manual for your humidifier. I felt that I had inconvenienced them to check the electrical outlet. This was not any emergency, I was shocked he had entered the room in my absence when the “DO NOT DISTURB” sign was on the door. It made it difficult for me to sleep when I was there.

His attitude of feeling like he was GOD that he could enter whenever he feels like it made me very uncomfortable It exacerbated my insomnia problems, not feeling safe there. I told him that, at my apartment, they have to give me 24 hours notice to enter. Adam stated, this is not an apartment, I can enter whether you have a Do not Disturb sign it or not. I couldnt understand why this rude/arrogant/I am GOD type of attitude.
I spent 8 months in Guatemala in 2007, in very cheap hotels at $5$10 p/night, and never was disrespected this way. I got better service from the General managers at the super cheap hotels in Guatemala! I got alot better view at the La Quinta Hotel on the 3rd floor. It makes me feel like I am in Oregon again. At the comfort suites hotel, I only had a view of some bushes, and the parking lot.

I tipped Adam before he had been rude, like he is GOD, and has the power to enter rooms whether you have a DO NOT DISTURB sign or not. That was because, he went out, and bought a rubber bath mat, when I told him I had almost slipped in the bathtub. He seemed to think it was a bother, but, said he would go get one, and just said hed keep it when I left for his other guests. So, I felt guilty for inconveniencing him, and even tipped him $5.00 for going to get the bathmat. I am not a rich person. My room was paid by Jon Berkeley Management Inc, a very professional company that I have never had a problem with treating me rudely as a tenant at Twin Pines Apartments in Davis. Some very major repairs were being done, that blocked the tenants access, hence the reason for my stay at Comfort Suites UC Davis.
I got great customer service from the front desk clerks, the breakfast bar attendant, and they do have great breakfasts. The housekeeping staff were all very pleasant, and I enjoyed any opportunity I had to practice my Spanish.
They also allowed me to bring in a hard of hearing telephone I had to use in my room. They accommodated an ambulance coming in to my room the previous Sat evening due to I was having some trouble breathing related to my Sleep Apnea.

This could have been a very pleasant stay. I am a senior citizen, I deserve to be treated better than this. I heard from some of the La Quinta Staff Hotel across the street, that they know he has done this to other guests as well, entering rooms without their permission.
I was exhausted, but, I stayed up to almost midnight, and by myself without a car carried five boxes, 2 suitcases, and several bags across the street to the better Hotel, La Quinta, and rented a room out of my own pocket, even though my room at Comfort Suite UC Davis was already being paid by JonBerkley, that was how angry I felt.

The front desk at La Quinta offered to move my vintage 1950s sewing machine over for me across the street as it was to heavy for me. As well as they provided a staff member to come with me while I checked out, as I wanted a witness that I had left the room clean, undamaged, and checked out before 11 AM. I recommend Comfort Suite, but, not the General Manager, Adam. Any further questions, or comments, email me at: redcrossalumni@yahoo.com
Maureen
cc: La Quinta Hotel Management —deals911

20120206 15:55:45   This comment is in response to the above comment submitted above by user deals911. Unfortunately the comments are incorrect and deliberately false.

The Comfort Suites staff and management did everything possible to accommodate this guest during her two week stay. On 02/01/2012 she informed our Front Desk staff that an outlet in her room did not work. She said she was going out and wanted it fixed she did not say to not enter the room (otherwise we would not have). Our maintenance staff found all outlets in the room to be in perfect working condition. When the guest returned I informed her that the outlet was fine. She commented that it may be broken, but also that it was new and she would return it. I suggested doublechecking the instructions first. Not once did she mention that she was unhappy at this point. It should also be noted that a few days prior she called from outside the hotel and asked me to go through her room and look for a safety deposit key she had lost. I informed her that that we only enter guest rooms for emergency or maintenance items. That request definitely contradicts her complaint.

Not once did I make any of the comments she stated. I honestly feel she is confused with an earlier conversation we had when she was inquiring about emergency procedures the day after she had called an ambulance. If she would have addressed the issue with myself or staff prior to moving I think both hotel and guest would have concluded that it was a misunderstanding. Its unfortunate and disappointing that after helping this guest for two weeks that this is the course of action she has chosen to take.

Adam Gough
General Manger
Comfort Suites —gmca642

A subsequent reply by deals911 was deleted because it contained personal attacks and potentially libelous remarks, but deals911 should feel free to leave another reply that does not include personal attacks.

20120207 21:38:51 nbsp Ive had many family members and friends stay here during the past few years. We even held a family reunion here which included several octogenarians and even one nonagenerian as well as multiple tots. And Ive had my elderly parents stay here several times. Each time, the staff has been terrific. Adam has been exceptionally helpful. So we are loyal patrons! Users/Chamoudah