text:
      subroutine dcsrch(stp,f,g,ftol,gtol,xtol,task,stpmin,stpmax,
     +                  isave,dsave)
      character*(*) task
      integer isave(2)
      double precision f, g, stp, ftol, gtol, xtol, stpmin, stpmax
      double precision dsave(13)
c     **********
c
c     Subroutine dcsrch
c
c     This subroutine finds a step that satisfies a sufficient
c     decrease condition and a curvature condition.
c
c     Each call of the subroutine updates an interval with
c     endpoints stx and sty. The interval is initially chosen
c     so that it contains a minimizer of the modified function
c
c           psi(stp) = f(stp) - f(0) - ftol*stp*f'(0).
c
c     If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
c     interval is chosen so that it contains a minimizer of f.
c
c     The algorithm is designed to find a step that satisfies
c     the sufficient decrease condition
c
c           f(stp) <= f(0) + ftol*stp*f'(0),
c
c     and the curvature condition
c
c           abs(f'(stp)) <= gtol*abs(f'(0)).
c
c     If ftol is less than gtol and if, for example, the function
c     is bounded below, then there is always a step which satisfies
c     both conditions.
c
c     If no step can be found that satisfies both conditions, then
c     the algorithm stops with a warning. In this case stp only
c     satisfies the sufficient decrease condition.
c
c     A typical invocation of dcsrch has the following outline:
c
c     Evaluate the function at stp = 0.0d0; store in f.
c     Evaluate the gradient at stp = 0.0d0; store in g.
c     Choose a starting step stp.
c
c     task = 'START'
c  10 continue
c        call dcsrch(stp,f,g,ftol,gtol,xtol,task,stpmin,stpmax,
c    +               isave,dsave)
c        if (task .eq. 'FG') then
c           Evaluate the function and the gradient at stp
c           go to 10
c           end if
c
c     NOTE: The user must not alter work arrays between calls.
c
c     The subroutine statement is
c
c       subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
c                         task,isave,dsave)
c     where
c
c       stp is a double precision variable.
c         On entry stp is the current estimate of a satisfactory
c            step. On initial entry, a positive initial estimate
c            must be provided.
c         On exit stp is the current estimate of a satisfactory step
c            if task = 'FG'. If task = 'CONV' then stp satisfies
c            the sufficient decrease and curvature condition.
c
c       f is a double precision variable.
c         On initial entry f is the value of the function at 0.
c            On subsequent entries f is the value of the
c            function at stp.
c         On exit f is the value of the function at stp.
c
c       g is a double precision variable.
c         On initial entry g is the derivative of the function at 0.
c            On subsequent entries g is the derivative of the
c            function at stp.
c         On exit g is the derivative of the function at stp.
c
c       ftol is a double precision variable.
c         On entry ftol specifies a nonnegative tolerance for the
c            sufficient decrease condition.
c         On exit ftol is unchanged.
c
c       gtol is a double precision variable.
c         On entry gtol specifies a nonnegative tolerance for the
c            curvature condition.
c         On exit gtol is unchanged.
c
c       xtol is a double precision variable.
c         On entry xtol specifies a nonnegative relative tolerance
c            for an acceptable step. The subroutine exits with a
c            warning if the relative difference between sty and stx
c            is less than xtol.
c         On exit xtol is unchanged.
c
c       task is a character variable of length at least 60.
c         On initial entry task must be set to 'START'.
c         On exit task indicates the required action:
c
c            If task(1:2) = 'FG' then evaluate the function and
c            derivative at stp and call dcsrch again.
c
c            If task(1:4) = 'CONV' then the search is successful.
c
c            If task(1:4) = 'WARN' then the subroutine is not able
c            to satisfy the convergence conditions. The exit value of
c            stp contains the best point found during the search.
c
c            If task(1:5) = 'ERROR' then there is an error in the
c            input arguments.
c
c         On exit with convergence, a warning or an error, the
c            variable task contains additional information.
c
c       stpmin is a double precision variable.
c         On entry stpmin is a nonnegative lower bound for the step.
c         On exit stpmin is unchanged.
c
c       stpmax is a double precision variable.
c         On entry stpmax is a nonnegative upper bound for the step.
c         On exit stpmax is unchanged.
c
c       isave is an integer work array of dimension 2.
c
c       dsave is a double precision work array of dimension 13.
c
c     Subprograms called
c
c       MINPACK-2 ... dcstep
c
c     MINPACK-1 Project. June 1983.
c     Argonne National Laboratory.
c     Jorge J. More' and David J. Thuente.
c
c     MINPACK-2 Project. November 1993.
c     Argonne National Laboratory and University of Minnesota.
c     Brett M. Averick, Richard G. Carter, and Jorge J. More'.
c
c     **********
      double precision zero, p5, p66
      parameter (zero=0.0d0,p5=0.5d0,p66=0.66d0)
      double precision xtrapl, xtrapu
      parameter (xtrapl=1.1d0,xtrapu=4.0d0)

      logical brackt
      integer stage
      double precision finit, ftest, fm, fx, fxm, fy, fym, ginit, gtest,
     +                 gm, gx, gxm, gy, gym, stx, sty, stmin, stmax,
     +                 width, width1

      external dcstep

c     Initialization block.

      if (task(1:5) .eq. 'START') then

c        Check the input arguments for errors.

         if (stp .lt. stpmin) task = 'ERROR: STP .LT. STPMIN'
         if (stp .gt. stpmax) task = 'ERROR: STP .GT. STPMAX'
         if (g .ge. zero) task = 'ERROR: INITIAL G .GE. ZERO'
         if (ftol .lt. zero) task = 'ERROR: FTOL .LT. ZERO'
         if (gtol .lt. zero) task = 'ERROR: GTOL .LT. ZERO'
         if (xtol .lt. zero) task = 'ERROR: XTOL .LT. ZERO'
         if (stpmin .lt. zero) task = 'ERROR: STPMIN .LT. ZERO'
         if (stpmax .lt. stpmin) task = 'ERROR: STPMAX .LT. STPMIN'

c        Exit if there are errors on input.

         if (task(1:5) .eq. 'ERROR') return

c        Initialize local variables.

         brackt = .false.
         stage = 1
         finit = f
         ginit = g
         gtest = ftol*ginit
         width = stpmax - stpmin
         width1 = width/p5

c        The variables stx, fx, gx contain the values of the step,
c        function, and derivative at the best step.
c        The variables sty, fy, gy contain the value of the step,
c        function, and derivative at sty.
c        The variables stp, f, g contain the values of the step,
c        function, and derivative at stp.

         stx = zero
         fx = finit
         gx = ginit
         sty = zero
         fy = finit
         gy = ginit
         stmin = zero
         stmax = stp + xtrapu*stp
         task = 'FG'

         go to 10

      else

c        Restore local variables.

         if (isave(1) .eq. 1) then
            brackt = .true.
         else
            brackt = .false.
         end if
         stage = isave(2)
         ginit = dsave(1)
         gtest = dsave(2)
         gx = dsave(3)
         gy = dsave(4)
         finit = dsave(5)
         fx = dsave(6)
         fy = dsave(7)
         stx = dsave(8)
         sty = dsave(9)
         stmin = dsave(10)
         stmax = dsave(11)
         width = dsave(12)
         width1 = dsave(13)

      end if

c     If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
c     algorithm enters the second stage.

      ftest = finit + stp*gtest
      if (stage .eq. 1 .and. f .le. ftest .and. g .ge. zero) stage = 2

c     Test for warnings.

      if (brackt .and. (stp .le. stmin .or. stp .ge. stmax))
     +    task = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS'
      if (brackt .and. stmax-stmin .le. xtol*stmax)
     +    task = 'WARNING: XTOL TEST SATISFIED'
      if (stp .eq. stpmax .and. f .le. ftest .and. g .le. gtest)
     +    task = 'WARNING: STP = STPMAX'
      if (stp .eq. stpmin .and. (f .gt. ftest .or. g .ge. gtest))
     +    task = 'WARNING: STP = STPMIN'

c     Test for convergence.

      if (f .le. ftest .and. abs(g) .le. gtol*(-ginit))
     +    task = 'CONVERGENCE'

c     Test for termination.

      if (task(1:4) .eq. 'WARN' .or. task(1:4) .eq. 'CONV') go to 10

c     A modified function is used to predict the step during the
c     first stage if a lower function value has been obtained but
c     the decrease is not sufficient.

      if (stage .eq. 1 .and. f .le. fx .and. f .gt. ftest) then

c        Define the modified function and derivative values.

         fm = f - stp*gtest
         fxm = fx - stx*gtest
         fym = fy - sty*gtest
         gm = g - gtest
         gxm = gx - gtest
         gym = gy - gtest

c        Call dcstep to update stx, sty, and to compute the new step.

         call dcstep(stx,fxm,gxm,sty,fym,gym,stp,fm,gm,brackt,stmin,
     +               stmax)

c        Reset the function and derivative values for f.

         fx = fxm + stx*gtest
         fy = fym + sty*gtest
         gx = gxm + gtest
         gy = gym + gtest

      else

c       Call dcstep to update stx, sty, and to compute the new step.

         call dcstep(stx,fx,gx,sty,fy,gy,stp,f,g,brackt,stmin,stmax)

      end if

c     Decide if a bisection step is needed.

      if (brackt) then
         if (abs(sty-stx) .ge. p66*width1) stp = stx + p5*(sty-stx)
         width1 = width
         width = abs(sty-stx)
      end if

c     Set the minimum and maximum steps allowed for stp.

      if (brackt) then
         stmin = min(stx,sty)
         stmax = max(stx,sty)
      else
         stmin = stp + xtrapl*(stp-stx)
         stmax = stp + xtrapu*(stp-stx)
      end if

c     Force the step to be within the bounds stpmax and stpmin.

      stp = max(stp,stpmin)
      stp = min(stp,stpmax)

c     If further progress is not possible, let stp be the best
c     point obtained during the search.

      if (brackt .and. (stp .le. stmin .or. stp .ge. stmax) .or.
     +    (brackt .and. stmax-stmin .le. xtol*stmax)) stp = stx

c     Obtain another function and derivative.

      task = 'FG'

   10 continue

c     Save local variables.

      if (brackt) then
         isave(1) = 1
      else
         isave(1) = 0
      end if
      isave(2) = stage
      dsave(1) = ginit
      dsave(2) = gtest
      dsave(3) = gx
      dsave(4) = gy
      dsave(5) = finit
      dsave(6) = fx
      dsave(7) = fy
      dsave(8) = stx
      dsave(9) = sty
      dsave(10) = stmin
      dsave(11) = stmax
      dsave(12) = width
      dsave(13) = width1

      end

max_stars_repo_name: lesserwhirls/scipy-cwt

max_issues_repo_name: lesserwhirls/scipy-cwt

max_forks_repo_name: lesserwhirls/scipy-cwt


text:
The UC Davis Department of Theatre and Dance is a part of the College of Letters and Science and is housed in Wright Hall. It offers one undergraduate Majors major (dramatic art, with a focus in either theater or dance), a Minors minor, an M.F.A. program (geared towards those pursuing a career in theater) and an interdisciplinary Ph.D. in Performance Studies.

The department includes over 120 undergraduate students majoring in dramatic art, 4045 graduate students, as well as 510 visiting professional artists and teachers each quarter. This mix of college students, midcareer professionals, and seasoned theater practitioners and scholars provides a dynamic environment for learning, experimentation, and growth, and a diverse, supportive community for all of its members.

 Undergraduate Program
