text:
! The kind of code only a scientist could love.
! TODO: The cape routine needs work to remove the GOTOs

!======================================================================
!
! !IROUTINE: TVIRTUAL -- Calculate virtual temperature (K)
!
! !DESCRIPTION:
!
!   This function returns a single value of virtual temperature in
!   K, given temperature in K and mixing ratio in kg/kg.  For an
!   array of virtual temperatures, use subroutine VIRTUAL_TEMP.
!
! !INPUT:
!    RATMIX - water vapor mixing ratio (kg/kg)
!    TEMP   - temperature (K)
!
! !OUTPUT:
!    TV     - Virtual temperature (K)
!

REAL(KIND=8) FUNCTION tvirtual(temp,ratmix)
    IMPLICIT NONE
    REAL(KIND=8),INTENT(IN) :: temp,ratmix
    REAL(KIND=8),PARAMETER :: EPS = .622D0

    tvirtual = temp*(EPS + ratmix)/(EPS*(1.D0 + ratmix))
    RETURN
END FUNCTION tvirtual

REAL(KIND=8) FUNCTION tonpsadiabat(thte,prs,PSADITHTE,PSADIPRS,PSADITMK,GAMMA,&
                                    throw_exception)
    IMPLICIT NONE
    EXTERNAL throw_exception
    REAL(KIND=8),INTENT(IN) :: thte
    REAL(KIND=8),INTENT(IN) :: prs
    REAL(KIND=8),DIMENSION(150),INTENT(IN) :: PSADITHTE
    REAL(KIND=8),DIMENSION(150),INTENT(IN) :: PSADIPRS
    REAL(KIND=8),DIMENSION(150,150),INTENT(IN) :: PSADITMK
    REAL(KIND=8),INTENT(IN) :: GAMMA

    REAL(KIND=8) :: fracjt
    REAL(KIND=8) :: fracjt2
    REAL(KIND=8) :: fracip
    REAL(KIND=8) :: fracip2

    INTEGER :: ip, ipch, jt, jtch

    !   This function gives the temperature (in K) on a moist adiabat
    !   (specified by thte in K) given pressure in hPa.  It uses a
    !   lookup table, with data that was generated by the Bolton (1980)
    !   formula for theta_e.

    !     First check if pressure is less than min pressure in lookup table.
    !     If it is, assume parcel is so dry that the given theta-e value can
    !     be interpretted as theta, and get temperature from the simple dry
    !     theta formula.

    IF (prs.LE.PSADIPRS(150)) THEN
      tonpsadiabat = thte * (prs/1000.D0)**GAMMA
      RETURN
    END IF

    !   Otherwise, look for the given thte/prs point in the lookup table.

    jt = -1
    DO jtch = 1,150 - 1
      IF (thte .GE. PSADITHTE(jtch) .AND. thte .LT. PSADITHTE(jtch+1)) THEN
          jt = jtch
          EXIT
          !GO TO 213
      END IF
    END DO

!    JT = -1
!213 CONTINUE
    ip = -1
    DO ipch = 1,150 - 1
      IF (prs.LE.PSADIPRS(ipch) .AND. prs.GT.PSADIPRS(ipch+1)) THEN
          ip = ipch
          EXIT
          !GO TO 215
      END IF
    END DO

!    IP = -1
!215 CONTINUE
    IF (jt.EQ.-1 .OR. ip.EQ.-1) THEN
        ! Need an exception here
        CALL throw_exception('capecalc3d: ','Outside of lookup table bounds. prs,thte=',prs,thte)
        !STOP
    END IF
    fracjt = (thte-PSADITHTE(jt)) / (PSADITHTE(jt+1)-PSADITHTE(jt))
    fracjt2 = 1.D0 - fracjt
    fracip = (PSADIPRS(ip)-prs) / (PSADIPRS(ip)-PSADIPRS(ip+1))
    fracip2 = 1.D0 - fracip
    IF (PSADITMK(ip,jt) .GT. 1D9 .OR. PSADITMK(ip+1,jt) .GT. 1D9 .OR. &
        PSADITMK(ip,jt+1) .GT. 1D9 .OR. PSADITMK(ip+1,jt+1) .GT. 1D9) THEN
        CALL throw_exception('capecalc3d: ','Tried to access missing temperature in lookup table.',&
        'Prs and Thte probably unreasonable. prs,thte=',prs,thte)
      !STOP
    END IF
    tonpsadiabat = fracip2*fracjt2*PSADITMK(ip,jt)+fracip*fracjt2*PSADITMK(ip+1,jt)+&
            fracip2*fracjt*PSADITMK(ip,jt+1)+fracip*fracjt*PSADITMK(ip+1,jt+1)

    RETURN
END FUNCTION tonpsadiabat

!     Historically, this routine calculated the pressure at full sigma
!     levels when RIP was specifically designed for MM4/MM5 output.
!     With the new generalized RIP (Feb '02), this routine is still
!     intended to calculate a set of pressure levels that bound the
!     layers represented by the vertical grid points, although no such
!     layer boundaries are assumed to be defined.  The routine simply
!     uses the midpoint between the pressures of the vertical grid
!     points as the bounding levels.  The array only contains mkzh
!     levels, so the pressure of the top of the uppermost layer is
!     actually excluded.  The kth value of pf is the lower bounding
!     pressure for the layer represented by kth data level.  At the
!     lower bounding level of the lowest model layer, it uses the
!     surface pressure, unless the data set is pressure-level data, in
!     which case it assumes the lower bounding pressure level is as far
!     below the lowest vertical level as the upper bounding pressure
!     level is above.
SUBROUTINE dpfcalc(prs,sfp,pf,miy,mjx,mkzh,ter_follow)

    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: prs
    REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) :: sfp
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(OUT) :: pf
    INTEGER,INTENT(IN) :: ter_follow,miy,mjx,mkzh

    INTEGER :: i,j,k

    !  do j=1,mjx-1  Artifact of MM5
    DO j = 1,mjx
    !  do i=1,miy-1  staggered grid
      DO i = 1,miy
          DO k = 1,mkzh
              IF (k .EQ. mkzh) THEN
    !  terrain-following data
                  IF (ter_follow .EQ. 1) THEN
                      pf(i,j,k) = sfp(i,j)
    !  pressure-level data
                  ELSE
                      pf(i,j,k) = .5D0 * (3.D0*prs(i,j,k) - prs(i,j,k-1))
                  END IF
              ELSE
                  pf(i,j,k) = .5D0 * (prs(i,j,k+1) + prs(i,j,k))
              END IF
          END DO
      END DO
    END DO

    RETURN
END SUBROUTINE dpfcalc

