text:
module mod_prescribed_mesh_motion_function
#include <messenger.h>
    use mod_kinds,                              only: rk, ik
    use type_prescribed_mesh_motion_function,   only: prescribed_mesh_motion_function_t
    use type_pmmfvector,                        only: pmmfvector_t

    !
    ! Import prescribed_mesh_motion_functions
    !
    use pmmf_static,                            only: static_pmmf
    use pmmf_sinusoidal,                        only: sinusoidal_pmmf
    use pmmf_sinusoidal_1d,                     only: sinusoidal_1d_pmmf
    use pmmf_sinusoidal_2d,                     only: sinusoidal_2d_pmmf
    use pmmf_sinusoidal_uniform_flow,           only: sinusoidal_uniform_flow_pmmf
    use pmmf_sinusoidal_convecting_vortex,      only: sinusoidal_convecting_vortex_pmmf
    use pmmf_isotropic_dilation,                only: isotropic_dilation_pmmf
    use pmmf_constant_motion_1d_x,              only: constant_motion_1d_x_pmmf
    use pmmf_hpaf_case1,                        only: hpaf_case1_pmmf
    use pmmf_hpaf_case2,                        only: hpaf_case2_pmmf
    use pmmf_hpaf_case3,                        only: hpaf_case3_pmmf
    use pmmf_hpaf_case1_blended,                only: hpaf_case1_blended_pmmf
    use pmmf_hpaf_case2_blended,                only: hpaf_case2_blended_pmmf
    use pmmf_hpaf_case3_blended,                only: hpaf_case3_blended_pmmf
    use pmmf_fic_case1,                         only: fic_case1_pmmf
    use pmmf_fic_case2,                         only: fic_case2_pmmf
    use pmmf_fic_case3,                         only: fic_case3_pmmf
    use pmmf_rigid_body_motion,                 only: rigid_body_motion_pmmf
    implicit none

    !
    ! Global vector of registered prescribed_mesh_motion_functions
    !
    type(pmmfvector_t)          :: registered_pmmfs
    logical                     :: initialized = .false.

contains

    !>  Register prescribed_mesh_motion_functions in a module vector.
    !!
    !!  This allows the available prescribed_mesh_motion_functions to be queried in the same way that they
    !!  are registered for allocation. Adapted from mod_functions
    !!
    !!  @author Eric Wolf
    !!  @date   3/15/2017
    !!
    !!
    !--------------------------------------------------------------------------------------------
    subroutine register_prescribed_mesh_motion_functions()
        integer :: npmmfs, ipmmf

        !
        ! Instantiate prescribed_mesh_motion_functions
        !
        type(static_pmmf)                               :: static
        type(sinusoidal_pmmf)                           :: sinusoidal
        type(sinusoidal_1d_pmmf)                        :: sinusoidal_1d
        type(sinusoidal_2d_pmmf)                        :: sinusoidal_2d
        type(sinusoidal_uniform_flow_pmmf)              :: sinusoidal_uniform_flow
        type(sinusoidal_convecting_vortex_pmmf)         :: sinusoidal_convecting_vortex
        type(isotropic_dilation_pmmf)                   :: isotropic_dilation
        type(constant_motion_1d_x_pmmf)                 :: constant_motion_1d_x
        type(hpaf_case1_pmmf)                           :: hpaf_case1
        type(hpaf_case2_pmmf)                           :: hpaf_case2
        type(hpaf_case3_pmmf)                           :: hpaf_case3
        type(hpaf_case1_blended_pmmf)                   :: hpaf_case1_blended
        type(hpaf_case2_blended_pmmf)                   :: hpaf_case2_blended
        type(hpaf_case3_blended_pmmf)                   :: hpaf_case3_blended
        type(fic_case1_pmmf)                            :: fic_case1
        type(fic_case2_pmmf)                            :: fic_case2
        type(fic_case3_pmmf)                            :: fic_case3
        type(rigid_body_motion_pmmf)                    :: rigid_body_motion

        if ( .not. initialized ) then
            !
            ! Register in global vector
            !
            call registered_pmmfs%push_back(static)
            call registered_pmmfs%push_back(sinusoidal)
            call registered_pmmfs%push_back(sinusoidal_1d)
            call registered_pmmfs%push_back(sinusoidal_2d)
            call registered_pmmfs%push_back(sinusoidal_uniform_flow)
            call registered_pmmfs%push_back(sinusoidal_convecting_vortex)
            call registered_pmmfs%push_back(isotropic_dilation)
            call registered_pmmfs%push_back(constant_motion_1d_x)
            call registered_pmmfs%push_back(hpaf_case1)
            call registered_pmmfs%push_back(hpaf_case2)
            call registered_pmmfs%push_back(hpaf_case3)
            call registered_pmmfs%push_back(hpaf_case1_blended)
            call registered_pmmfs%push_back(hpaf_case2_blended)
            call registered_pmmfs%push_back(hpaf_case3_blended)
            call registered_pmmfs%push_back(fic_case1)
            call registered_pmmfs%push_back(fic_case2)
            call registered_pmmfs%push_back(fic_case3)
            call registered_pmmfs%push_back(rigid_body_motion)

            !
            ! Initialize each boundary condition in set. Doesn't need modified.
            !
            npmmfs = registered_pmmfs%size()
            do ipmmf = 1,npmmfs
                call registered_pmmfs%data(ipmmf)%pmmf%init()
            end do

            !
            ! Confirm initialization
            !
            initialized = .true.

        end if

    end subroutine register_prescribed_mesh_motion_functions
    !********************************************************************************************

    !> Factory method for allocating concrete prescribed_mesh_motion_functions
    !!
    !!      - Allocate a concrete prescribed_mesh_motion_function_t based on the incoming string specification.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/9/2016
    !!
    !!  @param[in]      string  Character string used to select the appropriate boundary condition
    !!  @param[inout]   pmmf      Allocatable boundary condition
    !!
    !--------------------------------------------------------------
    subroutine create_prescribed_mesh_motion_function(pmmf,string)
        class(prescribed_mesh_motion_function_t),  allocatable,    intent(inout)   :: pmmf
        character(*),                       intent(in)      :: string

        integer(ik) :: ierr, pmmfindex

        if ( allocated(pmmf) ) then
            deallocate(pmmf)
        end if

        !
        ! Find equation set in 'registered_bcs' vector
        !
        pmmfindex = registered_pmmfs%index_by_name(trim(string))

        !
        ! Check prescribed_mesh_motion_function was found in 'registered_pmmfs'
        !
        if (pmmfindex == 0) call chidg_signal_one(FATAL,"create_prescribed_mesh_motion_function: prescribed_mesh_motion_function not recognized", trim(string))

        !
        ! Allocate conrete prescribed_mesh_motion_function_t instance
        !
        allocate(pmmf, source=registered_pmmfs%data(pmmfindex)%pmmf, stat=ierr)
        if (ierr /= 0) call chidg_signal(FATAL,"create_prescribed_mesh_motion_function: error allocating prescribed_mesh_motion_function from global vector.")

        !
        ! Check boundary condition was allocated
        !
        if ( .not. allocated(pmmf) ) call chidg_signal(FATAL,"create_prescribed_mesh_motion_function: error allocating concrete prescribed_mesh_motion_function.")

    end subroutine create_prescribed_mesh_motion_function
    !*****************************************************************************************

    !>  Print a list of the registered prescribed_mesh_motion_functions. Really just a utility for 'chidg edit' to
    !!  dynamically list the available 'prescribed_mesh_motion_function_t's.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   11/8/2016
    !!
    !!
    !------------------------------------------------------------------------------------------
    subroutine list_prescribed_mesh_motion_functions()
        integer                         :: npmmfs, ipmmf
        character(len=:),   allocatable :: pmmf_name

        npmmfs = registered_pmmfs%size()

        do ipmmf = 1,npmmfs

            pmmf_name = registered_pmmfs%data(ipmmf)%pmmf%get_name()
            call write_line(trim(pmmf_name))

        end do ! ipmmf

    end subroutine list_prescribed_mesh_motion_functions
    !******************************************************************************************

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   9/1/2016
    !!
    !!
    !!
    !------------------------------------------------------------------------------------------------------
    function check_pmmf_registered(state_string) result(state_found)
        character(len=*),   intent(in)  :: state_string

        integer(ik) :: state_index
        logical     :: state_found

        ! Find boundary condition string in 'registered_bcs' vector
        state_index = registered_pmmfs%index_by_name(trim(state_string))

        ! Set status of state_found
        state_found = (state_index /= 0)

    end function check_pmmf_registered
    !*******************************************************************************************************

end module mod_prescribed_mesh_motion_function

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
	SUBROUTINE SFL6HD  ( nlun, lun, iret )
C************************************************************************
C* SFL6HD								*
C*									*
C* This subroutine writes out the parameter name header.		*
C*									*
C* SFL6HD  ( NLUN, LUN, IRET )						*
C*									*
C* Input parameters:							*
C*	NLUN		INTEGER		Number of output devices 	*
C*	LUN  (NLUN)	INTEGER		LUNs 				*
C*									*
C* Output parameters:							*
C* 	IRET		INTEGER		Return code			*
C*					   0 = normal return		*
C**									*
C* Log:									*
C* I. Graffman/RDS	 7/87	GEMPAK4					*
C* M. desJardins/GSFC	 6/88	Leave space for 6 digit station		*
C* J. Whistler/SSAI	 2/91   Changed the header to incorporate new	*
C*				parameters				*
C* S. Jacobs/EAI	10/92	Changed horizontal spacing		*
C* K. Tyle/GSC		 1/97	Add a blank line at top			*
C************************************************************************
	INTEGER		lun (*)
C*
	CHARACTER	prmbuf*79
C*
	DATA 		prmbuf ( 1 : 49 )
     +		/ 'HH   STN  CLDL  CLDM  CLDH VSBY  WTHR  PMSL  TMP ' /
	DATA 		prmbuf ( 50 :  )
     +		/ 'DWP DDSS GST ALT  PTND PCPN SN' /
C------------------------------------------------------------------------
C*	Write header to requested luns.
C
	DO  ilun = 1, nlun
	    WRITE   ( lun ( ilun ), * ) ' '
	    WRITE   ( lun ( ilun ), 5000 )  prmbuf