The major in dramatic art provides a liberal arts education and preprofessional training in theater. Students develop their abilities through a wide range of courses combining practical and theoretical skills that emphasize critical thought and independence in a context of challenging cooperative work.  From the first year, majors may participate in every aspect of bringing a work to the stage: acting, dance, directing, playwriting, choreography, stage management and design. Technical practice in the design and creation of sets, costumes, sound and lighting accompanies a firm grounding in dramatic literature, dance and theatre history, performance theory and criticism. This combination of theory and practice prepares students for a broad range of career opportunities.

It is strongly recommended that all UCD students take DRA 010 at some point in their college careers; this is the nonmajors variation of Drama 20, which is the majors required Intro to Dramatic Art course. It is taught by graduate students. A very basic overview of theatrical styles and history is given, along with basics of acting and speech, with emphasis on actually getting up and doing things. If you have to have class at 8 am, this is the one.
 Graduate Program
The M.F.A. program offers a collaborative curriculum of interdisciplinary, teamtaught seminars that cover a full spectrum of performance practices, including acting, directing, choreography and design. Students in any area may elect to extend their expertise through advanced study in another discipline.

For information regarding academic programs, visit http://theatredance.ucdavis.edu.
 Granada ArtistsInResidence Program
The departments Granada ArtistsinResidence program is unique in American university theater, bringing prominent artists — directors, playwrights, choreographers, or filmmakers — to Davis during two academic quarters each year to teach and create a work for public performance. A special opportunity for graduate and undergraduate students to work closely with major theater and dance practitioners, it is a unique crosscultural experience.http://theatredance.ucdavis.edu/proartists/granada_history.aspx

 Department Buildings
Undergraduate students in the dramatic art major or minor (or even simply taking drama classes) can often be found in the Green Room.

The building houses not only a full proscenium theatre with raked seating (also known as the Main Stage), but also features a fullyoperational set shop, costume shop, and prop room. While the set shop is on the buildings ground level, the costume shop and prop room are in the basement.

 20102011 Theatre & Dance Season
For information about auditions, upcoming performances, visiting artists, and other events, sign up for the http://www.mailermailer.com/x?oid25079i Theatre and Dance eNews.

For more info, visit:
http://theatredance.ucdavis.edu/season/current.aspx

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


subroutine integrate_state(lo, hi, &
                           state   , s_l1, s_l2, s_l3, s_h1, s_h2, s_h3, &
                           diag_eos, d_l1, d_l2, d_l3, d_h1, d_h2, d_h3, &
                           a, half_dt, min_iter, max_iter) &
                           bind(C, name="integrate_state")

!
!   Calculates the sources to be added later on.
!
!   Parameters
!   ----------
!   lo : double array (3)
!       The low corner of the current box.
!   hi : double array (3)
!       The high corner of the current box.
!   state_* : double arrays
!       The state vars
!   diag_eos* : double arrays
!       Temp and Ne
!   src_* : double arrays
!       The source terms to be added to state (iterative approx.)
!   double array (3)
!       The low corner of the entire domain
!   a : double
!       The current a
!   half_dt : double
!       time step size, in Mpc km^-1 s ~ 10^12 yr.
!
!   Returns
!   -------
!   state : double array (dims) @todo
!       The state vars
!

    use amrex_fort_module, only : rt => amrex_real
    use meth_params_module, only : NVAR, heat_cool_type

    implicit none

    integer         , intent(in   ) :: lo(3), hi(3)
    integer         , intent(in   ) :: s_l1, s_l2, s_l3, s_h1, s_h2, s_h3
    integer         , intent(in   ) :: d_l1, d_l2, d_l3, d_h1, d_h2, d_h3
    real(rt), intent(inout) ::    state(s_l1:s_h1, s_l2:s_h2,s_l3:s_h3, NVAR)
    real(rt), intent(inout) :: diag_eos(d_l1:d_h1, d_l2:d_h2,d_l3:d_h3, 2)
    real(rt), intent(in   ) ::  a, half_dt
    integer         , intent(inout) :: min_iter, max_iter

    if (heat_cool_type .eq. 1) then
        call integrate_state_hc(lo, hi, state   , s_l1, s_l2, s_l3, s_h1, s_h2, s_h3, &
                                        diag_eos, d_l1, d_l2, d_l3, d_h1, d_h2, d_h3, &
                                a, half_dt, min_iter, max_iter)
    else if (heat_cool_type .eq. 3) then
        call integrate_state_vode(lo, hi, state   , s_l1, s_l2, s_l3, s_h1, s_h2, s_h3, &
                                          diag_eos, d_l1, d_l2, d_l3, d_h1, d_h2, d_h3, &
                                  a, half_dt, min_iter, max_iter)

    end if

end subroutine integrate_state


module mod
  integer :: i
  integer :: x(4)
  real, dimension(2,3) :: a
  real, allocatable, dimension(:,:) :: b
contains
  subroutine foo
    integer :: k
    k = 1
    a(1,2) = a(1,2)+3
  end subroutine foo
end module mod


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

!>    @brief linear system solver [x]:=[A^-1]x[b], dense matrix (quadratic)
!>    @param[in] NIJK system size
!>    @param[in] NRHS number of right sides/solutions
!>    @param[in] A system matrix
!>    @param[in] B rigth side
!>    @param[out] X solution
!#>   @details
      SUBROUTINE dense_solve(nijk,nrhs,a,b,x)
        use mod_linfos
        IMPLICIT NONE
        INTEGER nijk, nrhs, isnull
        INTEGER, ALLOCATABLE :: pivots(:)
        DOUBLE PRECISION, ALLOCATABLE :: tmp(:), tmpm(:,:)
        INTEGER error
        DOUBLE PRECISION a(nijk,nijk), b(nijk,nrhs)
        DOUBLE PRECISION x(nijk,nrhs), enough, enough2

        IF ((8.0D0*dble(2*nijk+nijk*nijk))>2.0D9) THEN
          WRITE(*,'(A,A,F9.4,A)') &
            'Error: array to large for direct dense', &
            ' matrix solver (', 8.0D0*dble(2*nijk+nijk*nijk)/1024.0D0 &
            **3, ' GBytes).'
          STOP
        END IF

        ALLOCATE(pivots(nijk))
        ALLOCATE(tmp(nijk))
        ALLOCATE(tmpm(nijk,nijk))
        error = 0
        CALL set_ival(nijk,0,pivots)

!     right side to [X]
        CALL dcopy(nijk*nrhs,b,1,x,1)
!     copy matrix
        CALL dcopy(nijk*nijk,a,1,tmpm,1)

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

!     solve dense matrix
        CALL dgesv(nijk,nrhs,a,nijk,pivots,x,nijk,error)

        IF (error>0) THEN
          WRITE(*,*) ' dense PLU: ', error, 'errors were found.'
        ELSE IF (linfos(4)>=1) THEN
!aw        write(*,*) ' dense PLU finished.'
        END IF

        IF ((linfos(4)>=1) .AND. (nrhs==1)) THEN
!       MVP
!       [tmp]==[r]:=[A]x[x]
          CALL dgemv('N',nijk,nijk,1.0D0,tmpm,nijk,x,1,0.0D0,tmp,1)

          enough = 0.0D0
          isnull = 0
          enough2 = 0.0D0
          CALL s_ddot(nijk,b,b,enough)
          CALL test_zero(enough,1,isnull)
          IF (isnull/=1) THEN
            CALL s_ddot(nijk,tmp,tmp,enough2)
            enough2 = abs(enough-enough2)/(enough)
          END IF

!       [r]:=[w]-[A]x[x]
          CALL dscal(nijk,-1.0D0,tmp,1)
          CALL daxpy(nijk,1.0D0,b,1,tmp,1)

          enough = 0.0D0
          CALL s_ddot(nijk,tmp,tmp,enough)
          WRITE(*,'(2(a,e12.5))') '   PLU: nrm2(R)', dsqrt(enough), &
            ' rel=', enough2
        END IF

!     restore matrix (only needed for debug)
        CALL dcopy(nijk*nijk,tmpm,1,a,1)

        DEALLOCATE(pivots)
        DEALLOCATE(tmp)
        DEALLOCATE(tmpm)

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

        RETURN
      END


text:
Dr. Carl K. Winter is a Professors professor of Food Science with the Cooperative Extension.  He is the Director of the UCD Food Safe Program.

He has written a lot of parody songs about food safety, such as Dont Get Sicky Wit It.  His musical stylings were heard on NPR, including the classic Fifty Ways to Eat Your Oysters.  All of his songs are available streaming from the http://foodsafe.ucdavis.edu/ Food Safe website.

Publications
  Heidinger, J.C., C.K. Winter and J.S. Cullor.  2009.  Quantatitive microbial risk assessment for Staphylococcus aureus and Staphylococcus Enterotoxin A in raw milk.  J. Food Protection 72(8): 16411653.
  Winter, C.K., A.M. Fraser, J.B. Gleason, S.K. Hovey, S.M. McCurdy, and O.S. Snider.  2009.  Food safety education using music parodies.  Journal of Food Science Education 8: 6267.
  Katz, J.M. and C.K. Winter.  2009.  Comparison of pesticide exposure from consumption of domestic and imported fruits and vegetables.  Food and Chemical Toxicology 47: 335338.
  Winter, C.K.  2009.  Organic food production and its influence on  naturally occurring toxins.  In: Kirchmann, H. and L. Bergstrom, eds. Organic Crop Production:  Ambitions and Limitations, Springer, New York pp. 231239.
  Winter, C.K.  2008.  Singing the songs of science:  Using musical parodies for microbiology education.  Focus on Microbiology Education 15(1): 1011.


text:
       subroutine impreg
c
c     + + + PURPOSE + + +
c
c     SR IMPREG is to modify the arrangment of flow regimes when
c     soil deposition occurs in an impoundment. This subroutine
c     was added as a new subroutine to the original WEPP code.
c
c     Called from: SR IMPMAI
c     Author(s):  S. Dun
c     Reference in User Guide:
c
c     Version:
c     Date recoded: 01/12/99
c     Recoded by: S.Dun
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
c
      include 'pmximp.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
c
c     + + + ARGUMENT DEFINITIONS + + +
c
c     + + + COMMON BLOCKS + + +
c
        include 'cimflg.inc'
c        flags for if a structure exists in an impondment
      include 'cimitf.inc'
c     write: a,b,c,d,e,ha,ht,hlm,a0,a1,a2,l0,l1,l2,qinf,isize
c
      include 'cimpnd.inc'
c        generic variables needed by the impoundment component
c
      include 'cimsed.inc'
c     write: d0,di,d100,sg,ndiv,ctd,cdp,hmin,hset,co,vs,cot
        include 'cimqot.inc'
c        here h is called
c
        include 'cimreg.inc'
c        coefficients used for rearranging the flow regimes.
        include 'cimds.inc'
c        Variables describing the feature of drop spillway.
        include 'cimcv1.inc'
c        Variables describing culvert 1's feature.
           include 'cimcv2.inc'
c        Variables describing culvert 2's feature.
        include 'cimrf.inc'
c        Variables describing rock-fill check dam's feature.
        include 'cimff.inc'
c        Variables describing filterfence,straw bales or trash barriers.
        include 'cimpr.inc'
c        Variables describing perforated riser.
c
        include 'cimcln.inc'
c        Flag for automatically cleanning.
        include 'cimsre.inc'
c        Sediment stage at last rearrangement time.
c
c     + + + LOCAL VARIABLES + + +
        integer regflag