!======================================================================
!
! !IROUTINE: capecalc3d -- Calculate CAPE and CIN
!
! !DESCRIPTION:
!
!   If i3dflag=1, this routine calculates CAPE and CIN (in m**2/s**2,
!   or J/kg) for every grid point in the entire 3D domain (treating
!   each grid point as a parcel).  If i3dflag=0, then it
!   calculates CAPE and CIN only for the parcel with max theta-e in
!   the column, (i.e. something akin to Colman's MCAPE).  By "parcel",
!   we mean a 500-m deep parcel, with actual temperature and moisture
!   averaged over that depth.
!
!   In the case of i3dflag=0,
!   CAPE and CIN are 2D fields that are placed in the k=mkzh slabs of
!   the cape and cin arrays.  Also, if i3dflag=0, LCL and LFC heights
!   are put in the k=mkzh-1 and k=mkzh-2 slabs of the cin array.
!

! Important!  The z-indexes must be arranged so that mkzh (max z-index) is the
! surface pressure.  So, pressure must be ordered in ascending order before
! calling this routine.  Other variables must be ordered the same (p,tk,q,z).

! Also, be advised that missing data values are not checked during the computation.
! Also also, Pressure must be hPa
!
SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
            PSADITHTE,PSADIPRS,PSADITMK,cmsg,i3dflag,ter_follow,&
            throw_exception,miy,mjx,mkzh)

    IMPLICIT NONE
    EXTERNAL throw_exception
    INTEGER,INTENT(IN) :: miy,mjx,mkzh,i3dflag,ter_follow
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: prs
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: tmk
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: qvp
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: ght
    REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) :: ter
    REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) ::sfp
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(INOUT) :: cape
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(INOUT) :: cin
    REAL(KIND=8),DIMENSION(150),INTENT(IN) :: PSADITHTE,PSADIPRS
    REAL(KIND=8),DIMENSION(150,150),INTENT(IN) :: PSADITMK
    REAL(KIND=8),INTENT(IN) :: cmsg

    ! local variables
    INTEGER :: i,j,k,ilcl,kel,kk,klcl,klev,klfc,kmax,kpar,kpar1,kpar2
    REAL(KIND=8) :: davg,ethmax,q,t,p,e,eth,tlcl,zlcl
    REAL(KIND=8) :: pavg,tvirtual,p1,p2,pp1,pp2,th,totthe,totqvp,totprs
    REAL(KIND=8) :: cpm,deltap,ethpari,gammam,ghtpari,qvppari,prspari,tmkpari
    REAL(KIND=8) :: facden,fac1,fac2,qvplift,tmklift,tvenv,tvlift,ghtlift
    REAL(KIND=8) :: eslift,tmkenv,qvpenv,tonpsadiabat
    REAL(KIND=8) :: benamin,dz,pup,pdn
    REAL(KIND=8),DIMENSION(150) :: buoy,zrel,benaccum
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh) :: prsf

    ! constants
    INTEGER,PARAMETER :: IUP = 6
    REAL(KIND=8),PARAMETER :: CELKEL = 273.15d0
    REAL(KIND=8),PARAMETER :: GRAV = 9.81d0
    ! hpa
    REAL(KIND=8),PARAMETER :: EZERO = 6.112d0
    REAL(KIND=8),PARAMETER :: ESLCON1 = 17.67d0
    REAL(KIND=8),PARAMETER :: ESLCON2 = 29.65d0
    REAL(KIND=8),PARAMETER :: EPS = 0.622d0
    ! j/k/kg
    REAL(KIND=8),PARAMETER :: RGAS = 287.04d0
    !  j/k/kg  note: not using bolton's value of 1005.7
    REAL(KIND=8),PARAMETER :: CP = 1004.d0
    REAL(KIND=8),PARAMETER :: GAMMA = RGAS/CP
    !  cp_moist=cp*(1.+cpmd*qvp)
    REAL(KIND=8),PARAMETER :: CPMD = .887d0
    !  rgas_moist=rgas*(1.+rgasmd*qvp)
    REAL(KIND=8),PARAMETER :: RGASMD = .608d0
    !  gamma_moist=gamma*(1.+gammamd*qvp)
    REAL(KIND=8),PARAMETER :: GAMMAMD = RGASMD - CPMD
    REAL(KIND=8),PARAMETER :: TLCLC1 = 2840.d0
    REAL(KIND=8),PARAMETER :: TLCLC2 = 3.5d0
    REAL(KIND=8),PARAMETER :: TLCLC3 = 4.805d0
    REAL(KIND=8),PARAMETER :: TLCLC4 = 55.d0
    !  k
    REAL(KIND=8),PARAMETER :: THTECON1 = 3376.d0
    REAL(KIND=8),PARAMETER :: THTECON2 = 2.54d0
    REAL(KIND=8),PARAMETER :: THTECON3 = .81d0

    ! To get rid of compiler warnings
    tmkpari = 0
    qvppari = 0
    klev = 0
    klcl = 0

    ! the comments were taken from a mark stoelinga email, 23 apr 2007,
    ! in response to a user getting the "outside of lookup table bounds"
    ! error message.

    ! tmkpari  - initial temperature of parcel, k
    !    values of 300 okay. (not sure how much from this you can stray.)

    ! prspari - initial pressure of parcel, hpa
    !    values of 980 okay. (not sure how much from this you can stray.)

    ! thtecon1, thtecon2, thtecon3
    !     these are all constants, the first in k and the other two have
    !     no units.  values of 3376, 2.54, and 0.81 were stated as being
    !     okay.

    ! tlcl - the temperature at the parcel's lifted condensation level, k
    !        should be a reasonable atmospheric temperature around 250-300 k
    !        (398 is "way too high")

    ! qvppari - the initial water vapor mixing ratio of the parcel,
    !           kg/kg (should range from 0.000 to 0.025)
    !

    !  calculated the pressure at full sigma levels (a set of pressure
    !  levels that bound the layers represented by the vertical grid points)

    CALL dpfcalc(prs,sfp,prsf,miy,mjx,mkzh,ter_follow)

    !  before looping, set lookup table for getting temperature on
    !  a pseudoadiabat.

    !call dlookup_table(psadithte,psadiprs,psaditmk,psafile)

    !   do j=1,mjx-1
    DO j = 1,mjx
    !   do i=1,miy-1
      DO i = 1,miy
          cape(i,j,1) = 0.d0
          cin(i,j,1) = 0.d0

          IF (i3dflag .EQ. 1) THEN
              kpar1 = 2
              kpar2 = mkzh
          ELSE
              ! find parcel with max theta-e in lowest 3 km agl.
              ethmax = -1.d0
              DO k = mkzh,1,-1
                  IF (ght(i,j,k)-ter(i,j) .LT. 3000.d0) THEN
                      q = MAX(qvp(i,j,k), 1.d-15)
                      t = tmk(i,j,k)
                      p = prs(i,j,k)
                      e = q*p/(EPS + q)
                      tlcl = TLCLC1/ (LOG(t**TLCLC2/e)-TLCLC3) + TLCLC4
                      eth = t * (1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))*&
                              EXP((THTECON1/tlcl - THTECON2)*q*(1.d0 + THTECON3*q))
                      IF (eth .GT. ethmax) THEN
                          klev = k
                          ethmax = eth
                      END IF
                  END IF
              END DO
              kpar1 = klev
              kpar2 = klev

              ! Establish average properties of that parcel
              ! (over depth of approximately davg meters)

              ! davg=.1
              davg = 500.d0
              pavg = davg*prs(i,j,kpar1)*&
                  GRAV/(RGAS*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1)))
              p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh))
              p1 = p2 - pavg
              totthe = 0.d0
              totqvp = 0.d0
              totprs = 0.d0
              DO k = mkzh,2,-1
                  IF (prsf(i,j,k) .LE. p1) GOTO 35
                  IF (prsf(i,j,k-1) .GE. p2) GOTO 34
                  p = prs(i,j,k)
                  pup = prsf(i,j,k)
                  pdn = prsf(i,j,k-1)
                  q = MAX(qvp(i,j,k),1.d-15)
                  th = tmk(i,j,k)*(1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))
                  pp1 = MAX(p1,pdn)
                  pp2 = MIN(p2,pup)
                  IF (pp2 .GT. pp1) THEN
                      deltap = pp2 - pp1
                      totqvp = totqvp + q*deltap
                      totthe = totthe + th*deltap
                      totprs = totprs + deltap
                  END IF
    34        CONTINUE
              END DO
    35        CONTINUE
              qvppari = totqvp/totprs
              tmkpari = (totthe/totprs)*&
                  (prs(i,j,kpar1)/1000.d0)**(GAMMA*(1.d0+GAMMAMD*qvp(i,j,kpar1)))
          END IF

          DO kpar = kpar1,kpar2

              ! Calculate temperature and moisture properties of parcel
              ! (note, qvppari and tmkpari already calculated above for 2d case.)

              IF (i3dflag .EQ. 1) THEN
                  qvppari = qvp(i,j,kpar)
                  tmkpari = tmk(i,j,kpar)
              END IF
              prspari = prs(i,j,kpar)
              ghtpari = ght(i,j,kpar)
              gammam = GAMMA * (1.d0 + GAMMAMD*qvppari)
              cpm = CP * (1.d0 + CPMD*qvppari)

              e = MAX(1.d-20,qvppari*prspari/(EPS + qvppari))
              tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4
              ethpari = tmkpari*(1000.d0/prspari)**(GAMMA*(1.d0 + GAMMAMD*qvppari))*&
                  EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.d0 + THTECON3*qvppari))
              zlcl = ghtpari + (tmkpari - tlcl)/(GRAV/cpm)

              ! Calculate buoyancy and relative height of lifted parcel at
              ! all levels, and store in bottom up arrays.  add a level at the lcl,
              ! and at all points where buoyancy is zero.
              !
              ! For arrays that go bottom to top
              kk = 0
              ilcl = 0

              IF (ghtpari .GE. zlcl) THEN
                  ! Initial parcel already saturated or supersaturated.
                  ilcl = 2
                  klcl = 1
              END IF

              DO k = kpar,1,-1
                  ! For arrays that go bottom to top
    33            kk = kk + 1

                  ! Model level is below lcl
                  IF (ght(i,j,k).lt.zlcl) THEN
                      qvplift = qvppari
                      tmklift = tmkpari - GRAV/cpm*(ght(i,j,k) - ghtpari)
                      tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k))
                      tvlift = tvirtual(tmklift, qvplift)
                      ghtlift = ght(i,j,k)
                  ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN
                      ! This model level and previous model level straddle the lcl,
                      ! so first create a new level in the bottom-up array, at the lcl.
                      tmklift = tlcl
                      qvplift = qvppari
                      facden = ght(i,j,k) - ght(i,j,k+1)
                      fac1 = (zlcl-ght(i,j,k+1))/facden
                      fac2 = (ght(i,j,k)-zlcl)/facden
                      tmkenv = tmk(i,j,k+1)*fac2 + tmk(i,j,k)*fac1
                      qvpenv = qvp(i,j,k+1)*fac2 + qvp(i,j,k)*fac1
                      tvenv = tvirtual(tmkenv, qvpenv)
                      tvlift = tvirtual(tmklift, qvplift)
                      ghtlift = zlcl
                      ilcl = 1
                  ELSE
                      tmklift = tonpsadiabat(ethpari, prs(i,j,k), PSADITHTE, PSADIPRS,&
                                             PSADITMK, GAMMA, throw_exception)
                      eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2))
                      qvplift = EPS*eslift/(prs(i,j,k) - eslift)
                      tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k))
                      tvlift = tvirtual(tmklift, qvplift)
                      ghtlift = ght(i,j,k)
                  END IF
                  !  Buoyancy
                  buoy(kk) = GRAV*(tvlift - tvenv)/tvenv
                  zrel(kk) = ghtlift - ghtpari

                  IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0d0)) THEN
                      ! Parcel ascent curve crosses sounding curve, so create a new level
                      ! in the bottom-up array at the crossing.
                      kk = kk + 1
                      buoy(kk) = buoy(kk-1)
                      zrel(kk) = zrel(kk-1)
                      buoy(kk-1) = 0.d0
                      zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/&
                          (buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2))
                  END IF

                  IF (ilcl .EQ. 1) THEN
                      klcl = kk
                      ilcl = 2
                      GOTO 33
                  END IF
              END DO

              kmax = kk
              IF (kmax .GT. 150) THEN
                  ! Need an exception here
                  CALL throw_exception('capecalc3d: kmax got too big. kmax=',kmax)
                  !STOP
              END IF

              ! If no lcl was found, set klcl to kmax.  it is probably not really
              ! at kmax, but this will make the rest of the routine behave
              ! properly.
              IF (ilcl .EQ. 0) klcl=kmax

              ! Get the accumulated buoyant energy from the parcel's starting
              ! point, at all levels up to the top level.
              benaccum(1) = 0.0d0
              benamin = 9d9
              DO k = 2,kmax
                  dz = zrel(k) - zrel(k-1)
                  benaccum(k) = benaccum(k-1) + .5d0*dz*(buoy(k-1) + buoy(k))
                  IF (benaccum(k) .LT. benamin) THEN
                      benamin = benaccum(k)
                  END IF
              END DO

              ! Determine equilibrium level (el), which we define as the highest
              ! level of non-negative buoyancy above the lcl. note, this may be
              ! the top level if the parcel is still buoyant there.

              DO k = kmax,klcl,-1
                  IF (buoy(k) .GE. 0.d0) THEN
                      ! k of equilibrium level
                      kel = k
                      GOTO 50
                  END IF
              END DO

              ! If we got through that loop, then there is no non-negative
              ! buoyancy above the lcl in the sounding.  in these situations,
              ! both cape and cin will be set to -0.1 j/kg. (see below about
              ! missing values in v6.1.0). also, where cape is
              ! non-zero, cape and cin will be set to a minimum of +0.1 j/kg, so
              ! that the zero contour in either the cin or cape fields will
              ! circumscribe regions of non-zero cape.

              ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return
              ! value of this function. at that time we decided to change -0.1
              ! to a more appropriate missing value, which is passed into this
              ! routine as cmsg.

              ! cape(i,j,kpar) = -0.1d0
              ! cin(i,j,kpar) = -0.1d0
              cape(i,j,kpar) = cmsg
              cin(i,j,kpar)  = cmsg
              klfc = kmax

              GOTO 102

    50        CONTINUE

              !   If there is an equilibrium level, then cape is positive.  we'll
              !   define the level of free convection (lfc) as the point below the
              !   el, but at or above the lcl, where accumulated buoyant energy is a
              !   minimum.  the net positive area (accumulated buoyant energy) from
              !   the lfc up to the el will be defined as the cape, and the net
              !   negative area (negative of accumulated buoyant energy) from the
              !   parcel starting point to the lfc will be defined as the convective
              !   inhibition (cin).

              !   First get the lfc according to the above definition.
              benamin = 9d9
              klfc = kmax
              DO k = klcl,kel
                  IF (benaccum(k) .LT. benamin) THEN
                      benamin = benaccum(k)
                      klfc = k
                  END IF
              END DO

              ! Now we can assign values to cape and cin

              cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1d0)
              cin(i,j,kpar) = MAX(-benamin, 0.1d0)

              ! cin is uninteresting when cape is small (< 100 j/kg), so set
              ! cin to -0.1 (see note about missing values in v6.1.0) in
              ! that case.

              ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return
              ! value of this function. at that time we decided to change -0.1
              ! to a more appropriate missing value, which is passed into this
              ! routine as cmsg.

              ! IF (cape(i,j,kpar).lt.100.d0) cin(i,j,kpar) = -0.1d0
              IF (cape(i,j,kpar) .LT. 100.d0) cin(i,j,kpar) = cmsg
    102   CONTINUE

          END DO

          IF (i3dflag .EQ. 0) THEN
              cape(i,j,mkzh) = cape(i,j,kpar1)
              cin(i,j,mkzh) = cin(i,j,kpar1)
    !  meters agl
              cin(i,j,mkzh-1) = zrel(klcl) + ghtpari - ter(i,j)
    !  meters agl
              cin(i,j,mkzh-2) = zrel(klfc) + ghtpari - ter(i,j)

          ENDIF
      END DO
    END DO

    RETURN
