!> @dir OrbitingParticles
!! @brief Contains files necessary for the Orbiting Particles test problem

!> @file problem.f90
!! @brief Main file for module Problem

!> @defgroup OrbitingParticles Orbiting Particle Module
!! @brief Module for setting up orbiting particles
!! @ingroup Modules

!> Orbiting Particle Module
!! @ingroup OrbitingParticles
MODULE Problem
  USE Outflows
  USE DataDeclarations
  USE ParticleDeclarations
  USE CoolingSrc
!  USE SourceDeclarations
!  USE Projections!n
  USE Fields
  USE Ambients!n
  USE Refinements
  USE ProcessingDeclarations
  IMPLICIT NONE
  SAVE
  PUBLIC ProblemModuleInit, ProblemGridInit, &
       ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
  INTEGER :: itracer2, itracer3, itracer4, itracer5
  REAL(KIND=qPREC) :: namb, njet, tamb, Rjet, tjet, vjet,torusalpha,torusbeta, tt, sigma, MDOT_wind, MDOT_jet, rho_jet, rho_wind, vwind, wind_temp, jet_width, open_angle, sigmaFACTOR, tf, balick_a, balick_aa, balick_b, balick_bb, balick_c
  REAL(KIND=qPREC), SAVE :: delta_mass, delta_mass_3D, local_mass, local_mass_3D
  LOGICAL :: lObliqueCone ! for turn on the 45 deg thing for specified taper, suggested by Bruce Balick 
  LOGICAL :: lPulsed, lAGBcenter !for M2-9
  REAL(KIND=qPrec) :: PulsedPeriod, Omega_Pulsed, orbitalSeparation, rStratified
  INTEGER :: PulseModule
  INTEGER :: StraTorus ! Setting stratified and torus, 1 for Martin's way, 2 for Baowei's way 
  INTEGER :: buff(0:MaxDepth) !number of extra cells of buffer around sink region to set err flags

  LOGICAL :: lCooling, rings, stratified, torus,lform
  TYPE(ParticleDef), POINTER :: Particle
  INTEGER :: outflowType=0
  INTEGER, PARAMETER :: collimated=1, conical=2, clump=3
  INTEGER :: oldframe = -1 
  REAL(KIND=qPrec) :: total_mass, total_mass_3D, total_delta_mass, total_delta_mass_3D, mass_flux, mass_flux_3D
  
  Private
     REAL(KIND=qPRec) :: WhiteDwarfPos(3), AGBstarPos(3)
CONTAINS

  !> Initializes module variables
   SUBROUTINE ProblemModuleInit()      
      REAL(KIND=qPREC) :: position(3)
      REAL(KIND=qPREC) :: vel(3)
      TYPE(RefinementDef), POINTER :: Refinement
      TYPE(AmbientDef), POINTER :: Ambient
      INTEGER :: i
      !	  TYPE(ProjectionDef), POINTER :: Projection
      NAMELIST /ProblemData/ outflowType, namb, njet, tamb, Rjet, &
           tjet, vjet, lCooling, rings, stratified, &
           torus,torusalpha,torusbeta, tt, sigma, MDOT_wind, MDOT_jet, &
           vwind, buff, wind_temp, jet_width, open_angle, tf,lObliqueCone,&
            lAGBcenter, orbitalSeparation, rStratified, lPulsed, PulsedPeriod, PulseModule, StraTorus, balick_a,&
           balick_aa, balick_b, balick_bb, balick_c

      delta_mass=0
      delta_mass_3d=0
      local_mass=0
      local_mass_3d=0
      OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
      READ(PROBLEM_DATA_HANDLE,NML=ProblemData)

      OPEN(UNIT=745, FILE='out/stdout.txt',STATUS="unknown")
      WRITE(745, *) "!======= Input Data from problem.data ======="
      WRITE(745, NML=ProblemData)
      WRITE(745, *) "!============================================"
      WRITE(745, *)
      WRITE(745, *)
      WRITE(745, *)
      WRITE(745, *) "!===== 3D Grid Mass and Mass Flow Rate: ======" 

      WhiteDwarfPos=(/0.0,0.0,0.0/)
      AGBstarPos=(/(orbitalSeparation*au/lScale),0d0,0d0/)

      print*, "AGBstarPos=", AGBstarPos

      !Pulsed frequency in C.U.
      Omega_Pulsed=2*Pi/(PulsedPeriod*yr/TimeScale) 

      CALL CreateAmbient(Ambient)
      CALL AddDiagnosticVar(CoolingStrength_Field)
      CALL AddTracer(itracer2, 'DM cooling')
      CALL AddTracer(itracer3, 'Clump tracer')
      CALL AddTracer(itracer4, 'Jet tracer')
      CALL AddTracer(itracer5, 'DensityFlag tracer')
      
      if (balick_a==0d0) balick_aa=0d0 ;  if (balick_b==0d0) balick_bb=0d0
      if (balick_a==0d0  .and.  balick_b==0d0) balick_c=0d0  

IF ( MPI_ID == 0 .and.&
     tt.ne.0d0   .and.&
     (balick_a*1d0**balick_aa + balick_b*1d0**balick_bb)*exp(-balick_c*1d0).ne.1d0 ) then
       print*,(balick_a*1d0**balick_aa + balick_b*1d0**balick_bb)*exp(-balick_c*1d0)
       print*,'' ; print*,''
       print*,'Use either tt (diff than 0) or balicks (diff than 0), but not both.'
       stop
END IF


IF (MPI_ID == 0 .and.&
    nDim==2     .and.&
	 lMHD        .and.&
	 iCylindrical.ne.2) THEN
       print*,'' ; print*,''
       print*,'ndim=2 & lMHD=.ture. MUST go with iCylindrical=2, please change and try again.'
       stop
END IF

IF (MPI_ID == 0 .and.&
    nDim==3     .and.&
	 iCylindrical.ne.0) THEN
       print*,'' ; print*,''
       print*,'3D sims MUST have iCylindrical=0, please change and try again.'
       stop
END IF

IF (MPI_ID == 0 .and.&
    nDim==3     .and.&
	 sigma.ne.0d0     ) then
       print*,'' ; print*,''
       print*,'ndim=3 & sigma is not 0. 3D MHD has not been testes. Please change signma to 0 and try again.'
       stop