c
      real c1wd,c1hitn,c1arn,c1hcv
        real c2wd,c2hitn,c2arn,c2hcv
        real rfh,ffh,hd
      real qpr(100), hpr(100), qb, hp, hpdel, qs, apr(5), hs,
     1     y, ko, ab, qw, qoo        ,hb,as
c
        real hl(15,30)
      integer  is, itf, ir(15), itr,i,j
c
c     + + + LOCAL DEFINITIONS + + +
c
c        regflag     - flag for determining whether or not to rearrange the flow
c                                  regimes of the impondment. 0 for not, 1 for yes.
c
c        c1arn                 - adjusted area of culvert inlet after sediment stage heigter
c                                  than inlet stage
c        c1hitn                - adjusted height of culvert inlet
c        c1wd                - culvert inlet's width
c        c1hcv                - stage of the bottom of culvert 1's inlet
c
c        c2arn                 - adjusted area of culvert inlet after sediment stage heigter
c                                  than inlet stage
c        c2hitn                - adjusted height of culvert inlet
c        c2wd                - culvert inlet's width
c        c2hcv                - stage of the bottom of culvert 1's inlet
c
c        rfh         - Stage of the bottom of rock-fill check dam's inlet
c
c        ffh         - Inlet bottom stage of filterfence,straw bales or trash barriers
c
c
c     hl(i,j)     - limiting lowest stage for outflow function i, flow regime j
c     i, j        - both are dummy do loop variables
c     ir          - initial flow regime dummy variable
c     is          - dummy do loop variable
c     itr         - flow regime indicator
c
c     + + + SAVES + + +
c
c     + + + SUBROUTINES CALLED + + +
c
c     + + + DATA INITIALIZATIONS + + +
c
c     + + + END SPECIFICATIONS + + +
c
c     If sediment stage haven't changed since last rearrangement,
c         this SR will not be used.
        if (abs(hmin(ipond)-sstage).lt.0.001) return
c
         regflag=0
c
c     IF the impoundment is cleaning automatically at this day
c        or the sediment stage is higher than the original inlet's
c     bottom stage of Drop Spillway, then reassign the coefficients
c     for corresponding structures.
c
        if(fids(ipond).ne.0) then
        if(hmin(ipond).gt.dhrs(ipond).or.atcln.eq.1) then
        regflag=1
        do 50 i=1,3
         dnfr(i,ipond) = 1
c
          do 40 j = 1, 10
            daf(i,j,ipond) = 0.0
            dbf(i,j,ipond) = 0.0
            dcf(i,j,ipond) = 0.0
            ddf(i,j,ipond) = 0.0
            def(i,j,ipond) = 0.0
            dhaf(i,j,ipond) = 0.0
            dhlf(i,j,ipond) = 0.0
            dhtf(i,j,ipond) = 0.0
   40     continue
   50        continue
c
        dbf(3,1,ipond)=10.0**8
c
        dcf(1,1,ipond) = 1.0
        dcf(2,1,ipond) = 1.0
        dcf(3,1,ipond) = 1.0
c
        dhaf(3,1,ipond) = -200.0
c
c     If the sediment stage is higher than the inlet's top stage of
c     the riser, then drop spillway cease to be effective.
c
       if (atcln.eq.1) then
c
            dcf(1,2,ipond) = 1.5
            dcf(2,2,ipond) = 0.5
            dcf(3,2,ipond) = 0.5
c
            dhaf(1,2,ipond) = dhrs(ipond)
            dhaf(2,2,ipond) = dhrs(ipond)
c
           if(fids(ipond).eq.1) then
            daf(3,2,ipond) = dhblot(ipond) + 0.6 * ddiabl(ipond)
c
            dbf(1,2,ipond) = dcoefw(ipond) * 22.0 / 7.0 * ddiars(ipond)
            dbf(2,2,ipond) = dcoefo(ipond) * 22.0 / 7.0 * ddiars(ipond)
     1           ** 2.0 / 4.0 * (2*32.2)** 0.5
            dbf(3,2,ipond) = 22.0 / 7.0 * ddiabl(ipond) ** 2.0 / 4.0 *
     1          (2*32.2) ** 0.5/ (1.0+dke(ipond)+dkb(ipond)+dkc(ipond)
     1                  *(dlbl(ipond)+dhrh(ipond))) ** 0.5
c
            dhaf(3,2,ipond) = dhrs(ipond) - (dhrh(ipond)+dsbl(ipond)*
     1                  dlbl(ipond)-0.6*ddiabl(ipond))
c
          else if (fids(ipond).eq.2) then
            daf(3,2,ipond) =dhblot(ipond) + 0.6 * ddiabl(ipond)
c
            dbf(1,2,ipond) = dcoefw(ipond) * 2.0 * (dlenrs(ipond)+
     1                 dwidrs(ipond))
            dbf(2,2,ipond) = dcoefo(ipond)*dlenrs(ipond)*dwidrs(ipond)*
     1                 (2*32.2) ** 0.5
            dbf(3,2,ipond) = 22.0/7.0*ddiabl(ipond)** 2.0/4.0 *(2*32.2)
     1          ** 0.5/ (1.0 +dke(ipond) +dkb(ipond) + dkc(ipond)*
     1                  (dlbl(ipond)+dhrh(ipond))) ** 0.5
c
            dhaf(3,2,ipond) = dhrs(ipond)-(dhrh(ipond)+dsbl(ipond)*
     1                  dlbl(ipond)-0.6*ddiabl(ipond))
c
          else if (fids(ipond).eq.3) then
            daf(3,2,ipond) = dhblot(ipond) + 0.6 * dhitbl(ipond)
c
            dbf(1,2,ipond) = dcoefw(ipond) * 2.0 * (dlenrs(ipond)+
     1                 dwidrs(ipond))
            dbf(2,2,ipond) = dcoefo(ipond)*dlenrs(ipond) *dwidrs(ipond)
     1                 * (2*32.2) ** 0.5
            dbf(3,2,ipond) = dhitbl(ipond) * dwdbl(ipond) * (2*32.2) **
     1         0.5 / (1.0+dke(ipond)+dkb(ipond)+dkc(ipond)*(dlbl(ipond)
     1                 +dhrh(ipond))) ** 0.5
c
            dhaf(3,2,ipond) = dhrs(ipond) - (dhrh(ipond)+dsbl(ipond)*
     1                  dlbl(ipond)-0.6*dhitbl(ipond))
c
          end if
c
          do 30 i = 1, 3
            dhlf(i,2,ipond) = dhaf(i,2,ipond)
            dhtf(i,2,ipond) = dhrs(ipond)
                  if(dhlf(i,2,ipond).gt.dhtf(i,2,ipond))
     1                dhtf(i,2,ipond)=dhlf(i,2,ipond)
            dnfr(i,ipond) = 2
   30                continue
          endif
        endif
        endif
c
c     IF the impoundment is cleaning automatically at this day
c        or the sediment stage is higher than the original inlet's
c     bottom stage of culvert 1, then reassign the coefficients
c     for corresponding structures.
c
c
        if(fcv1(ipond).ne.0) then
        if(hmin(ipond).gt.c1h(ipond).or.atcln.eq.1) then
        regflag=1
        do 60 i=4,6
         dnfr(i,ipond) = 1
c
          do 70 j = 1, 10
            daf(i,j,ipond) = 0.0
            dbf(i,j,ipond) = 0.0
            dcf(i,j,ipond) = 0.0
            ddf(i,j,ipond) = 0.0
            def(i,j,ipond) = 0.0
            dhaf(i,j,ipond) = 0.0
            dhlf(i,j,ipond) = 0.0
            dhtf(i,j,ipond) = 0.0
   70     continue
   60        continue
c
        daf(5,1,ipond) = 10.0**8
c
        dbf(4,1,ipond) = 1.0
        dbf(5,1,ipond) = 1.0
        dbf(6,1,ipond)=10.0**8
c
        dcf(4,1,ipond) = 1.0
        dcf(5,1,ipond) = 10.0**8
        dcf(6,1,ipond) = 1.0
c
        ddf(5,1,ipond) = 1.0
c
        dhaf(6,1,ipond) = -200.0
c
c     If the sediment stage is higher than the inlet's top stage of
c     culvert 1, then culvert 1 cease to be effective.
c
       if (hmin(ipond).lt.(c1h(ipond)+c1hit(ipond)).or.
     1         atcln.eq.1) then
c
          dnfr(4,ipond) = 2
          dnfr(5,ipond) = 2
          dnfr(6,ipond) = 2
c
c             calculate the changed height and area of culvert's inlet.
c             the inlet shape we consider as following is rectangular.
                   c1wd=c1ar(ipond)/c1hit(ipond)
                   c1hitn=(c1hit(ipond)-(hmin(ipond)-c1h(ipond)))
                   c1arn=c1hitn*c1wd
                   c1hcv=hmin(ipond)
c
                   if(atcln.eq.1) then
                        c1arn=c1ar(ipond)
                        c1hitn=c1hit(ipond)
                        c1hcv=c1h(ipond)
                   endif
c
          daf(4,2,ipond) = c1arn * c1hitn ** 0.5 * float(c1ncv(ipond))
          daf(5,2,ipond) = c1arn * c1hitn ** 0.5 * float(c1ncv(ipond))
          daf(6,2,ipond) = c1hot(ipond) + 0.6 * c1hit(ipond)
c
          dbf(4,2,ipond) = c1hitn * c1kus(ipond)
          dbf(5,2,ipond) = c1hitn
          dbf(6,2,ipond) = c1ar(ipond) * (2.*32.2) ** 0.5
     1    / (1.0+c1ke(ipond)+c1kb(ipond)+c1kc(ipond)*c1l(ipond)) ** 0.5
     1        * float(c1ncv(ipond))
c
          dcf(4,2,ipond) = 1 / c1mus(ipond)
          dcf(5,2,ipond) = 0.5 * c1s(ipond) - c1ys(ipond)
          dcf(6,2,ipond) = 0.5
c
          ddf(5,2,ipond) = c1cs(ipond)
c
          dhaf(4,2,ipond) = c1hcv
          dhaf(5,2,ipond) = c1hcv
          dhaf(6,2,ipond) = dhaf(6,2,ipond)
c
          dhlf(4,2,ipond) = c1hcv
          dhlf(5,2,ipond) = (dhaf(5,2,ipond)-dcf(5,2,ipond)
     1                                    *dbf(5,2,ipond)) + 0.0001
          dhlf(6,2,ipond) = dhaf(6,2,ipond)
c
          dhtf(4,2,ipond) = c1hcv
          dhtf(5,2,ipond) = c1h(ipond) + c1hit(ipond) - 0.0001
                 if(dhlf(5,2,ipond).gt.dhtf(5,2,ipond))
     1           dhtf(5,2,ipond)=dhlf(5,2,ipond)
          dhtf(6,2,ipond) = c1hcv
                      if(dhlf(6,2,ipond).gt.dhtf(6,2,ipond))
     1           dhtf(6,2,ipond)=dhlf(6,2,ipond)

          endif
        endif
        endif
c
c     IF the impoundment is cleaning automatically at this day
c        or the sediment stage is higher than the original inlet's
c     bottom stage of culvert 2, then reassign the coefficients
c     for corresponding structures.
c
        if(fcv2(ipond).ne.0) then
        if(hmin(ipond).gt.c2h(ipond).or.atcln.eq.1) then
        regflag=1
        do 80 i=7,9
         dnfr(i,ipond) = 1