END SUBROUTINE f_computecape

max_stars_repo_name: skysightsoaringweather/wrf-python

max_issues_repo_name: skysightsoaringweather/wrf-python

max_forks_repo_name: skysightsoaringweather/wrf-python


program example2

use nc4fortran, only : netcdf_file
implicit none (type, external)

character(:), allocatable :: filename
integer :: i32

type(netcdf_file) :: h

filename = 'nc4fortran_example2.nc'

call h%open(filename, action='w')
call h%write('x', 123)
call h%close()

call h%open(filename, action='r')
call h%read('x', i32)
if (i32 /= 123) error stop 'incorrect value read'

print *, 'OK: example 2'

end program


text:
../model/fir_2tapfilter_real.v
../model/fir_2tapfilter.v
../model/tx_driver_real.v
../model/alexander_pd.v
../model/dfe.v
../model/digital_lf.v
../model/dco.v
../model/bb_cdr.v
../../channel/model/channel.v
../../channel/model/channel_aug.v
../../ctle/model/ctle.v
../../phase_locked_loop/model/chgpmp.v
../../phase_locked_loop/model/div2.v
../../phase_locked_loop/model/lpf.v
../../phase_locked_loop/model/pfd.v
../../phase_locked_loop/model/pll2nd.v
../../phase_locked_loop/model/regulator.v
../../phase_locked_loop/model/ringosc.v
../../phase_locked_loop/model/vco.v
../../phase_locked_loop/model/txf_vco.v
../../phase_locked_loop/model/vcobuf.v

ext: f


text:
! This module is a dummy version of the module in "abaqus_utils_mod.f90".
! Having this allows the use of "use abaqus_utils_mod" (note that the ", only : " usage cannot be
! applied). The only change when using Abaqus or external compilation is that the present file
! should be included for the former, and "abaqus_utils_mod.f90" should be included for the latter.
module abaqus_utils_mod
implicit none

    contains

    ! Empty module

end module abaqus_utils_mod

max_stars_repo_name: KnutAM/AbaqusRolloverSimulation

max_issues_repo_name: KnutAM/AbaqusRolloverSimulation

max_forks_repo_name: KnutAM/AbaqusRolloverSimulation


!
!     --------------------------------------------------------------
!        T R B A K 1
!     --------------------------------------------------------------
!
!
      SUBROUTINE TRBAK1(NM, N, A, E, M, Z)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
!...Translated by Pacific-Sierra Research 77to90  4.3E  21:52:17  11/14/01
!...Switches:
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: NM
      INTEGER , INTENT(IN) :: N
      INTEGER , INTENT(IN) :: M
      REAL(DOUBLE) , INTENT(IN) :: A(NM,N)
      REAL(DOUBLE) , INTENT(IN) :: E(N)
      REAL(DOUBLE) , INTENT(INOUT) :: Z(NM,M)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: I, J, K, L
      REAL(DOUBLE) :: S
!-----------------------------------------------
!
!
!     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK1,
!     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
!     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
!
!     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
!     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
!     SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  TRED1.
!
!     ON INPUT:
!
!        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
!          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
!          DIMENSION STATEMENT;
!
!        N IS THE ORDER OF THE MATRIX;
!
!        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
!          FORMATIONS USED IN THE REDUCTION BY  TRED1
!          IN ITS STRICT LOWER TRIANGLE;
!
!        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
!          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY;
!
!        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED;
!
!        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
!          IN ITS FIRST M COLUMNS.
!
!     ON OUTPUT:
!
!        Z CONTAINS THE TRANSFORMED EIGENVECTORS
!          IN ITS FIRST M COLUMNS.
!
!     NOTE THAT TRBAK1 PRESERVES VECTOR EUCLIDEAN NORMS.
!
!     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
!     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
!
!     ------------------------------------------------------------------
!
      IF (M /= 0) THEN
         IF (N /= 1) THEN
!
            DO I = 2, N
               L = I - 1
               IF (E(I) == 0.0D0) CYCLE
!
               DO J = 1, M
