text:
C ###############################################################
C #                                                             #
C #                    THE VLIDORT  MODEL                       #
C #                                                             #
C #  Vectorized LInearized Discrete Ordinate Radiative Transfer #
C #  -          --         -        -        -         -        #
C #                                                             #
C ###############################################################

C ###############################################################
C #                                                             #
C #  Author :      Robert. J. D. Spurr                          #
C #                                                             #
C #  Address :      RT Solutions, Inc.                          #
C #            9 Channing Street                                #
C #             Cambridge, MA 02138, USA                        #
C #            Tel: (617) 492 1183                              #
C #                                                             #
C #  Email :      rtsolutions@verizon.net                       #
C #                                                             #
C #  Versions     :   2.0, 2.2, 2.3, 2.4, 2.4R, 2.4RT           #
C #  Release Date :   December 2005  (2.0)                      #
C #  Release Date :   March 2007     (2.2)                      #
C #  Release Date :   October 2007   (2.3)                      #
C #  Release Date :   December 2008  (2.4)                      #
C #  Release Date :   April/May 2009 (2.4R)                     #
C #  Release Date :   July 2009      (2.4RT)                    #
C #                                                             #
C #       NEW: TOTAL COLUMN JACOBIANS         (2.4)             #
C #       NEW: BPDF Land-surface KERNELS      (2.4R)            #
C #       NEW: Thermal Emission Treatment     (2.4RT)           #
C #                                                             #
C ###############################################################

C    #####################################################
C    #                                                   #
C    #   This Version of VLIDORT comes with a GNU-style  #
C    #   license. Please read the license carefully.     #
C    #                                                   #
C    #####################################################

C ###############################################################
C #                                                             #
C # Subroutines in this Module                                  #
C #                                                             #
C #       VLIDORT_INTEGRATED_OUTPUT (master)                    #
C #                                                             #
C #              QUADINTENS_LEVEL_UP                            #
C #              QUADINTENS_LEVEL_DN                            #
C #              QUADINTENS_OFFGRID_UP                          #
C #              QUADINTENS_OFFGRID_DN                          #
C #                                                             #
C ###############################################################

      SUBROUTINE VLIDORT_INTEGRATED_OUTPUT
     I  ( DO_INCLUDE_MVOUTPUT,
     I    DO_INCLUDE_DIRECTBEAM,
     I    DO_INCLUDE_THERMEMISS,
     I    FLUX_MULTIPLIER, IBEAM )

C  include files
C  -------------

C  include file of dimensions and numbers

      INCLUDE '../includes/VLIDORT.PARS'

C  Include files of input and bookkeeping variables

      INCLUDE '../includes/VLIDORT_INPUTS.VARS'
      INCLUDE '../includes/VLIDORT_BOOKKEEP.VARS'

C  include files of main model variables (input to this module)

      INCLUDE '../includes/VLIDORT_SETUPS.VARS'

C  include file of result variables (module output stored here)

      INCLUDE '../includes/VLIDORT_RESULTS.VARS'

C  Subroutine input arguments
C  --------------------------

C  Control inclusion

      LOGICAL          DO_INCLUDE_MVOUTPUT
      LOGICAL          DO_INCLUDE_DIRECTBEAM
      LOGICAL          DO_INCLUDE_THERMEMISS

C  PI number and flux multiplier

      INTEGER          IBEAM
      DOUBLE PRECISION FLUX_MULTIPLIER

C  local variables
C  ---------------

C  Local qudrature output (for debug)

      DOUBLE PRECISION QSTOKES_F
     &       ( MAX_USER_LEVELS, MAXSTREAMS, MAXSTOKES )

C  help variables

      INTEGER          I, IDIR, WDIR, UTA, UT, N, NLEVEL, O1
      DOUBLE PRECISION SUM_MI, SUM_FX, FTRANS
      DOUBLE PRECISION DIRECT_TRANS, DIRECT_FLUX, DIRECT_MEANI

C  mean intensity and flux
C  -----------------------

C  direction loop

      DO IDIR = 1, N_DIRECTIONS
        WDIR = WHICH_DIRECTIONS(IDIR)

C  Upwelling Stokes output at Quadrature angles

        IF ( WDIR .EQ. UPIDX ) THEN
          DO UTA = 1, N_USER_LEVELS
            NLEVEL = UTAU_LEVEL_MASK_UP(UTA)
            IF ( PARTLAYERS_OUTFLAG(UTA) ) THEN
              UT = PARTLAYERS_OUTINDEX(UTA)
              N  = PARTLAYERS_LAYERIDX(UT)
              CALL QUADINTENS_OFFGRID_UP ( N, UTA, UT, IBEAM,
     &            FLUX_MULTIPLIER, DO_INCLUDE_THERMEMISS, QSTOKES_F )
            ELSE
              CALL QUADINTENS_LEVEL_UP ( NLEVEL, UTA,
     &            FLUX_MULTIPLIER, QSTOKES_F )
            ENDIF
          ENDDO
        ENDIF

C  Downwelling Stokes output at Quadrature angles

        IF ( WDIR .EQ. DNIDX ) THEN
          DO UTA = 1, N_USER_LEVELS
            NLEVEL = UTAU_LEVEL_MASK_DN(UTA)
            IF ( PARTLAYERS_OUTFLAG(UTA) ) THEN
              UT = PARTLAYERS_OUTINDEX(UTA)
              N  = PARTLAYERS_LAYERIDX(UT)
              CALL QUADINTENS_OFFGRID_DN ( N, UTA, UT, IBEAM,
     &            FLUX_MULTIPLIER, DO_INCLUDE_THERMEMISS, QSTOKES_F )
            ELSE
              CALL QUADINTENS_LEVEL_DN ( NLEVEL, UTA,
     &            FLUX_MULTIPLIER, QSTOKES_F )
            ENDIF
          ENDDO
        ENDIF

C  Mean Intensity and Flux  output
C  -------------------------------

        IF ( DO_INCLUDE_MVOUTPUT ) THEN

C  Diffuse term output

          DO UTA = 1, N_USER_LEVELS
            DO O1 = 1, NSTOKES
              SUM_MI = ZERO
              SUM_FX = ZERO
              DO I = 1, NSTREAMS
                SUM_MI = SUM_MI + QUAD_WEIGHTS(I) * QSTOKES_F(UTA,I,O1)
                SUM_FX = SUM_FX + QUAD_STRMWTS(I) * QSTOKES_F(UTA,I,O1)
              ENDDO
              MEAN_STOKES(UTA,IBEAM,O1,WDIR) = SUM_MI * HALF
              FLUX_STOKES(UTA,IBEAM,O1,WDIR) = SUM_FX * PI2
           ENDDO
          ENDDO

C  nothing to do if no solar sources

          IF ( .NOT. DO_INCLUDE_DIRECTBEAM ) GO TO 455

C  For the downward direction, add the direct beam contributions

          IF ( WDIR .EQ. DNIDX ) THEN

C  loop over all the output optical depths

           DO UTA = 1, N_USER_LEVELS

C  For the offgrid values.......
C     .....Only contributions for layers above the PI cutoff

            IF ( PARTLAYERS_OUTFLAG(UTA) ) THEN
              UT = PARTLAYERS_OUTINDEX(UTA)
              N  = PARTLAYERS_LAYERIDX(UT)
              IF ( N .LE. LAYER_PIS_CUTOFF(IBEAM) ) THEN
                DIRECT_TRANS =
     &             INITIAL_TRANS(N,IBEAM) * T_UTDN_MUBAR(UT,IBEAM)
                DO O1 = 1, NSTOKES
                  FTRANS = FLUXVEC(O1) * DIRECT_TRANS
                  DIRECT_FLUX  = FTRANS * LOCAL_CSZA(N,IBEAM)
                  DIRECT_MEANI = FTRANS / PI4
                  MEAN_DIRECT(UTA,IBEAM,O1) = DIRECT_MEANI
                  FLUX_DIRECT(UTA,IBEAM,O1) = DIRECT_FLUX
                  MEAN_STOKES(UTA,IBEAM,O1,WDIR) =
     *                  MEAN_STOKES(UTA,IBEAM,O1,WDIR) + DIRECT_MEANI
                  FLUX_STOKES(UTA,IBEAM,O1,WDIR) =
     *                  FLUX_STOKES(UTA,IBEAM,O1,WDIR) + DIRECT_FLUX
                ENDDO
              ENDIF

C  For the on-grid values

             ELSE
              N = UTAU_LEVEL_MASK_DN(UTA)
              IF ( N .LE. LAYER_PIS_CUTOFF(IBEAM) ) THEN
                IF ( N .EQ. 0 ) THEN
                  DIRECT_TRANS = ONE
                ELSE
                  DIRECT_TRANS =
     &               INITIAL_TRANS(N,IBEAM)*T_DELT_MUBAR(N,IBEAM)
                ENDIF
                DO O1 = 1, NSTOKES
                  FTRANS = FLUXVEC(O1) * DIRECT_TRANS
                  DIRECT_FLUX  = FTRANS * LOCAL_CSZA(N,IBEAM)
                  DIRECT_MEANI = FTRANS / PI4
                  MEAN_DIRECT(UTA,IBEAM,O1) = DIRECT_MEANI
                  FLUX_DIRECT(UTA,IBEAM,O1) = DIRECT_FLUX
                  MEAN_STOKES(UTA,IBEAM,O1,WDIR) =
     *                  MEAN_STOKES(UTA,IBEAM,O1,WDIR) + DIRECT_MEANI
                  FLUX_STOKES(UTA,IBEAM,O1,WDIR) =
     *                  FLUX_STOKES(UTA,IBEAM,O1,WDIR) + DIRECT_FLUX
                ENDDO
              ENDIF
             ENDIF
           ENDDO

C  Finish downwelling direct contribution

          ENDIF

C  Finish MV output

        ENDIF

C  Continuation point for avoiding direct beam calculation

 455    CONTINUE

C  Finish directional loop

      ENDDO

C  Finish

      RETURN
      END

C

      SUBROUTINE QUADINTENS_LEVEL_UP
     I    ( NLEVEL, UTA, FLUX_MULTIPLIER,
     O      QSTOKES_F )

C  include files
C  -------------

C  include file of dimensions and numbers

      INCLUDE '../includes/VLIDORT.PARS'

C  Include files of input and bookkeeping variables

      INCLUDE '../includes/VLIDORT_INPUTS.VARS'
      INCLUDE '../includes/VLIDORT_BOOKKEEP.VARS'

C  include files of main model variables (input to this module)

      INCLUDE '../includes/VLIDORT_SETUPS.VARS'
      INCLUDE '../includes/VLIDORT_SOLUTION.VARS'
      INCLUDE '../includes/VLIDORT_THERMALSUP.VARS'
      INCLUDE '../includes/VLIDORT_REFLECTANCE.VARS'

C  Subroutine arguments
C  --------------------

C  inputs

      INTEGER          NLEVEL, UTA
      DOUBLE PRECISION FLUX_MULTIPLIER

C  output

      DOUBLE PRECISION QSTOKES_F
     &       ( MAX_USER_LEVELS, MAXSTREAMS, MAXSTOKES )

C  local variables
C  ---------------

      INTEGER          N, NL, I, I1, O1, K, KO1, K0, K1, K2
      DOUBLE PRECISION SPAR, SHOM, HOM1, HOM2, SHOM_R
      DOUBLE PRECISION SHOM_CR, HOM1CR, HOM2CR, TPROP, THELP
      DOUBLE PRECISION LXR, MXR, LXR_CR, LXR_CI, MXR_CR, MXR_CI

C  For those optical depths at layer boundaries
C  --------------------------------------------

C  This depends on the level mask - if this is 0 to NLAYERS - 1, then we are
C  looking at the intensity at the top of these layers. The
C  case where the level mask = NLAYERS is the upwelling intensity
C  at the bottom of the atmosphere (treated separately).

      NL = NLEVEL
      N  = NL + 1

C  Lowest level contributions
C  ==========================

C  Lowest level, thermal transmittance only

      IF ( NL .EQ. NLAYERS .and. DO_THERMAL_TRANSONLY ) THEN
        O1 = 1
        DO I = 1, NSTREAMS
          THELP = BOA_THTONLY_SOURCE(I)
          QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * THELP
        ENDDO
        RETURN
      ENDIF

C  For the lowest level, scattering solution
C  -----------------------------------------

      IF ( NL .EQ. NLAYERS ) THEN

C  Start main loops

        KO1 = K_REAL(NL) + 1
        DO I = 1, NSTREAMS
          I1 = I + NSTREAMS
          DO O1 = 1, NSTOKES

C  real homogeneous solutions

            SHOM_R = ZERO
            DO K = 1, K_REAL(NL)
              LXR = LCON(K,NL)*SOLA_XPOS(I1,O1,K,NL)
              MXR = MCON(K,NL)*SOLB_XNEG(I1,O1,K,NL)
              HOM1 = LXR * T_DELT_EIGEN(K,NL)
              HOM2 = MXR
              SHOM_R = SHOM_R + HOM1 + HOM2
            ENDDO