c
          do 90 j = 1, 10
            daf(i,j,ipond) = 0.0
            dbf(i,j,ipond) = 0.0
            dcf(i,j,ipond) = 0.0
            ddf(i,j,ipond) = 0.0
            def(i,j,ipond) = 0.0
            dhaf(i,j,ipond) = 0.0
            dhlf(i,j,ipond) = 0.0
            dhtf(i,j,ipond) = 0.0
   90     continue
   80        continue
c
        daf(8,1,ipond) = 10.0**8
c
        dbf(7,1,ipond) = 1.0
        dbf(8,1,ipond) = 1.0
        dbf(9,1,ipond)=10.0**8
c
        dcf(7,1,ipond) = 1.0
        dcf(8,1,ipond) = 10.0**8
        dcf(9,1,ipond) = 1.0
c
        ddf(8,1,ipond) = 1.0
c
        dhaf(9,1,ipond) = -200.0
c
c     If the sediment stage is higher than the inlet's top stage of
c     culvert 2, then culvert 2 cease to be effective.
c
       if (hmin(ipond).lt.(c2h(ipond)+c2hit(ipond)).or.
     1         atcln.eq.1) then
c
          dnfr(7,ipond) = 2
          dnfr(8,ipond) = 2
          dnfr(9,ipond) = 2
c
c          calculate the changed height and area of culvert's inlet.
c          the inlet shape we consider as following is rectangular.
                   c2wd=c2ar(ipond)/c2hit(ipond)
                   c2hitn=(c2hit(ipond)-(hmin(ipond)-c2h(ipond)))
                   c2arn=c2hitn*c2wd
                   c2hcv=hmin(ipond)
c
                   if(atcln.eq.1) then
                        c2arn=c2ar(ipond)
                        c2hitn=c2hit(ipond)
                        c2hcv=c2h(ipond)
                   endif
c
          daf(7,2,ipond) = c2arn * c2hitn ** 0.5 * float(c2ncv(ipond))
          daf(8,2,ipond) = c2arn * c2hitn ** 0.5 * float(c2ncv(ipond))
          daf(9,2,ipond) = c2hot(ipond) + 0.6 * c2hit(ipond)
c
          dbf(7,2,ipond) = c2hitn * c2kus(ipond)
          dbf(8,2,ipond) = c2hitn
          dbf(9,2,ipond) = c2ar(ipond) * (2.*32.2) ** 0.5
     1    / (1.0+c2ke(ipond)+c2kb(ipond)+c2kc(ipond)*c2l(ipond)) ** 0.5
     1        * float(c2ncv(ipond))
c
          dcf(7,2,ipond) = 1 / c2mus(ipond)
          dcf(8,2,ipond) = 0.5 * c2s(ipond) - c2ys(ipond)
          dcf(9,2,ipond) = 0.5
c
          ddf(8,2,ipond) = c2cs(ipond)
c
          dhaf(7,2,ipond) = c2hcv
          dhaf(8,2,ipond) = c2hcv
          dhaf(9,2,ipond) = dhaf(9,2,ipond)
c
          dhlf(7,2,ipond) = c2hcv
          dhlf(8,2,ipond) = (dhaf(8,2,ipond)-dcf(8,2,ipond)
     1                                    *dbf(8,2,ipond)) + 0.0001
          dhlf(9,2,ipond) = dhaf(9,2,ipond)
c
          dhtf(7,2,ipond) = c2hcv
          dhtf(8,2,ipond) = c2h(ipond) + c2hit(ipond) - 0.0001
                 if(dhlf(8,2,ipond).gt.dhtf(8,2,ipond))
     1           dhtf(8,2,ipond)=dhlf(8,2,ipond)
          dhtf(9,2,ipond) = c2hcv
                      if(dhlf(9,2,ipond).gt.dhtf(9,2,ipond))
     1           dhtf(9,2,ipond)=dhlf(9,2,ipond)

          endif
        endif
        endif
c
c     IF the impoundment is cleaning automatically at this day
c        or the sediment stage is higher than the original inlet's
c     bottom stage of the Rock-fill check dam, then reassign the
c     coefficients for corresponding structure.
c
        if(firf(ipond).ne.0) then
        if(hmin(ipond).gt.rhrf(ipond)*3.281.or.atcln.eq.1) then
        regflag=1
c
c        set coefficients for no flow regime
         dnfr(10,ipond) = 1
c
          do 100 j = 1, 10
            daf(10,j,ipond) = 0.0
            dbf(10,j,ipond) = 0.0
            dcf(10,j,ipond) = 0.0
            ddf(10,j,ipond) = 0.0
            def(10,j,ipond) = 0.0
            dhaf(10,j,ipond) = 0.0
            dhlf(10,j,ipond) = 0.0
            dhtf(10,j,ipond) = 0.0
  100     continue
c
        dbf(10,1,ipond) = 1.0
        dcf(10,1,ipond) = 1.0
c
c     If the sediment stage is higher than the inlet's top stage of
c     the rock-fill checkdam, only 2 flow regimes for this structure.
c        no flow and overtopping flow.
c
         if(hmin(ipond).ge.rhotrf(ipond)*3.281) then
                dnfr(10,ipond)=2
c
                dhaf(10,2,ipond)=hmin(ipond)
                dhtf(10,2,ipond)=dhaf(10,2,ipond)+0.001
                dhlf(10,2,ipond)=dhtf(10,2,ipond)
                dbf(10,2,ipond) = 1.0
                dcf(10,2,ipond) = 1.0
                ddf(10,2,ipond)=3.087*rwdrf(ipond)*3.281
                def(10,2,ipond)=dhaf(10,2,ipond)
          endif

c
       if (hmin(ipond).lt.(rhotrf(ipond))*3.281.or.atcln.eq.1) then
c
          dnfr(10,ipond) = 3
c
                   if(atcln.eq.1) then
                        rfh=rhrf(ipond)*3.281
                   else
                        rfh=hmin(ipond)
                   endif
c
                dhaf(10,2,ipond)=rfh
                dhtf(10,2,ipond)=dhaf(10,2,ipond)+0.001
                dhlf(10,2,ipond)=dhtf(10,2,ipond)
            daf(10,2,ipond) =rwdrf(ipond)*3.281**3.0
          dbf(10,2,ipond) = rlnrf(ipond)*rarf(ipond)*3.281
          dcf(10,2,ipond) = 1 / rbrf(ipond)
c
                dhaf(10,3,ipond)=rfh
                dhtf(10,3,ipond)=rhotrf(ipond)*3.281
                dhlf(10,3,ipond)=dhtf(10,3,ipond)
            daf(10,3,ipond) =daf(10,2,ipond)
          dbf(10,3,ipond) =dbf(10,2,ipond)
          dcf(10,3,ipond) =dcf(10,2,ipond)
                ddf(10,3,ipond)=3.087*rwdrf(ipond)*3.281
                def(10,3,ipond)=dhtf(10,3,ipond)
c
          endif
        endif
        endif
c
c     IF the impoundment is cleaning automatically at this day
c        or the sediment stage is higher than the original inlet's
c     bottom stage of filterfence or straw bales, then reassign
c     the coefficients for corresponding structure.
c
        if(fiff(ipond).ne.0) then
        if(hmin(ipond).gt.fhff(ipond).or.atcln.eq.1) then
        regflag=1
c
c        set coefficients for no flow regime
         dnfr(12,ipond) = 1
c
          do 110 j = 1, 10
            daf(12,j,ipond) = 0.0
            dbf(12,j,ipond) = 0.0
            dcf(12,j,ipond) = 0.0
            ddf(12,j,ipond) = 0.0
            def(12,j,ipond) = 0.0
            dhaf(12,j,ipond) = 0.0
            dhlf(12,j,ipond) = 0.0
            dhtf(12,j,ipond) = 0.0
  110     continue
c
c     If the sediment stage is higher than the inlet's top stage of
c     the filterfence or straw bales, only 2 flow regimes for this
c        structure.no flow and overtopping flow.
c
         if(hmin(ipond).ge.fhotff(ipond)) then
                dnfr(12,ipond)=2
c
                dhaf(12,2,ipond)=hmin(ipond)
                dhtf(12,2,ipond)=dhaf(12,2,ipond)+0.001
                dhlf(12,2,ipond)=dhtf(12,2,ipond)
                dbf(12,2,ipond)=3.087*fwdff(ipond)
          endif

c
       if (hmin(ipond).lt.(fhotff(ipond)).or.atcln.eq.1) then
c
          dnfr(12,ipond) = 3
c
                   if(atcln.eq.1) then
                        ffh=fhff(ipond)
                   else
                        ffh=hmin(ipond)
                   endif
c
                dhaf(12,2,ipond)=ffh
                dhtf(12,2,ipond)=dhaf(12,2,ipond)+0.001
                dhlf(12,2,ipond)=dhtf(12,2,ipond)
            daf(12,2,ipond) =fwdff(ipond)*fvsl(ipond)
c
                dhaf(12,3,ipond)=ffh
                dhtf(12,3,ipond)=fhotff(ipond)
                dhlf(12,3,ipond)=dhtf(12,3,ipond)
            daf(12,3,ipond) =daf(12,2,ipond)
                if(fiff(ipond).eq.1) then
c                choosing sharp crested weir for filterfence
                        dbf(12,3,ipond) =3.27*fwdff(ipond)
                dcf(12,3,ipond) =0.4/(fhotff(ipond)-fhff(ipond))
     1                                                *fwdff(ipond)
                else
c                choosing broad crested weir for straw bales
                        dbf(12,3,ipond)=3.087*fwdff(ipond)
                endif
                ddf(12,3,ipond)= fhotff(ipond)
c
          endif
        endif
        endif
c
c     IF the impoundment is cleaning automatically at this day
c        or the sediment stage is higher than the bottom stage of
c     the slots on a perforated riser, then reassign the
c     coefficients for corresponding structure.
c
        if(fipr(ipond).ne.0) then
        if(hmin(ipond).gt.phd(ipond).or.atcln.eq.1) then
        regflag=1
c
c        set coefficients for no flow regime
        do 120 i= 13,15
         dnfr(i,ipond) = 1
c
          do 125 j = 1, 10
            daf(i,j,ipond) = 0.0
            dbf(i,j,ipond) = 0.0
            dcf(i,j,ipond) = 0.0
            ddf(i,j,ipond) = 0.0
            def(i,j,ipond) = 0.0
            dhaf(i,j,ipond) = 0.0
            dhlf(i,j,ipond) = 0.0
            dhtf(i,j,ipond) = 0.0
  125      continue
  120   continue
c
        dbf(13,1,ipond) = 1.0
        dhaf(13,1,ipond) = -0.01
c
          dcf(15,1,ipond)= 1.0
c
c     If the sediment stage is higher than the inlet's top stage of
c     the perforated riser. This structure cease to be effective.
c
         if(hmin(ipond).le.phr(ipond).or.atcln.eq.1) then
                if(atcln.eq.1) then
                   hb=phb(ipond)
                   hs=phs(ipond)
                   hd=phd(ipond)
                   as=pas(ipond)
                else
                   hb=phb(ipond)+(hmin(ipond)-phd(ipond))
                   hd=hmin(ipond)
                   if(hmin(ipond).ge.(phd(ipond)+phs(ipond)-0.001)) then
                                hs=0.001
                                as=0.0
                        else
                                hs=phs(ipond)-(hmin(ipond)-phd(ipond))
                                as=pas(ipond)
                        endif
                endif
c
c         determining stage discharge relationship for unsubmerged
c         flow in riser
c
          ko = exp(-0.60721+0.329229*pdiab(ipond)/pdiar(ipond))
          ab = 22.0 / 7.0 * pdiab(ipond) ** 2 / 4.0