!
                  S = SUM(A(I,:L)*Z(:L,J))
!     :::::::::: DIVISOR BELOW IS NEGATIVE OF H FORMED IN TRED1.
!                DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ::::::::::
                  S = (S/A(I,L))/E(I)
!
                  Z(:L,J) = Z(:L,J) + S*A(I,:L)
!
               END DO
!
            END DO
!
         ENDIF
      ENDIF
      RETURN
!     :::::::::: LAST CARD OF TRBAK1 ::::::::::
      END SUBROUTINE TRBAK1


text:
! RUN: %python %S/test_errors.py %s %flang_fc1
! Check for semantic errors in NULLIFY statements

INTEGER, PARAMETER :: maxvalue=1024

Type dt
  Integer :: l = 3
End Type
Type t
  Type(dt) :: p
End Type

Type(t),Allocatable :: x(:)

Integer :: pi
Procedure(Real) :: prp

Allocate(x(3))
!ERROR: component in NULLIFY statement must have the POINTER attribute
Nullify(x(2)%p)

!ERROR: name in NULLIFY statement must have the POINTER attribute
Nullify(pi)

!ERROR: name in NULLIFY statement must have the POINTER attribute
Nullify(prp)

!ERROR: name in NULLIFY statement must be a variable or procedure pointer name
Nullify(maxvalue)

End Program

! Make sure that the compiler doesn't crash when NULLIFY is used in a context
! that has reported errors
module badNullify
  interface
    function ptrFun()
      integer, pointer :: ptrFun
    end function
  end interface
contains
  !ERROR: 'ptrfun' was not declared a separate module procedure
  module function ptrFun()
    integer, pointer :: ptrFun
    real :: realVar
    nullify(ptrFun)
    nullify(realVar)
  end function
end module

max_stars_repo_name: lxbndr/llvm-project

max_issues_repo_name: lxbndr/llvm-project

max_forks_repo_name: lxbndr/llvm-project


text:
Describe Users/OhHey here.

20110809 07:33:44 nbsp Welcome to the Wiki and welcome to Davis!  With starting so early are you going to King Hall? Users/TomGarberson

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE MPPOLA(DIST,AZIM,INC,TWIST)
      COMMON /MAP/MODEL(4,4),VIEW(4,4),PROJ(4,4),CPNEAR,CPFAR,VWPORT(4),
     *       MVP(4,4),VP(4,4),CPLINE(2,2,10),CPPLAN(2,3,10),PEYE(3),
     *       PLOOK(3),ETWIST,NCPLIN,NCPLAN,TMAT1(4,4),TMAT2(4,4),
     *       TMAT3(4,4),TVEC1(4),TVEC2(4),TVEC3(4),TVEC4(4),TARR1(32),
     *       TARR2(32),TARR3(32),TARR4(32),TARR5(32),TARR6(32),
     *       TARR7(32),TARR8(32)
      REAL MODEL,MVP
      REAL INC
      PARAMETER (DPR=57.2958)

      ANG = -AZIM/DPR
      CALL LDROTA('y',COS(ANG),SIN(ANG),TMAT1)
      ANG = -INC/DPR
      CALL LDROTA('x',COS(ANG),SIN(ANG),TMAT2)
      CALL MXMULT(4,TMAT1,TMAT2,TMAT3)
      ANG = -TWIST/DPR
      CALL LDROTA('z',COS(ANG),SIN(ANG),TMAT1)
      CALL MXMULT(4,TMAT3,TMAT1,TMAT2)
      CALL LDTRAN(0.,0.,-DIST,TMAT1)
      CALL MXMULT(4,TMAT2,TMAT1,VIEW)
      CALL MXMULT(4,VIEW,PROJ,VP)
      CALL MXMULT(4,MODEL,VP,MVP)
      RINC = INC/DPR
      RAZIM = AZIM/DPR
      PEYE(1) = DIST*COS(RINC)*SIN(RAZIM)
      PEYE(2) = -DIST*SIN(RINC)
      PEYE(3) = DIST*COS(RINC)*COS(RAZIM)
      PLOOK(1) = 0.
      PLOOK(2) = 0.
      PLOOK(3) = 0.
      ETWIST = TWIST
      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
c this file contains the following subroutines:
c
c readharmfile: reads a file with the harmonic expansion of the potential and density.
c all the relevant stuff ends up in a common block that only needs to be seen by
c the routines in this file (I think!). This routine must be called first.
c
c dens(r,z): density at point (r,z)
c
c densrpsi(r,psi): density as a function of r and potential psi.
c
c pot(r,z): potential at point (r,z)
c
c force(r,z,fr,fz): radial and z  force at (r,z)
c all the rest is (slightly nobbled in places) numerical recipes routines.
c

c
c calculate the force by interpolating the derived radial functions
c
c      subroutine force(lm,s,z,fs,fz,fs1,fz1)
      subroutine force(s,z,fs,fz,pot)
      parameter (pi=3.1415926535)
      parameter (pc0=0.282094792, pc2=0.630783131, pc4=0.846284375)
c constants cn are sqrt((2*l+1)/(4*pi))
      parameter (pc6=1.017107236, pc8=1.163106623)

      dimension pc(20), p(20), dp(20)
      common /potconstants/ apot(20,0:20000), frad(20,0:20000),
     +     dr, nr, lmax, potcor
      common /legendre/ plcon(0:40)
      common /gparameters/ a,b,c,v0,q,psi00,
     +                      psiouta, rho1a, sigbulge,
     +                      grmdisk, grdisk, gzdisk, goutdisk, gdrtrunc,
     +                      potcor1
      common /moreconstants/ v02, v03, rdisk2, diskconst, bulgea
      common /diskpars/ sigr0, disksr, nrdisk

      common /flags/ idiskflag, ibulgeflag, ihaloflag

c      write(*,*) 'Force entered with s,z=',s,z
c      write(*,*) 'lmax =',lmax
c      write(*,*) 'apot starts as:'
c      do ll=1,3
c         write(*,*) (apot(ll,i),i=1,10)
c      enddo
      r=sqrt(s*s+z*z)
      ihi=int(r/dr)+1
      if (ihi.lt.1) ihi=1
      if (ihi.gt.nr) ihi=nr
      r1=dr*(ihi-1)
      r2=dr*ihi
      redge = nr*dr
      t=(r-r1)/(r2-r1)
      ihim1 = ihi - 1
      tm1 = 1.0 - t
      if (r.eq.0.) then
         fs = 0.0
         fz = 0.0
      else
         costheta=z/r
         ct2 = costheta*costheta
         sintheta=s/r

         do l=0,lmax,2
            pc(l/2+1) = sqrt((2.0*l + 1)/(4.0*pi))
            p(l/2+1) = plgndr1(l,costheta)
			if( costheta .eq. 1.0 ) then
               dp(l/2+1) = 0.0
            else
               st2 = 1.0 - costheta*costheta
               dp(l/2+1) = l*(plgndr1(l-1, costheta) - costheta*p(l/2+1))/st2
            endif
         enddo
         do i=1,lmax/2+1
            p(i) = p(i)*pc(i)
            dp(i) = dp(i)*pc(i)
         enddo

         if( r .le. redge ) then
             fr = 0.0
             do i=1,lmax/2+1
                fr = fr + p(i)*(t*frad(i,ihi) + tm1*frad(i,ihim1))
             enddo
             fth = 0.0
             do i=2,lmax/2+1
                fth = fth - sintheta*dp(i)*(t*apot(i,ihi) + tm1*apot(i,ihim1))
             enddo
             pot = 0.0
             do i=1,lmax/2+1
                pot = pot + p(i)*(t*apot(i,ihi) + tm1*apot(i,ihim1))
             enddo
         else
! expressions below for fr and pot not correct if apot shifted!!
             fr = 0.0
             do i=1,lmax/2+1
                 l = 2*(i-1)
                 fr = fr - (l+1)*p(i)*apot(i,nr)/redge*(redge/r)**(l+2)
             enddo
             fth = 0.0
             do i=2,lmax/2+1
                 l = 2*(i-1)
                 fth = fth - sintheta*dp(i)*apot(i,nr)*(redge/r)**(l+1)
	     enddo
             pot = 0.0
             do i=1,lmax/2+1
                l = 2*(i-1)
                pot = pot + p(i)*apot(i,nr)*(redge/r)**(l+1)
             enddo
         endif
         if( idiskflag .eq. 1.OR.idiskflag.EQ.3 ) then
             pot = pot + appdiskpot(s,z)
         endif

      if( idiskflag .eq. 2 .or. idiskflag .eq. 3) then
       pot= pot + thickdiskpot(s,z)
!       fr=fr+thickdiskforce(r)
! fth (or rather fs and fz directly) tbd, see below
      endif

      if(ihaloflag.eq.2) then
       pot= pot + fixedpot(r)
       fr=fr-fixedforce(r)
! fixed force negative, fr positive
! fixed halo is spherical so no fth
      endif

         fs = -(sintheta*fr + costheta/r*fth)
         fz = -(costheta*fr - sintheta/r*fth)