C  complex homogeneous solutions

            SHOM_CR = ZERO
            DO K = 1, K_COMPLEX(NL)
              K0 = 2*K-2
              K1 = KO1 + K0
              K2 = K1  + 1
              LXR_CR =  LCON(K1,NL) * SOLA_XPOS(I1,O1,K1,NL) -
     &                  LCON(K2,NL) * SOLA_XPOS(I1,O1,K2,NL)
              LXR_CI =  LCON(K1,NL) * SOLA_XPOS(I1,O1,K2,NL) +
     &                  LCON(K2,NL) * SOLA_XPOS(I1,O1,K1,NL)
              MXR_CR =  MCON(K1,NL) * SOLB_XNEG(I1,O1,K1,NL) -
     &                  MCON(K2,NL) * SOLB_XNEG(I1,O1,K2,NL)
              HOM1CR =  LXR_CR*T_DELT_EIGEN(K1,NL)
     &                - LXR_CI*T_DELT_EIGEN(K2,NL)
              HOM2CR = MXR_CR
              SHOM_CR = SHOM_CR + HOM1CR + HOM2CR
            ENDDO

C  real part, add particular solution, complete result

            SHOM = SHOM_R + SHOM_CR
            SPAR = WLOWER(I1,O1,NL)
            QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * ( SPAR + SHOM )

C  Finish streams/stokes loops

          ENDDO
        ENDDO

C  End lowest level clause

      ENDIF

C  For other levels in the atmosphere
C  ==================================

C  For other levels, thermal transmittance only

      IF ( NL .NE. NLAYERS .and. DO_THERMAL_TRANSONLY ) THEN
        O1 = 1
        DO I = 1, NSTREAMS
          I1 = I + NSTREAMS
          THELP = BOA_THTONLY_SOURCE(I)
          DO K = NLAYERS, N, -1
            TPROP = T_WUPPER(I1,K)/QUAD_STREAMS(I)
            THELP = THELP * T_DELT_DISORDS(I,K) + TPROP
          ENDDO
          QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * THELP
        ENDDO
        RETURN
      ENDIF

C  For other levels, scattering solution
C  -------------------------------------

      IF ( NL .NE. NLAYERS ) THEN

C  start main loops

        KO1 = K_REAL(N) + 1
        DO I = 1, NSTREAMS
          I1 = I + NSTREAMS
          DO O1 = 1, NSTOKES

C  real homogeneous solutions

            SHOM_R = ZERO
            DO K = 1, K_REAL(N)
              LXR = LCON(K,N)*SOLA_XPOS(I1,O1,K,N)
              MXR = MCON(K,N)*SOLB_XNEG(I1,O1,K,N)
              HOM1 = LXR
              HOM2 = MXR * T_DELT_EIGEN(K,N)
              SHOM_R = SHOM_R + HOM1 + HOM2
            ENDDO

C  complex homogeneous solutions

            SHOM_CR = ZERO
            DO K = 1, K_COMPLEX(N)
              K0 = 2*K-2
              K1 = KO1 + K0
              K2 = K1  + 1
              LXR_CR =  LCON(K1,N) * SOLA_XPOS(I1,O1,K1,N) -
     &                  LCON(K2,N) * SOLA_XPOS(I1,O1,K2,N)
              MXR_CR =  MCON(K1,N) * SOLB_XNEG(I1,O1,K1,N) -
     &                  MCON(K2,N) * SOLB_XNEG(I1,O1,K2,N)
              MXR_CI =  MCON(K1,N) * SOLB_XNEG(I1,O1,K2,N) +
     &                  MCON(K2,N) * SOLB_XNEG(I1,O1,K1,N)
              HOM1CR = LXR_CR
              HOM2CR =  MXR_CR * T_DELT_EIGEN(K1,N)
     &                - MXR_CI * T_DELT_EIGEN(K2,N)
              SHOM_CR = SHOM_CR + HOM1CR + HOM2CR
            ENDDO

C  real part, add particular solution, complete result

            SHOM = SHOM_R + SHOM_CR
            SPAR = WUPPER(I1,O1,N)
            QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * ( SPAR + SHOM )

C  Finish streams/stokes loops

          ENDDO
        ENDDO
      ENDIF

C  Finish

      RETURN
      END

C

      SUBROUTINE QUADINTENS_LEVEL_DN
     I    ( NLEVEL, UTA, FLUX_MULTIPLIER,
     O      QSTOKES_F )

C  include files
C  -------------

C  include file of dimensions and numbers

      INCLUDE '../includes/VLIDORT.PARS'

C  Include files of input and bookkeeping variables

      INCLUDE '../includes/VLIDORT_INPUTS.VARS'
      INCLUDE '../includes/VLIDORT_BOOKKEEP.VARS'

C  include files of main model variables (input to this module)

      INCLUDE '../includes/VLIDORT_SETUPS.VARS'
      INCLUDE '../includes/VLIDORT_SOLUTION.VARS'
      INCLUDE '../includes/VLIDORT_THERMALSUP.VARS'
      INCLUDE '../includes/VLIDORT_REFLECTANCE.VARS'

C  Subroutine arguments
C  --------------------

C  inputs

      INTEGER          NLEVEL, UTA
      DOUBLE PRECISION FLUX_MULTIPLIER

C  output

      DOUBLE PRECISION QSTOKES_F
     &       ( MAX_USER_LEVELS, MAXSTREAMS, MAXSTOKES )

C  local variables
C  ---------------

      INTEGER          N, I, O1, K, KO1, K0, K1, K2
      DOUBLE PRECISION SPAR, SHOM, HOM1, HOM2, SHOM_R
      DOUBLE PRECISION SHOM_CR, HOM1CR, HOM2CR, TPROP, THELP
      DOUBLE PRECISION LXR, MXR, LXR_CR, LXR_CI, MXR_CR

C  For those optical depths at layer boundaries
C  --------------------------------------------

      N = NLEVEL

C  Downwelling radiation at TOA ( or N = 0 ) is zero
C      A check has been performed on this.

      IF ( NLEVEL .EQ. 0 ) THEN
        DO I = 1, NSTREAMS
          DO O1 = 1, NSTOKES
            QSTOKES_F(UTA,I,O1) = ZERO
          ENDDO
        ENDDO
        RETURN
      ENDIF

C  Other levels, Thermal transmittance-only solution

      IF ( NLEVEL .NE. 0 .and. DO_THERMAL_TRANSONLY ) THEN
        O1 = 1
        DO I = 1, NSTREAMS
          THELP = ZERO
          DO K = 1, N
            TPROP = T_WLOWER(I,K)/QUAD_STREAMS(I)
            THELP = THELP * T_DELT_DISORDS(I,K) + TPROP
          ENDDO
          QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * THELP
        ENDDO
        RETURN
      ENDIF

C  Other levels, scattering solution
C  ---------------------------------

      IF ( NLEVEL .NE. 0 ) THEN

C  start main loops

        KO1 = K_REAL(N) + 1
        DO I = 1, NSTREAMS
          DO O1 = 1, NSTOKES

C  real homogeneous solutions

            SHOM_R = ZERO
            DO K = 1, K_REAL(N)
              LXR  = LCON(K,N)*SOLA_XPOS(I,O1,K,N)
              MXR  = MCON(K,N)*SOLB_XNEG(I,O1,K,N)
              HOM1 = LXR * T_DELT_EIGEN(K,N)
              HOM2 = MXR
              SHOM_R = SHOM_R + HOM1 + HOM2
            ENDDO

C  complex homogeneous solutions

            SHOM_CR = ZERO
            DO K = 1, K_COMPLEX(N)
              K0 = 2*K-2
              K1 = KO1 + K0
              K2 = K1  + 1
              LXR_CR =  LCON(K1,N) * SOLA_XPOS(I,O1,K1,N) -
     &                  LCON(K2,N) * SOLA_XPOS(I,O1,K2,N)
              LXR_CI =  LCON(K1,N) * SOLA_XPOS(I,O1,K2,N) +
     &                  LCON(K2,N) * SOLA_XPOS(I,O1,K1,N)
              MXR_CR =  MCON(K1,N) * SOLB_XNEG(I,O1,K1,N) -
     &                  MCON(K2,N) * SOLB_XNEG(I,O1,K2,N)
              HOM1CR = LXR_CR*T_DELT_EIGEN(K1,N)
     &                -LXR_CI*T_DELT_EIGEN(K2,N)
              HOM2CR = MXR_CR
              SHOM_CR = SHOM_CR + HOM1CR + HOM2CR
            ENDDO

C  real part, add particular solution, complete result

            SHOM = SHOM_R + SHOM_CR
            SPAR = WLOWER(I,O1,N)
            QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * ( SPAR + SHOM )

C  Finish streams/stokes loops

          ENDDO
        ENDDO

      ENDIF

C  Finish

      RETURN
      END

C

      SUBROUTINE QUADINTENS_OFFGRID_UP
     I  ( N, UTA, UT, IBEAM, FLUX_MULTIPLIER, DO_INCLUDE_THERMEMISS,
     O    QSTOKES_F )

C  include files
C  -------------

C  include file of dimensions and numbers

      INCLUDE '../includes/VLIDORT.PARS'

C  Include files of input and bookkeeping variables

      INCLUDE '../includes/VLIDORT_INPUTS.VARS'
      INCLUDE '../includes/VLIDORT_BOOKKEEP.VARS'

C  include files of main model variables (input to this module)

      INCLUDE '../includes/VLIDORT_SETUPS.VARS'
      INCLUDE '../includes/VLIDORT_SOLUTION.VARS'
      INCLUDE '../includes/VLIDORT_THERMALSUP.VARS'
      INCLUDE '../includes/VLIDORT_MULTIPLIERS.VARS'
      INCLUDE '../includes/VLIDORT_REFLECTANCE.VARS'

C  Subroutine arguments
C  --------------------

C  inputs

      INTEGER          N, UTA, UT, IBEAM
      DOUBLE PRECISION FLUX_MULTIPLIER
      LOGICAL          DO_INCLUDE_THERMEMISS

C  output

      DOUBLE PRECISION QSTOKES_F
     &       ( MAX_USER_LEVELS, MAXSTREAMS, MAXSTOKES )

C  local variables
C  ---------------

      INTEGER          I, I1, O1, K, KO1, K0, K1, K2
      DOUBLE PRECISION SPAR, SHOM, HOM1, HOM2, SHOM_R
      DOUBLE PRECISION SHOM_CR, HOM1CR, HOM2CR, TPROP, THELP
      DOUBLE PRECISION LXR, MXR, LXR_CR, LXR_CI, MXR_CR, MXR_CI

C  Thermal Transmittance only
c  --------------------------

      IF ( DO_THERMAL_TRANSONLY ) THEN
        O1 = 1
        DO I = 1, NSTREAMS
          I1 = I + NSTREAMS
          THELP = BOA_THTONLY_SOURCE(I)
          DO K = NLAYERS, N + 1, -1
            TPROP = T_WUPPER(I1,K) / QUAD_STREAMS(I)
            THELP = THELP * T_DELT_DISORDS(I,K) + TPROP
          ENDDO
          TPROP =  UT_T_PARTIC(I1,UT) / QUAD_STREAMS(I)
          THELP = THELP * T_DISORDS_UTUP(I,UT) + TPROP
          QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * THELP
        ENDDO
        RETURN
      ENDIF

C  For those optical depths at off-grid levels
C  -------------------------------------------

C  Homogeneous solution

      KO1 = K_REAL(N) + 1
      DO I = 1, NSTREAMS
        I1 = I + NSTREAMS
        DO O1 = 1, NSTOKES

C  real homogeneous solutions

          SHOM_R = ZERO
          DO K = 1, K_REAL(N)
            LXR  = LCON(K,N)*SOLA_XPOS(I1,O1,K,N)
            MXR  = MCON(K,N)*SOLB_XNEG(I1,O1,K,N)
            HOM1 = LXR * T_UTDN_EIGEN(K,UT)
            HOM2 = MXR * T_UTUP_EIGEN(K,UT)
            SHOM_R = SHOM_R + HOM1 + HOM2
          ENDDO

C  complex homogeneous solutions

          SHOM_CR = ZERO
          DO K = 1, K_COMPLEX(N)
            K0 = 2*K-2
            K1 = KO1 + K0
            K2 = K1  + 1
            LXR_CR =  LCON(K1,N) * SOLA_XPOS(I1,O1,K1,N) -
     &                LCON(K2,N) * SOLA_XPOS(I1,O1,K2,N)
            LXR_CI =  LCON(K1,N) * SOLA_XPOS(I1,O1,K2,N) +
     &                LCON(K2,N) * SOLA_XPOS(I1,O1,K1,N)
            MXR_CR =  MCON(K1,N) * SOLB_XNEG(I1,O1,K1,N) -
     &                MCON(K2,N) * SOLB_XNEG(I1,O1,K2,N)
            MXR_CI =  MCON(K1,N) * SOLB_XNEG(I1,O1,K2,N) +
     &                MCON(K2,N) * SOLB_XNEG(I1,O1,K1,N)
            HOM1CR =    LXR_CR * T_UTDN_EIGEN(K1,UT)
     &                - LXR_CI * T_UTDN_EIGEN(K2,UT)
            HOM2CR =    MXR_CR * T_UTUP_EIGEN(K1,UT)
     &                - MXR_CI * T_UTUP_EIGEN(K2,UT)
            SHOM_CR = SHOM_CR + HOM1CR + HOM2CR
          ENDDO

C  real part

          SHOM = SHOM_R + SHOM_CR
          QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * SHOM

C  Finish streams/stokes loops

        ENDDO
      ENDDO

C  Add the thermal solution  (if flagged)

      IF ( DO_INCLUDE_THERMEMISS ) THEN
        O1 = 1
        DO I = 1, NSTREAMS
          I1 = I + NSTREAMS
          SPAR = UT_T_PARTIC(I1, UT)
           QSTOKES_F(UTA,I,O1) =
     &       QSTOKES_F(UTA,I,O1) + FLUX_MULTIPLIER * SPAR
        ENDDO
      ENDIF