c
          i = 1
          hpdel = 0.05
          hp = hpdel
          y = -hb
c
  240     qb = pcb(ipond) * ab * (64.4*(hb+y)) ** 0.5
c
          if (i.ge.100) then
            i = 1
            hpdel = hpdel * 2.0
            hp = hpdel
            y = -hb
            go to 240
          end if
c
          if (hp.lt.hs) then
c
            if (y.le.0.0) then
              qs = 2.0 / 3.0 * (pcs(ipond)*pas(ipond)/hs) * 64.4
     1                         ** 0.5 * hp ** 1.5
            else
              qs = (pcs(ipond)*pas(ipond)/hs) * 64.4 ** 0.5 *
     1             (y*(hp-y)**0.5+(2./3.)*(hp-y)**1.5)
            end if
c
          else
c
            if (hp.le.(phr(ipond)-hd)) then
c
              if (y.le.0.0) then
                qs = (2./3.) * (pcs(ipond)*pas(ipond)/hs) * 64.4 **
     1               0.5 * (hp**1.5-(hp-hs)**1.5)
c
              else if (y.le.hs) then
c
                qs = (pcs(ipond)*pas(ipond)/hs) * 64.4 ** 0.5 *
     1              (y*(hp-y)**0.5+(2./3.)*((hp-y)**1.5-(hp-hs)**1.5))
              else
                qs = (pcs(ipond)*pas(ipond)) * (64.4*(hp-y)) ** 0.5
              end if
c
            else
              qw = pcoefw(ipond) * 22.0 / 7.0 * pdiar(ipond)
     1                         * (hp-(phr(ipond)-hd)) ** 1.5
              qoo = pcoefo(ipond) * 22.0 / 7.0 * pdiar(ipond)
     1             ** 2 / 4.0 * (hp-(phr(ipond)-hd))** 0.5
c
              if (y.le.0.0) then
                qs = (2./3.) * (pcs(ipond)*pas(ipond)/hs) * 64.4 **
     1              0.5 * (hp**1.5-(hp-hs)**1.5) + min(qw,qoo)
c
              else if (y.le.hs) then
c
                qs = (pcs(ipond)*pas(ipond)/hs) * 64.4 ** 0.5 *
     1                          (y*(hp-y)**0.5+(2./3.)*(
     1              (hp-y)**1.5-(hp-hs)**1.5)) + min(qw,qoo)
              else
                qs = (pcs(ipond)*pas(ipond)) * (64.4*(hp-y)) ** 0.5 +
     1              min(qw,qoo)
              end if
c
            end if
c
          end if
c
          if (qb.lt.qs) then
            y = y + 0.0001
c
            if (y.ge.hp) then
              qpr(i) = qb
              hpr(i) = hp
              i = i + 1
              hp = hp + hpdel
              go to 240
c
            else if (y.gt.(phr(ipond)-hd)) then
c
              i = i - 1
              go to 250
            end if
c
            go to 240
c
          else
c
            qpr(i) = qs
            hpr(i) = hp
            i = i + 1
            hp = hp + hpdel
            go to 240
          end if
c
c         run regression routine on hpr and qpr to determine function
c         coefficients
c
  250     call impris(i,hpr,qpr,apr)
c
          daf(13,2,ipond) = 1.0
          daf(14,2,ipond) = pcb(ipond) * ab * 64.4 ** 0.5
c
          dbf(13,2,ipond) = apr(1)
          dbf(15,2,ipond) = 22.0 / 7.0 * pdiabl(ipond) ** 2.0 / 4.0 *
     1        (2*32.2) ** 0.5 /(1.0+pke(ipond)+pkb(ipond)+pkc(ipond)
     1                *plbl(ipond)+ko)
c
          dcf(13,2,ipond) = apr(2)
          dcf(15,2,ipond) = 0.5
c
          dhaf(13,2,ipond) = hd
          dhaf(14,2,ipond) = phd(ipond)-phb(ipond)
          dhaf(15,2,ipond) = phr(ipond) - (phrh(ipond)+psbl(ipond)*
     1                plbl(ipond)-0.6*pdiabl(ipond))
c
          dhlf(13,2,ipond) = hd
          dhlf(14,2,ipond) = phd(ipond)-phb(ipond)
          dhlf(15,2,ipond) = dhaf(15,2,ipond)
c
          dhtf(13,2,ipond) = hd
          dhtf(14,2,ipond) = hd
          dhtf(15,2,ipond) = hd
c
c         setting outflow functions for flow through a submerged riser
c         (flow regime #3)
c
          daf(13,3,ipond) = 10000000.0
          dbf(13,3,ipond) = 1.0
          dcf(13,3,ipond) = 1.0
c
          dhtf(13,3,ipond) = hpr(i)+hd
c
          dnfr(13,ipond) = 3
          dnfr(14,ipond) = 2
          dnfr(15,ipond) = 2
c
          endif
        endif
        endif

        if(regflag.eq.0) return
c
        sstage=hmin(ipond)
c
c     Rearrange the flow regimes.
        itr = 1
c
          do 10 j=1,mxrgm
              ht(j,ipond)=0.0
              hlm(j,ipond)=0.0
                do 20 i=1,mxstc
              a(i,j,ipond)=0.0
              b(i,j,ipond)=0.0
              c(i,j,ipond)=0.0
              d(i,j,ipond)=0.0
              e(i,j,ipond)=0.0
              ha(i,j,ipond)=0.0
20            continue
10            continue
c
        do 260 i = 1, 15
          ir(i) = 1
  260   continue
c
        ht(itr,ipond) = 0.0
c
  280   itf = 0
c
        do 290 is = 1, 15
c
          if (ht(itr,ipond).ge.dhtf(is,ir(is),ipond)) then
              a(is,itr,ipond) = daf(is,ir(is),ipond)
            b(is,itr,ipond) = dbf(is,ir(is),ipond)
            c(is,itr,ipond) = dcf(is,ir(is),ipond)
            d(is,itr,ipond) = ddf(is,ir(is),ipond)
            e(is,itr,ipond) = def(is,ir(is),ipond)
            ha(is,itr,ipond) = dhaf(is,ir(is),ipond)
            hl(is,itr) = dhlf(is,ir(is),ipond)
c
            if (ir(is).eq.(dnfr(is,ipond)+1)) then
              a(is,itr,ipond) = daf(is,(ir(is)-1),ipond)
              b(is,itr,ipond) = dbf(is,(ir(is)-1),ipond)
              c(is,itr,ipond) = dcf(is,(ir(is)-1),ipond)
              d(is,itr,ipond) = ddf(is,(ir(is)-1),ipond)
              e(is,itr,ipond) = def(is,(ir(is)-1),ipond)
              ha(is,itr,ipond) = dhaf(is,(ir(is)-1),ipond)
              hl(is,itr) = dhlf(is,(ir(is)-1),ipond)
            else
              ir(is) = ir(is) + 1
              itf = 1
            end if
c
          else
c
c           continue
c
            a(is,itr,ipond) = daf(is,(ir(is)-1),ipond)
            b(is,itr,ipond) = dbf(is,(ir(is)-1),ipond)
            c(is,itr,ipond) = dcf(is,(ir(is)-1),ipond)
            d(is,itr,ipond) = ddf(is,(ir(is)-1),ipond)
            e(is,itr,ipond) = def(is,(ir(is)-1),ipond)
            ha(is,itr,ipond) = dhaf(is,(ir(is)-1),ipond)
            hl(is,itr) = dhlf(is,(ir(is)-1),ipond)
          end if
c
  290   continue
c
        if (itf.eq.1) then
          itr = itr + 1
          ht(itr,ipond) = ht(itr-1,ipond) + 0.005
        else
          ht(itr,ipond) = ht(itr,ipond) + 0.005
        end if
c
        if (ht(itr,ipond).ge.(hot(ipond)+5.0)) then
          ht(itr,ipond) = hot(ipond) + 200.0
          nt(ipond) = itr
          go to 300
        else
          go to 280
        end if
c
c     determining the stage below which one of the outflow functions
c     blows up
  300   do 320 itr = 1, nt(ipond)
          hlm(itr,ipond) = 0.0
c
          do 310 is = 1, 15
            if (hl(is,itr).ge.hlm(itr,ipond)) then
              hlm(itr,ipond) = hl(is,itr)
            end if
  310     continue
c
  320   continue
c
        atcln=0
c
        return
        end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


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 chkwrfhdr (fl, cdfid)

!-------------------------------------------------------------------------------
! Name:     Check WRF Header
! Purpose:  Check WRF header variables from one WRF output file against the
!           "base" WRF output file used for this MCIP run to ensure that the
!           WRF output files are from the same simulation.
! Notes:    This routine is not thorough, but it should be enough to spot-check
!           key variables that would indicate a different WRF simulation.
!           This routine assumes that FL (input argument) is already opened.
! Revised:  15 May 2008  Original version.  (T. Otte)
!           11 May 2009  Correct bug in checking surface layer scheme in
!                        subsequent files.  (T. Otte)
!           25 Sep 2009  Removed netCDF file opening to prevent condition with
!                        too many open files for long simulations.  Added
!                        check on urban physics option and surface analysis
!                        nudging options.  Changed code to allow for GRID_FDDA
!                        to be greater than 1.  Corrected bug in checking to
!                        ensure that data are from same simulation (i.e.,
!                        restarted rather than reinitialized).  Corrected bug
!                        in checking observation nudging coefficient for
!                        temperature.  (T. Otte)
!           12 Feb 2010  Removed unused variables CDFID, DX, DY, N_TIMES, and
!                        VARID, and removed unused format statements 9600 and
!                        9700.  Changed RTOL to 1.0e-4.  (T. Otte)
!           18 Mar 2010  Added CDFID as an input argument.  Removed dependency
!                        on module WRF_NETCDF.  (T. Otte)
!           23 Aug 2011  Changed name of module FILE to FILES to avoid conflict
!                        with F90 protected intrinsic.  Updated netCDF commands
!                        to F90, and improved error handling.  (T. Otte)
!           07 Sep 2011  Updated disclaimer.  (T. Otte)
!-------------------------------------------------------------------------------

  USE files
  USE metinfo
  USE netcdf

  IMPLICIT NONE

  INTEGER,            INTENT(IN)    :: cdfid
  CHARACTER(LEN=80)                 :: cval
  CHARACTER(LEN=8)                  :: cval8
  CHARACTER(LEN=19)                 :: date_init
  INTEGER                           :: dimid
  CHARACTER(LEN=1024), INTENT(IN)    :: fl
  CHARACTER(LEN=1024)                :: fl1
  INTEGER                           :: ival
  CHARACTER(LEN=16),  PARAMETER     :: pname     = 'CHKWRFHDR'
  INTEGER                           :: rcode
  REAL,               PARAMETER     :: rtol      = 1.0e-4
  REAL                              :: rval

!-------------------------------------------------------------------------------
! Error-handling section.
!-------------------------------------------------------------------------------

  CHARACTER(LEN=1024), PARAMETER :: f9000 = "(/, 1x, 70('*'), &
    & /, 1x, '*** SUBROUTINE: ', a, &
    & /, 1x, '***   WRF FILES DO NOT SEEM TO BE FROM SAME DOMAIN', &
    & /, 1x, '***   VARIABLE = ', a)"

  CHARACTER(LEN=1024), PARAMETER :: f9100 = "( &
    & /, 1x, '***   FIRST FILE = ', a, &
    & /, 1x, '***   VALUE IN FIRST FILE = ', i4, &
    & /, 1x, '***   NEW FILE = ', a, &
    & /, 1x, '***   VALUE IN NEW FILE = ', i4, &
    & /, 1x, 70('*'))"

  CHARACTER(LEN=1024), PARAMETER :: f9200 = "( &
    & /, 1x, '***   FIRST FILE = ', a, &
    & /, 1x, '***   VALUE IN FIRST FILE = ', f13.3, &
    & /, 1x, '***   NEW FILE = ', a, &
    & /, 1x, '***   VALUE IN NEW FILE = ', f13.3, &
    & /, 1x, 70('*'))"

  CHARACTER(LEN=1024), PARAMETER :: f9300 = "( &
    & /, 1x, '***   FIRST FILE = ', a, &
    & /, 1x, '***   VALUE IN FIRST FILE = ', a, &
    & /, 1x, '***   NEW FILE = ', a, &
    & /, 1x, '***   VALUE IN NEW FILE = ', a, &
    & /, 1x, 70('*'))"

  CHARACTER(LEN=1024), PARAMETER :: f9400 = "(/, 1x, 70('*'), &
    & /, 1x, '*** SUBROUTINE: ', a, &
    & /, 1x, '***   ERROR RETRIEVING VARIABLE FROM WRF FILE', &
    & /, 1x, '***   VARIABLE = ', a, &
    & /, 1x, '***   NCF: ', a, &
    & /, 1x, 70('*'))"