!         write(*,*) 'Force before disk correction ',fr,fth
         if( idiskflag .eq. 1 .OR. idiskflag.EQ.3) then
            call appdiskforce(s,z,fsad,fzad)
c            write(*,*) 'Force of disk correction ',fsad,fzad
            fs = fs + fsad
            fz = fz + fzad
         endif

         if( idiskflag .eq. 2 .OR. idiskflag.EQ.3) then
!            print*,fs,fz,thickdiskrforce(s,z)
            fs = fs + thickdiskrforce(s,z)
            fz = fz + thickdiskzforce(s,z)
         endif
      endif
      return
      end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
MODULE inij

IMPLICIT NONE

CONTAINS

   SUBROUTINE initij
     USE glodata

     I_rms = 1.0
     J_rms = I_rms/(wind_thickness*wind_width)

     J_z(1:Nx+1,1:Ny+1) = 0.0

     J_z(N_core_x+N_core_wind_x+1:N_core_x+N_core_wind_x &
         +N_wind_x+1,N_dis_wind+1:N_dis_wind+1+N_wind_y) &
        = J_rms

   END SUBROUTINE initij

END MODULE inij

max_stars_repo_name: YudiXiao/Field-cannon-technique

max_issues_repo_name: YudiXiao/Field-cannon-technique

max_forks_repo_name: YudiXiao/Field-cannon-technique


text:
! MIT License
! Copyright (c) 2019 Sam Miller
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

module mod_local_evo_operator
#ifdef USE_OPENMP
  use omp_lib
#endif /* USE_OPENMP */
  use, intrinsic :: iso_fortran_env, only: &
    ik => int32, rk => real64, std_err => error_unit, std_out => output_unit
  use mod_globals, only: debug_print, print_recon_data, track_single_cell_cone, enable_debug_cone_output, &
                         du_midpoint_unit, lr_midpoint_unit, corner_unit, single_cell_cone_unit

  use, intrinsic :: ieee_arithmetic
  use mod_floating_point_utils, only: near_zero, equal, neumaier_sum
  use math_constants, only: pi, rad2deg
  use mod_input, only: input_t
  use mod_mach_cone_collection, only: mach_cone_collection_t

  use mod_grid_block_2d, only: grid_block_2d_t, C1, M1, C2, M2, C3, M3, C4, M4
  use mod_abstract_reconstruction, only: abstract_reconstruction_t

  implicit none

  private
  public :: local_evo_operator_t

  type :: local_evo_operator_t
    !< Local Evolution Operator (E0)
    !< This class provides the framework for any evolution operators. For example, the local evolution
    !< operator will inherit from this class. Most of the attributes are the same, just that the
    !< implementation varies between those who inherit this class

    class(grid_block_2d_t), pointer :: grid => null()
    !< pointer to the grid object

    real(rk), dimension(:, :, :), pointer :: reconstructed_rho => null()
    !< (point 1:8, i, j); pointer to reconstructed density
    real(rk), dimension(:, :, :), pointer :: reconstructed_u => null()
    !< (point 1:8, i, j); pointer to reconstructed x-velocity
    real(rk), dimension(:, :, :), pointer :: reconstructed_v => null()
    !< (point 1:8, i, j); pointer to reconstructed y-velocity
    real(rk), dimension(:, :, :), pointer :: reconstructed_p => null()
    !< (point 1:8, i, j); pointer to reconstructed pressure

    class(abstract_reconstruction_t), pointer :: reconstruction_operator => null()
    !< pointer to the R_Omega operator used to provide values at the P' location

    character(:), allocatable :: name  !< Name of the evolution operator
    real(rk) :: tau !< time increment to evolve
    integer(ik) :: iteration = 0
    real(rk) :: time = 0.0_rk
    real(rk) :: time_step = 0.0_rk
  contains
    procedure, public :: initialize
    procedure, public :: evolve
    procedure, public, nopass :: e0_operator
    final :: finalize
  endtype local_evo_operator_t