5000	    FORMAT  ( ' ', A )
	END DO
	iret = 0
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


    Subroutine nagf_blas_zddscl(n, d, incd, x, incx)

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     NAGF_BLAS_ZDDSCL performs the operation

!     x := diag( d )*x

!     Nag Fortran 77 O( n ) basic linear algebra routine.

!     -- Written on 12-February-1986.
!     Sven Hammarling, Nag Central Office.

!     .. Use Statements ..
      Use blas_interfaces, Only: zdscal
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Scalar Arguments ..
      Integer, Intent (In) :: incd, incx, n
!     .. Array Arguments ..
      Complex (Kind=dp), Intent (Inout) :: x(*)
      Real (Kind=dp), Intent (In) :: d(*)
!     .. Local Scalars ..
      Integer :: i, id, ix
!     .. Intrinsic Procedures ..
      Intrinsic :: abs, cmplx
!     .. Executable Statements ..
      Continue

      If (n>0) Then
        If ((incd==0) .And. (incx/=0)) Then
          Call zdscal(n=n, alpha=d(1), x=x, incx=abs(incx))
        Else If ((incd==incx) .And. (incd>0)) Then
          Do id = 1, 1 + (n-1)*incd, incd
            x(id) = cmplx(d(id), kind=dp)*x(id)
          End Do
        Else
          If (incx>=0) Then
            ix = 1
          Else
            ix = 1 - (n-1)*incx
          End If
          If (incd>0) Then
            Do id = 1, 1 + (n-1)*incd, incd
              x(ix) = cmplx(d(id), kind=dp)*x(ix)
              ix = ix + incx
            End Do
          Else
            id = 1 - (n-1)*incd
            Do i = 1, n
              x(ix) = cmplx(d(id), kind=dp)*x(ix)
              id = id + incd
              ix = ix + incx
            End Do
          End If
        End If
      End If

      Return

!     End of NAGF_BLAS_ZDDSCL.

    End Subroutine


text:
subroutine rp1(maxm,num_eqn,num_waves,num_aux,num_ghost,num_cells,ql,qr,auxl,auxr,fwave,s,amdq,apdq)

! Riemann solver for the 1D shallow water equations:

! (h)_t + (h*u)_x = 0
! (hu)_t + (h*u^2 + 1/2*grav*h^2)_x = -grav*h*(b)_x
!
! using the f-wave algorithm and Roe's approximate Riemann solver.

! waves:     2
! equations: 2

! Conserved quantities:
!       1 depth
!       2 momentum

! Auxiliary fields:
!       1 bathymetry

! The gravitational constant grav should be in the common block cparam.
!
! See http://www.clawpack.org/riemann.html for a detailed explanation
! of the Riemann solver API.

    implicit none

    integer, intent(in) :: maxm, num_eqn, num_waves, num_ghost, num_aux, num_cells
    real(kind=8), intent(in) :: ql(num_eqn, 1-num_ghost:maxm+num_ghost)
    real(kind=8), intent(in) :: qr(num_eqn, 1-num_ghost:maxm+num_ghost)
    real(kind=8), intent(in) :: auxl(num_aux, 1-num_ghost:maxm+num_ghost)
    real(kind=8), intent(in) :: auxr(num_aux, 1-num_ghost:maxm+num_ghost)
    real(kind=8), intent(out) :: s(num_waves, 1-num_ghost:maxm+num_ghost)
    real(kind=8), intent(out) :: fwave(num_eqn, num_waves, 1-num_ghost:maxm+num_ghost)
    real(kind=8), intent(out) :: amdq(num_eqn, 1-num_ghost:maxm+num_ghost)
    real(kind=8), intent(out) :: apdq(num_eqn, 1-num_ghost:maxm+num_ghost)

    real(kind=8) :: grav, dry_tolerance, sea_level
    common /cparam/ grav, dry_tolerance, sea_level

    real(kind=8) :: hl, ul, hr, ur, hbar, uhat, chat, bl, br, phil, phir
    real(kind=8) :: R(2,2), fluxdiff(2), beta(2)

    integer :: i, k

    amdq = 0.d0
    apdq = 0.d0

    do i=2-num_ghost,num_cells+num_ghost
        ! # Left states
        hl = qr(1, i - 1)
        if (hl > dry_tolerance) then
            ul = qr(2, i - 1) / hl
            phil = qr(1, i - 1) * ul**2 + 0.5d0 * grav * qr(1, i - 1)**2
        else
            ul = 0.d0
            phil = 0.d0
        end if
        bl = auxr(1, i - 1)

        ! # Right states
        hr = ql(1, i)
        if (hr > dry_tolerance) then
            ur = ql(2, i)/hr
            phir = ql(1, i) * ur**2 + 0.5d0 * grav * ql(1, i)**2
        else
            ur = 0.d0
            phir = 0.d0
        end if
        br = auxl(1, i)

        ! # Roe average states
        hbar = 0.5 * (hr + hl)
        uhat = (sqrt(hl) * ul + sqrt(hr) * ur) / (sqrt(hl) + sqrt(hr))
        chat = sqrt(grav * hbar)

        ! # Flux differences
        fluxdiff(1) = (hr * ur) - (hl * ul)
        fluxdiff(2) = phir - phil + grav * hbar * (br - bl)

        ! # Wave speeds
        s(1,i) = min(ul - sqrt(grav * hl), uhat - chat)
        s(2,i) = max(ur + sqrt(grav * hr), uhat + chat)

        ! # Right eigenvectors (column)
        R(1,1) = 1.d0
        R(2,1) = s(1, i)

        R(1,2) = 1.d0
        R(2,2) = s(2, i)

        ! Wave strengths
        beta(1) = (s(2, i) * fluxdiff(1) - fluxdiff(2)) / (s(2, i) - s(1, i))
        beta(2) = (fluxdiff(2) - s(1, i) * fluxdiff(1)) / (s(2, i) - s(1, i))

        ! # Flux waves
        do k=1,num_waves
            fwave(:,k,i) = beta(k) * R(:,k)
        enddo

        ! # Fluctuations
        do k=1, num_waves
            if (s(k, i) < 0.d0) then
                amdq(:, i) = amdq(:, i) + fwave(:, k, i)
            else if (s(k, i) > 0.d0) then
                apdq(:, i) = apdq(:, i) + fwave(:, k, i)
            else
                amdq(:, i) = amdq(:, i) + 0.5d0 * fwave(:, k, i)
                apdq(:, i) = apdq(:, i) + 0.5d0 * fwave(:, k, i)
            end if
        enddo
    enddo

end subroutine rp1

max_stars_repo_name: ranocha/riemann

max_issues_repo_name: ranocha/riemann

max_forks_repo_name: ranocha/riemann


text:
      program crc02
      integer n
      n = 17
      print *, n+25
      end


SUBROUTINE matrix_Rr (r1,r2 , r3)

! ----------------------------------------------------------------------
! SUBROUTINE: matrix_Rr.f90
! ----------------------------------------------------------------------
! Purpose:
!  Multiply 2 arrays : 3x3 multiply 3x1
! ----------------------------------------------------------------------
! Input arguments:
! - r1 : Array 1 with dimensions 3x3
! - r2 : Array 2 with dimensions 3x1
!
! Output arguments:
! - r3 : Array obtained from multiplying r1 * r2
! ----------------------------------------------------------------------
! Dr. Thomas Papanikolaou, Geoscience Australia            26 April 2016
! ----------------------------------------------------------------------

      USE mdl_precision
      IMPLICIT NONE

! ---------------------------------------------------------------------------
! Dummy arguments declaration
! ---------------------------------------------------------------------------
! IN
      REAL (KIND = prec_q), INTENT(IN) :: r1(3,3), r2(3)
! OUT
      REAL (KIND = prec_q), INTENT(OUT) :: r3(3)
! ---------------------------------------------------------------------------

      r3(1) = r1(1,1) * r2(1) + r1(1,2) * r2(2) + r1(1,3) * r2(3)
      r3(2) = r1(2,1) * r2(1) + r1(2,2) * r2(2) + r1(2,3) * r2(3)
      r3(3) = r1(3,1) * r2(1) + r1(3,2) * r2(2) + r1(3,3) * r2(3)

END


text:
!---The program is used to calculate k-distributions (32 GaussâChebyshev quadrature points)
!---Pressure: 0.1 bar ~ 80 bar
!---Temperature: recommended 300 K ~ 3000 K
!---Species: CO2 or H2O or CO or fv or the mixture of CO2&H2O&CO&fv
PROGRAM OneDimensionExact
!--------------------------------------------------------------------------------
use commonData, only : GasMixInfo
use MLP
use commonData, only : GasMixInfo
IMPLICIT NONE
!--------------------------------------------------------------------------------
  integer,parameter :: Nq=32
  real(8):: Tref ! reference temperature
  real(8),dimension(Nq) :: k ! k-values
  TYPE(GasMixInfo) :: locMix_f        ! Local mixture information
  OPEN(UNIT=32, FILE='results/k_new.dat', STATUS='REPLACE') ! new k-distributions
!--------------------------------------------------------------------------------
  locMix_f%p=1.0d0
  locMix_f%T=1500.0d0
  Tref=1500.0d0
  locMix_f%xCO2=0.1d0
  locMix_f%xH2O= 0.2d0
  locMix_f%xCO=0.05d0
  locMix_f%fv =1d-10
!~~~~~~~~~
  ! get k-values
  if (locMix_f%P<=1) then
    call get_k_mlp_1(locMix_f,Tref,k)
  else if (1<locMix_f%P .AND. locMix_f%P<10) then
    call get_k_mlp_2(locMix_f,Tref,k)
  else if (locMix_f%P>=10) then
    call get_k_mlp_3(locMix_f,Tref,k)
  end if
!~~~~~~~~~
  write(32,*) k
  write(*,*)k
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
END PROGRAM OneDimensionExact
!--------------------------------------------------------------------------------

max_stars_repo_path: Fortran codes/Main.f90

max_issues_repo_path: Fortran codes/Main.f90

max_forks_repo_path: Fortran codes/Main.f90


text:
      module SparseBLAS
! **********************************************************************
!     Author : C. Voemel
!
!     Date of last modification : 29.2.00
!
!     Description : SparseBLAS functions
! **********************************************************************

      use mod_info
      use mod_uscr
      use mod_usds
      use mod_usmv
      use mod_ussv
      use mod_usmm
      use mod_ussm

      end module SparseBLAS