END IF


      jet_width		 = Rjet 
      tf              = tf*Pi/180d0
      open_angle      = open_angle*Pi/180d0
      MDOT_wind       = MDOT_wind*(mSolar/yr)/& !from Msun/yr to gr/sec
								(mScale/TimeScale) !from real units to computational ones
      MDOT_jet        = MDOT_jet*(mSolar/yr)/& !ditto
								(mScale/TimeScale)

      if (vwind.ne.0d0) then
		   rho_wind= MDOT_wind/ &! in comp units
             (4d0*Pi*Rjet**2*(Vwind/VelScale)) !Rjet is in comp units already
          ! print *, "rho_wind=", rho_wind
          ! stop
      else
		   rho_wind= MDOT_wind/ &! in comp units
				 (4d0*Pi*Rjet**2*(1e6/VelScale)) !Rjet is in comp units already
          ! print *, "else rho_wind=", rho_wind
          ! stop
      end if
										 
      rho_jet= MDOT_jet/ &! ditto
           (4d0*Pi*Rjet**2*(Vjet/VelScale))

		sigmaFACTOR=dsqrt(sigma*4d0*Pi*njet/nscale*(vjet/velscale)**2d0)
		!sigmaFACTOR=dsqrt(sigma*4d0*Pi*rho_jet*(vjet/velscale)**2d0)

      CLOSE(PROBLEM_DATA_HANDLE)      

      CALL CreateRefinement(Refinement)
      CALL CreateShape(Refinement%Shape, RECTANGULAR_PRISM, (/Rjet,Rjet,Rjet/))
      Refinement%BufferCells=buff(0)

      !CALL CreateProjection(Projection)
      !  Projection%Field%id=CoolingStrength_Field
      !  Projection%Field%component=GASCOMP
      !  Projection%Field%name='Cooling_integrated_emiss'
      !  Projection%dim=3d0

   END SUBROUTINE ProblemModuleInit

  !> Applies initial conditions
  !! @param Info Info object
   SUBROUTINE ProblemGridInit(Info)
      TYPE(InfoDef) :: Info
      INTEGER :: i,j,k,l,m,ii,jj,kk
      INTEGER :: rmbc,zrmbc,level
      INTEGER :: mx, my, mz
      INTEGER :: iErr
      LOGICAL, SAVE :: once=.true.
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:)     :: q, aux
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: density
      REAL (KIND=qPrec), SAVE :: dxFIXED
      REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,r,r2,r_sc,dr,dt,r3,torus2,cos2theta,ex1,ex2,Ftheta,f,dec7, bmax
      LOGICAL, DIMENSION(:,:,:),ALLOCATABLE::  density_flag
      if (lRestart) return

      level=Info%level
      q=>Info%q
      aux=>Info%aux

      ! Calculating the number of ghost cells on each side of the grid.
      rmbc=levels(level)%gmbc(levels(level)%step)
      mx = Info%mX(1)
      my = Info%mX(2)
      mz = Info%mX(3)
      dx=levels(level)%dX
      dy=dx
      dz=dx
      xl=Info%xBounds(1,1)
      yl=Info%xBounds(2,1)
      zl=Info%xBounds(3,1)
      
      dt=levels(level)%dt
      SELECT CASE(nDim)
      CASE(2)
         zrmbc=0;mz=1;dz=0d0;zl=0d0
      CASE(3)
         zrmbc=rmbc
      END SELECT

      ALLOCATE(density_flag(1-rmbc:mx+rmbc, 1-rmbc:my+rmbc, 1-zrmbc:mz+zrmbc))

      ALLOCATE(density(1-rmbc:mx+rmbc, 1-rmbc:my+rmbc, 1-zrmbc:mz+zrmbc))
      density_flag=.true.
      density=Info%q(1-rmbc:mx+rmbc, 1-rmbc:my+rmbc, 1-zrmbc:mz+zrmbc, 1)


      SELECT CASE(StraTorus) ! 1 for Martin's version; 2 for Baowei/Bruce's version 
      CASE(1)
!!=====================================================================
!!
!! The following code was from Martin and used before Jan 6th 2015 
!!
!!======================================================================
      DO i=1-rmbc, mx+rmbc         ; x = (xl+(REAL(i,xPrec)-half)*dx)
         DO j=1-rmbc, my+rmbc      ; y = (yl+(REAL(j,xPrec)-half)*dy)
            DO k=1-zrmbc, mz+zrmbc ; z = (zl+(REAL(k,xPrec)-half)*dz)
               r =SQRT( x**2 + y**2 + z**2 ) ! from origin
               r2=SQRT( (y-2d0*Rjet)**2 + x**2 + z**2 ) ! from clump center
               dr =SQRT( dx**2 + dy**2 + dz**2 )

               if (outflowType==clump) r=r2 !b4 27jan'12

               !A M B I E N T 
               q(i,j,k,1) = namb/nScale !density

               if (r .le. Rjet) then
                  density_flag(i,j,k)=.false.
               end if

               if (stratified) then !27jan'12
                  if (outflowType /= clump) then
                     q(i,j,k,1) = q(i,j,k,1)/&
                          (r+1d0+.25d0*dr)**2
                  else
                     if (r.gt.Rjet) q(i,j,k,1) = q(i,j,k,1)/&
                          (r-Rjet+1d0+.25d0*dr)**2
                  end if
               end if


               !
!!! TORUS
               !IF (r > rjet) THEN
               !From Frnak & Mellema, 1994ApJ...430..800F:
               !produced a pole-to-equator density contrast
               torus2=1d0
               if ( torus ) then
                  if ( r.eq.0d0   .or.   (y-2d0*Rjet).eq.0d0) then
                     print*,'torus error';stop
                  else
                     !cos2theta=x**2/( x**2 + y**2 + z**2 )-z**2/x**2 !-------------A DIDNT WORK
                     IF (outflowType == clump) THEN
                        cos2theta=2d0*acos((y-2d0*Rjet)/r)!----------------------------B
                     ELSE 
                        cos2theta=2d0*acos((y-0d0*Rjet)/r)
                     END IF
                     !cos2theta=(x-2d0*Rjet)**2/r**2-z**2/(x-2d0*Rjet)**2!-----C DIDNT WORK
                  end if
                  ex1=-2d0*torusbeta
                  !ex2=torusbeta*cos2theta-torusbeta !---------------------------A DIDNT WORK
                  ex2=torusbeta*cos(cos2theta)-torusbeta!-----------------------B
                  !ex2=torusbeta*cos2theta-torusbeta !----------------------C DIDNT WORK
                  Ftheta=1d0-(  torusalpha*(dexp(ex2)-1d0)/(dexp(ex1)-1d0)  )
                  !