contains

  subroutine initialize(self, dt, grid_target, recon_operator_target)
    !< Constructor for the FVLEG operator
    class(local_evo_operator_t), intent(inout) :: self
    class(grid_block_2d_t), intent(in), target :: grid_target
    real(rk), intent(in) :: dt !< timestep
    class(abstract_reconstruction_t), intent(in), target :: recon_operator_target

    self%name = "FVLEG"
    self%time_step = dt
    self%grid => grid_target
    self%reconstruction_operator => recon_operator_target

  endsubroutine initialize

  subroutine finalize(self)
    !< Cleanup the operator type
    type(local_evo_operator_t), intent(inout) :: self

    call debug_print('Running local_evo_operator_t%finalize()', __FILE__, __LINE__)

    if(associated(self%grid)) nullify(self%grid)

    if(associated(self%reconstructed_rho)) nullify(self%reconstructed_rho)
    if(associated(self%reconstructed_u)) nullify(self%reconstructed_u)
    if(associated(self%reconstructed_v)) nullify(self%reconstructed_v)
    if(associated(self%reconstructed_p)) nullify(self%reconstructed_p)
    if(associated(self%reconstruction_operator)) nullify(self%reconstruction_operator)

    ! if (allocated(self%leftright_midpoint_neighbors)) deallocate(self%leftright_midpoint_neighbors)
    ! if (allocated(self%downup_midpoint_neighbors)) deallocate(self%downup_midpoint_neighbors)
    ! if (allocated(self%corner_neighbors)) deallocate(self%corner_neighbors)

    if(allocated(self%name)) deallocate(self%name)

  endsubroutine finalize

  subroutine evolve(self, dt, location, evolved_rho, evolved_u, evolved_v, evolved_p, lbounds, error_code)
    !< Create Mach cones and evolve the state at all of the left/right midpoints in the domain. Left/right midpoints
    !< are the midpoints defined by vectors pointing to the left and right.

    class(local_evo_operator_t), intent(inout) :: self
    integer(ik), dimension(2), intent(in) :: lbounds
    character(len=*), intent(in) :: location !< Mach cone location ['corner', 'left/right midpoint', or 'down/up midpoint']
    integer(ik), intent(out) :: error_code
    real(rk), intent(in) :: dt !< time step
    real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(out) :: evolved_rho
    real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(out) :: evolved_u
    real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(out) :: evolved_v
    real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(out) :: evolved_p

    ! Locals
    integer(ik) :: i, j
    integer(ik) :: ilo, ihi, jlo, jhi

    integer(ik), dimension(:, :, :, :), allocatable :: leftright_midpoint_neighbors
    !< ((i,j), (cell 1:2), i, j); (i,j) indices of the cells on either side of the left/right vector midpoint

    integer(ik), dimension(:, :, :, :), allocatable :: downup_midpoint_neighbors
    !< ((i,j), (cell 1:2), i, j); (i,j) indices of the cells on either side of the down/up vector midpoint

    integer(ik), dimension(:, :, :, :), allocatable :: corner_neighbors
    !< ((i,j), (cell 1:4), i, j); (i,j) indices of the cells surrounding the corner vector set

    type(mach_cone_collection_t) :: leftright_midpoint_mach_cones
    type(mach_cone_collection_t) :: downup_midpoint_mach_cones
    type(mach_cone_collection_t) :: corner_mach_cones

    real(rk) :: tau !< time increment

    real(rk), dimension(:, :, :), allocatable :: reconstructed_rho
    !< (point 1:8, i, j); the reconstructed density of the corners/midpoints with respect to each cell

    real(rk), dimension(:, :, :), allocatable :: reconstructed_u
    !< (point 1:8, i, j); the reconstructed x-velocity of the corners/midpoints with respect to each cell

    real(rk), dimension(:, :, :), allocatable :: reconstructed_v
    !< (point 1:8, i, j); the reconstructed y-velocity of the corners/midpoints with respect to each cell

    real(rk), dimension(:, :, :), allocatable :: reconstructed_p
    !< (point 1:8, i, j); the reconstructed pressure of the corners/midpoints with respect to each cell

    error_code = 0
    self%time_step = dt
    tau = dt / 2.0_rk
    ilo = lbound(evolved_rho, dim=1)
    ihi = ubound(evolved_rho, dim=1)
    jlo = lbound(evolved_rho, dim=2)
    jhi = ubound(evolved_rho, dim=2)

    select case(trim(location))
    case('corner')
      call debug_print('In local_evo_operator_t%evolve() -> location=corner', __FILE__, __LINE__)
      !       cell 4                   cell 3
      !      (i-1,j)                   (i,j)
      !  C4----M3----C3    P3   C4----M3----C3
      !  |            |    |    |            |
      !  M4    X4    M2    |    M4    X3    M2
      !  |            |    |    |            |
      !  C1----M1----C2    |    C1----M1----C2
      !                    |
      !  P4----------------O-----------------P2
      !                    |
      !  C4----M3----C3    |    C4----M3----C3
      !  |            |    |    |            |
      !  M4    X1    M2    |    M4    X2    M2
      !  |            |    |    |            |
      !  C1----M1----C2    P1   C1----M1----C2
      !      cell 1                  cell 2
      !     (i-1,j-1)               (i,j-1)

      ! For corners, the edge vectors go down (0-P1), right(O-P2), up(O-P3), right (O-P4)
      ! O is the corner point shared by the 4 cells
      if(.not. allocated(corner_neighbors)) then
        allocate(corner_neighbors(2, 4, ilo:ihi, jlo:jhi))
        do j = jlo, jhi
          do i = ilo, ihi
            corner_neighbors(:, 1, i, j) = [i - 1, j - 1] ! lower left
            corner_neighbors(:, 2, i, j) = [i, j - 1]     ! lower right
            corner_neighbors(:, 3, i, j) = [i, j]         ! upper right
            corner_neighbors(:, 4, i, j) = [i - 1, j]     ! upper left
          enddo
        enddo
      endif

      allocate(reconstructed_rho(4, ilo:ihi, jlo:jhi))
      allocate(reconstructed_u(4, ilo:ihi, jlo:jhi))
      allocate(reconstructed_v(4, ilo:ihi, jlo:jhi))
      allocate(reconstructed_p(4, ilo:ihi, jlo:jhi))

      !$omp parallel default(none) &
      !$omp shared(self) &
      !$omp shared(reconstructed_rho, reconstructed_u, reconstructed_v, reconstructed_p) &
      !$omp firstprivate(ilo,ihi,jlo,jhi) &
      !$omp private(i,j)
      !$omp do
      do j = jlo, jhi
        do i = ilo, ihi
          ! reconstructed_state indexing; ((rho, u ,v, p), point, node/midpoint, i, j)

          ! Cell 1: lower left cell -> corner is in the upper right (C3) of its parent cell
          reconstructed_rho(1, i, j) = self%reconstructed_rho(C3, i - 1, j - 1)
          reconstructed_u(1, i, j) = self%reconstructed_u(C3, i - 1, j - 1)
          reconstructed_v(1, i, j) = self%reconstructed_v(C3, i - 1, j - 1)
          reconstructed_p(1, i, j) = self%reconstructed_p(C3, i - 1, j - 1)

          ! Cell 2: lower right cell -> corner is in the upper left (C4) of its parent cell
          reconstructed_rho(2, i, j) = self%reconstructed_rho(C4, i, j - 1)
          reconstructed_u(2, i, j) = self%reconstructed_u(C4, i, j - 1)
          reconstructed_v(2, i, j) = self%reconstructed_v(C4, i, j - 1)
          reconstructed_p(2, i, j) = self%reconstructed_p(C4, i, j - 1)

          ! Cell 3: upper right cell-> corner is in the lower left (C1) of its parent cell
          reconstructed_rho(3, i, j) = self%reconstructed_rho(C1, i, j)
          reconstructed_u(3, i, j) = self%reconstructed_u(C1, i, j)
          reconstructed_v(3, i, j) = self%reconstructed_v(C1, i, j)
          reconstructed_p(3, i, j) = self%reconstructed_p(C1, i, j)

          ! Cell 4: upper left cell -> corner is in the lower right (C2) of its parent cell
          reconstructed_rho(4, i, j) = self%reconstructed_rho(C2, i - 1, j)
          reconstructed_u(4, i, j) = self%reconstructed_u(C2, i - 1, j)
          reconstructed_v(4, i, j) = self%reconstructed_v(C2, i - 1, j)
          reconstructed_p(4, i, j) = self%reconstructed_p(C2, i - 1, j)

        enddo
      enddo
      !$omp end do
      !$omp end parallel

      call corner_mach_cones%initialize(tau=tau, edge_vectors=self%grid%corner_edge_vectors, &
                                        reconstructed_rho=reconstructed_rho, &
                                        reconstructed_u=reconstructed_u, &
                                        reconstructed_v=reconstructed_v, &
                                        reconstructed_p=reconstructed_p, &
                                        cell_indices=corner_neighbors, &
                                        cone_location=trim(location), &
                                        lbounds=lbound(reconstructed_rho), &
                                        reconstruction_operator=self%reconstruction_operator)
      call self%e0_operator(cones=corner_mach_cones, &
                            rho=evolved_rho, u=evolved_u, v=evolved_v, p=evolved_p)

    case('left/right midpoint')
      call debug_print('In local_evo_operator_t%evolve() -> location=left/right midpoint', __FILE__, __LINE__)
      !       cell 1
      !      (i-1,j)
      !  C4----M3----C3
      !  |            |
      !  M4    X1     M2
      !  |            |
      !  C1----M1----C2
      !
      !  P1----O----P2  edge vector
      !
      !  C4----M3----C3
      !  |           |
      !  M4    X2    M2
      !  |           |
      !  C1----M1----C2
      !       cell 2
      !      (i,j-1)
      !
      ! For left/right midpoints, the edge vectors go left (O-P1) then right (O-P2).
      ! The neighboring cells are up (i, j+1) and down (i, j-1)

      if(.not. allocated(leftright_midpoint_neighbors)) then
        allocate(leftright_midpoint_neighbors(2, 2, ilo:ihi, jlo:jhi))
        do j = jlo, jhi
          do i = ilo, ihi
            leftright_midpoint_neighbors(:, 1, i, j) = [i, j]     ! above
            leftright_midpoint_neighbors(:, 2, i, j) = [i, j - 1] ! below
          enddo
        enddo
      endif

      allocate(reconstructed_rho(2, ilo:ihi, jlo:jhi))
      allocate(reconstructed_u(2, ilo:ihi, jlo:jhi))
      allocate(reconstructed_v(2, ilo:ihi, jlo:jhi))
      allocate(reconstructed_p(2, ilo:ihi, jlo:jhi))

      !$omp parallel default(none) &
      !$omp shared(self) &
      !$omp shared(reconstructed_rho, reconstructed_u, reconstructed_v, reconstructed_p) &
      !$omp firstprivate(ilo,ihi,jlo,jhi) &
      !$omp private(i,j)
      !$omp do
      do j = jlo, jhi
        do i = ilo, ihi
          ! Cell X1 -> use M1 from the cell above
          reconstructed_rho(1, i, j) = self%reconstructed_rho(M1, i, j)
          reconstructed_u(1, i, j) = self%reconstructed_u(M1, i, j)
          reconstructed_v(1, i, j) = self%reconstructed_v(M1, i, j)
          reconstructed_p(1, i, j) = self%reconstructed_p(M1, i, j)

          ! Cell X2 -> use M3 from the cell below
          reconstructed_rho(2, i, j) = self%reconstructed_rho(M3, i, j - 1)
          reconstructed_u(2, i, j) = self%reconstructed_u(M3, i, j - 1)
          reconstructed_v(2, i, j) = self%reconstructed_v(M3, i, j - 1)
          reconstructed_p(2, i, j) = self%reconstructed_p(M3, i, j - 1)
        enddo
      enddo
      !$omp end do
      !$omp end parallel

      call leftright_midpoint_mach_cones%initialize(tau=tau, edge_vectors=self%grid%leftright_midpoint_edge_vectors, &
                                                    reconstructed_rho=reconstructed_rho, &
                                                    reconstructed_u=reconstructed_u, &
                                                    reconstructed_v=reconstructed_v, &
                                                    reconstructed_p=reconstructed_p, &
                                                    cell_indices=leftright_midpoint_neighbors, &
                                                    cone_location=trim(location), &
                                                    lbounds=lbound(reconstructed_rho), &
                                                    reconstruction_operator=self%reconstruction_operator)
      call self%e0_operator(cones=leftright_midpoint_mach_cones, &
                            rho=evolved_rho, u=evolved_u, v=evolved_v, p=evolved_p)

    case('down/up midpoint')
      call debug_print('In local_evo_operator_t%evolve() -> location=down/up midpoint', __FILE__, __LINE__)
      !       cell 1          edge          cell 2
      !      (i-1,j)         vector          (i,j)
      !  C4----M3----C3       P2       C4----M3----C3
      !  |           |        |       |            |
      !  M4    X1    M2       O       M4    X2     M2
      !  |           |        |       |            |
      !  C1----M1----C2      P1       C1----M1----C2
      !
      ! For down/up midpoints, the edge vectors go down (O-P1) then up (O-P2).
      ! The neighboring cells are left (i, j) and right (i-1, j)
      if(.not. allocated(downup_midpoint_neighbors)) then
        allocate(downup_midpoint_neighbors(2, 2, ilo:ihi, jlo:jhi))
        do j = jlo, jhi
          do i = ilo, ihi
            downup_midpoint_neighbors(:, 1, i, j) = [i - 1, j] ! left
            downup_midpoint_neighbors(:, 2, i, j) = [i, j]     ! right
          enddo
        enddo
      endif

      allocate(reconstructed_rho(2, ilo:ihi, jlo:jhi))
      allocate(reconstructed_u(2, ilo:ihi, jlo:jhi))
      allocate(reconstructed_v(2, ilo:ihi, jlo:jhi))
      allocate(reconstructed_p(2, ilo:ihi, jlo:jhi))

      !$omp parallel default(none) &
      !$omp shared(self) &
      !$omp shared(reconstructed_rho, reconstructed_u, reconstructed_v, reconstructed_p) &
      !$omp firstprivate(ilo,ihi,jlo,jhi) &
      !$omp private(i,j)
      !$omp do
      do j = jlo, jhi
        do i = ilo, ihi
          ! Cell X1: cell to the left -> use M2
          reconstructed_rho(1, i, j) = self%reconstructed_rho(M2, i - 1, j)
          reconstructed_u(1, i, j) = self%reconstructed_u(M2, i - 1, j)
          reconstructed_v(1, i, j) = self%reconstructed_v(M2, i - 1, j)
          reconstructed_p(1, i, j) = self%reconstructed_p(M2, i - 1, j)

          ! Cell X2: cell to the right -> use M4
          reconstructed_rho(2, i, j) = self%reconstructed_rho(M4, i, j)
          reconstructed_u(2, i, j) = self%reconstructed_u(M4, i, j)
          reconstructed_v(2, i, j) = self%reconstructed_v(M4, i, j)
          reconstructed_p(2, i, j) = self%reconstructed_p(M4, i, j)
        enddo
      enddo
      !$omp end do
      !$omp end parallel

      call downup_midpoint_mach_cones%initialize(tau=tau, edge_vectors=self%grid%downup_midpoint_edge_vectors, &
                                                 reconstructed_rho=reconstructed_rho, &
                                                 reconstructed_u=reconstructed_u, &
                                                 reconstructed_v=reconstructed_v, &
                                                 reconstructed_p=reconstructed_p, &
                                                 cell_indices=downup_midpoint_neighbors, &
                                                 cone_location=trim(location), &
                                                 lbounds=lbound(reconstructed_rho), &
                                                 reconstruction_operator=self%reconstruction_operator)
      call self%e0_operator(cones=downup_midpoint_mach_cones, &
                            rho=evolved_rho, u=evolved_u, v=evolved_v, p=evolved_p)

    case default
      error stop "Unsupported location in local_evo_operator_t%evolve()"
    endselect

    if(allocated(reconstructed_rho)) deallocate(reconstructed_rho)
    if(allocated(reconstructed_u)) deallocate(reconstructed_u)
    if(allocated(reconstructed_v)) deallocate(reconstructed_v)
    if(allocated(reconstructed_p)) deallocate(reconstructed_p)
  endsubroutine evolve

  subroutine e0_operator(cones, rho, u, v, p)
    type(mach_cone_collection_t), intent(in) :: cones
    real(rk), dimension(:, :), intent(inout) :: rho
    real(rk), dimension(:, :), intent(inout) :: u
    real(rk), dimension(:, :), intent(inout) :: v
    real(rk), dimension(:, :), intent(inout) :: p

    real(rk), allocatable, dimension(:, :) :: rho_a_tilde
    integer(ik) :: i, j, c
    real(rk) :: t1, t2, t3
    real(rk), dimension(3) :: t_tot

    call debug_print('Running local_evo_operator_t%e0_operator()', __FILE__, __LINE__)

    allocate(rho_a_tilde, mold=cones%reference_density)

    associate(a_tilde => cones%reference_sound_speed, &
              dtheta => cones%dtheta, &
              sin_dtheta => cones%sin_dtheta, &
              cos_dtheta => cones%cos_dtheta, &
              sin_d2theta => cones%sin_d2theta, &
              cos_d2theta => cones%cos_d2theta, &
              rho_p_prime => cones%density_p_prime, &
              pressure_p_prime => cones%pressure_p_prime, &
              u_i => cones%u, &
              v_i => cones%v, &
              p_i => cones%p)

      !$omp parallel default(none) &
      !$omp shared(cones, p, rho, u, v, rho_a_tilde) &
      !$omp private(i,j) &
      !$omp private(t1, t2, t3, t_tot)
      !$omp do
      do j = 1, cones%nj
        do i = 1, cones%ni
          rho_a_tilde(i, j) = cones%reference_density(i, j) * cones%reference_sound_speed(i, j)
        enddo
      enddo
      !$omp end do

      !$omp do
      do j = 1, cones%nj
        do i = 1, cones%ni
          t1 = neumaier_sum(p_i(:, i, j) * dtheta(:, i, j))
          t2 = neumaier_sum(-rho_a_tilde(i, j) * u_i(:, i, j) * sin_dtheta(:, i, j))
          t3 = neumaier_sum(rho_a_tilde(i, j) * v_i(:, i, j) * cos_dtheta(:, i, j))
          t_tot = [t1, t2, t3]
          p(i, j) = neumaier_sum(t_tot) / (2.0_rk * pi)
        enddo
      enddo
      !$omp end do

      !$omp do
      do j = 1, cones%nj
        do i = 1, cones%ni
          t_tot = [rho_p_prime(i, j),(p(i, j) / a_tilde(i, j)**2), -pressure_p_prime(i, j) / a_tilde(i, j)**2]
          rho(i, j) = neumaier_sum(t_tot)
        enddo
      enddo
      !$omp end do

      !$omp do
      do j = 1, cones%nj
        do i = 1, cones%ni
          t1 = neumaier_sum(-p_i(:, i, j) * sin_dtheta(:, i, j)) / rho_a_tilde(i, j)
          if(abs(t1) < 1e-5_rk * a_tilde(i, j)) t1 = 0.0_rk

          t2 = neumaier_sum(u_i(:, i, j) * ((dtheta(:, i, j) / 2.0_rk) + (sin_d2theta(:, i, j) / 4.0_rk)))
          if(abs(t2) < 1e-5_rk * a_tilde(i, j)) t2 = 0.0_rk

          t3 = neumaier_sum(-v_i(:, i, j) * cos_d2theta(:, i, j)) / 4.0_rk
          if(abs(t3) < 1e-5_rk * a_tilde(i, j)) t3 = 0.0_rk

          t_tot = [t1, t2, t3]
          u(i, j) = neumaier_sum(t_tot) / pi
        enddo
      enddo
      !$omp end do

      !$omp do
      do j = 1, cones%nj
        do i = 1, cones%ni
          t1 = neumaier_sum(p_i(:, i, j) * cos_dtheta(:, i, j)) / rho_a_tilde(i, j)
          if(abs(t1) < 1e-5_rk * a_tilde(i, j)) t1 = 0.0_rk

          t2 = neumaier_sum(-u_i(:, i, j) * cos_d2theta(:, i, j)) / 4.0_rk
          if(abs(t2) < 1e-5_rk * a_tilde(i, j)) t2 = 0.0_rk

          t3 = neumaier_sum(v_i(:, i, j) * ((dtheta(:, i, j) / 2.0_rk) - (sin_d2theta(:, i, j) / 4.0_rk)))
          if(abs(t3) < 1e-5_rk * a_tilde(i, j)) t3 = 0.0_rk

          t_tot = [t1, t2, t3]
          v(i, j) = neumaier_sum(t_tot) / pi
        enddo
      enddo
      !$omp end do
      !$omp end parallel
    endassociate
    deallocate(rho_a_tilde)
  endsubroutine e0_operator