max_stars_repo_name: NetlibRepository/lapack-www

max_issues_repo_name: NetlibRepository/lapack-www

max_forks_repo_name: NetlibRepository/lapack-www


! Copyright 2019-2020 Free Software Foundation, Inc.
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.

program pointers

  type :: two
    integer, allocatable :: ivla1 (:)
    integer, allocatable :: ivla2 (:, :)
  end type two

  logical, target :: logv
  complex, target :: comv
  character, target :: charv
  character (len=3), target :: chara
  integer, target :: intv
  integer, target, dimension (10,2) :: inta
  real, target :: realv
  type(two), target :: twov

  logical, pointer :: logp
  complex, pointer :: comp
  character, pointer :: charp
  character (len=3), pointer :: charap
  integer, pointer :: intp
  integer, pointer, dimension (:,:) :: intap
  real, pointer :: realp
  type(two), pointer :: twop

  nullify (logp)
  nullify (comp)
  nullify (charp)
  nullify (charap)
  nullify (intp)
  nullify (intap)
  nullify (realp)
  nullify (twop)

  logp => logv    ! Before pointer assignment
  comp => comv
  charp => charv
  charap => chara
  intp => intv
  intap => inta
  realp => realv
  twop => twov

  logv = associated(logp)     ! Before value assignment
  comv = cmplx(1,2)
  charv = "a"
  chara = "abc"
  intv = 10
  inta(:,:) = 1
  inta(3,1) = 3
  realv = 3.14

  allocate (twov%ivla1(3))
  allocate (twov%ivla2(2,2))
  twov%ivla1(1) = 11
  twov%ivla1(2) = 12
  twov%ivla1(3) = 13
  twov%ivla2(1,1) = 211
  twov%ivla2(2,1) = 221
  twov%ivla2(1,2) = 212
  twov%ivla2(2,2) = 222

  intv = intv + 1 ! After value assignment

end program pointers


text:
program sbcalibrand

  use globals
  use gdal
  use mapio
  use slezanbear
  use geo

  implicit none

  type dataentry
    real(dp) :: lat, lng
    real(dp) :: sky, err
    real(dp) :: model
  end type

  real(dp) :: llatsrc, llngsrc, ulatsrc, ulngsrc
  type(modelpar) :: par
  type(dataentry), dimension(1024), target :: dat_all
  type(dataentry), dimension(:), pointer :: dat
  integer :: ndata, i, j

  character(128), dimension(9) :: fndem

  real(sp), dimension(:,:), allocatable :: mapi, maph
  type(geotransform) :: gti, gth
  integer :: errno
  real(dp) :: sky(4), hobs

  character(*), parameter :: fni = 'SVDNB_npp_20150101-20151231_75N060W_vcm-orm-ntl_v10_c201701311200.avg_rade9.tif'
  character(*), parameter :: datafmt = '(F9.4,1X,F9.4,3X,F5.2,1X,F4.2,3X,F5.2)'

  terrain_attenuation = .false.

  ndata = 0
  open(33, file = 'data.txt', action = 'read')
  do i = 1,size(dat_all)
    read(33, *, iostat = errno) dat_all(i) % lat, dat_all(i) % lng, &
      & dat_all(i) % sky, dat_all(i) % err
    if (errno .ne. 0) exit
    ndata = ndata + 1
  end do
  close(33)

  dat => dat_all(1:ndata)

  llat = minval(dat % lat)
  ulat = maxval(dat % lat)
  llng = minval(dat % lng)
  ulng = maxval(dat % lng)

  fndem(1) = "eudem_dem_5deg_n45e010.tif"
  fndem(2) = "eudem_dem_5deg_n45e015.tif"
  fndem(3) = "eudem_dem_5deg_n45e020.tif"
  fndem(4) = "eudem_dem_5deg_n50e010.tif"
  fndem(5) = "eudem_dem_5deg_n50e015.tif"
  fndem(6) = "eudem_dem_5deg_n50e020.tif"
  fndem(7) = "eudem_dem_5deg_n55e010.tif"
  fndem(8) = "eudem_dem_5deg_n55e015.tif"
  fndem(9) = "eudem_dem_5deg_n55e020.tif"

  write (0,'("dmax = ",F0.2," km")') dmax / 1e3
  call estimate_data_limits(llat, llng, ulat, ulng, dmax * 1.05, &
  &       llatsrc, llngsrc, ulatsrc, ulngsrc)
  print '("data limits = ",F6.2,1X,F7.2)', llatsrc, llngsrc, ulatsrc, ulngsrc

  call GDALAllRegister

  call map_read_section(fni, llatsrc, llngsrc, ulatsrc, ulngsrc, mapi, gti)
  write (0, '("satellite raster size = ", 2I6)') size(mapi,1), size(mapi,2)

  read_dem: block
    integer nox, noy

    call map_gen_size(llatsrc, llngsrc, ulatsrc, ulngsrc, elev_grid_meters, &
    &              nox, noy, gth)

    allocate( maph(nox,noy) )
    write (0,'("dem raster size = ",2I6)') size(maph,1), size(maph,2)

    call map_read_tiles(fndem, maph, gth)
  end block read_dem

  call random_seed()

  printcalib: block
    type(source), dimension(:,:), allocatable :: src
    real :: a(5)
    real(dp) :: err(ndata)

    allocate( src(size(mapi,1), size(mapi,2)) )
    call mksources(mapi, gti, maph, gth, src)

    open(34, file = 'calib.txt', action = 'write')
    write (34,'(A8,A7,A8,A6,A7,A10,A10,A7)') 'SCTR', 'ABS', 'HSCAL', 'SIDE', &
    & 'BG', 'CHISQ', 'DEVITY', 'DEVMAG'

    iterate_parameters: do j = 1,16807
      call random_number(a)

      par % alpha =  M(a(1), 3.0, 24.0) * 1e-6
      par % relabs = M(a(2), 0.00, 12.00)
      par % hscale = M(a(3), 2.0, 12.0) * 1e3
      par % beta =   M(a(4), 0.0, 0.4)
      par % skybg =  mag2ity(M(a(5), 22.1, 21.6))

      !$omp parallel do private(sky,hobs)
      iterate_datpoints: do i = 1, ndata
        call onepoint(par, dat(i) % lat, dat(i) % lng, &
        &       src, maph, gth, sky, hobs)
        dat(i) % model = sky(4)
      end do iterate_datpoints
      !$omp end parallel do

      write (34,'(F8.3,F7.3,F8.1,F6.3,F7.3,ES10.3,ES10.3,F7.4)') &
        & par % alpha * 1e6, par % relabs, &
        & par % hscale, par % beta, (par % skybg) * 1e7, &
        &      sum((mag2ity(dat % sky) - dat % model)**2  / dat % model), &
        & sqrt(sum((mag2ity(dat % sky) - dat % model)**2) / ndata), &
        & sqrt(sum((dat % sky - ity2mag(dat % model))**2) / ndata)
      flush(34)
    end do iterate_parameters

    close(34)

    deallocate(src)
  end block printcalib

  deallocate(mapi, maph)

contains

  elemental real(dp) function m(a, l, u) result(b)
    real, intent(in) :: a
    real, intent(in) :: l, u
    b = l + a * (u - l)
  end function

end program sbcalibrand

max_stars_repo_name: gronki/slezanbear

max_issues_repo_name: gronki/slezanbear

max_forks_repo_name: gronki/slezanbear