!!!
                  !1)
                  !torus2=Rjet**2/r**2/Ftheta 
                  !2)
                  torus2=1d0/Ftheta
!!!
                  !
               end if
               		!print*,'tours2=',torus2
               q(i,j,k,1)=q(i,j,k,1)*torus2
!!!
               !																				

!!!!!!!!!!RINGS!!!!!!!!!!!!!!
               f=0.5d0

               dec7=2.5d0/1.5d0  
               do ii=1,1000
                  f=f+exp(-(6d0*(r-dec7))**2)
                  dec7=dec7+1.5d0/1.5d0
               end do

               !dec7=3d0/1.5d0  
               !do ii=1,1000
               !   f=f+exp(-(10d0*(r-dec7))**2)
               !	dec7=dec7+.8d0/1.5d0
               !end do
               f=2d0*f
               if (rings) q(i,j,k,1) = q(i,j,k,1)*f
!!!!!!!!!!RINGS!!!!!!!!!!!!!!
               !END IF

               q(i,j,k,imom(1:nDim)) = 0d0
               q(i,j,k,iE) = q(i,j,k,1)*tamb/TempScale/gamma1
               q(i,j,k,itracer3)=0d0 ; q(i,j,k,itracer4)=0d0
!            q(i,j,k,itracer2)=GetCoolingStrength(q(i,j,k,:),lform)*pScale/TimeScale !erg/cc/s

!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
! C L U M P:
               if (outflowType == clump .AND. r2.le.Rjet) then !Martin version
          !     if (outflowType == clump .AND. r2.lt.Rjet) then !Baowei 1.2.2015
			
               !			q(i,j,k,itracer3)=1d0 !martin's original set up for clump tracer
                  q(i,j,k,1) = namb/nScale+&
                       njet/nscale*(1d0-(r2/Rjet)**2) ![cu]
               
     			q(i,j,k,itracer3)=q(i,j,k,1) !set the clump tracer as density

                  q(i,j,k,3) = vjet*q(i,j,k,1)/velscale ![cu]
                  q(i,j,k,iE) = q(i,j,k,1)*tjet/TempScale/gamma1 +& !internal enegy     [cu]
						   half*q(i,j,k,3)**2/q(i,j,k,1)
                  !if (ndim.eq.2   .and.   iCylindrical.ne.2) then 
                  !   q(i,j,k,iE) =q(i,j,k,iE) + &
					!		   half*SUM(q(i,j,k,2:3)**2)/q(i,j,k,1)            !kinetic energy
                  !else
                  !   q(i,j,k,iE) =q(i,j,k,iE) + &
                  !        half*SUM(q(i,j,k,2:4)**2)/q(i,j,k,1)            !kinetic energy
                  !end if

               end if!r<Rjet
               
            END DO
         END DO
      END DO

      aux(:,:,:,:)=0d0 
      q(:,:,:,iBx:iBz) = 0d0


      if (outflowType==clump .and. lMHD) then
		if(nDim==3 .and. sigma.ne.0d0 .and. MPI_ID == 0) THEN
         print*,''
         print*,'3D magnetic clump not implemented yet';stop
         END IF

			if (once) then ; once=.false. ; dxFIXED=dx/(2d0**3d0)
			end if

           DO i=1-rmbc, mx+rmbc  ; x = (xl+(REAL(i,xPrec)-half)*dx) 
           DO j=1-rmbc, my+rmbc  ; y = (yl+(REAL(j,xPrec)-half)*dy)
           DO k=1-zrmbc,mz+zrmbc ; z = (zl+(REAL(k,xPrec)-half)*dz)
				 r=SQRT( (y-2d0*Rjet)**2 + x**2 + z**2 ) ! from clump center
             if (r>Rjet) cycle

         	bmax=        ((3e12/lScale)/dxFIXED   )**2d0*&
         	   ((dxFIXED     /(3e12/lScale))-1d0)

              if (nDim==2) then
         		  q(i,j,k,iBz) = ((3e12/lScale)/sqrt(x**2d0))**2d0*&
         		     ((sqrt(x**2d0)/(3e12/lScale))-1d0)*& !B
         			  sigmaFACTOR/bmax
              else if (nDim==3) then
!NOT IMPLEMENTED IN TERMA OF AUX YET
         		  q(i,j,k,iBx) = ((3e12/lScale)/sqrt(x**2d0))**2d0*&
         		     ((sqrt(x**2d0)/(3e12/lScale))-1d0)*& !B
         			  sigmaFACTOR/bmax
         	
         		  q(i,j,k,iBz) = ((3e12/lScale)/sqrt(z**2d0))**2d0*&
         		     ((sqrt(z**2d0)/(3e12/lScale))-1d0)*& !B
         			  sigmaFACTOR/bmax
         	
				  end if
              if (q(i,j,k,iBz)>1.e8   .or.   q(i,j,k,iBx)>1.e8) then
         	     print*,'bx or bz problems, problem.data';stop
         	  end if


!                       !!!A, based in http://iopscience.iop.org/1538-4357/469/2/L127/pdf/5447.pdf
!!!                    r =SQRT( x**2 + y**2 + z**2 )
!!!						if (r.le.Rjet) cycle
!!!                   dr =SQRT( dx**2 + dy**2 + dz**2 )
!!!								bmax= 		 ((3e12/lScale)/dr)**2d0*&
!!!								             ((dr/(3e12/lScale))-1d0)
!!!                        q(i,j,k,iBz) = ((3e12/lScale)/r)**2d0*&
!!!								             ((r/(3e12/lScale))-1d0)*&
!!!								   sigmaFACTOR/bmax

!!!							if (dabs(x)==dx) q(i,j,k,iBz)=0d0

           q(i,j,k,iE) = q(i,j,k,iE)+ half*SUM(q(i,j,k,iBx:iBz)**2)
													
                 END DO
              END DO
           END DO
     END if!mhd=t

!!==================================================================
!!
!! The above code was from Martine and used before Jan 6th 2015
!!
!!===================================================================
     CASE(2)