!-------------------------------------------------------------------------------
! Check NX, NY, and NZ.
!-------------------------------------------------------------------------------

  fl1 = file_mm(1)

  rcode = nf90_get_att (cdfid, nf90_global, 'WEST-EAST_GRID_DIMENSION', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_nx ) THEN
      WRITE (*,f9000) TRIM(pname), 'WEST-EAST_GRID_DIMENSION'
      WRITE (*,f9100) TRIM(fl1), met_nx, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'WEST-EAST_GRID_DIMENSION',  &
                    TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'SOUTH-NORTH_GRID_DIMENSION', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_ny ) THEN
      WRITE (*,f9000) TRIM(pname), 'SOUTH-NORTH_GRID_DIMENSION'
      WRITE (*,f9100) TRIM(fl1), met_ny, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'SOUTH-NORTH_GRID_DIMENSION',  &
                    TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'BOTTOM-TOP_GRID_DIMENSION', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival-1 /= met_nz ) THEN
      WRITE (*,f9000) TRIM(pname), 'BOTTOM-TOP_GRID_DIMENSION'
      WRITE (*,f9100) TRIM(fl1), met_nz, TRIM(fl), ival-1
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'BOTTOM-TOP_GRID_DIMENSION',  &
                    TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

!-------------------------------------------------------------------------------
! Check domain attributes.
!-------------------------------------------------------------------------------

  rcode = nf90_get_att (cdfid, nf90_global, 'TITLE', cval)
  IF ( rcode == nf90_noerr ) THEN
    cval8 = '        '
    IF ( cval(18:18) == "V" ) THEN
      cval8(1:2) = cval(18:19)
    ENDIF
    IF ( cval(20:20) == '.' ) THEN
      cval8(3:4) = cval(20:21)
    ENDIF
    IF ( cval(22:22) == '.' ) THEN
      cval8(5:6) = cval(22:23)
    ENDIF
    IF ( cval(24:24) == '.' ) THEN
      cval8(7:8) = cval(24:25)
    ENDIF
    IF ( cval8 /= met_release ) THEN
      WRITE (*,f9000) TRIM(pname), 'WRF VERSION'
      WRITE (*,f9300) TRIM(fl1), met_release, TRIM(fl), cval8
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'TITLE', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'DX', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_resoln ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'DX'
      WRITE (*,f9200) TRIM(fl1), met_resoln, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'DX', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'DY', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_resoln ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'DY'
      WRITE (*,f9200) TRIM(fl1), met_resoln, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'DY', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'MAP_PROJ', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_mapproj ) THEN
      WRITE (*,f9000) TRIM(pname), 'MAP_PROJ'
      WRITE (*,f9100) TRIM(fl1), met_mapproj, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'MAP_PROJ', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'STAND_LON', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_proj_clon ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'STAND_LON'
      WRITE (*,f9200) TRIM(fl1), met_proj_clon, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'STAND_LON', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'MOAD_CEN_LAT', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_proj_clat ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'MOAD_CEN_LAT'
      WRITE (*,f9200) TRIM(fl1), met_proj_clat, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'MOAD_CEN_LAT', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'CEN_LON', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_x_centd ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'CEN_LON'
      WRITE (*,f9200) TRIM(fl1), met_x_centd, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'CEN_LON', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'CEN_LAT', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_y_centd ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'CEN_LAT'
      WRITE (*,f9200) TRIM(fl1), met_y_centd, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'CEN_LAT', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'TRUELAT1', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_tru1 ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'TRUELAT1'
      WRITE (*,f9200) TRIM(fl1), met_tru1, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'TRUELAT1', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'TRUELAT2', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_tru2 ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'TRUELAT2'
      WRITE (*,f9200) TRIM(fl1), met_tru2, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'TRUELAT2', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

!-------------------------------------------------------------------------------
! Extract model run options.
!-------------------------------------------------------------------------------

  rcode = nf90_get_att (cdfid, nf90_global, 'MMINLU', cval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( cval(1:4) /= met_lu_src(1:4) ) THEN
      WRITE (*,f9000) TRIM(pname), 'MMINLU'
      WRITE (*,f9300) TRIM(fl1), met_lu_src(1:4), TRIM(fl), cval(1:4)
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'MMINLU', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'ISWATER', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_lu_water ) THEN
      WRITE (*,f9000) TRIM(pname), 'ISWATER'
      WRITE (*,f9100) TRIM(fl1), met_lu_water, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'ISWATER', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_inq_dimid (cdfid, 'soil_layers_stag', dimid)
  IF ( rcode /= nf90_noerr ) THEN
    WRITE (*,f9400) TRIM(pname), 'ID for soil_layers_stag',  &
                    TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF
  rcode = nf90_inquire_dimension (cdfid, dimid, len=ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_ns ) THEN
      WRITE (*,f9000) TRIM(pname), 'soil_layers_stag'
      WRITE (*,f9100) TRIM(fl1), met_ns, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'value for soil_layers_stag',  &
                    TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'RA_LW_PHYSICS', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_lw_rad ) THEN
      WRITE (*,f9000) TRIM(pname), 'RA_LW_PHYSICS'
      WRITE (*,f9100) TRIM(fl1), met_lw_rad, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'RA_LW_PHYSICS', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'RA_SW_PHYSICS', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_sw_rad ) THEN
      WRITE (*,f9000) TRIM(pname), 'RA_SW_PHYSICS'
      WRITE (*,f9100) TRIM(fl1), met_sw_rad, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'RA_SW_PHYSICS', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'CU_PHYSICS', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_cumulus ) THEN
      WRITE (*,f9000) TRIM(pname), 'CU_PHYSICS'
      WRITE (*,f9100) TRIM(fl1), met_cumulus, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'CU_PHYSICS', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'MP_PHYSICS', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_expl_moist ) THEN
      WRITE (*,f9000) TRIM(pname), 'MP_PHYSICS'
      WRITE (*,f9100) TRIM(fl1), met_expl_moist, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'MP_PHYSICS', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'BL_PBL_PHYSICS', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_pbl ) THEN
      WRITE (*,f9000) TRIM(pname), 'BL_PBL_PHYSICS'
      WRITE (*,f9100) TRIM(fl1), met_pbl, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'BL_PBL_PHYSICS', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'SF_SFCLAY_PHYSICS', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_sfc_lay ) THEN
      WRITE (*,f9000) TRIM(pname), 'SF_SFCLAY_PHYSICS'
      WRITE (*,f9100) TRIM(fl1), met_sfc_lay, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'SF_SFCLAY_PHYSICS', TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'SF_SURFACE_PHYSICS', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_soil_lsm ) THEN
      WRITE (*,f9000) TRIM(pname), 'SF_SURFACE_PHYSICS'
      WRITE (*,f9100) TRIM(fl1), met_soil_lsm, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'SF_SURFACE_PHYSICS',  &
                    TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

  IF ( TRIM(met_release) >= 'V3.1' ) THEN
    rcode = nf90_get_att (cdfid, nf90_global, 'SF_URBAN_PHYSICS', ival)
    IF ( rcode == nf90_noerr ) THEN
      IF ( ival /= met_urban_phys ) THEN
        WRITE (*,f9000) TRIM(pname), 'SF_URBAN_PHYSICS'
        WRITE (*,f9100) TRIM(fl1), met_urban_phys, TRIM(fl), ival
        CALL graceful_stop (pname)
      ENDIF
    ELSE
      WRITE (*,f9400) TRIM(pname), 'SF_URBAN_PHYSICS',  &
                      TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

!-------------------------------------------------------------------------------
! Extract WRF start date and time information.
!-------------------------------------------------------------------------------

  rcode = nf90_get_att (cdfid, nf90_global, 'SIMULATION_START_DATE', date_init)
  IF ( rcode == nf90_noerr ) THEN
    date_init(11:11) = "-"  ! change from "_" to "-" for consistency
    IF ( date_init /= met_startdate(1:19) ) THEN
      WRITE (*,f9000) TRIM(pname), 'SIMULATION_START_DATE'
      WRITE (*,f9300) TRIM(fl1), met_startdate(1:19), TRIM(fl),  date_init
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    WRITE (*,f9400) TRIM(pname), 'SIMULATION_START_DATE',  &
                    TRIM(nf90_strerror(rcode))
    CALL graceful_stop (pname)
  ENDIF