text:
subroutine varden()

  use variables
  use network
  use burner_module, only: burner_init
  use actual_rhs_module, only: actual_rhs_init
  use geometry
  use base_state_module
  use base_io_module
  use aux_data_module
  use estdt_module
  use firstdt_module
  use advance_timestep_module
  use box_util_module
  use fabio_module
  use checkpoint_module
  use make_plotfile_module, only: MAX_FILENAME_LEN, get_plot_names, make_plotfile, &
                                  make_filename
  use restart_module
  use probin_module
  use runtime_init_module
  use bl_constants_module
  use average_module
  use make_grav_module
  use sponge_module
  use make_beta0_module
  use define_bc_module
  use fill_3d_module
  use eos_module
  use conductivity_module
  use divu_iter_module
  use initial_proj_module
  use make_gamma_module
  use enforce_HSE_module
  use rhoh_vs_t_module
  use initialize_module
  use init_particles_module
  use make_new_grids_module
  use regrid_module
  use make_eta_module
  use diag_module, only: flush_diag, diag_finalize
  use omp_module
  use check_cutoff_module, only: check_cutoff_values
  use time_module, only: time
  use particle_module
  use cputime_module, only: start_cputime_clock
  use simple_log_module, only: simple_log_init, simple_log_finalize
  use pert_form_module, only: put_in_pert_form
  use model_parser_module, only: model_parser_finalize

  implicit none

  integer    :: init_step,istep
  integer    :: istep_divu_iter,istep_init_iter
  integer    :: i,n,r,numcell,dm,nlevs
  integer    :: last_plt_written,last_chk_written
  real(dp_t) :: smin,smax,smaxold,nuclear_dt_scalefac
  real(dp_t) :: umin,umax,vmin,vmax,wmin,wmax
  real(dp_t) :: dt,dtold

  type(ml_layout) :: mla
  type(bc_tower)  :: the_bc_tower

  real(dp_t)  , pointer     :: dx(:,:)

  ! initialized below. not checkpointed.
  type(multifab), allocatable :: unew(:)
  type(multifab), allocatable :: snew(:)
  type(multifab), allocatable :: S_cc_new(:)
  type(multifab), allocatable :: normal(:)
  type(multifab), allocatable :: sponge(:)
  type(multifab), allocatable :: nodalrhs(:)

  !!!!!!!!!!!!!
  ! these are pointers because they need to be allocated and built within
  !   another function (in initialize.f90)

  ! these are in the checkpoint, and need to be passed into regrid
  type(multifab), pointer :: uold(:)
  type(multifab), pointer :: sold(:)
  type(multifab), pointer :: gpi(:)
  type(multifab), pointer :: dSdt(:)
  type(multifab), pointer :: S_cc_old(:)

  ! these are in the checkpoint only so the plotfile after checkpoint
  ! matches.  They are not persistent, and do not need to be passed into regrid
  type(multifab), pointer :: pi(:)
  type(multifab), pointer :: rho_omegadot2(:)
  type(multifab), pointer :: rho_Hnuc2(:)
  type(multifab), pointer :: rho_Hext(:)
  type(multifab), pointer :: thermal2(:)

  type(multifab), pointer :: tag_mf(:)
  type(multifab), pointer :: chkdata(:)

  character(len=MAX_FILENAME_LEN) :: plot_file_name, check_file_name

  integer :: npartdata
  integer, allocatable :: index_partdata(:)
  character(len=16), allocatable :: names_partdata(:)

  real(dp_t), parameter :: SMALL = 1.d-13
  real(dp_t)            :: runtime1, runtime2

  logical :: init_mode

  real(dp_t), pointer :: beta0_old(:,:)
  real(dp_t), pointer :: beta0_new(:,:)
  real(dp_t), pointer :: gamma1bar_old(:,:)
  real(dp_t), pointer :: gamma1bar_new(:,:)
  real(dp_t), pointer :: s0_init(:,:,:)
  real(dp_t), pointer :: rho0_old(:,:)
  real(dp_t), pointer :: rhoh0_old(:,:)
  real(dp_t), pointer :: rho0_new(:,:)
  real(dp_t), pointer :: rhoh0_new(:,:)
  real(dp_t), pointer :: p0_init(:,:)
  real(dp_t), pointer :: p0_old(:,:)
  real(dp_t), pointer :: p0_new(:,:)
  real(dp_t), pointer :: w0(:,:)
  real(dp_t), pointer :: etarho_ec(:,:)
  real(dp_t), pointer :: etarho_cc(:,:)
  real(dp_t), pointer :: psi(:,:)
  real(dp_t), pointer :: tempbar(:,:)
  real(dp_t), pointer :: tempbar_init(:,:)
  real(dp_t), pointer :: grav_cell_old(:,:)
  real(dp_t), pointer :: grav_cell_new(:,:)

  real(dp_t), allocatable :: psi_temp(:,:)
  real(dp_t), allocatable :: etarho_cc_temp(:,:)
  real(dp_t), allocatable :: rho0_temp(:,:)
  real(dp_t), allocatable :: etarho_ec_temp(:,:)
  real(dp_t), allocatable :: w0_temp(:,:)
  real(dp_t), allocatable :: tempbar_init_temp(:,:)

  logical :: dump_plotfile, dump_checkpoint, abort_maestro
  real(dp_t) :: write_pf_time

  type(particle_container) :: particles
  integer :: numparticles

  logical :: have_overview

  type(plot_t) :: plt(2)  ! 1 = main plotfile; 2 = mini plotfile
  integer :: ip

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! initialization
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! keep track of cputime
  call start_cputime_clock()

  last_plt_written = -1
  last_chk_written = -1

  call runtime_init()
  call init_spherical()
  call init_center()
  call init_rotation()

  call init_variables()

  ! initialize the microphysics modules
  call network_init()
  call actual_rhs_init()
  call burner_init()
  call eos_init(small_temp=small_temp,small_dens=small_dens)
  call conductivity_init()

  ! setup the plotfiles
  call init_plot_variables(plt(1), plot_int, plot_deltat, plot_base_name)
  call get_plot_names(plt(1))

  call init_miniplot_variables(plt(2), mini_plot_int, mini_plot_deltat, mini_plot_base_name)
  call get_plot_names(plt(2))

  ! particle initialization
  call particle_setverbose(.true.)

  ! a simple logging facility for writting both to the screen and a file
  call simple_log_init()

  ! collect the indices into the s multifab and the associated variable
  ! names for the data we will store along with each particle.
  npartdata = 2 + nspec
  allocate(index_partdata(npartdata))
  allocate(names_partdata(npartdata))

  index_partdata(1) = rho_comp
  names_partdata(1) = "density"

  index_partdata(2) = temp_comp
  names_partdata(2) = "temperature"

  do n = 1, nspec
     index_partdata(2+n) = spec_comp -1 + n
     names_partdata(2+n) = spec_names(n)
  enddo

  if (restart >= 0) then

     call initialize_from_restart(mla,restart,dt,pmask,dx,uold,sold,gpi,pi, &
                                  dSdt,S_cc_old, &
                                  rho_omegadot2,rho_Hnuc2,rho_Hext,thermal2,the_bc_tower, &
                                  beta0_old,beta0_new,gamma1bar_old,gamma1bar_new, &
                                  s0_init,rho0_old,rhoh0_old,rho0_new,rhoh0_new,p0_init, &
                                  p0_old,p0_new,w0,etarho_ec,etarho_cc,psi, &
                                  tempbar,tempbar_init,grav_cell_old,grav_cell_new)

     check_file_name = make_filename(check_base_name, restart)

     if (use_particles) then
        call particle_container_restart(particles,check_file_name,&
                                        mla,dx,prob_lo)
     endif

  else if (test_set /= '') then

     if(use_particles) call build(particles)

     call initialize_with_fixed_grids(mla,dt,pmask,dx,uold,sold,gpi,pi,dSdt, &
                                      S_cc_old, &
                                      rho_omegadot2,rho_Hnuc2,rho_Hext,thermal2, &
                                      the_bc_tower, &
                                      beta0_old,beta0_new,gamma1bar_old, &
                                      gamma1bar_new,s0_init,rho0_old,rhoh0_old, &
                                      rho0_new,rhoh0_new,p0_init,p0_old,p0_new,w0, &
                                      etarho_ec,etarho_cc,psi,tempbar,tempbar_init, &
                                      grav_cell_old,grav_cell_new)

  else

     if (use_particles) call build(particles)

     call initialize_with_adaptive_grids(mla,dt,pmask,dx,uold,sold,gpi,pi,dSdt, &
                                         S_cc_old, &
                                         rho_omegadot2,rho_Hnuc2,rho_Hext,thermal2, &
                                         the_bc_tower, &
                                         beta0_old,beta0_new,gamma1bar_old, &
                                         gamma1bar_new,s0_init,rho0_old,rhoh0_old, &
                                         rho0_new,rhoh0_new,p0_init,p0_old,p0_new,w0, &
                                         etarho_ec,etarho_cc,psi,tempbar,tempbar_init, &
                                         grav_cell_old,grav_cell_new)

  end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! error checking
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  dm = dm_in
  nlevs = mla%nlevel

  ! check to make sure dimensionality is consistent in the inputs file
  if (dm .ne. get_dim(mla%mba)) then
     call bl_error('dm_in not properly set in inputs file')
  end if

  ! check to make sure spherical is only used for 3d
  if (spherical .eq. 1 .and. dm .ne. 3) then
     call bl_error("spherical = 1 and dm != 3")
  end if

  ! check to make sure our grid is square -- the solvers assume this
  if (dm == 2) then
     if (abs(dx(1,1) - dx(1,2)) > SMALL) then
        call bl_error('zones must be square')
     end if
  else if (dm == 3) then
     if (abs(dx(1,1) - dx(1,2)) > SMALL .OR. abs(dx(1,1) - dx(1,3)) > SMALL) then
        call bl_error('zones must be square')
     end if
  end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warn the user if we are burning and heating - this might not be what we want
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (do_burning .and. do_heating .and. parallel_IOProcessor()) then
     print *, 'WARNING: attempting to apply heating and a reaction network'
     print *, '         make sure this is what you want to do, otherwise check'
     print *, '         the do_heating and do_burning flag settings.'
  endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! print processor and grid info
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  if (parallel_IOProcessor()) then
     print *, 'number of MPI processes = ', parallel_nprocs()
     print *, 'number of threads       = ', omp_get_max_threads()
     print *, ' '
     print *, 'number of dimensions    = ', dm
     do n = 1, nlevs
        print *, 'level: ', n
        print *, '   number of boxes = ', nboxes(S_cc_old(n)%la)
        print *, '   maximum zones   = ', (extent(mla%mba%pd(n),i),i=1,dm)
     end do
     print *, ' '
  end if

  if (verbose .ge. 1) then
     do n=1,nlevs
        numcell = multifab_volume(S_cc_old(n),.false.)
        if (parallel_IOProcessor()) then
           print*,"Number of valid cells at level        ",n,numcell
        end if
        numcell = multifab_volume(S_cc_old(n),.true.)
        if (parallel_IOProcessor()) then
           print*,"Number of valid + ghost cells at level",n,numcell
        end if
     end do
     if (parallel_IOProcessor()) then
        print*,""
     end if
  end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Initialize all remaining arrays and multifabs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  if (spherical .eq. 0) then
     allocate(       psi_temp(max_levs,0:nr_fine-1))
     allocate( etarho_cc_temp(max_levs,0:nr_fine-1))
     allocate(      rho0_temp(max_levs,0:nr_fine-1))
     allocate( etarho_ec_temp(max_levs,0:nr_fine))
     allocate(        w0_temp(max_levs,0:nr_fine))
     allocate(tempbar_init_temp(max_levs,0:nr_fine-1))
  else
     allocate(       psi_temp(1,0:nr_fine-1))
     allocate( etarho_cc_temp(1,0:nr_fine-1))
     allocate(      rho0_temp(1,0:nr_fine-1))
     allocate( etarho_ec_temp(1,0:nr_fine))
     allocate(        w0_temp(1,0:nr_fine))
     allocate(tempbar_init_temp(1,0:nr_fine-1))
  end if

  allocate(unew(nlevs),snew(nlevs),S_cc_new(nlevs))
  allocate(normal(nlevs),sponge(nlevs),nodalrhs(nlevs))
  allocate(tag_mf(nlevs))

  do n = 1,nlevs
     call multifab_build(    unew(n), mla%la(n),    dm, nghost(uold(n)))
     call multifab_build(    snew(n), mla%la(n), nscal, nghost(sold(n)))
     call multifab_build(S_cc_new(n), mla%la(n),     1, 0)
     if (dm .eq. 3) then
        call multifab_build(normal(n), mla%la(n),    dm, 1)
     end if
     call multifab_build(  sponge(n), mla%la(n),     1, 0)
     call multifab_build( nodalrhs(n), mla%la(n),     1, 0, nodal)
     call multifab_build(  tag_mf(n), mla%la(n), 1, 0)

     call setval(      unew(n), ZERO, all=.true.)
     call setval(      snew(n), ZERO, all=.true.)
     call setval(  S_cc_new(n), ZERO, all=.true.)
     call setval(    sponge(n), ONE,  all=.true.)
     call setval(   nodalrhs(n), ZERO, all=.true.)
     call setval(    tag_mf(n), ZERO, all=.true.)
  end do

  ! Create normal now that we have defined center and dx
  call make_normal(normal,dx)

  if (do_sponge) then
     call init_sponge(rho0_old(1,:))
  end if

  call make_grav_cell(grav_cell_old,rho0_old)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Print out some diagnostics and warnings about cutoff / sponge parameter
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  call check_cutoff_values()