endmodule mod_local_evo_operator

max_issues_repo_name: smillerc/cato

max_forks_repo_name: smillerc/cato


text:
      subroutine impday(iyr,sdate)
c
c     + + + PURPOSE + + +
c
c     SR IMPDAY outputs and sums daily impoundment variables.
c
c     Called from: SR WSHDRV
c     Author(s): Jim Ascough II
c     Reference in User Guide:
c
c     Version:
c     Date recoded:
c     Recoded by:
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
c
      include 'pmximp.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
c
      integer iyr, sdate
c
c     + + + ARGUMENT DEFINITIONS + + +
c
c     iyr   -
c     sdate -
c
c     + + + COMMON BLOCKS + + +
c
      include 'cflags.inc'
c     read: ipdout
c
      include 'cimday.inc'
      include 'cimeos.inc'
      include 'cimitf.inc'
      include 'cimmon.inc'
      include 'cimpnd.inc'
      include 'cimsed.inc'
      include 'cimyrs.inc'
c
c     + + + LOCAL VARIABLES + + +
c
      character*4 otchar
      character*6 filchr
c
c     + + + LOCAL DEFINITIONS + + +
c
c     otchar - overtopping character string
c     filchr - sediment filling character string
c
c     + + + SAVES + + +
c
c     + + + SUBROUTINES CALLED + + +
c
c     + + + DATA INITIALIZATIONS + + +
c
c     + + + END SPECIFICATIONS + + +
c
c
c     begin multiple impoundment loop
c
      do 10 ipond = 1, npond
c
        otchar = '    '
        filchr = '      '
c
c       calculate daily total sediment
c
        totin(ipond) = ciin(ipond) * vi(ipond) * 0.4536
        ret(ipond) = totin(ipond) - cout(ipond)
c
        if (ret(ipond).lt.0.0) ret(ipond) = 0.0