!!=====================================================================
!!
!! The following code was from Baowei, updated from Martin's version 
!! and according to Bruce's request
!!  Jan 6th 2015 
!!
!!======================================================================

      DO i=1-rmbc, mx+rmbc         ; x = (xl+(REAL(i,xPrec)-half)*dx)
         DO j=1-rmbc, my+rmbc      ; y = (yl+(REAL(j,xPrec)-half)*dy)
            DO k=1-zrmbc, mz+zrmbc ; z = (zl+(REAL(k,xPrec)-half)*dz)
               r =SQRT( x**2 + y**2 + z**2 ) ! from origin
               r2=SQRT( (y-2d0*Rjet)**2 + x**2 + z**2 ) ! from clump center
               dr =SQRT( dx**2 + dy**2 + dz**2 )

               !! if the Abient centering at the AGB star
               if (lAGBcenter) then
                   r_sc=SQRT((x-AGBstarPos(1))**2 + (y-AGBstarPos(2))**2 + (z-AGBstarPos(3))**2)
               end if
 
               if (outflowType==clump) r=r2 !b4 27jan'12

               !A M B I E N T 
               q(i,j,k,1) = namb/nScale !density
               q(i,j,k,itracer3)=0d0 ; q(i,j,k,itracer4)=0d0
               q(i,j,k,itracer5)=0.0

               if(r .le. Rjet) then
                  density_flag(i,j,k)=.false.
               end if


               if (stratified) then !27jan'12
                  if (outflowType /= clump) then
                    !!========= Martin's version ==========
                     !q(i,j,k,1) = q(i,j,k,1)/&
                     !     (r+1d0+.25d0*dr)**2
                    !!=======================================
                    if(lAGBcenter) then
                       r=r_sc
                       if(r.gt.rStratified) then
                         q(i,j,k,1) = (namb/nScale)*((rStratified/r)**2)
                       else
                         q(i,j,k,1) = (namb/nScale)
                       end if
                    else  
                      if(r.gt.Rjet) q(i,j,k,1) = (namb/nScale)*((Rjet/r)**2)
                    end if
                  else
                     !========= Martin's version =========
                     if (r.gt.Rjet) then
                         q(i,j,k,1) = q(i,j,k,1)/&
                          (r-Rjet+1d0+.25d0*dr)**2
                     else
                         density_flag(i,j,k) = .false.
                         q(i,j,k,itracer3)=q(i,j,k,1)
                     end if
                     !====================================

                  end if
               end if


               !
!!! TORUS
               !IF (r > rjet) THEN
               !From Frnak & Mellema, 1994ApJ...430..800F:
               !produced a pole-to-equator density contrast
               torus2=1d0
               if ( torus ) then
                  if ( r.eq.0d0   .or.   (y-2d0*Rjet).eq.0d0) then
                     print*,'torus error';stop
                  else
                     !cos2theta=x**2/( x**2 + y**2 + z**2 )-z**2/x**2 !-------------A DIDNT WORK
                     IF (outflowType == clump) THEN
                        cos2theta=2d0*acos((y-2d0*Rjet)/r)!----------------------------B
                     ELSE 
                        cos2theta=2d0*acos((y-0d0*Rjet)/r)
                     END IF
                     !cos2theta=(x-2d0*Rjet)**2/r**2-z**2/(x-2d0*Rjet)**2!-----C DIDNT WORK
                  end if
                  ex1=-2d0*torusbeta
                  !ex2=torusbeta*cos2theta-torusbeta !---------------------------A DIDNT WORK
                  ex2=torusbeta*cos(cos2theta)-torusbeta!-----------------------B
                  !ex2=torusbeta*cos2theta-torusbeta !----------------------C DIDNT WORK
                  Ftheta=1d0-(  torusalpha*(dexp(ex2)-1d0)/(dexp(ex1)-1d0)  )
                  !
!!!
                  !1)
                  !torus2=Rjet**2/r**2/Ftheta 
                  !2)
                  torus2=1d0/Ftheta
!!!
                !!========= Baowei's version ========
                 if (r.ge.Rjet) then
                     q(i,j,k,1) = (namb/nScale)*((Rjet/r)**2)
                    ! print*, "torus(i,j,k)=", i,j,k, torus2
                 else
                     density_flag(i,j,k)=.false.
                     if(outflowType == clump) then
                        q(i,j,k,itracer3) = q(i,j,k,1)
                     else
                        q(i,j,k,itracer4) = q(i,j,k,1) 
                     end if
                 end if
                !===================================
                  !if (stratified) torus2 = 1.0 + torus2 ! if both stratified and torus are true add the two density together according to Bruce; changed back to Martin's way on Jan 14th 2015
                  !
               q(i,j,k,1)=q(i,j,k,1)*torus2

               end if !if torus

!!!!!!!!!!RINGS!!!!!!!!!!!!!!
               f=0.5d0

               dec7=2.5d0/1.5d0  
               do ii=1,1000
                  f=f+exp(-(6d0*(r-dec7))**2)
                  dec7=dec7+1.5d0/1.5d0
               end do

               !dec7=3d0/1.5d0  
               !do ii=1,1000
               !   f=f+exp(-(10d0*(r-dec7))**2)
               !	dec7=dec7+.8d0/1.5d0
               !end do
               f=2d0*f
               if (rings) q(i,j,k,1) = q(i,j,k,1)*f
!!!!!!!!!!RINGS!!!!!!!!!!!!!!
               !END IF

               q(i,j,k,imom(1:nDim)) = 0d0
               q(i,j,k,iE) = q(i,j,k,1)*tamb/TempScale/gamma1
!            q(i,j,k,itracer2)=GetCoolingStrength(q(i,j,k,:),lform)*pScale/TimeScale !erg/cc/s