C  Finished if no solar terms

      IF ( .NOT. DO_SOLAR_SOURCES ) RETURN

C  Add the solar particular solutions
C     Green function solution is still a Placeholder....!

      IF ( DO_CLASSICAL_SOLUTION ) THEN
        DO I = 1, NSTREAMS
          I1 = I + NSTREAMS
          DO O1 = 1, NSTOKES
            SPAR = WUPPER(I1,O1,N) * T_UTDN_MUBAR(UT,IBEAM)
            QSTOKES_F(UTA,I,O1) =
     &        QSTOKES_F(UTA,I,O1) + FLUX_MULTIPLIER * SPAR
          ENDDO
        ENDDO
      ELSE
C                P L A C E H O L D E R
      ENDIF

C  Finish

      RETURN
      END

C

      SUBROUTINE QUADINTENS_OFFGRID_DN
     I  ( N, UTA, UT, IBEAM, FLUX_MULTIPLIER, DO_INCLUDE_THERMEMISS,
     O    QSTOKES_F )

C  include files
C  -------------

C  include file of dimensions and numbers

      INCLUDE '../includes/VLIDORT.PARS'

C  Include files of input and bookkeeping variables

      INCLUDE '../includes/VLIDORT_INPUTS.VARS'
      INCLUDE '../includes/VLIDORT_BOOKKEEP.VARS'

C  include files of main model variables (input to this module)

      INCLUDE '../includes/VLIDORT_SETUPS.VARS'
      INCLUDE '../includes/VLIDORT_SOLUTION.VARS'
      INCLUDE '../includes/VLIDORT_THERMALSUP.VARS'
      INCLUDE '../includes/VLIDORT_MULTIPLIERS.VARS'

C  Subroutine arguments
C  --------------------

C  inputs

      INTEGER          N, UTA, UT, IBEAM
      DOUBLE PRECISION FLUX_MULTIPLIER
      LOGICAL          DO_INCLUDE_THERMEMISS

C  output

      DOUBLE PRECISION QSTOKES_F
     &       ( MAX_USER_LEVELS, MAXSTREAMS, MAXSTOKES )

C  local variables
C  ---------------

      INTEGER          I, O1, K, KO1, K0, K1, K2
      DOUBLE PRECISION SPAR, SHOM, HOM1, HOM2, SHOM_R
      DOUBLE PRECISION SHOM_CR, HOM1CR, HOM2CR, TPROP, THELP
      DOUBLE PRECISION LXR, MXR, LXR_CR, LXR_CI, MXR_CR, MXR_CI

C  Thermal Transmittance only
c  --------------------------

      IF ( DO_THERMAL_TRANSONLY ) THEN
        O1 = 1
        DO I = 1, NSTREAMS
          THELP = ZERO
          DO K = 1, N - 1
            TPROP = T_WLOWER(I,K) / QUAD_STREAMS(I)
            THELP = THELP*T_DELT_DISORDS(I,K) + TPROP
          ENDDO
          TPROP = UT_T_PARTIC(I,UT) / QUAD_STREAMS(I)
          THELP = THELP * T_DISORDS_UTDN(I,UT) + TPROP
          QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * THELP
        ENDDO
        RETURN
      ENDIF

C  For those optical depths at off-grid levels
C  -------------------------------------------

C  Homogeneous solution

      KO1 = K_REAL(N) + 1
      DO I = 1, NSTREAMS
        DO O1 = 1, NSTOKES

C  real homogeneous solutions

          SHOM_R = ZERO
          DO K = 1, K_REAL(N)
            LXR  = LCON(K,N)*SOLA_XPOS(I,O1,K,N)
            MXR  = MCON(K,N)*SOLB_XNEG(I,O1,K,N)
            HOM1 = LXR * T_UTDN_EIGEN(K,UT)
            HOM2 = MXR * T_UTUP_EIGEN(K,UT)
            SHOM_R = SHOM_R + HOM1 + HOM2
          ENDDO

C  complex homogeneous solutions

          SHOM_CR = ZERO
          DO K = 1, K_COMPLEX(N)
            K0 = 2*K-2
            K1 = KO1 + K0
            K2 = K1  + 1
            LXR_CR =  LCON(K1,N) * SOLA_XPOS(I,O1,K1,N) -
     &                LCON(K2,N) * SOLA_XPOS(I,O1,K2,N)
            LXR_CI =  LCON(K1,N) * SOLA_XPOS(I,O1,K2,N) +
     &                LCON(K2,N) * SOLA_XPOS(I,O1,K1,N)
            MXR_CR =  MCON(K1,N) * SOLB_XNEG(I,O1,K1,N) -
     &                MCON(K2,N) * SOLB_XNEG(I,O1,K2,N)
            MXR_CI =  MCON(K1,N) * SOLB_XNEG(I,O1,K2,N) +
     &                MCON(K2,N) * SOLB_XNEG(I,O1,K1,N)
            HOM1CR =    LXR_CR * T_UTDN_EIGEN(K1,UT)
     &                - LXR_CI * T_UTDN_EIGEN(K2,UT)
            HOM2CR =    MXR_CR * T_UTUP_EIGEN(K1,UT)
     &                - MXR_CI * T_UTUP_EIGEN(K2,UT)
            SHOM_CR = SHOM_CR + HOM1CR + HOM2CR
          ENDDO

C  real part

          SHOM = SHOM_R + SHOM_CR
          QSTOKES_F(UTA,I,O1) = FLUX_MULTIPLIER * SHOM

C  Finish streams/stokes loops

        ENDDO
      ENDDO

C  Add the thermal solution  (if flagged)

      IF ( DO_INCLUDE_THERMEMISS ) THEN
        O1 = 1
        DO I = 1, NSTREAMS
          SPAR = UT_T_PARTIC(I,UT)
           QSTOKES_F(UTA,I,O1) =
     &       QSTOKES_F(UTA,I,O1) + FLUX_MULTIPLIER * SPAR
        ENDDO
      ENDIF

C  Finished if no solar terms

      IF ( .NOT. DO_SOLAR_SOURCES ) RETURN

C  Add the solar particular solutions
C     Green function solution is still a Placeholder....!

      IF ( DO_CLASSICAL_SOLUTION ) THEN
        DO I = 1, NSTREAMS
          DO O1 = 1, NSTOKES
            SPAR = WUPPER(I,O1,N) * T_UTDN_MUBAR(UT,IBEAM)
            QSTOKES_F(UTA,I,O1) =
     &        QSTOKES_F(UTA,I,O1) + FLUX_MULTIPLIER * SPAR
          ENDDO
        ENDDO
      ELSE
C                P L A C E H O L D E R
      ENDIF

C  Finish

      RETURN
      END

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
  subroutine tinvu (sDCM, tin, cmin, tn, cmn, sintin, costin, sinfin, &
       & cosfin, sintn, costn, sinfn, cosfn, v, u, tin1)

! ======================================================================
!
!     Calculation of tin1, v, & u
!     tin1 is kinetic energy of the cascade particle in a frame in
!         which the nuclear partner particle is at rest.
!     u is the total energy available in the center-of-momentum frame.
!     v is the velocity of the center-of-momentum frame of the
!         interacting particles with respect to the nucleus.
!
!     The four variables ....in correspond to the direction of the
!     momentum of the cascade particle; those ....n to the direction
!     of the target particle (or two-particle state, for absorption)
!     in the lab.
!
!     Called by: ABSORP POINTE
!
!   CEM95 written by S. G. Mashnik
!
!   Edited by A. J. Sierk,  LANL  T-2  February-March, 1996.
!   Edited by AJS  July, 1997.
!   Edited by A. J. Sierk, LANL T-16  October, 2003.
!   Edited by AJS, LANL T-2, December, 2011.
!   Edited by LMK, XCP-3, July 2013 (included error protection)
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64
    use standardDCMParams, only: one, two

    implicit none
    class(StandardDCM), intent(inout) :: sDCM
    real(real64), intent(in   ) :: tin      ! Kinetic and rest energies of particle 1
    real(real64), intent(in   ) :: cmin
    real(real64), intent(in   ) :: tn       ! Kinetic and rest energies of particle 2
    real(real64), intent(in   ) :: cmn
    real(real64), intent(in   ) :: sintin   ! Angles of one particle 1
    real(real64), intent(in   ) :: costin
    real(real64), intent(in   ) :: sinfin
    real(real64), intent(in   ) :: cosfin
    real(real64), intent(in   ) :: sintn    ! Angles of one particle 2
    real(real64), intent(in   ) :: costn
    real(real64), intent(in   ) :: sinfn
    real(real64), intent(in   ) :: cosfn
    real(real64), intent(  out) :: v(3)     ! Velocity in CP system of interacting particles
    real(real64), intent(  out) :: u        ! Energy available in CP system
    real(real64), intent(  out) :: tin1     ! Kinetic energy of cascade partner

    real(real64) :: ein, en, fac, pin, pinx, piny, pinz, pn, pnx, pny, &
         & pnz, temp, v2

! ======================================================================

!   Momentum and total energy of cascade particle in lab frame:
    pin = sqrt(abs(tin*(tin + two*cmin)))
    pinx = pin*sintin*cosfin
    piny = pin*sintin*sinfin
    pinz = pin*costin
    ein = tin + cmin
!   Momentum and total energy of Fermi sea particle (or of pair of
!   nucleons for pion absorption) in lab rframe:
    pn = sqrt(abs(tn*(tn + two*cmn)))
    pnx = pn*sintn*cosfn
    pny = pn*sintn*sinfn
    pnz = pn*costn
    en = tn + cmn
!   Total energy in lab frame.
    fac = ein + en
!   Velocity of the center-of-momentum system of the 2 (or 3)
!   particles in the lab frame.
    if (fac < div0Lim .and. fac > -div0Lim) then
       fac = div0Lim
       write(sDCM%io%message,1000) "65-67"
       call sDCM%io%print(4, 3, sDCM%io%message)
    end if
    v(1) = (pinx + pnx)/fac
    v(2) = (piny + pny)/fac
    v(3) = (pinz + pnz)/fac
    v2 = v(1)**2 + v(2)**2 + v(3)**2
!   Lorentz transformation of total energy into C.M. frame.
    u = fac*sqrt(abs(one - v2))
!   KINETIC energy of projectile in target rest frame
    temp = two*cmn
    if (temp < div0Lim .and. temp > -div0Lim) then
       temp = div0Lim
       write(sDCM%io%message, 1000) "77"
       call sDCM%io%print(4, 3, sDCM%io%message)
    end if
    tin1 = (u**2 - (cmin + cmn)**2)/(temp)
    return

! ======================================================================
1000 format("Divide by zero error prevented in 'tinvu.f90', line(s) ", A)
! ======================================================================
  end subroutine tinvu

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


text:
      SUBROUTINE fixaray(pexp, qexp, mexp)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE Vname0
      USE Vname1
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(in) :: mexp
      REAL(rprec), INTENT(in) :: pexp, qexp
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: l, nn0, n, m
C-----------------------------------------------
!
!       This routine stores toroidal and poloidal mode number arrays
!
!       MPOL = NUMBER OF POLOIDAL MODES USED IN CURVE-FIT
!       NPHI = NUMBER OF EQUALLY-SPACED PHI PLANES PER FIELD
!              PERIOD IN WHICH THE CURVE-FITTING IS PERFORMED
!       IMPORTANT: NPHI MUST BE EVEN FOR NON-SYMMETRIC SYSTEMS
!
!       MPNT = NUMBER OF R AND Z MODES APPEARING IN FIT
!       NTHETA = NUMBER OF THETA POINTS IN EACH PHI PLANE
!       N2   = TOTAL NUMBER OF RESIDUALS IN FSQ (PER PHI PLANE)
!            == 2 (m=0 modes) + 2*mrho (rhoc, rhos) + ntheta
!       NFP  = NUMBER OF TOROIDAL FIELD PERIODS (IRRELEVANT)
!
      mrho = mpol-1
      nphi2 = nphi/2
      IF (nphi2 .eq. 0) nphi2 = 1
      n2 = 2*(mrho+1) + ntheta

      l = 0
      ntor = MAX0(1,nphi-1)

      IF (nphi .eq. 2) ntor = 2

      nn0 = 1 - (ntor + 1)/2
      DO n = 1, ntor
         nn(n) = (nn0 + (n - 1))*nfp
      END DO
      DO m = 1, mpol
         mm(m) = m - 1
         DO n = 1, ntor
            IF (mm(m).ne.0 .or. nn(n).ge.0) THEN
               l = l + 1
               m1(l) = mm(m)
               n1(l) = nn(n)
            ENDIF
         END DO
      END DO
      mpnt = l
      dnorm = 2._dp/ntheta

      DO m = 0, mpol
         dm1(m) = m
         IF (m .eq. 0) THEN
            xmpq(m,1) = zero
            xmpq(m,2) = zero
         ELSE
            xmpq(m,1) = dm1(m)**(pexp+qexp)
            xmpq(m,2) = dm1(m)**(pexp)
         ENDIF
      END DO

      t1m(0) = 0; t2m(0) = 0
      t1m(1) = one
      DO m = 2, mrho
        t1m(m) = (REAL(m-1, rprec)/m)**mexp
      END DO
      DO m = 1,mrho-1
        t2m(m) = (REAL(m+1, rprec)/m)**mexp
      END DO
!      t2m(mrho-1) = zero
      t2m(mrho) = zero

!!SPH031214: BETTER FIT IF AXIS IS FIXED (mimics polar m=1 coupling)
      t2m(0) = t1m(2)

      END SUBROUTINE fixaray

ext: f