1010 format(1x, "restarted from step ", i7, 2x, "MPI tasks = ", i6, 2x, "OMP threads = ", i3)

  ! write the maestro-overview.out file
  if (restart < 0) then
     call write_job_info("", mla%mba, the_bc_tower, 0.0d0)
  else
     if (parallel_IOProcessor()) then
        inquire(file="maestro-overview.out", exist=have_overview)
        if (have_overview) then
           open(unit=99, file="maestro-overview.out", form = "formatted", &
                action="write", status="old", position="append")
           write(99,1010) restart, parallel_nprocs(), omp_get_max_threads()
           close(unit=99)
        endif
     endif
  endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Initialization
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  if (restart < 0) then

     !----------------------------------------------------------------------
     ! Do an initial projection with omegadot = 0 and rho_Hext = 0
     !----------------------------------------------------------------------

     call make_gamma1bar(mla,sold,gamma1bar_old,p0_old,dx)

     call make_beta0(beta0_old,rho0_old,p0_old,gamma1bar_old,grav_cell_old)

     if(do_initial_projection) then
        call initial_proj(uold,sold,pi,gpi,S_cc_old,normal,nodalrhs,thermal2, &
                          beta0_old,p0_old,gamma1bar_old,dx,the_bc_tower,mla)
     end if

     !----------------------------------------------------------------------
     ! Compute the initial time step
     !----------------------------------------------------------------------

     call firstdt(mla,the_bc_tower%bc_tower_array,uold,gpi,sold,S_cc_old, &
                  rho0_old,p0_old,grav_cell_old,gamma1bar_old,dx,cflfac,dt)

     !------------------------------------------------------------------------
     ! do the initial iterations to define a velocity field consistent with S
     !------------------------------------------------------------------------
111 format(78('-'))

     if ( parallel_IOProcessor() ) then
        print *, ' '
        write (*,111)
        print *, 'DOING',init_divu_iter ,'INITIAL DIVU ITERATIONS'
        write (*,111)
        print *, ' '
     end if

     do istep_divu_iter=1,init_divu_iter

        call divu_iter(istep_divu_iter,uold,sold,pi,gpi,thermal2, &
                       S_cc_old,normal,nodalrhs,dSdt,beta0_old,rho0_old,p0_old, &
                       gamma1bar_old,tempbar_init,w0,grav_cell_old,dx,dt,the_bc_tower,mla)

     end do

  else if (restart >= 0) then

     if ( plot_int > 0 ) then

        !------------------------------------------------------------------------
        ! write out a plotfile
        !------------------------------------------------------------------------

        ! we generate this only for the initial plotfile
        if (do_sponge) then
           call make_sponge(sponge,dx,dt,mla)
        end if

        do ip = 1, 2
           if (ip == 2 .and. plt(ip)%plot_int < 0 .and. plt(ip)%plot_dt < 0.0) cycle

           plot_file_name = make_filename(plt(ip)%base_name, restart)

           call make_plotfile(plt(ip),plot_file_name,mla,uold,sold,pi,gpi,rho_omegadot2, &
                              rho_Hnuc2,rho_Hext, &
                              thermal2,S_cc_old,sponge,mla%mba,dx, &
                              the_bc_tower,w0,rho0_old,rhoh0_old,p0_old, &
                              tempbar,gamma1bar_old,etarho_cc, &
                              normal,dt,particles,write_pf_time)

           call write_base_state(restart, plot_file_name, &
                                 rho0_old, rhoh0_old, p0_old, gamma1bar_old, &
                                 w0, etarho_ec, etarho_cc, &
                                 beta0_old, psi, tempbar, tempbar_init, prob_lo(dm))

           call write_job_info(plot_file_name, mla%mba, the_bc_tower, write_pf_time)
        enddo

     end if

  end if

  if (restart < 0) then
     istep = 0
  else
     istep = restart
  end if

  if (stop_time >= 0.d0 .and. time+dt > stop_time) then
     dt = stop_time - time
     if (parallel_IOProcessor()) then
        print*, "Stop time limits dt =",dt
     end if
  end if

  dtold = dt

  if (restart <= 0) then
     ! initialize S_cc_new to the S_cc_old for the initialization or if
     ! restarting from checkpoint0
     do n = 1,nlevs
        call multifab_copy_c(S_cc_new(n),1,S_cc_old(n),1,1)
     end do
  end if

  if (restart < 0) then

     !------------------------------------------------------------------------
     ! Begin the initial pressure iterations
     !------------------------------------------------------------------------

     if (do_sponge) then
        call make_sponge(sponge,dx,dt,mla)
     end if

     if (init_iter > 0) then
        if (parallel_IOProcessor() .and. verbose .ge. 1) then
           print *,''
           write (*,111)
           print*,'DOING',init_iter,'INITIAL PRESSURE ITERATIONS'
           write (*,111)
           print*,''
        end if

        !----------------------------------------------------------------------
        ! main initial iteration loop
        !----------------------------------------------------------------------

        do istep_init_iter = 1,init_iter

           ! Advance a single timestep at all levels.
           init_mode = .true.

           runtime1 = parallel_wtime()

           call advance_timestep(init_mode,mla,uold,sold,unew,snew,gpi,pi,normal, &
                                 rho0_old,rhoh0_old,rho0_new,rhoh0_new,p0_old,p0_new, &
                                 tempbar,gamma1bar_old,gamma1bar_new,w0,rho_omegadot2,rho_Hnuc2, &
                                 rho_Hext,thermal2, &
                                 beta0_old,beta0_new,grav_cell_old,grav_cell_new, &
                                 dx,dt,dtold, &
                                 the_bc_tower,dSdt,S_cc_old,S_cc_new,etarho_ec, &
                                 etarho_cc,psi,sponge,nodalrhs,tempbar_init,particles)

           ! copy pi from snew to sold
           do n=1,nlevs
              call multifab_copy_c(sold(n),pi_comp,snew(n),pi_comp,1,0)
           end do

           runtime2 = parallel_wtime() - runtime1
           call parallel_reduce(runtime1, runtime2, MPI_MAX, proc=parallel_IOProcessorNode())
           if (parallel_IOProcessor()) then
              print*,'Time to advance timestep: ',runtime1,' seconds'
           end if

           call print_and_reset_fab_byte_spread()

        end do ! end do istep_init_iter = 1,init_iter

     end if ! end if (init_iter > 0)

     ! initialize any passively-advected particles
     if (use_particles) then
        call init_particles(particles,sold,rho0_old,rhoh0_old,p0_old,tempbar, &
                            mla,dx,1)

        numparticles = particle_global_numparticles(particles)

        if ( parallel_IOProcessor()) then
           print *,""
           print *,"number of particles initialized = ", numparticles
           print *,""
        endif
     endif

     if ( chk_int > 0 ) then

        !-----------------------------------------------------------------------
        ! write a checkpoint file
        !-----------------------------------------------------------------------

        allocate(chkdata(nlevs))
        do n = 1,nlevs
           call multifab_build(chkdata(n), mla%la(n), 2*dm+nscal, 0)
           call multifab_copy_c(chkdata(n),1                ,uold(n), 1,dm)
           call multifab_copy_c(chkdata(n),rho_comp+dm      ,sold(n), 1,nscal)
           call multifab_copy_c(chkdata(n),rho_comp+dm+nscal,gpi(n),1,dm)
        end do

        check_file_name = make_filename(check_base_name, istep)

        call checkpoint_write(check_file_name, chkdata, &
                              pi, dSdt, S_cc_old, &
                              rho_omegadot2, rho_Hnuc2, rho_Hext, thermal2, &
                              mla%mba%rr, dt)

        call write_base_state(istep, check_file_name, &
                              rho0_old, rhoh0_old, p0_old, gamma1bar_old, &
                              w0, etarho_ec, etarho_cc, &
                              beta0_old, psi, tempbar, tempbar_init, prob_lo(dm))

        call write_aux_data(istep, check_file_name)

        last_chk_written = istep

        do n = 1,nlevs
           call destroy(chkdata(n))
        end do
        deallocate(chkdata)

        if (use_particles) then
           call particle_container_checkpoint(particles,check_file_name,mla)
        endif
     end if

     !-----------------------------------------------------------------------
     ! write a plotfile
     !-----------------------------------------------------------------------

     do ip = 1, 2
        if (ip == 2 .and. plt(ip)%plot_int < 0 .and. plt(ip)%plot_dt < 0.0) cycle

        if ( plt(ip)%plot_int > 0 .or. plt(ip)%plot_dt > ZERO) then

           plot_file_name = make_filename(plt(ip)%base_name, istep)

           call make_plotfile(plt(ip),plot_file_name,mla,uold,sold,pi,gpi,rho_omegadot2, &
                              rho_Hnuc2,rho_Hext, &
                              thermal2,S_cc_old,sponge,mla%mba,dx, &
                              the_bc_tower,w0,rho0_old,rhoh0_old,p0_old, &
                              tempbar,gamma1bar_old,etarho_cc, &
                              normal,dt,particles,write_pf_time)

           call write_base_state(istep, plot_file_name, &
                                 rho0_old, rhoh0_old, p0_old, gamma1bar_old, &
                                 w0, etarho_ec, etarho_cc, &
                                 beta0_old, psi, tempbar, tempbar_init, prob_lo(dm))

           call write_job_info(plot_file_name, mla%mba, the_bc_tower, write_pf_time)
           last_plt_written = istep
        end if

     enddo

  end if ! end if (restart < 0)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Main evolution loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  if (restart < 0) then
     init_step = 1
  else
     init_step = restart+1
  end if

  if ( parallel_IOProcessor()) then
     print*,""
     print*,"BEGIN MAIN EVOLUTION LOOP WITH dt =",dt
     print*,""
  end if

  ! we'd eventually like to read this in from checkpoint if we are restarting
  nuclear_dt_scalefac = 1.d0

  if ( (max_step >= init_step) .and. (time < stop_time .or. stop_time < 0.d0) ) then

     do istep = init_step,max_step

        if (drive_initial_convection .and. istep >= stop_initial_convection) &
             drive_initial_convection = .false.

        if (drive_initial_convection .and. istep >= stop_initial_convection) &
             drive_initial_convection = .false.

        if ( verbose .ge. 1 ) then
           if ( parallel_IOProcessor() ) then
              print *, 'MEMORY STATS AT START OF TIMESTEP ', istep
              print*, ' '
           end if
           call print(multifab_mem_stats(),    "    multifab")
           call print(fab_mem_stats(),         "         fab")
           call print(boxarray_mem_stats(),    "    boxarray")
           call print(layout_mem_stats(),      "      layout")
           call print(boxassoc_mem_stats(),    "    boxassoc")
           call print(fgassoc_mem_stats(),     "     fgassoc")
           call print(syncassoc_mem_stats(),   "   syncassoc")
           call print(copyassoc_mem_stats(),   "   copyassoc")
           call print(fluxassoc_mem_stats(),   "   fluxassoc")
           if ( parallel_IOProcessor() ) print*, ''

           do n = 1,nlevs
              umax = multifab_max_c(uold(n),1)
              umin = multifab_min_c(uold(n),1)
              if (dm.eq.1) then
                 if ( parallel_IOProcessor()) then
                    write(6,1001) n,umax
                    write(6,1011) n,umin
                 end if
              else if (dm .eq. 2) then
                 vmax = multifab_max_c(uold(n),2)
                 vmin = multifab_min_c(uold(n),2)
                 if ( parallel_IOProcessor()) then
                    write(6,1002) n,umax,vmax
                    write(6,1012) n,umin,vmin
                 end if

              else if (dm .eq. 3) then
                 vmax = multifab_max_c(uold(n),2)
                 vmin = multifab_min_c(uold(n),2)
                 wmax = multifab_max_c(uold(n),3)
                 wmin = multifab_min_c(uold(n),3)
                 if ( parallel_IOProcessor()) then
                    write(6,1003) n,umax,vmax,wmax
                    write(6,1013) n,umin,vmin,wmin
                 end if
              end if
           end do
        end if