!!!!!!!!!!
!!!!!!!!!!
!!!!!!!!!!
! C L U M P:
               if (outflowType == clump .AND. r2.lt.Rjet) then !Martin version
          !     if (outflowType == clump .AND. r2.lt.Rjet) then !Baowei 1.2.2015
          !			q(i,j,k,itracer3)=1d0 !martin's original set up for clump tracer
                  q(i,j,k,1) = namb/nScale+&
                       njet/nscale*(1d0-(r2/Rjet)**2) ![cu]
                  q(i,j,k,itracer3)=q(i,j,k,1) !set the clump tracer as density

                  q(i,j,k,3) = vjet*q(i,j,k,1)/velscale ![cu]

                  density_flag(i,j,k)=.false.

                  q(i,j,k,iE) = q(i,j,k,1)*tjet/TempScale/gamma1 +& !internal enegy     [cu]
                             half*q(i,j,k,3)**2/q(i,j,k,1)
                  !if (ndim.eq.2   .and.   iCylindrical.ne.2) then 
                  !   q(i,j,k,iE) =q(i,j,k,iE) + &
                  !		   half*SUM(q(i,j,k,2:3)**2)/q(i,j,k,1)            !kinetic energy
                  !else
                  !   q(i,j,k,iE) =q(i,j,k,iE) + &
                  !        half*SUM(q(i,j,k,2:4)**2)/q(i,j,k,1)            !kinetic energy
                  !end if

               end if!r<Rjet
              

          !! add tracers, Baowei Liu 04/16/15 
         ! if(outflowType == collimated) q(i,j,k,itracer4)=q(i,j,k,1) ! Jet tracer
         ! if(outflowType == clump) q(i,j,k,itracer3)=q(i,j,k,1) ! clump tracer
         !  if(i .eq. 0 .AND. j.eq.0 .AND. k.eq.0) print*, "jettracer=",q(i,j,k,itracer4)  
           if(density_flag(i,j,k)) q(i,j,k,itracer5)=1.0
 
            END DO
         END DO
      END DO

      aux(:,:,:,:)=0d0 
      q(:,:,:,iBx:iBz) = 0d0


      if (outflowType==clump .and. lMHD) then
      if(nDim==3 .and. sigma.ne.0d0 .and. MPI_ID == 0) THEN
       print*,''
       print*,'3D magnetic clump not implemented yet';stop
       END IF
      
      if (once) then ; once=.false. ; dxFIXED=dx/(2d0**3d0)
      end if
           DO i=1-rmbc, mx+rmbc  ; x = (xl+(REAL(i,xPrec)-half)*dx) 
           DO j=1-rmbc, my+rmbc  ; y = (yl+(REAL(j,xPrec)-half)*dy)
           DO k=1-zrmbc,mz+zrmbc ; z = (zl+(REAL(k,xPrec)-half)*dz)
              r=SQRT( (y-2d0*Rjet)**2 + x**2 + z**2 ) ! from clump center
             if (r>=Rjet) cycle
              bmax=        ((3e12/lScale)/dxFIXED   )**2d0*&
                           ((dxFIXED     /(3e12/lScale))-1d0)

              if (nDim==2) then
                 q(i,j,k,iBz) = ((3e12/lScale)/sqrt(x**2d0))**2d0*&
                   ((sqrt(x**2d0)/(3e12/lScale))-1d0)*& !B
                   sigmaFACTOR/bmax
              else if (nDim==3) then
!NOT IMPLEMENTED IN TERMA OF AUX YET
               q(i,j,k,iBx) = ((3e12/lScale)/sqrt(x**2d0))**2d0*&
                  ((sqrt(x**2d0)/(3e12/lScale))-1d0)*& !B
                   sigmaFACTOR/bmax
          
               q(i,j,k,iBz) = ((3e12/lScale)/sqrt(z**2d0))**2d0*&
                  ((sqrt(z**2d0)/(3e12/lScale))-1d0)*& !B
                 sigmaFACTOR/bmax
          
            end if
            if (q(i,j,k,iBz)>1.e8   .or.   q(i,j,k,iBx)>1.e8) then
               print*,'bx or bz problems, problem.data';stop
            end if


           q(i,j,k,iE) = q(i,j,k,iE)+ half*SUM(q(i,j,k,iBx:iBz)**2)

          !! add tracers, Baowei Liu 04/16/15 
   !       if(outflowType == collimated) q(i,j,k,itracer4)=q(i,j,k,1) ! Jet tracer
   !       if(outflowType == clump) q(i,j,k,itracer3)=q(i,j,k,1) ! clump tracer
                 END DO
              END DO
           END DO
     END if!mhd=t
     