!-------------------------------------------------------------------------------
! Determine FDDA options.
!-------------------------------------------------------------------------------

  rcode = nf90_get_att (cdfid, nf90_global, 'GRID_FDDA', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_fdda_3dan ) THEN
      WRITE (*,f9000) TRIM(pname), 'GRID_FDDA'
      WRITE (*,f9100) TRIM(fl1), met_fdda_3dan, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( TRIM(met_release) >= 'V2.2' ) THEN
      WRITE (*,f9400) TRIM(pname), 'GRID_FDDA', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'GUV', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_fdda_gv3d ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'GUV'
      WRITE (*,f9200) TRIM(fl1), met_fdda_gv3d, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( ( TRIM(met_release) >= 'V2.2' ) .AND. ( met_fdda_3dan == 1 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'GUV', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ELSE IF ( ( TRIM(met_release) >= 'V3.1' ) .AND.  &
              ( met_fdda_3dan == 2 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'GUV', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'GT', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_fdda_gt3d ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'GT'
      WRITE (*,f9200) TRIM(fl1), met_fdda_gt3d, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( ( TRIM(met_release) >= 'V2.2' ) .AND. ( met_fdda_3dan == 1 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'GT', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ELSE IF ( ( TRIM(met_release) >= 'V3.1' ) .AND. ( met_fdda_3dan == 2 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'GT', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'GQ', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_fdda_gq3d ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'GQ'
      WRITE (*,f9200) TRIM(fl1), met_fdda_gq3d, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( ( TRIM(met_release) >= 'V2.2' ) .AND. ( met_fdda_3dan == 1 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'GQ', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'GRID_SFDDA', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_fdda_sfan ) THEN
      WRITE (*,f9000) TRIM(pname), 'GRID_SFDDA'
      WRITE (*,f9100) TRIM(fl1), met_fdda_sfan, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( TRIM(met_release) >= 'V3.1' ) THEN
      WRITE (*,f9400) TRIM(pname), 'GRID_SFDDA', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'GUV_SFC', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_fdda_gvsfc ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'GUV_SFC'
      WRITE (*,f9200) TRIM(fl1), met_fdda_gvsfc, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( ( TRIM(met_release) >= 'V3.1' ) .AND. ( met_fdda_sfan == 1 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'GUV_SFC', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'GT_SFC', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_fdda_gtsfc ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'GT_SFC'
      WRITE (*,f9200) TRIM(fl1), met_fdda_gtsfc, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( ( TRIM(met_release) >= 'V3.1' ) .AND. ( met_fdda_sfan == 1 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'GT_SFC', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'GQ_SFC', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_fdda_gqsfc ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'GQ_SFC'
      WRITE (*,f9200) TRIM(fl1), met_fdda_gqsfc, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( ( TRIM(met_release) >= 'V3.1' ) .AND. ( met_fdda_sfan == 1 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'GQ_SFC', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'OBS_NUDGE_OPT', ival)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ival /= met_fdda_obs ) THEN
      WRITE (*,f9000) TRIM(pname), 'OBS_NUDGE_OPT'
      WRITE (*,f9100) TRIM(fl1), met_fdda_obs, TRIM(fl), ival
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( TRIM(met_release) >= 'V2.2' ) THEN
      WRITE (*,f9400) TRIM(pname), 'OBS_NUDGE_OPT', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'OBS_COEF_WIND', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_fdda_giv ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'OBS_COEF_WIND'
      WRITE (*,f9200) TRIM(fl1), met_fdda_giv, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( ( TRIM(met_release) >= 'V2.2' ) .AND. ( met_fdda_obs == 1 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'OBS_COEF_WIND', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'OBS_COEF_TEMP', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_fdda_git ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'OBS_COEF_TEMP'
      WRITE (*,f9200) TRIM(fl1), met_fdda_git, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( ( TRIM(met_release) >= 'V2.2' ) .AND. ( met_fdda_obs == 1 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'OBS_COEF_TEMP', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

  rcode = nf90_get_att (cdfid, nf90_global, 'OBS_COEF_MOIS', rval)
  IF ( rcode == nf90_noerr ) THEN
    IF ( ABS( rval - met_fdda_giq ) > rtol ) THEN
      WRITE (*,f9000) TRIM(pname), 'OBS_COEF_MOIS'
      WRITE (*,f9200) TRIM(fl1), met_fdda_giq, TRIM(fl), rval
      CALL graceful_stop (pname)
    ENDIF
  ELSE
    IF ( ( TRIM(met_release) >= 'V2.2' ) .AND. ( met_fdda_obs == 1 ) ) THEN
      WRITE (*,f9400) TRIM(pname), 'OBS_COEF_MOIS', TRIM(nf90_strerror(rcode))
      CALL graceful_stop (pname)
    ENDIF
  ENDIF

END SUBROUTINE chkwrfhdr

max_stars_repo_name: ClimateIT/CMAQ

max_issues_repo_name: ClimateIT/CMAQ

max_forks_repo_name: ClimateIT/CMAQ


text:
!==============================================================================!
  subroutine Info_Mod_Bulk_Start()
!------------------------------------------------------------------------------!
!   Essentially creates a box in which iteration residuls will be written.     !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
! use pro_mod
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  integer :: n  ! inner iteration
!-----------------------------------[Locals]-----------------------------------!
  integer :: i
!==============================================================================!

  ! Create a frame
  do i = 1, L_LINE, L_LINE-1
    bulk_info % line_lead (i:i) = '#'
    bulk_info % lines(1)  (i:i) = '#'
    bulk_info % line_sep  (i:i) = '#'
    bulk_info % lines(2)  (i:i) = '#'
    bulk_info % lines(3)  (i:i) = '#'
    bulk_info % line_trail(i:i) = '#'
  end do

  do i = 2, L_LINE-1
    bulk_info % line_lead (i:i) = '='
    bulk_info % line_sep  (i:i) = '-'
    bulk_info % line_trail(i:i) = '-'
  end do

  ! Create separators
  bulk_info % line_lead(L_LINE/2+1:L_LINE/2+1) = '+'
  bulk_info % lines(1) (L_LINE/2+1:L_LINE/2+1) = '|'
  bulk_info % line_sep (L_LINE/2+1:L_LINE/2+1) = '+'

  bulk_info % line_sep  (27:27) = '+'
  bulk_info % lines(2)  (27:27) = '|'
  bulk_info % lines(3)  (27:27) = '|'
  bulk_info % line_trail(27:27) = '+'

  bulk_info % line_sep  (54:54) = '+'
  bulk_info % lines(2)  (54:54) = '|'
  bulk_info % lines(3)  (54:54) = '|'
  bulk_info % line_trail(54:54) = '+'

  end subroutine

max_stars_repo_name: mhadziabdic/T-FlowS

max_issues_repo_name: mhadziabdic/T-FlowS

max_forks_repo_name: mhadziabdic/T-FlowS


text:
      program kwn !Ni-7.5at%Al-8.5at%Cr at 600C
C...define parameters
      implicit double precision (a-h, o-z),integer(i-n)
C...parameters about growth, rate nucleation rate etc.
      double precision Gmgp,Q,DGM,DGV,GNV,DM,RC,REFF
      double precision nig,alg,crg,nigp,algp,crgp,nige,alge,crge
      double precision gnimu,galmu,gcrmu
      double precision rate,zeld,belta
C...number of new gamma prime and total gamma prime in 100^3 nm cell
      double precision ngpcell,tgpcell,tgpsum
C...iteration times
      integer id,jd,kd,ld,md
      integer ka,la
      integer idi
C...sum and average of effective diffusion distance
      double precision EDSUM,EDAVG
C...result of evolution N f NiRi R alau algp crau crgp
      double precision RES(0:1000000,1:8)   !!!change time interval
C...EVR array i.e. radius array in 100nm^3 cell
      double precision EVR(1,100000)        !!!change array size
C...EVO array i.e. X Y Z R sigmar(effective diffusion distance)
      double precision EVO(1:5,100000)      !!!change array size
C...define the X,Y,Z,DIST,COMB,COMS to calculate the xyz array
      double precision X,Y,Z,COMB,COMS
      double precision DIST(1,100000)       !!!change array size
C...define the Input and Output array of matlab
      double precision ainp(100000,3),bout(100000,1)    !!!notice column and row
C...critical nucleus compositions
      double precision cnal,cncr
C...define parameters about gfile
      parameter(nwg=800000,nwp=5000000)
      dimension iwsg(nwg),iwse(nwp)
      character*256 tcpath,tmppath
      integer iwsg,iwse
      character*32 gfile
      logical sg2err
C...define parameters about matlab interface
      integer*8,external::engOpen,engClose,mxCreateDoubleMatrix
      integer*8,external::mxGetPr,engPutVariable,engGetVariable
      integer*8,external::engEvalString
      integer*8 ep,T
      integer status

C...define constants
      QNI=288310D0                    !J/mol
      QAL=281590D0                    !J/mol
      QCR=281150D0                    !J/mol
      DNI=9.64D-22                    !M**2/s Ni in Ni at 873K
      DAl=8.31D-21                    !M**2/s Al in Ni at 873K
      DCR=4.39D-21                    !M**2/s Cr in Ni at 873K
      PI=3.141592654D0
      SFE=0.0285D0                    !J/M**2   surface energy  !!!change
      BOL=1.3806505D-23				!JOL/K
      DMS=3.0D27                      !density of nucleation site  !!!change
 	VM=6.59D-6			            !M**3/mol   molar volume of gamma prime
      Va=1.09D-29                     !M**3   volume per atom in the matrix
      alattice=3.4D-10                !M lattice constant of gamma prime
      Temp=873.15D0                   !K
C...define initial values
      tevo=0D0
      tgpcell=0D0   !!!change initial number of gamma prime in 100^3 nm cell
C...al and cr in full austenite at time 0  !!!change intial label and value
      RES(0,5)=0.075D0
      RES(0,7)=0.085D0
C...al and cr in gamma prime at time 0
      RES(0,6)=0D0
      RES(0,8)=0D0
C...composition in critical nucleus
      cnal=0.168D0
      cncr=0.081D0
C...N,f,sum(Ri),Ravg at time 0
      RES(0,1)=0D0   !!!change initial number density
      RES(0,2)=0D0
      RES(0,3)=0D0
      RES(0,4)=0D0   !!!change initial average size
C...define the initial value X Y Z R for calculating the DIST
      EVO(1,0)=0.0D0  !X
      EVO(2,0)=0.0D0  !Y
      EVO(3,0)=0.0D0  !Z
      EVO(4,0)=0.0D0  !R
      EVO(5,0)=0.0D0  !sigma r effective diffusion distance

C...comments
      write(*,1)
1     format('Ni-Al-Cr at 600C evolution')
C...open print file
      open(2,file='C:\Users\Public\Documents\Thermo-Calc\2021a\SDK\TQ\
     &Windows\fortran\KWN\kwnresult.txt')
      open(3,file='C:\Users\Public\Documents\Thermo-Calc\2021a\SDK\TQ\
     &Windows\fortran\KWN\kwnnr.txt')

C...open gfile
C...gfile
      gfile='NIALCR.GES5'
C...initiate the workspace
      tcpath=' '
      tmppath=' '
      call tqini3(tcpath,tmppath,nwg,nwp,iwsg,iwse)
C...read the thermodynamic data file which was created by using
C...the GES module inside the Thermo-Calc software package
      call tqrfil(gfile,iwsg,iwse)
C...open matlab interface and set working path
C...open matlab
      ep=engOpen('')
      if (ep == 0) then
          write(*,*) 'can''t start matlab engine'
          stop
      end if
C...set path
      if (engEvalString(ep," addpath('C:\Program Files\MATLAB\R2017b\
     &toolbox\mpt\mpt')") /= 0) then
          write(*,*) 'engEvalString failed'
          stop
      end if

C...evolution(time interval changes)
      do 10 i=1,1000000                        !!!change time interval 1000000s
          if (i.LE.1000000D0) then             !!!change time interval 1000000s
              tevo=1.0D0*i                     !!!change time interval 1s
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...get new nucleus density reff and etc.(subrountine)
C...get component indexes
              call tqgsci(indexni,'NI',iwsg,iwse)
              call tqgsci(indexal,'AL',iwsg,iwse)
              call tqgsci(indexcr,'CR',iwsg,iwse)
C...get phase indexes
              call tqgpi(indexg,'FCC_L12',iwsg,iwse)
              call tqgpi(indexgp,'FCC_L12#2',iwsg,iwse)
C...set T/K
              call tqssu('T','K',iwsg,iwse)
C...set condition of equilibrium
              call tqsetc('T',-1,-1,Temp,icont,iwsg,iwse)
              call tqsetc('N',-1,-1,1.00D0,iconn,iwsg,iwse)
              call tqsetc('P',-1,-1,101325.0D0,iconp,iwsg,iwse)
              call tqsetc('X',-1,indexal,RES(i-1,5),iconal,iwsg,iwse)   !!!change initial value
              call tqsetc('X',-1,indexcr,RES(i-1,7),iconcr,iwsg,iwse)   !!!change initial value
C...calculate equilibrium
              call tqce(' ',0,0,0.0D0, iwsg,iwse)
              if(sg2err(ierr)) print *, 'Calculation failed!'
C...get mole fraction of gamma prime
              call tqget1('X',indexgp,indexni,nigp,iwsg,iwse)
              call tqget1('X',indexgp,indexal,algp,iwsg,iwse)
              call tqget1('X',indexgp,indexcr,crgp,iwsg,iwse)
C...get mole fraction of gamma after precipitation
              call tqget1('X',indexg,indexni,nige,iwsg,iwse)
              call tqget1('X',indexg,indexal,alge,iwsg,iwse)
              call tqget1('X',indexg,indexcr,crge,iwsg,iwse)
C...get Gm of fcc_l12
              Gmgp=tqggm(indexgp,iwsg,iwse)
C...calculate deltaG* and Qeff
              call tqcsp(indexg,'ENTERED',1.0D0,iwsg,iwse)
              call tqcsp(indexgp,'SUSPENDED',0.0D0,iwsg,iwse)
C...calculate equilibrium
              call tqce(' ',0,0,0.0D0,iwsg,iwse)
              if(sg2err(ierr)) print *, 'Calculation failed!'
C...get mole fraction of elements in fcc_a1
              call tqget1('X',-1,indexni,nig,iwsg,iwse)
              call tqget1('X',-1,indexal,alg,iwsg,iwse)
              call tqget1('X',-1,indexcr,crg,iwsg,iwse)
C...get Qeff
              Q=QNI+alg*QAL+crg*QCR
C...get chemical potential of elements in fcc_a1
              call tqget1('MU',-1,indexni,gnimu,iwsg,iwse)
              call tqget1('MU',-1,indexal,galmu,iwsg,iwse)
              call tqget1('MU',-1,indexcr,gcrmu,iwsg,iwse)
C...calculate deltaGm deltaGv
              DGM=nigp*gnimu+algp*galmu+crgp*gcrmu-Gmgp
              DGV=DGM/VM
C...get deltaG*
              GNV=16*PI*(SFE**3)/DGV/DGV/3
C...get critical nucleation radium and reff
              RC=2*SFE/DGV
              REFF=RC+0.5*((BOL*Temp/PI/SFE)**0.5)
C...get nucleation rate
              zeld=Va*DGV*DGV/8/PI/((SFE*SFE*SFE*BOL*Temp)**0.5)     !zeldovich factor
C...rate of attachment of solute atoms
              belta=4*PI*RC*RC/alattice**4/
     &(1/DNI/nig+1/DAL/alg+1/DCR/crg)
C...nucleation rate
              DM=DMS*(1-RES(i-1,2))*zeld*belta*exp(-GNV/BOL/Temp)
C...rate of growth
              rate=8.314D0*Temp*((nigp-nige)**2/nige/DNI+(algp-alge)**2/
     &alge/DAL+(crgp-crge)**2/crge/DCR)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...get the EVR array in 100nm^3 cell
C...get the nulceation rate
              ngpcell=ANINT(1.0D-21*DM)   !!! new added number   !!!change time interval 1s
              tgpcell=tgpcell+ngpcell
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C...get the EVR array
              do 11 j=1,tgpcell
                  if (j .LE. tgpcell-ngpcell) then
                      EVR(1,j)=EVR(1,j)
                  else
                      EVR(1,j)=REFF
                  end if
11            continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...get the random coordination XYZ and R array (i.e. EVO array)
              do 12  ka=1,tgpcell
C...get the radius
                  EVO(4,ka)=EVR(1,ka)
C...generate random coordiantion
121               call random_number(X)
                  call random_number(Y)
                  call random_number(Z)
C...check whether every former point satisfies D/2 .GL. Rbigger & D/2 .LT. Rsmaller
                  do 122 la=1,ka
                      DIST(1,la)=sqrt((1D-7*X-EVO(1,la-1))**2
     &+(1D-7*Y-EVO(2,la-1))**2+(1D-7*Z-EVO(3,la-1))**2)
                      COMB=0.5D0*DIST(1,la)-EVO(4,la)
                      COMS=0.5D0*DIST(1,la)-EVO(4,la-1)
                      if(COMB .GT. 0.0D0 .AND. COMS .GT. 0.0D0) then
                          go to 122
                      else
                          go to 121
                      end if
122               continue
C...make sure the coordination and get the array ainp
                  EVO(1,ka)=1D-7*X
                  EVO(2,ka)=1D-7*Y
                  EVO(3,ka)=1D-7*Z
                  ainp(ka,1)=1D-7*X
                  ainp(ka,2)=1D-7*Y
                  ainp(ka,3)=1D-7*Z
12            continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...set the rest of ainp as 0
              do 13 l=tgpcell+1,100000   !!!change size of ainp array
                  ainp(l,1)=0.0D0
                  ainp(l,2)=0.0D0
                  ainp(l,3)=0.0D0
13            continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...calculate effective diffusion distance by voronoi in matlab
C...get the ainp array (x,y,z) in matlab from Fortran
              T=mxCreateDoubleMatrix(100000,3,0)     !!!change size of ainp array
              call mxCopyReal8ToPtr(ainp,mxGetPr(T),100000*3)   !!!change size of ainp array
              status=engPutVariable(ep,'T',T)
              if (status /= 0) then
                  write(*,*) 'engPutVariable failed'
                  stop
              end if
C...give the element, which is not equeal to 0 in array A, to array B
              if (engEvalString(ep,"A=1E9*T;for m=1:100000;if(A(m,1)>0);        !!!change size of ainp array
     &B(m,1)=A(m,1);end;if(A(m,2)>0);B(m,2)=A(m,2);end;
     &if(A(m,3)>0);B(m,3)=A(m,3);end;end") /= 0) then
                  write(*,*) 'engEvalString failed'
                  stop
              end if
C...draw the voronoi construction by B array
              if (engEvalString(ep,"Options.plot=1;v=[0 0 0;0 0 100;
     &0 100 0;0 100 100;100 0 0;100 0 100;100 100 0;100 100 100];
     &P=polytope(v);Options.pbound=P;
     &Pn=mpt_voronoi(B,Options)") /= 0) then
                  write(*,*) 'engEvalString failed'
                  stop
              end if
C...calculate effective diffusion distance
              if (engEvalString(ep,"for n=1:size(B,1);V=extreme(Pn(n));
     &for o=1:size(V,1);
     &D(1,o)=norm([B(n,1),B(n,2),B(n,3)]-[V(o,1),V(o,2),V(o,3)]);
     &end;E(n,1)=sum(D)/size(V,1);clear D;clear V;end") /= 0) then
                  write(*,*) 'engEvalString failed'
                  stop
              end if
C...calculate E and bout array
              if (engEvalString(ep,"for p=1:100000;if(p>size(B,1));    !!!change size of ainp array
     &E(p,1)=0;end;end") /= 0) then
                  write(*,*) 'engEvalString failed'
                  stop
              end if
C...get the result in Fortran from matlab
              status = engGetVariable(ep,"E")
              call mxCopyPtrToReal8(mxGetPr(status),bout,100000*1)     !!!change size of ainp array
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C...get the EVO(5,id) array effective diffusion distance
              do 14 id=1,tgpcell
                  EVO(5,id)=1.0D-9*bout(id,1)-EVO(4,id)
14            continue
C...calculate average of diffusion distance
              EDSUM=0.0D0
              do 141 idi=1,tgpcell
                  EDSUM=EDSUM+EVO(5,idi)
141           continue
              EDAVG=EDSUM/tgpcell
C...clear all
              if (engEvalString(ep,"clc;clear all;close all") /= 0) then
                  write(*,*) 'engEvalString failed'
                  stop
              end if
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...update the EVO by growth rate equation (EVO(4,jd))
              do 15 jd=1,tgpcell
                  EVO(4,jd)=EVO(4,jd)+
     &1.0D0*(DGM-2*SFE*VM/EVO(4,jd))/rate/EVO(5,jd)
C...if radius is less than 0, set EVO(4,jd) as 0
                  if (EVO(4,jd).LT.1.0D-10) then
                      EVO(4,jd)=0D0
                  end if
15            continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C... define initial values of RES
              RES(i,2)=0
              RES(i,3)=0
C...calculate the RES array
              do 16 kd=1,tgpcell
C...sum f (not unity)
                  RES(i,2)=RES(i,2)+4*PI*EVO(4,kd)*EVO(4,kd)*EVO(4,kd)/3
C...sum (Ri)
                  RES(i,3)=RES(i,3)+EVO(4,kd)
16            continue
C...sum N in 1m^3
              RES(i,1)=1.0D21*tgpcell
C...Ravg
              RES(i,4)=RES(i,3)/tgpcell
C...sum f in 1m3
              RES(i,2)=1.0D21*RES(i,2)
C...aveal in gamma prime
              RES(i,6)=RES(i-1,6)*RES(i-1,2)/RES(i,2)+
     &4*PI*cnal*(REFF**3)*1.0D0*DM/RES(i,2)/3+
     &algp*(RES(i,2)-RES(i-1,2)-4*PI*(REFF**3)*1.0D0*DM/3)/RES(i,2)
C...aual
              RES(i,5)=(RES(0,5)-RES(i,6)*RES(i,2))/(1-RES(i,2))
C...avecr in gamma prime
              RES(i,8)=RES(i-1,8)*RES(i-1,2)/RES(i,2)+
     &4*PI*cncr*(REFF**3)*1.0D0*DM/RES(i,2)/3+
     &crgp*(RES(i,2)-RES(i-1,2)-4*PI*(REFF**3)*1.0D0*DM/3)/RES(i,2)
C...aucr
              RES(i,7)=(RES(0,7)-RES(i,8)*RES(i,2))/(1-RES(i,2))
              print 22,tevo,RES(i,1),RES(i,2),RES(i,4),REFF,DM,DGM,GNV,
     &RES(i,5),RES(i,7),RES(i,6),RES(i,8),EDAVG
              write(2,22),tevo,RES(i,1),RES(i,2),RES(i,4),REFF,DM,DGM,
     &GNV,RES(i,5),RES(i,7),RES(i,6),RES(i,8),EDAVG

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...update the EVR array and tgpcell
C...define initial value
              tgpsum=0.0D0
              do 17 ld=1,tgpcell
                  if (EVO(4,ld) .NE. 0D0) then
                      tgpsum=tgpsum+1
                      EVR(1,tgpsum)=EVO(4,ld)
                  end if
17            continue
C...get the new number in 100nm^3
              tgpcell=tgpsum
C...print EVR
              write(3,33),tevo,(EVR(1,md),md=1,tgpcell)
33            format(100001ES11.3)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...end if of iteration 10
          end if
10    continue

C...print content
22    format(13ES11.3)

      pause
      end program kwn

max_stars_repo_path: Voronoi based KWN model.for

max_issues_repo_path: Voronoi based KWN model.for

max_forks_repo_path: Voronoi based KWN model.for


text:
  function new_StandardDCMResults(clientProgeny) result(results)

! ==============================================================================
!
! This function returns a standard DCM results data type (stores progeny array
! into results type)
!
!
! Written by CMJ, XCP-3, 12/2018
!
! ==============================================================================

    implicit none
    type(sDCMProgeny), intent(in   ), dimension(:), target :: clientProgeny
    type(StandardDCMResults)                               :: results

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

    ! Point to the desired progeny array
    results%progenyBnk => clientProgeny

    ! Set progeny information
    results%maxProgeny = size( results%progenyBnk )
    results%maxProgenyM3 = results%maxProgeny - 3

    ! Flag that results object was properly initialized
    results%initialized = .TRUE.

    return
! ==============================================================================
  end function new_StandardDCMResults

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