1001 format('... level / max : old vels   ',i2,2x,e17.10)
1011 format('... level / min : old vels   ',i2,2x,e17.10)
1002 format('... level / max : old vels   ',i2,2x,e17.10,2x,e17.10)
1012 format('... level / min : old vels   ',i2,2x,e17.10,2x,e17.10)
1003 format('... level / max : old vels   ',i2,2x,e17.10,2x,e17.10,2x,e17.10)
1013 format('... level / min : old vels   ',i2,2x,e17.10,2x,e17.10,2x,e17.10)

        !---------------------------------------------------------------------
        ! regrid
        !---------------------------------------------------------------------
        if (max_levs > 1 .and. regrid_int > 0 .and. (mod(istep-1,regrid_int) .eq. 0)) then

           ! Regrid psi, etarho_cc, etarho_ec, and w0.
           ! We do not regrid these in spherical since the base state array only
           ! contains 1 level of refinement.
           ! We do not regrid these if evolve_base_state=F since they are
           ! identically zero, set this way in initialize.
           if (spherical .eq. 0) then
              if (evolve_base_state) then

                 ! copy the coarsest level only, interpolate to all
                 ! the other levels and then copy the valid data from
                 ! the old arrays onto the new this must be done
                 ! before we call init_multilevel or else we lose
                 ! track of where the old valid data was

                 ! copy the coarsest level of the real arrays into the
                 ! temp arrays
                 psi_temp(1,:)        = psi(1,:)
                 etarho_cc_temp(1,:)  = etarho_cc(1,:)
                 etarho_ec_temp(1,:)  = etarho_ec(1,:)
                 w0_temp(1,:)         = w0(1,:)

                 ! piecewise linear interpolation to fill the cc temp arrays
                 do n=2,max_levs
                    do r=0,nr(n)-1
                       if (r .eq. 0 .or. r .eq. nr(n)-1) then
                          psi_temp(n,r) = psi_temp(n-1,r/2)
                          etarho_cc_temp(n,r)  = etarho_cc_temp(n-1,r/2)
                       else
                          if (mod(r,2) .eq. 0) then
                             psi_temp(n,r) = 0.75d0*psi_temp(n-1,r/2) &
                                  + 0.25d0*psi_temp(n-1,r/2-1)
                             etarho_cc_temp(n,r) = 0.75d0*etarho_cc_temp(n-1,r/2) &
                                  + 0.25d0*etarho_cc_temp(n-1,r/2-1)
                          else
                             psi_temp(n,r) = 0.75d0*psi_temp(n-1,r/2) &
                                  + 0.25d0*psi_temp(n-1,r/2+1)
                             etarho_cc_temp(n,r) = 0.75d0*etarho_cc_temp(n-1,r/2) &
                                  + 0.25d0*etarho_cc_temp(n-1,r/2+1)
                          end if
                       end if
                    end do
                 end do

                 ! piecewise linear interpolation to fill the edge-centered temp arrays
                 do n=2,max_levs
                    do r=0,nr(n)
                       if (mod(r,2) .eq. 0) then
                          etarho_ec_temp(n,r) = etarho_ec_temp(n-1,r/2)
                          w0_temp(n,r) = w0_temp(n-1,r/2)
                       else
                          etarho_ec_temp(n,r) = &
                               HALF*(etarho_ec_temp(n-1,r/2)+etarho_ec_temp(n-1,r/2+1))
                          w0_temp(n,r) = HALF*(w0_temp(n-1,r/2)+w0_temp(n-1,r/2+1))
                       end if
                    end do
                 end do

                 ! copy valid data into temp
                 do n=2,nlevs_radial
                    do i=1,numdisjointchunks(n)
                       do r=r_start_coord(n,i),r_end_coord(n,i)
                          psi_temp(n,r)        = psi(n,r)
                          etarho_cc_temp(n,r)  = etarho_cc(n,r)
                       end do
                    end do
                 end do
                 do n=2,nlevs_radial
                    do i=1,numdisjointchunks(n)
                       do r=r_start_coord(n,i),r_end_coord(n,i)+1
                          etarho_ec_temp(n,r) = etarho_ec(n,r)
                          w0_temp(n,r)        = w0(n,r)
                       end do
                    end do
                 end do

                 ! copy temp array back into the real thing
                 psi        = psi_temp
                 etarho_cc  = etarho_cc_temp
                 etarho_ec  = etarho_ec_temp
                 w0         = w0_temp

              else
                 ! evolve_base_state == F and spherical == 0

                 ! Here we want to fill in the rho0 array so there is
                 ! valid data in any new grid locations that are created
                 ! during the regrid.

                 ! copy the coarsest level of the real arrays into the
                 ! temp arrays
                 rho0_temp(1,:) = rho0_old(1,:)

                 ! piecewise linear interpolation to fill the cc temp arrays
                 do n=2,max_levs
                    do r=0,nr(n)-1
                       if (r .eq. 0 .or. r .eq. nr(n)-1) then
                          rho0_temp(n,r) = rho0_temp(n-1,r/2)
                       else
                          if (mod(r,2) .eq. 0) then
                             rho0_temp(n,r) = 0.75d0*rho0_temp(n-1,r/2) &
                                  + 0.25d0*rho0_temp(n-1,r/2-1)
                          else
                             rho0_temp(n,r) = 0.75d0*rho0_temp(n-1,r/2) &
                                  + 0.25d0*rho0_temp(n-1,r/2+1)
                          end if
                       end if
                    end do
                 end do

                 ! copy valid data into temp -- this way we haven't
                 ! overwritten any of the original information.
                 do n=2,nlevs_radial
                    do i=1,numdisjointchunks(n)
                       do r=r_start_coord(n,i),r_end_coord(n,i)
                          rho0_temp(n,r) = rho0_old(n,r)
                       end do
                    end do
                 end do

                 ! don't copy rho0_temp back into the rho0 array just
                 ! yet.  Wait until we regrid

              endif

              ! regardless of evolve_base_state, if new grids were
              ! created, we need to initialize tempbar_init there, in
              ! case drive_initial_convection = T

              ! copy the coarsest level of the real arrays into the
              ! temp arrays
              tempbar_init_temp(1,:)  = tempbar_init(1,:)

              ! piecewise linear interpolation to fill the cc temp arrays
              do n=2,max_levs
                 do r=0,nr(n)-1
                    if (r .eq. 0 .or. r .eq. nr(n)-1) then
                       tempbar_init_temp(n,r) = tempbar_init_temp(n-1,r/2)
                    else
                       if (mod(r,2) .eq. 0) then
                          tempbar_init_temp(n,r) = 0.75d0*tempbar_init_temp(n-1,r/2) &
                               + 0.25d0*tempbar_init_temp(n-1,r/2-1)
                       else
                          tempbar_init_temp(n,r) = 0.75d0*tempbar_init_temp(n-1,r/2) &
                               + 0.25d0*tempbar_init_temp(n-1,r/2+1)
                       end if
                    end if
                 end do
              end do

              ! copy valid data into temp
              do n=2,nlevs_radial
                 do i=1,numdisjointchunks(n)
                    do r=r_start_coord(n,i),r_end_coord(n,i)
                       tempbar_init_temp(n,r) = tempbar_init(n,r)
                    end do
                 end do
              end do

              ! copy temp array back into the real thing
              tempbar_init = tempbar_init_temp

           end if ! end regridding of base state

           ! we can pass as an auxillary tagging quantity either the
           ! energy generate rate (per volume) or tpert
           if (use_tpert_in_tagging) then
              do n = 1, nlevs
                 call multifab_copy_c(tag_mf(n), 1, sold(n), temp_comp, 1)
              enddo

              call put_in_pert_form(mla,tag_mf,tempbar,dx,1, &
                                    foextrap_comp,.true., &
                                    the_bc_tower%bc_tower_array)

           else
              ! figure out if we are tagging off of heating or rxns
              if (do_heating) then
                 do n = 1, nlevs
                    call multifab_copy_c(tag_mf(n), 1, rho_Hext(n), 1, 1)
                 enddo
              else
                 do n = 1, nlevs
                    call multifab_copy_c(tag_mf(n), 1, rho_Hnuc2(n), 1, 1)
                 enddo
              endif
           endif

           do n=1,nlevs
              call multifab_destroy(unew(n))
              call multifab_destroy(snew(n))
              call multifab_destroy(S_cc_new(n))
              if (dm .eq. 3) then
                 call multifab_destroy(normal(n))
              end if
              call multifab_destroy(sponge(n))
              call multifab_destroy(nodalrhs(n))

              call multifab_destroy(pi(n))
              call multifab_destroy(rho_omegadot2(n))
              call multifab_destroy(rho_Hnuc2(n))
              call multifab_destroy(rho_Hext(n))
              call multifab_destroy(thermal2(n))
           end do

           ! create new grids and fill in data on those grids

           call regrid(istep,mla,uold,sold,gpi,dSdt,S_cc_old,dx,the_bc_tower, &
                       rho0_old,rhoh0_old,.false.,tag_mf)

           ! nlevs is local so we need to reset it
           nlevs = mla%nlevel

           if (nlevs .ne. max_levs) then
              call bl_error('varden.f90: nlevs .ne. max_levs not supported yet')
           end if

           call init_multilevel(sold)

           do n = 1,nlevs
              call multifab_build(         unew(n), mla%la(n),    dm, nghost(uold(n)))
              call multifab_build(         snew(n), mla%la(n), nscal, nghost(sold(n)))
              call multifab_build(     S_cc_new(n), mla%la(n),     1, 0)
              if (dm .eq. 3) then
                 call multifab_build(normal(n), mla%la(n),    dm, 1)
              end if
              call multifab_build(       sponge(n), mla%la(n),     1, 0)
              call multifab_build(      nodalrhs(n), mla%la(n),     1, 0, nodal)

              call multifab_build(           pi(n), mla%la(n),     1, 0, nodal)
              call multifab_build(rho_omegadot2(n), mla%la(n), nspec, 0)
              call multifab_build(    rho_Hnuc2(n), mla%la(n),     1, 0)
              call multifab_build(     rho_Hext(n), mla%la(n),     1, 0)
              call multifab_build(     thermal2(n), mla%la(n),     1, 1)

              call setval(      unew(n), ZERO, all=.true.)
              call setval(      snew(n), ZERO, all=.true.)
              call setval(  S_cc_new(n), ZERO, all=.true.)
              call setval(    sponge(n), ONE,  all=.true.)
              call setval(   nodalrhs(n), ZERO, all=.true.)
           end do

           ! Create normal now that we have defined center and dx
           call make_normal(normal,dx)

           if (evolve_base_state) then
              ! force rho0 to be the average of rho
              call average(mla,sold,rho0_old,dx,rho_comp)

           else

              if (spherical .eq. 0) then
                 ! copy the old base state density with piecewise linear
                 ! interpolated data in the new positions -- this is
                 ! only necessary for evolve_base_state = F and
                 ! spherical = 0.
                 rho0_old = rho0_temp

                 ! zero out any data where there is no corresponding full
                 ! state array
                 do n=2,nlevs_radial
                    do i=1,numdisjointchunks(n)
                       if (i .eq. numdisjointchunks(n)) then
                          do r=r_end_coord(n,i)+1,nr(n)-1
                             rho0_old(n,r) = 0.d0
                          end do
                       else
                          do r=r_end_coord(n,i)+1,r_start_coord(n,i+1)-1
                             rho0_old(n,r) = 0.d0
                          end do
                       end if
                    end do
                 end do

              endif

           endif

           ! recompute p0 based on the new rho0
           call compute_cutoff_coords(rho0_old)

           call make_grav_cell(grav_cell_old,rho0_old)

           ! enforce HSE
           call enforce_HSE(rho0_old,p0_old,grav_cell_old)

           if (use_tfromp) then
              ! compute full state T = T(rho,p0,X)
              call makeTfromRhoP(sold,p0_old,mla,the_bc_tower%bc_tower_array,dx)
           else
              ! compute full state T = T(rho,h,X)
              call makeTfromRhoH(sold,p0_old,mla,the_bc_tower%bc_tower_array,dx)
           end if

           ! force tempbar to be the average of temp
           call average(mla,sold,tempbar,dx,temp_comp)

           ! gamma1bar needs to be recomputed
           call make_gamma1bar(mla,sold,gamma1bar_old,p0_old,dx)

           ! beta0_old needs to be recomputed
           call make_beta0(beta0_old,rho0_old,p0_old,gamma1bar_old,grav_cell_old)

           ! redistribute the particles to their new processor locations
           if (use_particles) call redistribute(particles,mla,dx,prob_lo)

        end if ! end regridding

        !---------------------------------------------------------------------
        ! get the new timestep
        !---------------------------------------------------------------------
        dtold = dt

        if (istep > 1) then

           call estdt(mla,the_bc_tower,uold,sold,gpi,S_cc_old,dSdt, &
                      w0,rho0_old,p0_old,gamma1bar_old,grav_cell_old,dx,cflfac,dt)

           if (parallel_IOProcessor() .and. verbose .ge. 1) then
              print*,''
              print*,"Call to estdt at beginning of step",istep
              print*,"gives dt =",dt
           end if

           if (nuclear_dt_scalefac .lt. 1.d0) then
              dt = nuclear_dt_scalefac*dtold
              if (parallel_IOProcessor()) then
                 print*, "Applying nuclear_dt_fac; new dt=",dt
              end if
           end if

           if(dt .gt. max_dt_growth*dtold) then
              dt = max_dt_growth*dtold
              if (parallel_IOProcessor() .and. verbose .ge. 1) then
                 print*,'dt_growth factor limits the new dt =',dt
              end if
           end if

           if(dt .gt. max_dt) then
              if (parallel_IOProcessor() .and. verbose .ge. 1) then
                 print*,'max_dt limits the new dt =',max_dt
              end if
              dt = max_dt
           end if

           if(fixed_dt .ne. -1.0d0) then
              dt = fixed_dt
              if (parallel_IOProcessor()) then
                 print*, "Setting fixed dt =",dt
              end if
           end if

           if (parallel_IOProcessor() .and. verbose .ge. 1) then
              print*,''
           end if

        end if

        if (stop_time >= 0.d0) then
           if (time+dt > stop_time) then
              dt = stop_time - time
              if (parallel_IOProcessor()) then
                 print*, "Stop time limits dt =",dt
              end if
           end if
        end if

        !---------------------------------------------------------------------
        ! Advance a single timestep at all levels.
        !---------------------------------------------------------------------
        init_mode = .false.
        if (do_sponge) then
           call init_sponge(rho0_old(1,:))
           call make_sponge(sponge,dx,dt,mla)
        end if
        runtime1 = parallel_wtime()

        call advance_timestep(init_mode,mla,uold,sold,unew,snew,gpi,pi,normal,rho0_old, &
                              rhoh0_old,rho0_new,rhoh0_new,p0_old,p0_new,tempbar,gamma1bar_old, &
                              gamma1bar_new,w0,rho_omegadot2,rho_Hnuc2,rho_Hext,thermal2, &
                              beta0_old,beta0_new, &
                              grav_cell_old,grav_cell_new, &
                              dx,dt,dtold,the_bc_tower,dSdt,S_cc_old, &
                              S_cc_new,etarho_ec,etarho_cc,psi,sponge,nodalrhs,tempbar_init, &
                              particles)

        ! limit the timestep if the temperature is changing too rapidly
        if (nuclear_dt_fac .gt. 0.d0) then
           smaxold = 0.d0
           smax    = 0.d0
           do n=1,nlevs
              smaxold = max(smaxold,multifab_norm_inf_c(sold(n),temp_comp,1))
              smax    = max(smax   ,multifab_norm_inf_c(snew(n),temp_comp,1))
           end do
           if (smax .gt. smaxold) then
              nuclear_dt_scalefac = nuclear_dt_fac*( smaxold / (smax-smaxold) )
           else
              nuclear_dt_scalefac = 1.d0
           end if
        end if

        runtime2 = parallel_wtime() - runtime1
        call parallel_reduce(runtime1, runtime2, MPI_MAX, proc = parallel_IOProcessorNode())
        if (parallel_IOProcessor()) print*,'Time to advance timestep: ',runtime1,' seconds'

        call print_and_reset_fab_byte_spread()

        time = time + dt

        if ( verbose .ge. 1 ) then

           if ( parallel_IOProcessor() ) then
              print *, 'MEMORY STATS AT END OF TIMESTEP ', istep
              print*, ' '
           end if
           call print(multifab_mem_stats(),    "    multifab")
           call print(fab_mem_stats(),         "         fab")
           call print(boxarray_mem_stats(),    "    boxarray")
           call print(layout_mem_stats(),      "      layout")
           call print(boxassoc_mem_stats(),    "    boxassoc")
           call print(fgassoc_mem_stats(),     "     fgassoc")
           call print(syncassoc_mem_stats(),   "   syncassoc")
           call print(copyassoc_mem_stats(),   "   copyassoc")
           call print(fluxassoc_mem_stats(),   "   fluxassoc")
           if ( parallel_IOProcessor() ) print*, ''

           do n = 1,nlevs
              if (parallel_IOProcessor()) write(6,1100) n

              smin = multifab_min_c(unew(n),1)
              smax = multifab_max_c(unew(n),1)
              if (parallel_IOProcessor()) write(6,1101) smin,smax

              if (dm .ge. 2) then
                 smin = multifab_min_c(unew(n),2)
                 smax = multifab_max_c(unew(n),2)
                 if (parallel_IOProcessor()) write(6,1102) smin,smax
              end if

              if (dm .eq. 3) then
                 smin = multifab_min_c(unew(n),3)
                 smax = multifab_max_c(unew(n),3)
                 if (parallel_IOProcessor()) write(6,1103) smin,smax
              end if

              smin = multifab_min_c(gpi(n),1)
              smax = multifab_max_c(gpi(n),1)
              if (parallel_IOProcessor()) write(6,1104) smin,smax

              if (dm .ge. 2) then
                 smin = multifab_min_c(gpi(n),2)
                 smax = multifab_max_c(gpi(n),2)
                 if (parallel_IOProcessor()) write(6,1105) smin,smax
              end if

              if (dm .eq. 3) then
                 smin = multifab_min_c(gpi(n),3)
                 smax = multifab_max_c(gpi(n),3)
                 if (parallel_IOProcessor()) write(6,1106) smin,smax
              end if

              if (parallel_IOProcessor()) write(6,1107)
           end do
        end if