max_stars_repo_name: ORNL-Fusion/DESCUR

max_issues_repo_name: ORNL-Fusion/DESCUR

max_forks_repo_name: ORNL-Fusion/DESCUR


text:
! ###########################################################
! #                                                         #
! #                    THE LIDORT FAMILY                    #
! #                                                         #
! #      (LInearized Discrete Ordinate Radiative Transfer)  #
! #        --           -            -        -        -    #
! #                                                         #
! ###########################################################

! ###########################################################
! #                                                         #
! #  Author :      Robert. J. D. Spurr                      #
! #                                                         #
! #  Address :     RT Solutions, Inc.                       #
! #                9 Channing Street                        #
! #                Cambridge, MA 02138, USA                 #
! #                                                         #
! #  Tel:          (617) 492 1183                           #
! #  Email :        rtsolutions@verizon.net                 #
! #                                                         #
! #  This Version :   3.7 F90                               #
! #  Release Date :   June 2014                             #
! #                                                         #
! #       NEW: THERMAL SUPPLEMENT INCLUDED     (3.2)        #
! #       NEW: OUTGOING SPHERICITY CORRECTION  (3.2)        #
! #       NEW: TOTAL COLUMN JACOBIANS          (3.3)        #
! #       VLIDORT COMPATIBILITY                (3.4)        #
! #       THREADED/OPTIMIZED F90 code          (3.5)        #
! #       EXTERNAL SS / NEW I/O STRUCTURES     (3.6)        #
! #                                                         #
! #       Surface-leaving, BRDF Albedo-scaling     (3.7)    #
! #       Taylor series, BBF Jacobians, ThreadSafe (3.7)    #
! #                                                         #
! ###########################################################

!    #####################################################
!    #                                                   #
!    #   This Version of LIDORT comes with a GNU-style   #
!    #   license. Please read the license carefully.     #
!    #                                                   #
!    #####################################################

      MODULE BRDF_Sup_Mod

!  BRDF standard i/o type definitions

      USE BRDF_Sup_Inputs_def
      USE BRDF_Sup_Outputs_def

!  BRDF main standard file

      USE BRDF_Sup_Masters_M

      END MODULE BRDF_Sup_Mod


text:
        FUNCTION PR_AWNM ( wnum, twnm, vwnm, pprb )