c
        clin(ipond) = pcl(ipond) * totin(ipond)
        slin(ipond) = psl(ipond) * totin(ipond)
        sain(ipond) = psa(ipond) * totin(ipond)
        sdin(ipond) = psd(ipond) * totin(ipond)
        lain(ipond) = pla(ipond) * totin(ipond)
c
        clret(ipond) = clin(ipond) - clout(ipond)
        slret(ipond) = slin(ipond) - slout(ipond)
        saret(ipond) = sain(ipond) - saout(ipond)
        sdret(ipond) = sdin(ipond) - sdout(ipond)
        laret(ipond) = lain(ipond) - laout(ipond)
c
c       calculate sediment balance remaining in impoundment
c
        totint(ipond) = totint(ipond) + totin(ipond)
        coutt(ipond) = coutt(ipond) + cout(ipond)
        rett(ipond) = totint(ipond) - coutt(ipond)
c
c       convert variables for printout
c
c       inflow volume from ft^3 to m^3
c
        vi(ipond) = vi(ipond) / 3.281 ** 3
c
c       inflow and outflow rates from ft^3/s to m^3/s
c
        qiin(ipond) = qiin(ipond) / 3.281 ** 3
        qomx(ipond) = qomx(ipond) / 3.281 ** 3
c
c       stages from ft to m
c
        hmax(ipond) = hmax(ipond) / 3.281
        hmin(ipond) = hmin(ipond) / 3.281
c
c       concentration from lbs/ft^3 to mg/l
c
        ciin(ipond) = ciin(ipond) / 0.00006241778465
c
        if (hmax(ipond).gt.(hot(ipond)/3.281)) then
c
          otflg(ipond) = 1
          otyr(ipond) = otyr(ipond) + 1
          oteos(ipond) = oteos(ipond) + 1
c
          if (hmin(ipond).gt.hfull(ipond)) then
c
            if (filflg(ipond).ne.1) then
              filflg(ipond) = 1
              filday(ipond) = sdate
              filyr(ipond) = iyr
            end if
c
c           case 1: impoundment overtopped and filled with sediment
c
            otchar = '>Hot'
            filchr = '>Hfull'
          else
c
c           case 2: impoundment overtopped but not filled with sediment
c
            otchar = '>Hot'
          end if
c
        else
c
          if (hmin(ipond).gt.hfull(ipond)) then
c
            if (filflg(ipond).ne.1) then
              filflg(ipond) = 1
              filday(ipond) = sdate
              filyr(ipond) = iyr
            end if
c
c           case 3: impoundment not overtopped but filled with sediment
c
            filchr = '>Hfull'
          else
c
c         case 4: impoundment not overtopped and not filled with sediment
c
c         otchar, filchr initialized above
c
          end if
        end if
c
        if ((npond.eq.1).and.(ipdout.eq.1)) then
c
          write (57,1000) iyr, sdate, qomx(ipond), volo(ipond),
     1        qiin(ipond), vi(ipond), hmax(ipond), otchar,
     1        cpeak(ipond), ca(ipond), ciin(ipond), hmin(ipond), filchr
c
          write (56,1100) iyr, sdate, totin(ipond), cout(ipond),
     1        ret(ipond), rett(ipond), clin(ipond), clout(ipond),
     1        clret(ipond), slin(ipond), slout(ipond), slret(ipond),
     1        sain(ipond), saout(ipond), saret(ipond), sdin(ipond),
     1        sdout(ipond), sdret(ipond), lain(ipond), laout(ipond),
     1        laret(ipond)
c
        end if
c
        hmin(ipond) = hmin(ipond) * 3.281
c
c       monthly variables
c
        vom(ipond) = vom(ipond) + volo(ipond)
        vim(ipond) = vim(ipond) + vi(ipond)
c
        if (qiin(ipond).gt.qiinm(ipond)) qiinm(ipond) = qiin(ipond)
        if (qomx(ipond).gt.qomxm(ipond)) qomxm(ipond) = qomx(ipond)
        if (hmax(ipond).gt.hmaxm(ipond)) hmaxm(ipond) = hmax(ipond)
c
        totinm(ipond) = totinm(ipond) + totin(ipond)
c
        clinm(ipond) = clinm(ipond) + clin(ipond)
        slinm(ipond) = slinm(ipond) + slin(ipond)
        sainm(ipond) = sainm(ipond) + sain(ipond)
        sdinm(ipond) = sdinm(ipond) + sdin(ipond)
        lainm(ipond) = lainm(ipond) + lain(ipond)
c
        coutm(ipond) = coutm(ipond) + cout(ipond)
        cloutm(ipond) = cloutm(ipond) + clout(ipond)
        sloutm(ipond) = sloutm(ipond) + slout(ipond)
        saoutm(ipond) = saoutm(ipond) + saout(ipond)
        sdoutm(ipond) = sdoutm(ipond) + sdout(ipond)
        laoutm(ipond) = laoutm(ipond) + laout(ipond)
c
c       yearly variables
c
        voy(ipond) = voy(ipond) + volo(ipond)
        viy(ipond) = viy(ipond) + vi(ipond)
c
        if (qiin(ipond).gt.qiiny(ipond)) qiiny(ipond) = qiin(ipond)
        if (qomx(ipond).gt.qomxy(ipond)) qomxy(ipond) = qomx(ipond)
        if (hmax(ipond).gt.hmaxy(ipond)) hmaxy(ipond) = hmax(ipond)
c
        totiny(ipond) = totiny(ipond) + totin(ipond)
c
        cliny(ipond) = cliny(ipond) + clin(ipond)
        sliny(ipond) = sliny(ipond) + slin(ipond)
        sainy(ipond) = sainy(ipond) + sain(ipond)
        sdiny(ipond) = sdiny(ipond) + sdin(ipond)
        lainy(ipond) = lainy(ipond) + lain(ipond)
c
        couty(ipond) = couty(ipond) + cout(ipond)
        clouty(ipond) = clouty(ipond) + clout(ipond)
        slouty(ipond) = slouty(ipond) + slout(ipond)
        saouty(ipond) = saouty(ipond) + saout(ipond)
        sdouty(ipond) = sdouty(ipond) + sdout(ipond)
        laouty(ipond) = laouty(ipond) + laout(ipond)
c
c       end of simulation variables
c
        voe(ipond) = voe(ipond) + volo(ipond)
        vie(ipond) = vie(ipond) + vi(ipond)
c
        if (qiin(ipond).gt.qiine(ipond)) qiine(ipond) = qiin(ipond)
        if (qomx(ipond).gt.qomxe(ipond)) qomxe(ipond) = qomx(ipond)
        if (hmax(ipond).gt.hmaxe(ipond)) hmaxe(ipond) = hmax(ipond)
c
        totine(ipond) = totine(ipond) + totin(ipond)
c
        cline(ipond) = cline(ipond) + clin(ipond)
        sline(ipond) = sline(ipond) + slin(ipond)
        saine(ipond) = saine(ipond) + sain(ipond)
        sdine(ipond) = sdine(ipond) + sdin(ipond)
        laine(ipond) = laine(ipond) + lain(ipond)
c
        coute(ipond) = coute(ipond) + cout(ipond)
        cloute(ipond) = cloute(ipond) + clout(ipond)
        sloute(ipond) = sloute(ipond) + slout(ipond)
        saoute(ipond) = saoute(ipond) + saout(ipond)
        sdoute(ipond) = sdoute(ipond) + sdout(ipond)
        laoute(ipond) = laoute(ipond) + laout(ipond)
c
   10 continue
c
      return
c
 1000 format (i4,2x,i3,1x,f10.6,f10.2,f10.6,f10.2,f6.2,2x,a4,3f8.0,2x,
     1    'N/A',1x,f6.2,2x,a6)
 1100 format (i4,2x,i3,1x,19(f12.2,1x))
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
      PROGRAM LIFE
c     Main game
      INTEGER plane(100,100)
      INTEGER x, y, psx, psy
      INTEGER alive

      alive = 1

c     Clears game of all life
      CALL NUKE(plane)

c     PostScript
      PRINT *, '%!PS'
      x = 1
      y = 1
      psx = 100
      psy = 100
      DO WHILE (y .LT. 101)
         psx = 0
         psy = 0
         DO WHILE (x .LT. 101)
            plane (x,y) = 0
            PRINT *, 5, 5, 'translate'
            PRINT *, 1, alive, alive, 'csquare'
            x = x + 1
            psx = psx + 10
            IF (alive .EQ. 1) THEN
               count = 1
               DO WHILE (count .LT. 10001)
                  PRINT *, '000000'
                  count = count + 1
               ENDDO
            ENDIF
         ENDDO
         x = 1
         psx = 100
         y = y + 1
         psy = psy + 10
      ENDDO

      PRINT *, 'stroke'
      PRINT *, 'showpage'

      END PROGRAM LIFE

max_stars_repo_name: asnowmang/GameofLife

max_issues_repo_name: asnowmang/GameofLife

max_forks_repo_name: asnowmang/GameofLife