1100    format('At level ',i3)
1101    format('... min/max : x-velocity       ',e17.10,2x,e17.10)
1102    format('... min/max : y-velocity       ',e17.10,2x,e17.10)
1103    format('... min/max : z-velocity       ',e17.10,2x,e17.10)
1104    format('... min/max : gpix              ',e17.10,2x,e17.10)
1105    format('... min/max : gpiy              ',e17.10,2x,e17.10)
1106    format('... min/max : gpiz              ',e17.10,2x,e17.10)
1107    format(' ')

        if (parallel_IOProcessor()) then
           write(6,1000) istep,time,dt
           write(6,*)
        end if

        !---------------------------------------------------------------------
        ! save the old velocity, scalars and source terms for the next step
        !---------------------------------------------------------------------

        do n = 1,nlevs
           call multifab_copy_c(    uold(n),1,    unew(n),1,dm,   nghost(uold(n)))
           call multifab_copy_c(    sold(n),1,    snew(n),1,nscal,nghost(uold(n)))
           call multifab_copy_c(S_cc_old(n),1,S_cc_new(n),1,    1)
        end do

        ! update beta0, grav_cell, gamma1bar
        beta0_old = beta0_new
        grav_cell_old = grav_cell_new
        gamma1bar_old = gamma1bar_new

        ! Copy the base state
        rho0_old = rho0_new
        rhoh0_old = rhoh0_new
        p0_old = p0_new

        !---------------------------------------------------------------------
        ! output
        !---------------------------------------------------------------------

        ! output any particle information
        if (use_particles) then
           if (store_particle_vels) then
              call timestamp(particles, 'timestamp', sold, index_partdata, &
                             names_partdata, time, uold)
           else
              call timestamp(particles, 'timestamp', sold, index_partdata, &
                             names_partdata, time)
           endif
        endif

        ! if the file .dump_checkpoint exists in our output directory, then
        ! automatically dump a plotfile
        inquire(file=".dump_checkpoint", exist=dump_checkpoint)

        if (chk_int > 0 .or. dump_checkpoint) then
           if (mod(istep,chk_int) .eq. 0 .or. dump_checkpoint) then

              ! write out any buffered diagnostic information
              call flush_diag()

              allocate(chkdata(nlevs))
              do n = 1,nlevs
                 call multifab_build(chkdata(n), mla%la(n), 2*dm+nscal, 0)
                 call multifab_copy_c(chkdata(n),1                ,unew(n),1,dm)
                 call multifab_copy_c(chkdata(n),rho_comp+dm      ,snew(n),1,nscal)
                 call multifab_copy_c(chkdata(n),rho_comp+dm+nscal,gpi(n) ,1,dm)
              end do

              check_file_name = make_filename(check_base_name, istep)

              call checkpoint_write(check_file_name, chkdata, &
                                    pi, dSdt, S_cc_new, &
                                    rho_omegadot2, rho_Hnuc2, rho_Hext, &
                                    thermal2, mla%mba%rr, &
                                    dt)

              call write_base_state(istep, check_file_name, &
                                    rho0_new, rhoh0_new, p0_new, gamma1bar_new, &
                                    w0, etarho_ec, etarho_cc, &
                                    beta0_old, psi, tempbar, tempbar_init, prob_lo(dm))

              call write_aux_data(istep, check_file_name)

              last_chk_written = istep

              do n = 1,nlevs
                 call destroy(chkdata(n))
              end do
              deallocate(chkdata)

              if (use_particles) then
                 call particle_container_checkpoint(particles,check_file_name,mla)
              endif

           end if
        end if

        ! if the file .dump_plotfile exists in our output directory, then
        ! automatically dump a plotfile
        inquire(file=".dump_plotfile", exist=dump_plotfile)

        do ip = 1, 2
           if (ip == 2 .and. plt(ip)%plot_int < 0 .and. plt(ip)%plot_dt < 0.0) cycle

           if ( plt(ip)%plot_int > 0 .or. plt(ip)%plot_dt > ZERO .or. dump_plotfile) then
              if ( (plt(ip)%plot_int > 0 .and. mod(istep,plt(ip)%plot_int) .eq. 0) .or. &
                   (plt(ip)%plot_dt > ZERO .and. &
                   mod(time - dt,plt(ip)%plot_dt) > mod(time,plt(ip)%plot_dt)) .or. &
                   dump_plotfile) then

                 plot_file_name = make_filename(plt(ip)%base_name, istep)

                 call make_plotfile(plt(ip),plot_file_name,mla,unew,snew,pi,gpi,rho_omegadot2, &
                                    rho_Hnuc2,rho_Hext, &
                                    thermal2,S_cc_new,sponge,mla%mba,dx, &
                                    the_bc_tower,w0,rho0_new,rhoh0_new,p0_new, &
                                    tempbar,gamma1bar_new,etarho_cc, &
                                    normal,dt,particles,write_pf_time)

                 call write_base_state(istep, plot_file_name, &
                                       rho0_new, rhoh0_new, p0_new, gamma1bar_new, &
                                       w0, etarho_ec, etarho_cc, &
                                       beta0_old, psi, tempbar, tempbar_init, prob_lo(dm))

                 call write_job_info(plot_file_name, mla%mba, the_bc_tower, write_pf_time)
                 last_plt_written = istep
              end if
           end if
        enddo

        ! if the file .abort_maestro exists in our output directory, then
        ! automatically end the run.  This has the effect of also dumping
        ! a final checkpoint file.
        inquire(file=".abort_maestro", exist=abort_maestro)
        if (abort_maestro) exit

        ! have we reached the stop time?
        if (stop_time >= 0.d0) then
           if (time >= stop_time) goto 999
        end if

     end do  ! end main evolution loop