!!=====================================================================
!!
!! The above code was from Baowei, updated from Martin's version 
!! and according to Bruce's request
!!  Jan 6th 2015 
!!
!!======================================================================

    END SELECT

    delta_mass=0.0
    delta_mass_3D=0.0
    CALL ProblemBeforeStep(Info)
    DO i=1,mx
       r=info%xbounds(1,1)+(real(i)-half)*dx
       IF (level < MaxLevel) THEN
          local_mass_3D=local_mass_3D+sum(Info%q(i,1:my,1:mz,1), density_flag(i,1:my,1:mz) .AND. Info%childmask(i,1:my,1:mz) == 0) * 2*pi*r*dx**2
       ELSE
          local_mass_3D=local_mass_3D+sum(Info%q(i,1:my,1:mz,1), density_flag(i,1:my,1:mz)) * 2*pi*r*dx**2
       END IF
    END DO

   
    DEALLOCATE(density)

  END SUBROUTINE ProblemGridInit


  !> Applies Boundary conditions
  !! @param Info Info object
   SUBROUTINE ProblemBeforeStep(Info)
      TYPE(InfoDef) :: Info
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:)     :: q, aux
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: density
      INTEGER :: i,j,k,ii,jj,kk,mx,my,mz,rmbc,zrmbc,level,iErr,mm
      REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,time, r,dt, dr,velFactor,torus2,cos2theta,ex1,ex2,Ftheta, drho
      REAL(KIND=qPREC) :: qnew(NrHydroVars), qwind(NrHydroVars), qjet(NrHydroVars), v(3), rhoFact, fact, fact2, f, timef, bmax, theta
      REAL(KIND=qPREC) :: mom_flow=1d0
      REAL(KIND=qPREC) :: MODULATION
      LOGICAL, SAVE :: once=.true.
      REAL(KIND=xprec), SAVE :: dxFIXED=.5d0

      level=Info%level
      dt=levels(level)%dt
      time=levels(level)%tnow
      ! if (time.eq.dt) return 
      !we don't want to overwrite the initial conditions
      !and for the clump case we don't want to use this routine

      if (outflowType==clump) return

      !INITIAL ACCELERATION OF JET/WIND
      IF (time < tt) THEN
         timef = sin(Pi/2d0*time/tt)
      ELSE
         timef=1d0
      END IF

      !time-variable momentum flow requested by bbalick.
      mom_flow=1d0 
      if (balick_a.ne.0d0 .or. balick_b.ne.0d0) &
        mom_flow=(balick_a*time**balick_aa + balick_b*time**balick_bb)*&
              exp(-balick_c*time)

      q=>Info%q
      aux=>Info%aux

      !Calculating the number of ghost cells on each side of the grid.
      rmbc=levels(level)%gmbc(levels(level)%step)
      mx = Info%mX(1)
      my = Info%mX(2)
      mz = Info%mX(3)
      dx=levels(level)%dX
      dy=dx
      dz=dx
      if (once) then
         once=.false.
         dxFIXED=dx/(2d0**3d0)
      end if
      xl=Info%xBounds(1,1)
      yl=Info%xBounds(2,1)
      zl=Info%xBounds(3,1)
      SELECT CASE(nDim)
      CASE(2)
         zrmbc=0;mz=1;dz=0d0;zl=0d0
      CASE(3)
         zrmbc=rmbc
      END SELECT

      ALLOCATE(density(1-rmbc:mx+rmbc, 1-rmbc:my+rmbc, 1-zrmbc:mz+zrmbc))
      density=Info%q(1-rmbc:mx+rmbc, 1-rmbc:my+rmbc, 1-zrmbc:mz+zrmbc, 1)

      DO i=1-rmbc, mx+rmbc  ; x = (xl+(REAL(i,xPrec)-half)*dx)
      DO j=1-rmbc, my+rmbc  ; y = (yl+(REAL(j,xPrec)-half)*dy)
      DO k=1-zrmbc,mz+zrmbc ; z = (zl+(REAL(k,xPrec)-half)*dz)

         r =DSQRT( x**2 + y**2 + z**2 )
         dr =DSQRT( dx**2 + dy**2+ dz**2 )

         IF (y>0d0 .and. r.le.Rjet) THEN !20may2014 Martin version, before 1.2.2015
       !  !IF (y>0d0 .and. r.lt.Rjet) THEN !Baowei 1.2.2015
       !  !IF (x>0d0 .and. y>0d0 .and. r.le.Rjet) THEN !b 4 28feb2014

           CALL ConvertTotalToInternalEnergy(q(i,j,k,:))               
           qwind=0d0
           qjet=0d0

           fact=1d0
           if (cos(open_angle) .le. abs(y)/r) fact=0d0
           if (outflowType == collimated) fact=0d0
           qwind(1)=fact*namb/nScale/&
                   (r-Rjet+1d0+.25d0*dr)**2
           !!! TORUS
           torus2=1d0
           if ( torus ) then
              if ( r.eq.0d0   .or.   (y-2d0*Rjet).eq.0d0) then
                 print*,'torus error';stop
              else
                 cos2theta=2d0*acos((y-0d0*Rjet)/r)
              end if
              ex1=-2d0*torusbeta
              ex2=torusbeta*cos(cos2theta)-torusbeta
              Ftheta=1d0-(  torusalpha*(dexp(ex2)-1d0)/(dexp(ex1)-1d0)  )
              torus2=1d0/Ftheta
           end if
           qwind(1)=qwind(1)*torus2
           v(:)=0d0
           qwind(imom(1:nDim))=v(1:nDim)*qwind(1)
           qwind(iE)=qwind(1)*wind_temp/TempScale*gamma7
           !where gamma7=1/gamma

           IF (lMHD   .and.    .false.) THEN
              qwind(iBx)=0d0; qwind(iBy)=0d0; aux(i,j,k,1:2)=0d0
              bmax=        ((3e12/lScale)/dxFIXED   )**2d0*&
                 ((dxFIXED     /(3e12/lScale))-1d0)
              
              qwind(iBz) = ((3e12/lScale)/sqrt(x**2+z**2d0))**2d0*&
                         ((sqrt(x**2d0+z**2d0)/(3e12/lScale))-1d0)*& !B
                         sigmaFACTOR/bmax

              if (qwind(iBz)>1e8) then
                print*,'bz problems, problem.data, line 347';stop
              end if
           ELSE
              qwind(iBx)=0d0; qwind(iBy)=0d0; qwind(iBz)=0d0; aux(i,j,k,1:2)=0d0
           END IF!mhd=t

           !! Jet tracer, by Baowei 4.20.15
           qwind(itracer4)=qwind(1)
           !======== J E T=========:
           IF (outflowType == collimated) then
              !q(i,j,k,itracer4)=1d0  !old tracer

              !q(i,j,k,itracer4)=q(i,j,k,1) ! new tracer by Baowei
              fact=1d0 !10 mar 2014
              !fact=exp(-(x**2+z**2)/jet_width**2) ! b 4 10 mar 2014
              qjet(1)=njet/nScale*fact
              v=vjet/velscale*(/0d0,1d0,0d0/)/fact*timef !7 Nov 2014 get rid the Rjet for vjet
              ! v=vjet/velscale*(/0d0,1d0,0d0/)/Rjet*fact*timef !10 mar 2014
              !v=vjet/velscale*(/0d0,y,0d0/)/Rjet*fact*timef !b 4 10 mar 2014
              qjet(imom(1:nDim))=v(1:nDim)*qjet(1)*& !ramp up velocity
                               mom_flow !5 may 2014, time dependent mom flux requested by bbalick. 
              qjet(iE)=qjet(1)*tjet/TempScale*gamma7

              !========= W I N D ===========:							
           ELSEIF (outflowType == conical) THEN
              !q(i,j,k,itracer4)=1d0  ! old tracer
              !q(i,j,k,itracer4)=q(i,j,k,1) !new tracer by Baowei
              if (y.ne.0d0) then 
                  theta=datan( sqrt(x**2+z**2)/y)
              else ; theta=datan( sqrt(x**2+z**2)/dxFIXED)
              end if

              fact =1d0!20 feb 2014
              fact2=1d0!4mra2014
              !F)4mar2014
              if (cos(open_angle) .gt. abs(y)/r) fact2=0d0
            
              !E) NO: very close to spherical, 28 feb 2014
              !if (open_angle .gt. abs(acos(y/r))) &
              !fact=exp( -abs(acos(y/r))+open_angle )
            
              !D) MAYBE: good one, 27feb 2014
              if (cos(open_angle) .gt. abs(y)/r) &
                fact=1.00d0*&!28feb2014
                !0.50d0*&!28feb2014
                 !0.25d0*&!28feb2014
                 (Pi*half-abs(acos(y/r)))/&
                     (Pi*half-open_angle)
                             
               !C)
               !if (cos(open_angle) .gt. abs(y)/r) fact=.01d0
               !B)
               !if (cos(open_angle) .gt. abs(y)/r) & 	!24feb2014
               !fact=0.5d0*abs(y)/r!24feb2014
               !A)
               !if (  atan(sqrt(x**2+z**2)/abs(y)).gt.open_angle  ) & 	!20feb2014
                                          !fact=0.2d0*
             
              !fact=exp(-(atan(sqrt(x**2+z**2)/y**2)/open_angle)**2) !b4 feb20.2014
              qjet(1)=njet/nscale*fact2
              v=vjet/velscale*(/x,y,z/)/Rjet*fact*timef

              !tempering, a la Lee & Sahai, 2003
              if (tf.ne.0d0) then
                 ! modified by Baowei, Nov 20 2014: add lObliqueCone and 45 deg part
                 if(lObliqueCone) then
                 !============= Bruce's suggesttion with 45 deg  ==============
                  ! print*,"45 deg theta=", theta
                   qjet(1)=qjet(1)*max(0.001d0,dexp(-((theta-0.25*PI)/tf)**2d0))
                   v      =v      *max(0.001d0,dexp(-((theta-0.25*PI)/tf)**2d0))
                 !========================================================
                 else
                  ! print*,"not 45 deg theta=", theta
                 !============= Martin's version ==============
                   qjet(1)=qjet(1)*max(0.001d0,dexp(-(theta/tf)**2d0))
                   v      =v      *max(0.001d0,dexp(-(theta/tf)**2d0))
                 !========================================
                 end if

                 if(lPulsed) then
                   !MODULATION=abs(cos(Omega_Pulsed*time))**PulseModule
                   MODULATION=abs(cos(Omega_Pulsed*time))

                   if(PulseModule > 4) PulseModule=4
                   MODULATION=MODULATION**PulseModule
                   qjet(1)=qjet(1)*MODULATION
                   v = v*MODULATION 
                 end if
              end if
              qjet(imom(1:nDim))=v(1:nDim)*qjet(1)*& !ramp up velocity
                       mom_flow !5 may 2014, time dependent mom flux requested by bbalick. 
              qjet(iE)=qjet(1)*tjet/TempScale*gamma7
              !where gamma7=1/gamma
           ELSE
              qjet(1:iE)=0d0
           END IF

           IF (lMHD) THEN
              qjet(iBx)=0d0; qjet(iBy)=0d0; aux(i,j,k,1:2)=0d0
              if (nDim==3) aux(i,j,k,3)=0d0
              bmax=    ((3e12/lScale)/dxFIXED   )**2d0*&
                 ((dxFIXED/(3e12/lScale))-1d0)

              if (nDim==2) then
                 qjet(iBz) = ((3e12/lScale)/sqrt(x**2d0))**2d0*&
                    ((sqrt(x**2d0)/(3e12/lScale))-1d0)*& !B
                     sigmaFACTOR/bmax
                 !if (sqrt(x**2d0)<.8d0*Rjet) & 
                 !  qjet(iBz) =qjet(iBz)*&
                 !    sin(sqrt(x**2d0)/(.8d0*Rjet)*Pi*half)!21ap14
                 if (sqrt(x**2d0)<4d0*dxFIXED) then
                    qjet(iBz) =qjet(iBz)*.25d0!18mar14
                 end if
              else if (nDim==3) then
                 qjet(iBx) = ((3e12/lScale)/sqrt(x**2d0))**2d0*&
                           ((sqrt(x**2d0)/(3e12/lScale))-1d0)*& !B
                           sigmaFACTOR/bmax
                 qjet(iBz) = ((3e12/lScale)/sqrt(z**2d0))**2d0*&
                          ((sqrt(z**2d0)/(3e12/lScale))-1d0)*& !B
                          sigmaFACTOR/bmax
              end if
              if (qjet(iBz)>1.e8   .or.   qjet(iBx)>1.e8) then
                 print*,'bx or bz problems, problem.data';stop
              end if
           ELSE
              qjet(iBx)=0d0; qjet(iBy)=0d0; qjet(iBz)=0d0; aux(i,j,k,1:2)=0d0
              if (nDim==3) aux(i,j,k,3)=0d0
           END IF

           !! jet tracer, Baowei, 4.20.15
           qjet(itracer4)=qjet(1)
           qnew=qwind+qjet
           q(i,j,k,1:NrHydroVars)=q(i,j,k,1:NrHydroVars)+(qnew(1:NrHydroVars)-q(i,j,k,1:NrHydroVars))
           CALL ConvertInternaltoTotalEnergy(q(i,j,k,:))
           
           !END IF!r<rjet
          END IF !y>0 or r>rjet
      END DO
      END DO
      END DO

      !!write(*,*) 'a', info%level, MPI_ID, delta_mass, delta_mass_3D
      DO i=1,mx
         r=info%xbounds(1,1)+(real(i)-half)*dx
         IF (level < MaxLevel) THEN
            delta_mass=delta_mass+sum(Info%q(i,1:my,1:mz,1)-density(i,1:my,1:mz), Info%childmask(i,1:my,1:mz)== 0) * dx**2
            delta_mass_3D=delta_mass_3D+sum(Info%q(i,1:my,1:mz,1)-density(i,1:my,1:mz), Info%childmask(i,1:my,1:mz) == 0) * 2*pi*r*dx**2
         ELSE
            delta_mass=delta_mass+sum(Info%q(i,1:my,1:mz,1)-density(i,1:my,1:mz)) * dx**2
            delta_mass_3D=delta_mass_3D+sum(Info%q(i,1:my,1:mz,1)-density(i,1:my,1:mz)) * 2*pi*r*dx**2
         END IF
      END DO
      !!write(*,*) 'b', info%level, MPI_ID, delta_mass, delta_mass_3D
   END SUBROUTINE ProblemBeforeStep

  !> Could be used to update grids pre-output
  !! @param Info Info Object
  SUBROUTINE ProblemAfterStep(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemAfterStep

  !> Could be used to set force refinement
  !! @param Info Info object
  SUBROUTINE ProblemSetErrFlag(Info)
     TYPE(InfoDef) :: Info
     REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:)     :: q
     INTEGER :: i,j,k,mx,my,mz,rmbc,zrmbc,level
     REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz
     
     return

     if (.not.outflowType   .or.   outflowType.eq.3) return
     q=>Info%q
     level=Info%level
     rmbc=levels(level)%gmbc(levels(level)%step)
     mx = Info%mX(1)
     my = Info%mX(2)
     mz = Info%mX(3)
     dx=levels(level)%dX
     dy=dx;dz=dx
     xl=Info%xBounds(1,1)
     yl=Info%xBounds(2,1)
     zl=Info%xBounds(3,1) 
     IF (nDim.eq.2) then 
        mz=1;dz=0;zl=0
     END IF

     DO j=1, my ; y = (yl+(REAL(j,xPrec)-half)*dy)
        DO k=1,mz ; z = (zl+(REAL(k,xPrec)-half)*dz) 
           DO i=1, mx ; x = (xl+(REAL(i,xPrec)-half)*dx)
              if (sqrt(x**2+y**2+z**2).lt.Rjet*1.15d0) Info%ErrFlag(i,j,k)=1
           END DO
        END DO
     END DO
  END SUBROUTINE ProblemSetErrFlag

  SUBROUTINE ProblemBeforeGlobalStep(n)
     INTEGER :: n 
     INTEGER :: ierr
     INTEGER :: massPrintTime=0 !! only need to print the initial mass once.
     if (outflowType==clump) return

     if (n == 0 .AND. levels(n)%tnow==start_time) THEN
        delta_mass=0.0
        delta_mass_3D=0.0
     END if

     IF (n == 0) THEN
        IF (CURRENT_frame == oldframe) RETURN
        oldframe=current_frame
        CALL MPI_ALLREDUCE(local_mass, total_mass, 1, MPI_DOUBLE_PRECISION, MPI_SUM, levels(n)%MPI_COMM, ierr)
        CALL MPI_ALLREDUCE(local_mass_3D, total_mass_3D, 1, MPI_DOUBLE_PRECISION, MPI_SUM, levels(n)%MPI_COMM, ierr)
        CALL MPI_ALLREDUCE(delta_mass, total_delta_mass,  1, MPI_DOUBLE_PRECISION, MPI_SUM, levels(n)%MPI_COMM, ierr)
        CALL MPI_ALLREDUCE(delta_mass_3D, total_delta_mass_3D,  1, MPI_DOUBLE_PRECISION, MPI_SUM, levels(n)%MPI_COMM, ierr)
    
        !! compute the mass flow rate (in solar mass per year)  
        mass_flux=2.0*total_delta_mass/levels(n)%tnow*(mscale/timeScale)/(mSolar/yr)
        mass_flux_3D=2.0*total_delta_mass_3D/levels(n)%tnow*(mscale/timeScale)/(mSolar/yr)

        total_mass_3D=total_mass_3D + total_delta_mass_3D
         
        if(MPI_ID == 0 .AND. levels(n)%tnow .ge. start_time) then
         ! if(levels(n)%tnow == start_time) then !print the initial grid mass only
         !   if( GxBounds(2,1) == 0) then !if only the top half calculated, double the values
         !     print *, "the initial 2D grid mass (in solar mass) is: ", 2.0*total_mass/(mSolar/mscale)
         !     !print *, "the 2.5D mass flow rate (in solar masses per year) is: ", 2.0*mass_flux 
         !     print *, "the initial 3D grid mass (in solar mass) is: ", 2.0*total_delta_mass_3D/(mSolar/mscale)
         !     !print *, "the equivalent 3D mass flow rate (in solar masses per year) is: ", 2.0*mass_flux_3D
         !   else
         !     print *, "the initial 2D grid mass (in solar mass) is: ", total_mass/(mSolar/mscale)
         !    ! print *, "the 2.5D mass flow rate (in solar masses per year) is: ", mass_flux 
         !     print *, "the initial 3D grid mass (in solar mass) is: ", total_delta_mass_3D/(mSolar/mscale)
         !    ! print *, "the equivalent 3D mass flow rate (in solar masses per year) is: ", mass_flux_3D
         !   end if

         ! else !print the flow rate only

         if( GxBounds(2,1) == 0) then !if only the top half calculated, double the values
            if (massPrintTime ==0) then
               print *
               write(745,*)
               !print *, "dx=", levels(n)%dx
               !print *, "total_mass=", total_mass
               !print *, "mSolar=", mSolar
               !print *, "mscale=", mscale
               !print *, "mSolar CU=", mSolar/mScale
               
               !print *, "the 2.5D grid mass (in solar mass) is: ", 2.0*total_mass/(mSolar/mscale)
               !print *, "the equivalent 3D grid mass (in solar mass) is: ", 2.0*total_mass_3D/(mSolar/mscale)
               write(745,*) "Time(y)   Total mass on the grid (solar masses)  3D mass flow rate(solar mass per year)"
               write(745,*) levels(n)%tnow*8250.0,"  ", 2.0*total_mass_3D/(mSolar/mscale), "   ", 0.0
               massPrintTime = 1
             else if (massPrintTime == 1 .AND. levels(n)%tnow .gt. start_time) then
                write(745,*)
                write(745,*) levels(n)%tnow*8250.0,"  ", 2.0*total_mass_3D/(mSolar/mscale), 2.0*mass_flux_3D
             end if

             !print *, "the 2D mass flow rate (in solar masses per year) is: ", 2.0*mass_flux 
             !print *, "current time in cu: ",levels(n)%tnow, "the 3D mass flow rate (in solar masses per year) is: ", 2.0*mass_flux_3D
             !write(745,*) "current time in cu: ",levels(n)%tnow, "the 3D mass flow rate (in solar masses per year) is: ", 2.0*mass_flux_3D
             !print *
           else
             !print *
             !if (massPrintTime ==0) then
             !print *, "the 2.5D grid mass (in solar mass) is: ", total_mass/(mSolar/mscale)
             !print *, "the 2D mass flow rate (in solar masses per year) is: ", mass_flux 
             !  massPrintTime = 1
             !end if

             !print *, "the equivalent 3D grid mass (in solar mass) is: ", total_mass_3D/(mSolar/mscale)
             !write(745,*) "the equivalent 3D grid mass (in solar mass) is: ", total_mass_3D/(mSolar/mscale)
             !print *, "the 3D mass flow rate (in solar masses per year) is: ", mass_flux_3D
             !write(745,*)
             write(745,*)
             write(745,*) levels(n)%tnow*8250.0, "  ",total_mass_3D/(mSolar/mscale), "  ", mass_flux_3D
             !write(745,*) "the 3D mass flow rate (in solar masses per year) is: ", mass_flux_3D
             !print *
           end if
           !end if !!start_time or not
       end if !!MPI_ID
     END IF
  END SUBROUTINE ProblemBeforeGlobalStep

END MODULE Problem