C************************************************************************
C* PR_AWNM                                                              *
C*                                                                      *
C* This function computes AWNM, the numeric weather code for prevailing/*
C* temporary/probability/vicinity weather, with probability 30 weather  *
C* taking precedence over vicinity weather, vicinity weather taking     *
C* precedence over temporary weather, and temporary weather taking      *
C* precedence over prevailing weather.                                  *
C*                                                                      *
C* REAL PR_AWNM  ( WNUM, TWNM, VWNM, PPRB )                             *
C*                                                                      *
C* Input parameters:                                                    *
C*	WNUM		REAL	Prevailing numeric weather code         *
C*	TWNM 		REAL    Tempo./probability numeric weather code *
C*      VWNM		REAL	Vicinity numeric weather code           *
C*	PPRB		REAL	Probability for fcst change indicator   *
C*                                                                      *
C* Output parameters:                                                   *
C*      PR_AWNM         REAL    All weather numeric weather code        *
C**                                                                     *
C* Log:                                                                 *
C* D. Kidwell/NCEP	 5/03                                           *
C************************************************************************
        INCLUDE   	'GEMPRM.PRM'
	INCLUDE	  	'ERMISS.FNC'
C------------------------------------------------------------------------
	PR_AWNM = PR_TPWN ( twnm, vwnm, pprb )
	IF ( ERMISS ( PR_AWNM ) )  PR_AWNM = wnum
C*
        RETURN
        END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      SUBROUTINE MFFTC6(C,FAC)
*
*   PURPOSE:
*       ELEMENTARY GENTLEMAN-SANDE RADIX 5 STEP APPLIED TO A VECTOR-OF
*       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
*       SEE REF.[1] FOR NOTATIONS.
*       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDM, WHICH CONTROLS
*       ITS OPERATION THROUGH THE MFFTPA COMMON
*
*   DUMMY ARGUMENTS :
*
*   C   ARRAY BEING FOURIER  TRANSFORMED
*   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      INTEGER NUSTEP
      DOUBLE COMPLEX C(0:NUSTEP-1,0:4),FAC(0:*)
      DOUBLE COMPLEX T1,T2,T3,T4,T5,F1,F2,F3,F4
      REAL *8 SIN72,RAD5D4,S36D72
      PARAMETER (
     $ SIN72 =  9.51056516295153572116439333E-1,
     $ RAD5D4 =  5.59016994374947424102293417E-1,
     $ S36D72 =  6.18033988749894848204586834E-1 )
*..  MU > 0
            MUF=LX
            DO 200 MU=MSTEP,MLIM,MSTEP
              F1=CONJG(FAC(MUF))
              F2=CONJG(FAC(MUF*2))
              F3=CONJG(FAC(MUF*3))
              F4=CONJG(FAC(MUF*4))
              DO 190 LAM=MU,MU+LLIM,LSTEP
                DO 180 IV=LAM,LAM+IVLIM,IVS
                DO 180 I=IV,IV+ILIM,IES
                  T1=C(I,1)+C(I,4)
                  T2=C(I,2)+C(I,3)
                  T3=(C(I,1)-C(I,4))*SIN72
                  T4=(C(I,2)-C(I,3))*SIN72
                  T5=T1+T2
                  T1=RAD5D4*(T1-T2)
                  T2=C(I,0)-0.25*T5
                  C(I,0)=C(I,0)+T5
                  T5=T2+T1
                  T2=T2-T1
                  T1=T3+S36D72*T4
                  T3=S36D72*T3-T4
                  C(I,1)=(T5-DCMPLX(-DIMAG(T1),DREAL(T1)))*F1
                  C(I,4)=(T5+DCMPLX(-DIMAG(T1),DREAL(T1)))*F4
                  C(I,2)=(T2-DCMPLX(-DIMAG(T3),DREAL(T3)))*F2
                  C(I,3)=(T2+DCMPLX(-DIMAG(T3),DREAL(T3)))*F3
  180           CONTINUE
  190         CONTINUE
              MUF=MUF+LX
200         CONTINUE
            DO 193 LAM=0,LLIM,LSTEP
              DO 183 IV=LAM,LAM+IVLIM,IVS
              DO 183 I=IV,IV+ILIM,IES
                T1=C(I,1)+C(I,4)
                T2=C(I,2)+C(I,3)
                T3=(C(I,1)-C(I,4))*SIN72
                T4=(C(I,2)-C(I,3))*SIN72
                T5=T1+T2
                T1=RAD5D4*(T1-T2)
                T2=C(I,0)-0.25*T5
                C(I,0)=C(I,0)+T5
                T5=T2+T1
                T2=T2-T1
                T1=T3+S36D72*T4
                T3=S36D72*T3-T4
                C(I,1)=T5-DCMPLX(-DIMAG(T1),DREAL(T1))
                C(I,4)=T5+DCMPLX(-DIMAG(T1),DREAL(T1))
                C(I,2)=T2-DCMPLX(-DIMAG(T3),DREAL(T3))
                C(I,3)=T2+DCMPLX(-DIMAG(T3),DREAL(T3))
183           CONTINUE
193         CONTINUE
      END

max_stars_repo_name: cea-trust-platform/trust-code

max_issues_repo_name: pledac/trust-code

max_forks_repo_name: pledac/trust-code


text:
module SWESetupModule
!------------------------------------------------------------------------------
! Lagrangian Particle / Panel Method - Spherical Model
!------------------------------------------------------------------------------
!
!> @author
!> Peter Bosler, Department of Mathematics, University of Michigan
!
!> @defgroup SWEVorticityAndDivergence SWE vorticity and divergence
!> Provides data structures and methods for defining test cases for the shallow water equations
!
!
! DESCRIPTION:
!> @file
!> Provides data structures and methods for defining test cases for the shallow water equations!
!
!------------------------------------------------------------------------------
use NumberKindsModule
use LoggerModule
use SphereGeomModule
use ParticlesModule
use PanelsModule
use SphereMeshModule

implicit none

private
public SWESetup
public New, Delete
public InitNonlinearSteadyStateTC2, SetNonlinearSteadyStateTC2OnMesh, TC2_NREAL, TC2_NINT
public InitGravityWave, SetGravityWaveOnMesh, GRAVITY_WAVE_NINT, GRAVITY_WAVE_NREAL

!
!----------------
! Types and module constants
!----------------
!
type SWESetup
	integer(kint), pointer :: integers(:) => null()	!	integers required by test case definitions
	real(kreal), pointer :: reals(:) => null() 		! 	real numbers required by test case definitions
end type

integer(kint), parameter :: TC2_NREAL = 3, &
							TC2_NINT = 0, &
							GRAVITY_WAVE_NINT = 0, &
							GRAVITY_WAVE_NREAL = 5

!
!----------------
! Logging
!----------------
!
logical(klog), save :: logInit = .FALSE.
type(Logger) :: log
character(len=28), save :: logKey = 'SWESetup'
integer(kint), parameter :: logLevel = DEBUG_LOGGING_LEVEL
character(len=128) :: logString
character(len=24) :: formatString
!
!----------------
! Interfaces
!----------------
!
interface New
	module procedure NewPrivate
end interface

interface Delete
	module procedure DeletePrivate
end interface

contains

!
!----------------
! Standard methods : Constructor / Destructor
!----------------
!
subroutine NewPrivate(self, nInt, nReal)
	type(SWESetup), intent(out) :: self
	integer(kint), intent(in) :: nInt, nReal

	if (.NOT. logInit) call InitLogger(log,procRank)

	if ( nInt > 0 ) then
		allocate(self%integers(nInt))
		self%integers = 0
	else
		nullify(self%integers)
	endif
	if ( nReal > 0 ) then
		allocate(self%reals(nReal))
		self%reals = 0.0_kreal
	else
		nullify(self%reals)
	endif
end subroutine

subroutine DeletePrivate(self)
	type(SWESetup), intent(inout) :: self
	if ( associated( self%integers)) deallocate(self%integers)
	if ( associated(self%reals)) deallocate(self%reals)
end subroutine

!
!----------------
! Public functions
!----------------
!
subroutine InitNonlinearSteadyStateTC2(self, u0, rotAlpha, h0)
	type(SWESetup), intent(inout) :: self
	real(kreal), intent(in) :: u0, rotAlpha, h0

	if ( rotAlpha /= 0.0_kreal ) then
		call LogMessage(log, WARNING_LOGGING_LEVEL,logKey,' WARNING : rotated coriolis not implemented in timestepping yet.')
	endif

	self%reals(1) = u0
	self%reals(2) = rotAlpha
	self%reals(3) = h0
end subroutine

subroutine InitGravityWave(self, lat0, lon0, h0, hp, beta)
! set parameters to define a Gaussian bump in the h field.
!	 center at (lon0, lat0)
!	 h0 = background height
!	 hp = perturbation height
!	 beta = Gaussian shape parameter
	type(SWESetup), intent(inout) :: self
	real(kreal), intent(in) :: lat0, lon0, h0, hp, beta
	self%reals(1) = lat0
	self%reals(2) = lon0
	self%reals(3) = h0
	self%reals(4) = hp
	self%reals(5) = beta
end subroutine

subroutine SetGravityWaveOnMesh(aMesh, gWave)
	type(SphereMesh), intent(inout) :: aMesh
	type(SWESetup), intent(in) :: gWave
	!
	type(Particles), pointer :: aParticles
	type(Panels), pointer :: aPanels
	integer(kint) :: j
	real(kreal) :: xyzCent(3), h0, hp, beta

	if ( aMesh%problemKind /= SWE_SOLVER ) then
		call LogMessage(log,ERROR_LOGGING_LEVEL,logKey,' problemKind must be SWE_SOLVER.')
		return
	endif

	aParticles => aMesh%particles
	aPanels => aMesh%panels

	xyzCent = EARTH_RADIUS * [ cos(gwave%reals(1))*cos(gwave%reals(2)), cos(gwave%reals(1))*sin(gwave%reals(2)), sin(gwave%reals(1)) ]
	h0 = gwave%reals(3)
	hp = gwave%reals(4)
	beta = gwave%reals(5)

	do j=1,aParticles%N
		aParticles%relVort(j) = 0.0_kreal
		aParticles%div(j) = 0.0_kreal
		aParticles%h(j) = GravityWaveHeight(aParticles%x(:,j), xyzCent, h0, hp, beta)
		aParticles%potVort(j) = 0.0_kreal
	enddo
	do j=1,aPanels%N
		if (.NOT. aPanels%hasChildren(j)) then
			aPanels%relVort(j) = 0.0_kreal
			aPanels%div(j) = 0.0_kreal
			aPanels%h(j) = GravityWaveHeight(aPanels%x(:,j),xyzCent, h0, hp, beta)
			aPanels%potVort(j) = 0.0_kreal
		endif
	enddo

end subroutine

subroutine SetNonlinearSteadyStateTC2OnMesh(aMesh, testCase2)
	type(SphereMesh), intent(inout) :: aMesh
	type(SWESetup), intent(in) :: testCase2
	!
	type(Particles), pointer :: aParticles
	type(Panels), pointer :: aPanels
	integer(kint) :: j
	real(kreal) :: u0, rotAlpha, h0, lat, lon, f

	if ( aMesh%problemKind /= SWE_SOLVER ) then
		call LogMessage(log,ERROR_LOGGING_LEVEL,logKey,' problemKind must be SWE_SOLVER.')
		return
	endif

	aParticles => aMesh%particles
	aPanels => aMesh%panels

	u0 = testCase2%reals(1)
	rotAlpha = testCase2%reals(2)
	h0 = testCase2%reals(3)

	do j=1, aParticles%N
!		lat = Latitude(aParticles%x0(:,j))
!		lon = Longitude(aParticles%x0(:,j))
!		f = 2.0_kreal*OMEGA*(-cos(lon)*cos(lat)*sin(rotAlpha) + sin(lat)*cos(rotAlpha))
!		aParticles%relVort(j) = 2.0_kreal*u0/EARTH_RADIUS * ( -cos(lon)*cos(lat)*sin(rotAlpha) + sin(lat)*cos(rotAlpha))
!		aParticles%div(j) = 0.0_kreal
!		aParticles%h(j) = GRAV*h0 - (EARTH_RADIUS*OMEGA*u0 + 0.5_kreal*u0*u0)*&
!			(-cos(lon)*cos(lat)*sin(rotAlpha) + sin(lat)*cos(rotAlpha))*&
!			(-cos(lon)*cos(lat)*sin(rotAlpha) + sin(lat)*cos(rotAlpha))
!		aParticles%potVort(j) = ( aParticles%relVort(j) + f ) / aParticles%h(j)
		aParticles%relVort(j) = 2.0_kreal*u0/EARTH_RADIUS * &
			( -sin(rotAlpha)*aParticles%x0(1,j)/EARTH_RADIUS + cos(rotAlpha)*aParticles%x0(3,j)/EARTH_RADIUS )
		aParticles%div(j) = 0.0_kreal
		aParticles%h(j) = GRAV*h0 - (EARTH_RADIUS*OMEGA*u0 + u0*u0/2.0_kreal) * &
			( -sin(rotAlpha)*aParticles%x0(1,j)/EARTH_RADIUS + cos(rotAlpha)*aParticles%x0(3,j)/EARTH_RADIUS)*&
			( -sin(rotAlpha)*aParticles%x0(1,j)/EARTH_RADIUS + cos(rotAlpha)*aParticles%x0(3,j)/EARTH_RADIUS)
		aParticles%potVort(j) = (aParticles%relVort(j) + 2.0_kreal*OMEGA*&
			( -sin(rotAlpha)*aParticles%x0(1,j)/EARTH_RADIUS + cos(rotAlpha)*aParticles%x0(3,j)/EARTH_RADIUS))/&
			aParticles%h(j)
	enddo
	do j=1, aPanels%N
		if (.NOT. aPanels%hasChildren(j) ) then
			aPanels%relVort(j) = 2.0_kreal*u0/EARTH_RADIUS * &
				( -sin(rotAlpha)*aPanels%x0(1,j)/EARTH_RADIUS + cos(rotAlpha)*aPanels%x0(3,j)/EARTH_RADIUS )
			aPanels%div(j) = 0.0_kreal
			aPanels%h(j) = GRAV*h0 - (EARTH_RADIUS*OMEGA*u0 + u0*u0/2.0_kreal) * &
				( -sin(rotAlpha)*aPanels%x0(1,j)/EARTH_RADIUS + cos(rotAlpha)*aPanels%x0(3,j)/EARTH_RADIUS)*&
				( -sin(rotAlpha)*aPanels%x0(1,j)/EARTH_RADIUS + cos(rotAlpha)*aPanels%x0(3,j)/EARTH_RADIUS)
			aPanels%potVort(j) = (aPanels%relVort(j) + 2.0_kreal*OMEGA*&
				( -sin(rotAlpha)*aPanels%x0(1,j)/EARTH_RADIUS + cos(rotAlpha)*aPanels%x0(3,j)/EARTH_RADIUS))/&
				aPanels%h(j)
		else
			aPanels%relVort(j) = 0.0_kreal
			aPanels%div(j) = 0.0_kreal
			aPanels%h(j) = 0.0_kreal
			aPanels%potVort(j) = 0.0_kreal
		endif
	enddo

end subroutine

!
!----------------
! Module methods : module- or type-specific private functions
!----------------
!
function GeneralGaussian(xyz,xyzCent, bb, maxVal)
	real(kreal) :: GeneralGaussian
	real(kreal), intent(in) :: xyz(3), xyzCent(3), bb, maxVal
	GeneralGaussian = maxVal * exp( -2.0_kreal*bb*bb/EARTH_RADIUS/EARTH_RADIUS *( sum((xyz-xyzCent)*(xyz-xyzCent))))
end function

function GravityWaveHeight(xyz, xyzCent, h0, hp, beta)
	real(kreal) :: GravityWaveHeight
	real(kreal), intent(in) :: xyz(3), xyzCent(3), h0, hp, beta
	GravityWaveHeight = h0 + GeneralGaussian(xyz,xyzCent,beta,hp)
end function

subroutine InitLogger(aLog,rank)
	type(Logger), intent(inout) :: aLog
	integer(kint), intent(in) :: rank
	write(logKey,'(A,A,I0.2,A)') trim(logKey),'_',rank,' : '
	call New(aLog,logLevel)
	logInit = .TRUE.
end subroutine

end module

max_stars_repo_name: pbosler/LPPM

max_issues_repo_name: pbosler/LPPM

max_forks_repo_name: pbosler/LPPM


      subroutine basin_ls_pest_output

      use output_ls_pesticide_module
      use plant_module
      use plant_data_module
      use time_module
      use basin_module
      use output_landscape_module
      use constituent_mass_module
      use hydrograph_module, only : sp_ob, sp_ob1, ob

      implicit none

      integer :: ipest
      integer :: ls
      integer :: iob
      real :: const

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine outputs HRU variables on daily, monthly and annual time steps

      iob = sp_ob1%hru

      !! print balance for each pesticide
      do ipest = 1, cs_db%num_pests
        bpestb_d%pest(ipest) = pestbz

          do ls = 1, sp_ob%hru
            bpestb_d%pest(ipest) = bpestb_d%pest(ipest) + hpestb_d(ls)%pest(ipest)
          end do

       bpestb_m%pest(ipest) = bpestb_m%pest(ipest) + bpestb_d%pest(ipest)

      !! daily print
        if (pco%day_print == "y" .and. pco%int_day_cur == pco%int_day) then
          if (pco%pest%d == "y") then
             write (2864,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", ob(iob)%name, cs_db%pests(ipest), bpestb_d%pest(ipest)   !! pesticide balance
             if (pco%csvout == "y") then
               write (2868,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", ob(iob)%name, cs_db%pests(ipest), bpestb_d%pest(ipest)
             end if
          end if
        end if
        !! zero daily output

        !! check end of month
        if (time%end_mo == 1) then
          bpestb_y%pest(ipest) = bpestb_y%pest(ipest) + bpestb_m%pest(ipest)
          const = float (ndays(time%mo + 1) - ndays(time%mo))
          bpestb_m%pest(ipest) = bpestb_m%pest(ipest) // const

          !! monthly print
           if (pco%pest%m == "y") then
             write (2865,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", ob(iob)%name, cs_db%pests(ipest), bpestb_m%pest(ipest)
               if (pco%csvout == "y") then
                 write (2869,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", ob(iob)%name, cs_db%pests(ipest), bpestb_m%pest(ipest)
               end if
           end if

          bpestb_m%pest(ipest) = pestbz
        end if

        !! check end of year
        if (time%end_yr == 1) then
          bpestb_a%pest(ipest) = bpestb_a%pest(ipest) + bpestb_y%pest(ipest)
          const = time%day_end_yr
          bpestb_y%pest(ipest) = bpestb_y%pest(ipest) // const

          !! yearly print
           if (time%end_yr == 1 .and. pco%pest%y == "y") then
             write (2866,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", ob(iob)%name, cs_db%pests(ipest), bpestb_y%pest(ipest)
               if (pco%csvout == "y") then
                 write (2870,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", ob(iob)%name, cs_db%pests(ipest), bpestb_y%pest(ipest)
               end if
           end if

        end if

!!!!! average annual print
         if (time%end_sim == 1 .and. pco%pest%a == "y") then
           bpestb_a%pest(ipest) = bpestb_a%pest(ipest) / time%yrs_prt
           bpestb_a%pest(ipest) = bpestb_a%pest(ipest) // time%days_prt
           write (2867,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", ob(iob)%name, cs_db%pests(ipest), bpestb_a%pest(ipest)
           if (pco%csvout == "y") then
             write (2871,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", ob(iob)%name, cs_db%pests(ipest), bpestb_a%pest(ipest)
           end if
           bpestb_a%pest(ipest) = pestbz
         end if

      end do    !pesticide loop
      return

100   format (4i6,2a,2x,2a,13e12.4)

      end subroutine basin_ls_pest_output


text:
      subroutine plasticctr(iflg)
!***********************************************************************
! Copyright 2011 Los Alamos National Security, LLC  All rights reserved
! Unless otherwise indicated,  this information has been authored by an
! employee or employees of the Los Alamos National Security, LLC (LANS),
! operator of the  Los  Alamos National  Laboratory  under Contract  No.
! DE-AC52-06NA25396  with  the U. S. Department  of  Energy.  The  U. S.
! Government   has   rights  to  use,  reproduce,  and  distribute  this
! information.  The  public may copy  and  use this  information without
! charge, provided that this  Notice and any statement of authorship are
! reproduced on all copies.  Neither  the  Government nor LANS makes any
! warranty,   express   or   implied,   or   assumes  any  liability  or
! responsibility for the use of this information.
!***********************************************************************
!
! Driver routine for plastic computations using FEHM algorithm
!
! Author : Sai Rapaka
!

      use comai
      use comsi

      implicit none

      integer iflg
      integer i,j,k,iterm,jterm
      logical null1
      character*10 macro1
      double precision, dimension(1:3, 1:3) :: tmp

      if(iflg.eq.initPlastic) then

         if(icnl.eq.0) then
           allocate(row(1:3,1:3,1:9))
           allocate(col(1:3,1:3,1:9))
           tmp(1,1) = 1; tmp(1,2) = 4; tmp(1,3) = 6;
           tmp(2,1) = 4; tmp(2,2) = 2; tmp(2,3) = 5;
           tmp(3,1) = 6; tmp(3,2) = 5; tmp(3,3) = 3;
           do i=1,3
             do j=1,3
               do k=1,9
                 iterm = floor((k-0.5)/3) + 1
                 jterm = k - (iterm-1)*3
                 row(i,j,k) = tmp(i,jterm)
                 col(i,j,k) = tmp(iterm,j)
!                 write(iout,*) i,j,k,row(i,j,k), col(i,j,k)
               enddo
             enddo
           enddo
         else
             allocate(row(1:2,1:2,1:4))
             allocate(col(1:2,1:2,1:4))
         endif

      elseif (iflg.eq.assemblePlastic) then
!        if(icnl.eq.0) then
!          call geneq_plastic_3D()
!        endif
      endif

      end subroutine plasticctr

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
module mod1
implicit none
save
integer, parameter :: precision=15, range=307
integer, parameter :: dp = selected_real_kind(precision, range)
real(dp), parameter :: gravity=9.81d0
real(dp), parameter :: pi=3.1415926539
real(dp), parameter :: RU=8314d0
real(dp), parameter :: zero=0._dp, one=1._dp
real(dp), parameter :: cd=1.1, rhob=1.225,tamb=300d0
real(dp), parameter :: mwair=28.96,surfrocket=pi/4
! assuminng a 1.1 drag coefficient and

real(dp):: cp,cv,g,rgas,mw,vol=one,dia,cf,id,od,length,rref,rhos,psipa,pref
real(dp):: db=zero,dt,tmax,Tflame
real(dp):: thrust=zero, area, r, n, surf,mdotgen,mdotout,edotgen,edotout,energy
real(dp):: mdotos=zero, edotos, texit, dsigng,pamb,p,t
real(dp):: mcham,echam,time=zero,propmass=zero,drag=zero,netthrust=zero
integer nsteps,i
real(dp):: accel=zero, vel=zero, altitude=zero, rocketmass=zero
real(dp), allocatable :: output(:,:)
real(dp) den ! air density
end module

subroutine propwt ! calculate weight of propellent
  use mod1
  implicit none
  propmass=pi/4*(od**2-id**2)*length*rhos
  rocketmass=0.15*propmass ! assume 85% propellant loading and 15% extra wt of rocket
end subroutine

subroutine burnrate
  use mod1
  implicit none
  r=rref*(p/pref)**n ! calculate burn rate
  db=db+r*dt  ! calculate incremental burn distance
end subroutine

subroutine calcsurf
  ! cylinder burning from id outward and from both ends along the length
  use mod1
  implicit none

  surf=pi*(id+2.0d0*db)*(length-2.0d0*db)+pi*(od**2.0d0-(id+2.0*db)**2.0d0)*0.5

  if(id+2d0*db.gt.od.or.db.gt.length/2d0) THEN
     surf=0d0  ! we hit the wall and burned out
     r=0  ! turn off burn rate so burn distance stops increasing
   endif

vol=vol+r*surf*dt ! increment the interior volume of the chamber a little
end subroutine

subroutine calmdotgen
  use mod1
  implicit none
  mdotgen=rhos*r*surf
  edotgen=mdotgen*cp*Tflame
end subroutine

subroutine massflow
   USE mod1
   implicit none
   REAL (8)::mdtx,engyx
   REAL (8)::tx,gx,rx,px,cpx,pcrit,facx,term1,term2,pratio,cstar,ax,hx
   REAL (8):: p1,p2

   mdotos=0.
   edotos=0.  ! initially set them to zero prior to running this loop

     p1=p
     p2=pamb
     ax=area
     IF(p1.GT.p2) THEN
        dsigng=1
        tx=t
        gx=g
        rx=rgas
        px=p
        cpx=cp
        hx=cp*t
        pratio=p1/p2
     else
        dsigng=-1
        tx=tamb
        gx=g
        rx=rgas
        px=pamb
        cpx=cp
        hx=cp*tamb
        pratio=p2/p1
    end if

    pcrit=(2./(gx+1.))**(gx/(gx-1.))
    IF((1./pratio).LT.pcrit) then
        ! choked flow
        cstar=sqrt((1./gx)*((gx+1.)/2.)**((gx+1.)/(gx-1.))*rx*tx)
        mdtx=px*ax/cstar
    else
        ! unchoked flow
      facx=pratio**((gx-1.)/gx)
      term1=SQRT(gx*rx*tx/facx)
      term2=SQRT((facx-1.)/(gx-1.))
      mdtx=SQRT(2.)*px/pratio/rx/tx*facx*term1*term2*ax
    end if
    engyx=mdtx*hx  ! reformulate based on enthalpy of the chamber
    mdotos=mdtx*dsigng ! exiting mass flow (could be negative "dsigng")
    edotos=engyx*dsigng ! exiting enthalpy
end subroutine

subroutine addmass
    use mod1
    implicit none
    mcham=mcham+(mdotgen-mdotos)*dt
    echam=echam+(edotgen-edotos)*dt
end subroutine

subroutine calct
    use mod1
    implicit none
    t=echam/mcham/cv
end subroutine

subroutine calcp
    use mod1
    implicit none
    p=mcham*rgas*t/vol
end subroutine

subroutine calcthrust
    use mod1
    implicit none
    thrust=(p-pamb)*area*cf ! correction to thrust (actual vs vacuum thrust)
    den=rhob*exp(-gravity*mwair*altitude/RU/tamb)
    drag=-cd*0.5*den*vel*abs(vel)*surfrocket

    netthrust=thrust+drag
end subroutine

subroutine height
  use mod1
  implicit none
  propmass=propmass-mdotgen*dt ! incremental change in propellant mass
  accel=netthrust/(propmass+rocketmass+mcham)-gravity
  vel=vel+accel*dt
  altitude=altitude+vel*dt
end subroutine

!!  Main program

function legacy_rocket( &
    dt_, &
    t_max_, &
    c_p_, &
    MW_, &
    temperature_, &
    pressure_, &
    T_flame_, &
    r_ref_, &
    n_, &
    id_, &
    od_, &
    length_, &
    rho_solid_, &
    dia_, &
    C_f_)
  !! this is a basic program of a single stage
  !! rocket motor flowing out of a nozzle, assuming
  !! a thrust coefficient and ignoring the complexities of
  !! what happens to thrust at low pressures, i.e. shock in the nozzle

use mod1
implicit none

real(dp), intent(in) :: dt_, t_max_
real(dp), intent(in) :: c_p_, MW_
real(dp), intent(in) :: temperature_, pressure_
real(dp), intent(in) :: T_flame_, r_ref_, n_
real(dp), intent(in) :: id_, od_, length_, rho_solid_
real(dp), intent(in) :: dia_, C_f_
real(dp), allocatable :: legacy_rocket(:,:)

dt   = dt_
tmax = t_max_
cp = c_p_
mw = MW_
t = temperature_
p = pressure_
Tflame = T_flame_
rref   = r_ref_
n      = n_
id     = id_
od     = od_
length = length_
rhos   = rho_solid_
dia = dia_
cf  = C_f_

!  propellent grain is a cylinder burning radially outward and axially inward from one end.
!  the other end is considered inhibited.
! outer diameter is inhibited because this is a cast propellent: it was poured
! into the tube/chamber and only the inner diameter burns when ignited.

  ! propellant burn rate information
  psipa=6894.76d0 ! pascals per psi (constant)
  pref=3000d0*psipa ! reference pressure (constant)

  nsteps=nint(tmax/dt) ! number of time steps

! preallocate an output file for simulation infomration
  allocate(output(0:nsteps,11))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11

!! now begin calculating and initializing
! gas variables
  rgas=ru/mw
  cv=cp-rgas
  g=cp/cv

  area=pi/4d0*dia**2.0d0 ! nozzle area

  pamb=101325d0 ! atmospheric pressure

!  calculate initial mass and energy in the chamber
  mcham=p*vol/rgas/t  ! use ideal gas law to determine mass in chamber
  echam=mcham*cv*t ! initial internal energy in chamber

  output(0,:) = [time,p,t,mdotos,thrust,drag,netthrust,vol,accel,vel,altitude]

  call propwt
  do i=1,nsteps
   call burnrate
   call calcsurf
   call calmdotgen  ! [mdot,engy,dsign]= massflow(p1,pamb,t1,tamb,cp,cp,rgas,rgas,g,g,area)
   call massflow
   call addmass
   call calct
   call calcp
   call calcthrust
   call height
   time=time+dt
   output(i,:)=[time,p,t,mdotos,thrust,drag,netthrust,vol,accel,vel,altitude]

  enddo

  legacy_rocket = output

end function legacy_rocket

max_stars_repo_name: jbrzensk/rocket-science

max_issues_repo_name: jbrzensk/rocket-science


text:
!------------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in           !
!  continuous development by various groups and is based on information        !
!  from these groups: Federal Government employees, contractors working        !
!  within a United States Government contract, and non-Federal sources         !
!  including research institutions.  These groups give the Government          !
!  permission to use, prepare derivative works of, and distribute copies       !
!  of their work in the CMAQ system to the public and to permit others         !
!  to do so.  The United States Environmental Protection Agency                !
!  therefore grants similar permission to use the CMAQ system software,        !
!  but users are requested to provide copies of derivative works or            !
!  products designed to operate in the CMAQ system to the United States        !
!  Government without restrictions as to use by others.  Software              !
!  that is used with the CMAQ system but distributed under the GNU             !
!  General Public License or the GNU Lesser General Public License is          !
!  subject to their copyright restrictions.                                    !
!------------------------------------------------------------------------------!

SUBROUTINE getpblht (c, r, wspd)

!-------------------------------------------------------------------------------
! Name:     Get PBL Height
! Purpose:  Calculates PBL height from bulk Richardson number.
! Notes:    This routine is used when PBL from input meteorology is 0.0.
!           Bulk Richardson number calculation taken from Byun's pblpkg.F
!           subroutine.  PBL height algorithm from Alapaty's VMM for MM5.
! Revised:  19 Sep 2001  Original version.  (T. Otte)
!           08 Jan 2002  Adjusted minimum PBL height to be height of lowest
!                        mid-layer rather than height of lowest full layer.
!                        (T. Otte)
!           12 Jul 2004  Added provision to specify PBL height if KTOP is 1.
!                        Changed local array allocation to occur only on
!                        initial call to subroutine to avoid memory
!                        fragmentation.  (T. Otte)
!           09 Apr 2007  Added IMPLICIT NONE.  Added explicit declarations for
!                        E_AERK, QSURF, and DTEMP.  Modified calculations of
!                        DTMPV, DTEMP, and RIBK to remove dependencies on
!                        function CVMGZ.  (T. Otte)
!           01 Sep 2011  Replaced statement functions with external routines.
!                        (T. Otte)
!           07 Sep 2011  Updated disclaimer.  (T. Otte)
!-------------------------------------------------------------------------------

  USE mcipparm
  USE const
  USE xvars

  IMPLICIT NONE

  REAL                         :: apbl
  INTEGER,       INTENT(IN)    :: c
  REAL                         :: densd
  REAL                         :: densw
  REAL                         :: dtemp
  REAL,          PARAMETER     :: dtempmax   =  5.0 ! max sfc_T - air_T
  REAL,          PARAMETER     :: dtempmin   = -3.0 ! min sfc_t - air_T
  REAL                         :: dtmpv
  REAL,          EXTERNAL      :: e_aerk
  REAL                         :: fract1
  REAL                         :: fract2
  INTEGER                      :: k
  INTEGER                      :: ktop
  REAL                         :: psurf
  REAL,          EXTERNAL      :: ptemp
  REAL                         :: qratio
  REAL                         :: qsurf
  INTEGER,       INTENT(IN)    :: r
  REAL                         :: rhumid
  REAL,    SAVE, ALLOCATABLE   :: ribk       ( : )
  REAL                         :: tsurf
  REAL                         :: vpress
  REAL                         :: vptemp
  REAL,    SAVE, ALLOCATABLE   :: vptempk    ( : )
  REAL                         :: vptsurf
  REAL                         :: vsat
  REAL,          EXTERNAL      :: vtemp
  REAL,          INTENT(IN)    :: wspd       ( : )
  REAL                         :: zpbl

!-------------------------------------------------------------------------------
! Allocate necessary arrays.
!-------------------------------------------------------------------------------

  IF ( .NOT. ALLOCATED ( ribk    ) ) ALLOCATE ( ribk    ( metlay ) )
  IF ( .NOT. ALLOCATED ( vptempk ) ) ALLOCATE ( vptempk ( metlay ) )

!-------------------------------------------------------------------------------
! Compute difference of virtual temperature.  Algorithm taken from Byun's
! pblpkg.F.
!-------------------------------------------------------------------------------

  ! Saturation vapor pressure over water [mb]

  ! Estimating water vapor at surface with TSURF = 0.5*(TEMPM(1)+TEMPG)
  ! assuming it has the same relative humidity as the layer 1 air,
  ! but at the reference temperature

  vsat    = e_aerk( xtempm(c,r,1) - stdtemp )
  vpress  = xdenswm(c,r,1) * rwvap * xtempm(c,r,1)

  rhumid  = vpress / vsat

  psurf   = xpresm(c,r,1) + grav*xdensam(c,r,1)*x3htm(c,r,1)  ! [Pa]
  tsurf   = 0.5 * ( xtempm(c,r,1) + xtempg(c,r) )
  vsat    = e_aerk( tsurf - stdtemp )                         ! [Pa]

  vpress  = vsat * rhumid
  densw   = vpress / ( rwvap * tsurf )
  densd   = xdensaf(c,r,0)

  qsurf   = densw / densd

  ! Virtual temperature of ground and air above ground

  vptsurf = vtemp( tsurf, qsurf )
  vptsurf = ptemp( vptsurf, psurf )

  qratio  = xdenswm(c,r,1) / ( xdensam(c,r,1) - xdenswm(c,r,1) )

  vptemp  = vtemp( xtempm(c,r,1), qratio )
  vptemp  = ptemp( vptemp, xpresm(c,r,1) )

  dtmpv   = vptsurf - vptemp
  IF ( dtmpv < dtempmin ) dtmpv = dtempmin
  IF ( dtmpv > dtempmax ) dtmpv = dtempmax
  IF ( dtmpv == 0.0 ) THEN
    dtmpv = 1.0e-10
  ENDIF

!-------------------------------------------------------------------------------
! Calculate bulk Richardson number.  Algorithm taken from Byun's pblpkg.F.
!-------------------------------------------------------------------------------

  ribloop: DO k = 1, metlay

    vptempk(k) = vtemp( xtempm(c,r,k), qratio )
    vptempk(k) = ptemp( vptempk(k), xpresm(c,r,k) )

    dtemp = - ( dtmpv + vptempk(1) - vptempk(k) )
    IF ( dtemp == 0.0 ) THEN
      dtemp = 1.0e-10
    ENDIF

    ribk(k) = grav * x3htm(c,r,k) * dtemp / ( tsurf * wspd(k)**2 )
    IF ( ribk(k) == 0.0 ) THEN
      ribk(k) = -1.0e-10
    ENDIF

    IF ( ribk(k) >= 0.25 ) THEN
      ktop = k
      EXIT ribloop
    ENDIF

  ENDDO ribloop

!-------------------------------------------------------------------------------
! Calculate PBL height.  Algorithm adapted from Alapaty's VMM for Blackadar
! in MM5.
!-------------------------------------------------------------------------------

  IF ( ktop /= 1 ) THEN
    fract1 = ( 0.25 - ribk(ktop-1) ) / ( ribk(ktop) - ribk(ktop-1) )
    fract2 = 1.0 - fract1
    apbl   = x3htm(c,r,ktop) * fract1
    zpbl   = apbl + ( x3htm(c,r,ktop) * fract2 )
  ELSE
    zpbl   = x3htm(c,r,ktop)
  ENDIF

  IF ( zpbl < x3htm(c,r,1) ) THEN
    xpbl(c,r) = x3htm(c,r,1)
  ELSE
    xpbl(c,r) = zpbl
  ENDIF

!-------------------------------------------------------------------------------
! Deallocate arrays.
!-------------------------------------------------------------------------------

! DEALLOCATE ( ribk    )  ! commented out to avoid memory fragmentation
! DEALLOCATE ( vptempk )  ! commented out to avoid memory fragmentation

END SUBROUTINE getpblht

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


text:
!! R755 forall-triplet-spec
!    is index-name = subscript : subscript [ : stride]

10 forall (ii = 1:10, jj = 1:10, kk=4:8)
11 end forall
forall (ii = 1:nn)
end forall
forall (ii = 1:nn:5)
end forall
forall (ii = 1:30:N, jj=K:5:3)
end forall

end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:58:24 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE COMPUTE_ATTENUATION_COEFFS__genmod
          INTERFACE
            SUBROUTINE COMPUTE_ATTENUATION_COEFFS(N,QREF,F0,F_MIN,F_MAX,&
     &TAU_EPSILON,TAU_SIGMA)
              INTEGER(KIND=4), INTENT(IN) :: N
              REAL(KIND=8), INTENT(IN) :: QREF
              REAL(KIND=8), INTENT(IN) :: F0
              REAL(KIND=8), INTENT(IN) :: F_MIN
              REAL(KIND=8), INTENT(IN) :: F_MAX
              REAL(KIND=8), INTENT(OUT) :: TAU_EPSILON(1:N)
              REAL(KIND=8), INTENT(OUT) :: TAU_SIGMA(1:N)
            END SUBROUTINE COMPUTE_ATTENUATION_COEFFS
          END INTERFACE
        END MODULE COMPUTE_ATTENUATION_COEFFS__genmod


text:
! <rlsrcfn.f90 - A component of the City-scale
!                 Chemistry Transport Model EPISDOE-CityChem>
!*****************************************************************************!
!*
!* EPISODE - An urban-scale air quality model
!* ==========================================
!* Copyright (C) 2018  NILU - Norwegian Institute for Air Research
!*                     Instituttveien 18
!*                     PO Box 100
!*                     NO-2027 Kjeller
!*                     Norway
!*
!*                     Contact persons: Gabriela Sousa Santos - gss@nilu.no
!*                                      Paul Hamer - pdh@nilu.no
!*
!* Unless explicitly acquired and licensed from Licensor under another license,
!* the contents of this file are subject to the Reciprocal Public License ("RPL")
!* Version 1.5, https://opensource.org/licenses/RPL-1.5 or subsequent versions as
!* allowed by the RPL, and You may not copy or use this file in either source code
!* or executable form, except in compliance with the terms and conditions of the RPL.
!*
!* All software distributed under the RPL is provided strictly on an "AS IS" basis,
!* WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, AND LICENSOR HEREBY
!* DISCLAIMS ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF
!* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT, OR NON-INFRINGEMENT.
!* See the RPL for specific language governing rights and limitations under the RPL.
!*
!* ==========================================
!* The dispersion model EPISODE (Grønskei et. al., 1993; Larssen et al., 1994;
!* Walker et al., 1992, 1999; Slørdal et al., 2003, 2008) is an Eulerian grid model
!* with embedded subgrid models for calculations of pollutant concentrations resulting
!* from different types of sources (area-, line- and point sources). EPISODE solves the
!* time dependent advection/-diffusion equation on a 3 dimensional grid.
!* Finite difference numerical methods are applied to integrate the solution forward in time.
!* It also includes extensions as the implementation of a simplified EMEP photochemistry
!* scheme for urban areas (Walker et al. 2004) and a scheme for Secondary Organic Aerosol
!* implemented by Håvard Slørdal
!*
!* Grønskei, K.E., Walker, S.E., Gram, F. (1993) Evaluation of a model for hourly spatial
!*    concentrations distributions. Atmos. Environ., 27B, 105-120.
!* Larssen, S., Grønskei, K.E., Gram, F., Hagen, L.O., Walker, S.E. (1994) Verification of
!*    urban scale time-dependent dispersion model with sub-grid elements in Oslo, Norway.
!*    In: Air poll. modelling and its appl. X. New York, Plenum Press.
!* Slørdal, L.H., McInnes, H., Krognes, T. (2008): The Air Quality Information System AirQUIS.
!*    Info. Techn. Environ. Eng., 1, 40-47, 2008.
!* Slørdal, L.H., Walker, S.-E., Solberg, S. (2003) The Urban Air Dispersion Model EPISODE
!*    applied in AirQUIS. Technical Description. NILU TR 12/2003. ISBN 82-425-1522-0.
!* Walker, S.E., Grønskei, K.E. (1992) Spredningsberegninger for on-line overvåking i Grenland.
!*    Programbeskrivelse og brukerveiledning. Lillestrøm,
!*    Norwegian Institute for Air Research (NILU OR 55/92).
!* Walker, S.E., Slørdal, L.H., Guerreiro, C., Gram, F., Grønskei, K.E. (1999) Air pollution
!*    exposure monitoring and estimation. Part II. Model evaluation and population exposure.
!*    J. Environ. Monit, 1, 321-326.
!* Walker, S.-E., Solberg, S., Denby, B. (2003) Development and implementation of a
!*    simplified EMEP photochemistry scheme for urban areas in EPISODE. NILU TR 13/2013.
!*    ISBN 82-425-1524-7
!*****************************************************************************!

      subroutine RLSRCFN

! *** The subroutine reads main line sources filenames.
! ----------------------------------------------------------------------------------
! Based on:
! Version Episode 5.5 (May29, 2012) prepared for BB-Stand-Alone
! Original source code of EPISODE by Sam-Erik Walker (NILU)
!
! Sam-Erik Walker
! Norwegian Institute for Air Research (NILU)
! Instituttveien 18 P.O. Box 100
! N-2027 Kjeller, NORWAY
! Tel: +47 63898000 Fax: +47 63898050
! E-mail: sam-erik.walker@nilu.no
!
! ----------------------------------------------------------------------------------
! ----------------------------------------------------------------------------------
! REVISION HISTORY
!
!           2016  M. Karl: error handling for missing line source file
!           2017  M. Karl: reading LSRCSCALE and LSRCCANYON
!
! ----------------------------------------------------------------------------------

      use mod_util
      use mod_main
      use mod_site
      use mod_time
      use mod_mete
      use mod_conc
      use mod_lsrc

      implicit none

! *** Local variables

      REAL QLRMINV
      REAL QLRDV(2)
      INTEGER I
      INTEGER IQL
      character(len=256) :: TXTSTR
      LOGICAL LEOF
      INTEGER IC

! *** QLRMINV - Line source minimum influence distance
! *** QLRDV   - Line source distance to receptor points
! *** I       - Index
! *** IQL     - Line source index
! *** TXTSTR  - Textstring
! *** LEOF    - If end of file then true else false

! *** Read line sources subgrid model add to result indicator:
      IF (MAINFE) THEN
          CALL GETDAT(MAINUN,TXTSTR,LEOF)
          IF (TXTSTR(1:1) .NE. '<') THEN
              READ (TXTSTR,*) LSRCAI       ! Default = 1
          ENDIF
          print *,'rlsrcfn: add line sources: ',LSRCAI
      ENDIF

      IF (LSRCAI .LT. 0) LSRCAI = 0
      IF (LSRCAI .GT. 1) LSRCAI = 1

      IF (MAINFE) WRITE (MESSUN,2000) LSRCAI

! *** Read line sources static data filename:
      IF (MAINFE) THEN
          CALL GETDAT(MAINUN,TXTSTR,LEOF)
          print *,'rlsrcfn: static line src file: ',TXTSTR
          IF (TXTSTR(1:1) .NE. '<') THEN
              READ (TXTSTR,1000) LSRSFN
          ELSE

! ***        Read minimum influence distance and distances to receptor
! ***        points on each side of each road. These data are not
! ***        coming from AirQUIS.

              IF (NQL .GT. 0) THEN
                  IF (.NOT. ALLOCATED(QLRMIN)) ALLOCATE (QLRMIN(NQL))
                  IF (.NOT. ALLOCATED(QLRD))   ALLOCATE (QLRD(2*NQL))
              ENDIF

              READ (TXTSTR(2:256),*) QLRMINV,QLRDV(1),QLRDV(2)
              DO IQL = 1,NQL
                  QLRMIN(IQL)     = QLRMINV
                  QLRD(2*IQL - 1) = QLRDV(1)
                  QLRD(2*IQL    ) = QLRDV(2)
              ENDDO

          ENDIF

!MSK      ELSE
!MSK! ***     If we get data from AirQUIS:
!MSK          TXTSTR(1:1) = '<'
!MSK          CALL GETFNV(TXTSTR,LSRSFN,LSRSFE,1,1,LSRSFV)
      ENDIF

!      IF (MESSFE) WRITE (MESSUN,2010)
!     &            (LSRSFN(I:I),I=1,INDEX(LSRSFN,' ') - 1)

      IF (.NOT. ALLOCATED(LSRVFN)) ALLOCATE(LSRVFN(NC))
      IF (.NOT. ALLOCATED(LSRVUN)) ALLOCATE(LSRVUN(NC))
      IF (.NOT. ALLOCATED(LSRVFE)) ALLOCATE(LSRVFE(NC))
      IF (.NOT. ALLOCATED(LSRVFV)) ALLOCATE(LSRVFV(NC))
      IF (.NOT. ALLOCATED(LSRVFM)) ALLOCATE(LSRVFM(NC))

! *** Read line sources variable data filename:
!MSK There is a variable emission file for each compound
      DO 110 IC = 1,NC

! *** Read line source filenames or values:
        IF (MAINFE) THEN
          CALL GETDAT(MAINUN,TXTSTR,LEOF)
          print *,'rlsrcfn: line sources: ',TXTSTR
!MSK start
          if(TXTSTR(1:1)==' ') then
             LSRVFE(IC) = .false.
             LSRVFV(IC) = MISS
             print *,'no line src file for ',ic
          else
            CALL GETFNV(TXTSTR,LSRVFN(IC),LSRVFE(IC),1,1,LSRVFV(IC))
          endif
!MSK end
        ELSE
!MSK start
          IF ( LSRVFN(IC)(1:1) .EQ. ' ') THEN
            TXTSTR(1:1) = '<'
          ELSE
            TXTSTR = LSRVFN(IC)
          ENDIF
!MSK end
          CALL GETFNV(TXTSTR,LSRVFN(IC),LSRVFE(IC),1,1,LSRVFV(IC))
        ENDIF

  110 CONTINUE

!          CALL GETDAT(MAINUN,TXTSTR,LEOF)
!          print *,'rlsrcfn: dynamic line src file: ',TXTSTR
!          IF (TXTSTR(1:1) .NE. '<') THEN
!              READ (TXTSTR,1000) LSRVFN
!          ENDIF
!      ELSE
!          TXTSTR(1:1) = '<'
!          CALL GETFNV(TXTSTR,LSRVFN,LSRVFE,1,1,LSRVFV)
!      ENDIF
!MSK end

!MSK start
! *** Read emission scaling factor for subgrid scale line source model:
      IF (MAINFE) THEN
          CALL GETDAT(MAINUN,TXTSTR,LEOF)
          IF (TXTSTR(1:1) .NE. '<') THEN
              READ (TXTSTR,*) LSRCSCALE
          ENDIF
          print *,'rlsrcfn:  emission scaling factor: ',LSRCSCALE
      ENDIF

      IF (MESSFE) WRITE (MESSUN,2040) LSRCSCALE
!MSK end

!MSK start
! *** Read street canyon option for subgrid scale line source model:
      IF (MAINFE) THEN
          CALL GETDAT(MAINUN,TXTSTR,LEOF)
          IF (TXTSTR(1:1) .NE. '<') THEN
              READ (TXTSTR,*) LSRCCANYON
          ENDIF
          print *,'rlsrcfn:  street canyon option: ',LSRCCANYON
      ENDIF

      IF (MESSFE) WRITE (MESSUN,2050) LSRCCANYON
!MSK end

! *** Read minimum windspeed for subgrid scale line source model:
      IF (MAINFE) THEN
          CALL GETDAT(MAINUN,TXTSTR,LEOF)
          IF (TXTSTR(1:1) .NE. '<') THEN
              READ (TXTSTR,*) LSRCFFMIN
          ENDIF
          print *,'rlsrcfn:  min wind speed: ',LSRCFFMIN
      ENDIF

      IF (MESSFE) WRITE (MESSUN,2030) LSRCFFMIN

      RETURN

 1000 format (A256)

 2000 format ('RLSRCFN: Line sources sub grid model add to',       &
             ' indicator = ',I3)
 2010 format ('RLSRCFN: Line sources static   data file = ',256A1)
 2020 format ('RLSRCFN: Line sources variable data file = ',256A1)
 2030 format ('RLSRCFN: Minimum windspeed for line source subgrid',       &
             ' model = ',F7.2,' m/s')
!MSK start
 2040 format ('RLSRCFN: Emission scaling for line source subgrid',       &
             ' model = ',F7.2)
 2050 format ('RLSRCFN: Line sources sub grid model street canyon',       &
             ' option = ',I3)
!MSK end

! *** End of subroutine RLSRCFN

      end subroutine rlsrcfn

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


text:
      subroutine movfin
c
c     ******************************************************************
c     *                                                                *
c     *   move to a finer mesh                                         *
c     *                                                                *
c     ******************************************************************
c
      use dims
      use dimsc
c
c     ******************************************************************
c
      use flo_var
      use mesh_var
      use solv_var
      use mg_var
c
c     ******************************************************************
c
      use flo_param
      use solv_param
      use geo_param
      use mg_param
c
c     ******************************************************************
c
      implicit none
c
c     ******************************************************************
c
c     local variables
c
c     ******************************************************************
c
      integer  :: i,j,k
      integer  :: nc,ile,ite
      integer  :: j2
c
c     ******************************************************************
c
      real     :: rqq
c
c     ******************************************************************
c
c     set the mesh dimensions to those of the fine mesh
c
      nx        = nx  +nx
      ny        = ny  +ny
      il        = nx  +1
      jl        = ny  +1
      ie        = nx  +2
      je        = ny  +2
      ib        = nx  +3
      jb        = ny  +3
c
c     set the limits of the profile
c
      ite       = .50000005*xte*(il -1)
      ile       = il/2  +1
      itl       = ile  -ite
      itu       = ile  +ite
c
c     set the coarse mesh dimensions
c
      iil       = nx/2  +1
      jjl       = ny/2  +1
      iie       = nx/2  +2
      jje       = ny/2  +2
      iib       = nx/2  +3
      jjb       = ny/2  +3
c
c     reset the pointers to those of the fine mesh
c
      nc        = (iib  +1)*(jjb  +1)
      k1        = k1  -nc
      kx        = kx  -2*nc
      kw        = kw  -4*nc
c
c     exchange variables between the coarse and fine meshes
c
      call shift (ww(kw),ww1(kw),wwr(kw),xx(kx),xxc(kx),volc(k1),
     .            rlvc(k1),revc(k1))
c
c     interpolate the corrections to the next finer mesh
c
      call addw  (ww(kw),ww1(kw))
c
c     calculate the pressure
c
      j2        = jl
      if (mode.eq.0) j2 = je
      do j=2,j2
      do i=1,ie
         rqq       = .5*(w(i,j,2)**2  +w(i,j,3)**2)/w(i,j,1)
         p(i,j)    = (gamma  -1.)*dim(w(i,j,4),rqq)
      end do
      end do

      if (mode.le.0) return
c
c     flag the wall and far field boundary points
c
      call bound
c
c     update the surface pressure
c
      call bcwall
c
c     option to update the far field
c
      if (fbc.gt.0.) call bcfar
c
c     set extended values of the flow variables
c     in the halo surrounding the mesh
c
      call halo

      return

      end

max_stars_repo_name: AlexT-L/RANS

max_issues_repo_name: AlexT-L/RANS

max_forks_repo_name: AlexT-L/RANS


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

      SUBROUTINE CRCSIZ (MAXNP, X, Y, NINT, N, XEND, YEND, XCEN, YCEN,
     &   THETA1, THETA2, TANG, AA, BB, ERR, TEST, XNOLD, YNOLD, NXKOLD,
     &   LINKEG, LISTEG, BMESUR, MLINK, NPNOLD, NPEOLD, NNXK, REMESH,
     &   REXMIN, REXMAX, REYMIN, REYMAX, IDIVIS, SIZMIN, EMAX, EMIN,
     &   GRAPH)
C***********************************************************************

C  SUBROUTINE CRCSIZ = GETS INTERVALS ON AN ARC LINE BASED ON ERROR
C                      SIZE

C***********************************************************************

      DIMENSION X (MAXNP), Y (MAXNP)

      DIMENSION XNOLD(NPNOLD), YNOLD(NPNOLD), NXKOLD(NNXK, NPEOLD)
      DIMENSION LINKEG(2, MLINK), LISTEG(4 * NPEOLD), BMESUR(NPNOLD)

      LOGICAL GRAPH, REMESH, TEST, ERR, SGRAPH, MOVED

      IF (GRAPH) THEN
         SGRAPH = .TRUE.
      ELSE
         SGRAPH = .FALSE.
      ENDIF

      ITERAT = 100
      EPS = .01

      DELANG = 0.
      ANGNOW = 0.
      INTNOW = 0
      IF (GRAPH) THEN
         CALL SYMBOL (1, X(1), Y(1), 'DIAMND')
         CALL PLTFLU
      ENDIF
  100 CONTINUE
      INTNOW = INTNOW + 1
      IF ( ((TANG .GT. 0.) .AND.
     &   (THETA1 + ANGNOW + (DELANG * 1.3) .GT. TANG))
     &   .OR. ((TANG .LT. 0.) .AND.
     &   (THETA1 + ANGNOW + (DELANG * 1.3) .LT. TANG))
     &   ) THEN

C  THE END OF THE LINE (OR CLOSE ENOUGH) HAS BEEN REACHED

C  IF WE ARE TESTING OR THE INTERVALS MATCH, THEN SIMPLY FINISH THE
C  LINE.

         IF ((TEST) .OR. (INTNOW .EQ. NINT)) THEN
            NINT = INTNOW
            N = NINT + 1
            X(N) = XEND
            Y(N) = YEND
            IF (GRAPH) THEN
               CALL SYMBOL (1, X(INTNOW), Y(INTNOW), 'DIAMND')
               CALL MPD2VC (1, X(INTNOW), Y(INTNOW),
     &            X(N), Y(N))
               CALL SYMBOL (1, X(N), Y(N), 'DIAMND')
               CALL PLTFLU
            ENDIF
         ELSE

C  OTHERWISE, MAKE SURE THE INTERVALS ARE ALRIGHT AND ADD THE EXTRA ONE

            EPS = .001
            IF (INTNOW + 1 .NE. NINT) THEN
               CALL MESAGE ('** PROBLEMS WITH INTNOW '//
     &            'IN PLINE **')
               ERR = .TRUE.
               GOTO 160
            ENDIF
            ANG = THETA1 + ANGNOW +
     &         ((TANG - (THETA1 + ANGNOW)) * .5)
            RADIUS = BB * EXP (AA * ANG)
            X (INTNOW + 1) = XCEN + COS (ANG) * RADIUS
            Y (INTNOW + 1) = YCEN + SIN (ANG) * RADIUS
            N = NINT + 1
            X(N) = XEND
            Y(N) = YEND
            IF (GRAPH) THEN
               CALL SYMBOL (1, X(INTNOW + 1), Y(INTNOW + 1),
     &            'DIAMND')
               CALL SYMBOL (1, X(N), Y(N),
     &            'DIAMND')
               CALL MPD2VC (1, X(INTNOW), Y(INTNOW),
     &            X(INTNOW+1), Y(INTNOW+1))
               CALL MPD2VC (1, X(INTNOW+1), Y(INTNOW+1),
     &            X(N), Y(N))
               CALL PLTFLU
            ENDIF
         ENDIF
      ELSE

C  NOT TO THE END YET

         CALL GETSIZ (XNOLD, YNOLD, NXKOLD, LINKEG, LISTEG,
     &      BMESUR, MLINK, NPNOLD, NPEOLD, NNXK, REMESH,
     &      REXMIN, REXMAX, REYMIN, REYMAX, IDIVIS, SIZMIN, EMAX,
     &      EMIN, X(INTNOW), Y(INTNOW), S1)
         DELANG = S1 / (BB * EXP (AA * (THETA1 + ANGNOW)))
         IF (TANG .LT. 0.) DELANG = - DELANG
         ANG1 = ANGNOW + DELANG
         RAD1 = BB * EXP (AA * (THETA1 + ANG1))
         XNEW1 = XCEN + COS (THETA1 + ANG1) * RAD1
         YNEW1 = YCEN + SIN (THETA1 + ANG1) * RAD1
         CALL GETSIZ (XNOLD, YNOLD, NXKOLD, LINKEG, LISTEG,
     &      BMESUR, MLINK, NPNOLD, NPEOLD, NNXK, REMESH,
     &      REXMIN, REXMAX, REYMIN, REYMAX, IDIVIS, SIZMIN, EMAX,
     &      EMIN, XNEW1, YNEW1, S2)
         DELANG = ((S1 + S2) * .5) / (BB * EXP (AA * (THETA1 + ANGNOW)))
         IF (TANG .LT. 0.) DELANG = - DELANG
         ANG1 = ANGNOW + DELANG
         RAD1 = BB * EXP (AA * (THETA1 + ANG1))
         XNEW1 = XCEN + COS (THETA1 + ANG1) * RAD1
         YNEW1 = YCEN + SIN (THETA1 + ANG1) * RAD1
         CALL GETSIZ (XNOLD, YNOLD, NXKOLD, LINKEG, LISTEG,
     &      BMESUR, MLINK, NPNOLD, NPEOLD, NNXK, REMESH,
     &      REXMIN, REXMAX, REYMIN, REYMAX, IDIVIS, SIZMIN, EMAX,
     &      EMIN, XNEW1, YNEW1, S3)
         SIZNOW = (((S1 + S2) * .5) + S3) * .5
         DELANG = SIZNOW / (BB * EXP (AA * (THETA1 + ANGNOW)))
         IF (TANG .LT. 0.) DELANG = - DELANG
         ANGNOW = ANGNOW + DELANG
         RADIUS = BB * EXP (AA * (THETA1 + ANGNOW))
         X (INTNOW + 1) = XCEN + COS (THETA1 + ANGNOW) * RADIUS
         Y (INTNOW + 1) = YCEN + SIN (THETA1 + ANGNOW) * RADIUS
         IF (GRAPH) THEN
            CALL SYMBOL (1, X(INTNOW + 1), Y(INTNOW + 1),
     &         'DIAMND')
            CALL MPD2VC (1, X(INTNOW), Y(INTNOW),
     &         X(INTNOW+1), Y(INTNOW+1))
            CALL PLTFLU
         ENDIF
         GOTO 100
      ENDIF

C  ERASE THE LINES FOR SMOOTHING IF NEEDED

      IF ((.NOT. SGRAPH) .AND. (GRAPH)) THEN
         DO 110 J = 2, NINT
            CALL LCOLOR ('BLACK')
            CALL MPD2VC (1, X(J), Y(J), X(J+1), Y(J+1))
            CALL MPD2VC (1, X(J), Y(J), X(J-1), Y(J-1))
            CALL SYMBOL (1, X(J), Y(J), 'DIAMND')
            CALL LCOLOR ('WHITE')
            CALL PLTFLU
  110    CONTINUE
      ENDIF

C  NOW SMOOTH THE NODES ALONG THE LINE

      DO 130 I = 1, ITERAT
         MOVED = .FALSE.
         ANGNOW = 0.
         DO 120 J = 2, NINT
            CALL GETSIZ (XNOLD, YNOLD, NXKOLD, LINKEG,
     &         LISTEG, BMESUR, MLINK, NPNOLD, NPEOLD, NNXK,
     &         REMESH, REXMIN, REXMAX, REYMIN, REYMAX,
     &         IDIVIS, SIZMIN, EMAX, EMIN, X(J-1), Y(J-1), SIZE1)
            CALL GETSIZ (XNOLD, YNOLD, NXKOLD, LINKEG,
     &         LISTEG, BMESUR, MLINK, NPNOLD, NPEOLD, NNXK,
     &         REMESH, REXMIN, REXMAX, REYMIN, REYMAX,
     &         IDIVIS, SIZMIN, EMAX, EMIN, X(J), Y(J), SIZE2)
            CALL GETSIZ (XNOLD, YNOLD, NXKOLD, LINKEG,
     &         LISTEG, BMESUR, MLINK, NPNOLD, NPEOLD, NNXK,
     &         REMESH, REXMIN, REXMAX, REYMIN, REYMAX,
     &         IDIVIS, SIZMIN, EMAX, EMIN, X(J+1), Y(J+1), SIZE3)
            DIST1 = SQRT ( ((X(J-1) - X(J)) **2) +
     &         ((Y(J-1) - Y(J)) **2) )
            DIST2 = SQRT ( ((X(J) - X(J+1)) **2) +
     &         ((Y(J) - Y(J+1)) **2) )
            DTOTAL = DIST1 + DIST2
            RATIO = DIST1 / DTOTAL
            DRATIO = ((SIZE1 + SIZE2) * .5) /
     &         ( ((SIZE1 + SIZE2) * .5) +
     &         ((SIZE2 + SIZE3) * .5) )
            TRATIO = (RATIO + DRATIO) * .5
            IF (SGRAPH) THEN
               CALL LCOLOR ('BLACK')
               CALL MPD2VC (1, X(J), Y(J), X(J+1), Y(J+1))
               CALL MPD2VC (1, X(J), Y(J), X(J-1), Y(J-1))
               CALL SYMBOL (1, X(J), Y(J), 'DIAMND')
               CALL LCOLOR ('WHITE')
               CALL PLTFLU
            ENDIF
            SIZNOW = DTOTAL * TRATIO
            DELANG = SIZNOW / (BB * EXP (AA * (THETA1 + ANGNOW)))
            IF (TANG .LT. 0.) DELANG = -DELANG
            ANGNOW = ANGNOW + DELANG
            ANG = THETA1 + ANGNOW
            RADIUS = BB * EXP (AA * ANG)
            X (J) = XCEN + COS (ANG) * RADIUS
            Y (J) = YCEN + SIN (ANG) * RADIUS
            IF (SGRAPH) THEN
               CALL MPD2VC (1, X(J), Y(J), X(J+1), Y(J+1))
               CALL MPD2VC (1, X(J), Y(J), X(J-1), Y(J-1))
               CALL SYMBOL (1, X(J), Y(J), 'DIAMND')
               CALL PLTFLU
            ENDIF
            DX1 = DIST1 / (.5 * (SIZE1 + SIZE2))
            DX2 = DIST2 / (.5 * (SIZE2 + SIZE3))
            IF (J .EQ. 2) THEN
               DXMAX = AMAX1 (DX1, DX2)
            ELSE
               DXMAX = AMAX1 (DXMAX, DX1, DX2)
            ENDIF
            DT = ABS((TRATIO * DTOTAL) - DIST1)
            IF (DT/DTOTAL .GT. EPS) MOVED = .TRUE.
  120    CONTINUE
         IF (.NOT. MOVED) GOTO 140
  130 CONTINUE
  140 CONTINUE
      IF ((.NOT. SGRAPH) .AND. (GRAPH)) THEN
         DO 150 J = 2, NINT
            CALL MPD2VC (1, X(J), Y(J), X(J+1), Y(J+1))
            CALL MPD2VC (1, X(J), Y(J), X(J-1), Y(J-1))
            CALL SYMBOL (1, X(J), Y(J), 'DIAMND')
            CALL PLTFLU
  150    CONTINUE
      ENDIF

  160 CONTINUE

      RETURN

      END

max_stars_repo_name: utkarshayachit/seacas

max_issues_repo_name: utkarshayachit/seacas

max_forks_repo_name: utkarshayachit/seacas