999  continue
     if (istep > max_step) istep = max_step

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! write the final checkpoint and plotfile
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

1000 format('STEP = ',i7,1x,' TIME = ',es16.10,1x,' DT = ',es16.10)

     if ( chk_int > 0 .and. last_chk_written .ne. istep ) then
        !       This writes a checkpoint file.

        ! write out any buffered diagnostic information
        call flush_diag()

        allocate(chkdata(nlevs))
        do n = 1,nlevs
           call multifab_build(chkdata(n), mla%la(n), 2*dm+nscal, 0)
           call multifab_copy_c(chkdata(n),1,unew(n),1,dm)
           call multifab_copy_c(chkdata(n),rho_comp+dm,snew(n),1,nscal)
           call multifab_copy_c(chkdata(n),rho_comp+dm+nscal,gpi(n),1,dm)
        end do

        check_file_name = make_filename(check_base_name, istep)

        call checkpoint_write(check_file_name, chkdata, &
                              pi, dSdt, S_cc_new, &
                              rho_omegadot2, rho_Hnuc2, rho_Hext, thermal2, &
                              mla%mba%rr, dt)

        call write_base_state(istep, check_file_name, &
                              rho0_new, rhoh0_new, p0_new, gamma1bar_new, &
                              w0, etarho_ec, etarho_cc, &
                              beta0_old, psi, tempbar, tempbar_init, prob_lo(dm))

        call write_aux_data(istep, check_file_name)

        do n = 1,nlevs
           call destroy(chkdata(n))
        end do
        deallocate(chkdata)

        if (use_particles) then
           call particle_container_checkpoint(particles,check_file_name,mla)
        endif

     end if

     do ip = 1, 2
        if (ip == 2 .and. plt(ip)%plot_int < 0 .and. plt(ip)%plot_dt < 0.0) cycle

        if ( ( plt(ip)%plot_int > 0 .or. &
             (plt(ip)%plot_dt > ZERO .and. &
             mod(time - dt,plt(ip)%plot_dt) > mod(time,plt(ip)%plot_dt)) ) .and. &
             last_plt_written .ne. istep ) then

           plot_file_name = make_filename(plt(ip)%base_name, istep)

           call make_plotfile(plt(ip),plot_file_name,mla,unew,snew,pi,gpi,rho_omegadot2, &
                              rho_Hnuc2,rho_Hext, &
                              thermal2,S_cc_new,sponge,mla%mba,dx, &
                              the_bc_tower,w0,rho0_new,rhoh0_new,p0_new, &
                              tempbar,gamma1bar_new,etarho_cc, &
                              normal,dt,particles,write_pf_time)

           call write_base_state(istep, plot_file_name, &
                                 rho0_new, rhoh0_new, p0_new, gamma1bar_new, &
                                 w0, etarho_ec, etarho_cc, &
                                 beta0_old, psi, tempbar, tempbar_init, prob_lo(dm))

           call write_job_info(plot_file_name, mla%mba, the_bc_tower, write_pf_time)
        end if
     enddo
  endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! clean-up
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  do n=1,nlevs
     call destroy(unew(n))
     call destroy(snew(n))
     call destroy(S_cc_new(n))
     if(dm .eq. 3) then
        call destroy(normal(n))
     end if
     call destroy(sponge(n))
     call destroy(nodalrhs(n))

     call destroy(uold(n))
     call destroy(sold(n))
     call destroy(gpi(n))
     call destroy(dSdt(n))
     call destroy(S_cc_old(n))

     call destroy(pi(n))
     call destroy(rho_omegadot2(n))
     call destroy(rho_Hnuc2(n))
     call destroy(rho_Hext(n))
     call destroy(thermal2(n))

     call destroy(tag_mf(n))
  end do

  call destroy(mla)

  call bc_tower_destroy(the_bc_tower)

  call destroy(particles)

  call destroy_geometry()

  call eos_finalize()
  call network_finalize()

  call model_parser_finalize()
  call simple_log_finalize()

  call diag_finalize()

  call runtime_close()

  deallocate(uold,sold,pi,gpi,dSdt,S_cc_old,S_cc_new,rho_omegadot2)
  deallocate(rho_Hnuc2,rho_Hext,thermal2,tag_mf)
  deallocate(dx,beta0_old,beta0_new,gamma1bar_old,gamma1bar_new,s0_init,rho0_old)
  deallocate(rhoh0_old,rho0_new,rhoh0_new,p0_init,p0_old,p0_new,w0,etarho_ec,etarho_cc)
  deallocate(psi,tempbar,tempbar_init,grav_cell_old,grav_cell_new)

end subroutine varden

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO