/home/hemv_g/srf_a/RAMSES_CODE_CORI/ramses06062019_grackle311Cosmo/patch/mom2/courant_fine.f90 subroutine courant_fine(ilevel) use amr_commons use hydro_commons use poisson_commons use mpi_mod implicit none #ifndef WITHOUTMPI integer::info real(kind=8),dimension(3)::comm_buffin,comm_buffout #endif integer::ilevel !---------------------------------------------------------------------- ! Using the Courant-Friedrich-Levy stability condition, ! ! this routine computes the maximum allowed time-step. ! !---------------------------------------------------------------------- integer::i,ivar,idim,ind,ncache,igrid,iskip integer::nleaf,ngrid,nx_loc integer,dimension(1:nvector),save::ind_grid,ind_cell,ind_leaf real(dp)::dt_lev,dx,vol,scale real(kind=8)::mass_loc,ekin_loc,eint_loc,dt_loc real(kind=8)::mass_all,ekin_all,eint_all,dt_all real(dp),dimension(1:nvector,1:nvar),save::uu real(dp),dimension(1:nvector,1:ndim),save::gg real(dp),dimension(1:nvector)::pp if(numbtot(1,ilevel)==0)return if(verbose)write(*,111)ilevel mass_all=0.0d0; mass_loc=0.0d0 ekin_all=0.0d0; ekin_loc=0.0d0 eint_all=0.0d0; eint_loc=0.0d0 dt_all=dtnew(ilevel); dt_loc=dt_all ! Mesh spacing at that level nx_loc=icoarse_max-icoarse_min+1 scale=boxlen/dble(nx_loc) dx=0.5D0**ilevel*scale vol=dx**ndim ! Loop over active grids by vector sweeps ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Loop over cells do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=ind_grid(i)+iskip end do ! Gather leaf cells nleaf=0 do i=1,ngrid if(son(ind_cell(i))==0)then nleaf=nleaf+1 ind_leaf(nleaf)=ind_cell(i) end if end do ! Gather hydro variables do ivar=1,nvar do i=1,nleaf uu(i,ivar)=uold(ind_leaf(i),ivar) end do end do ! Gather gravitational acceleration gg=0.0d0 if(poisson)then do idim=1,ndim do i=1,nleaf gg(i,idim)=f(ind_leaf(i),idim) end do end do end if ! Compute total mass do i=1,nleaf mass_loc=mass_loc+uu(i,1)*vol end do ! Compute total energy do i=1,nleaf ekin_loc=ekin_loc+uu(i,ndim+2)*vol end do ! Compute total internal energy do i=1,nleaf eint_loc=eint_loc+uu(i,ndim+2)*vol end do do ivar=1,ndim do i=1,nleaf eint_loc=eint_loc-0.5d0*uu(i,1+ivar)**2/max(uu(i,1),smallr)*vol end do end do #if NENER>0 do ivar=1,nener do i=1,nleaf eint_loc=eint_loc-uu(i,ndim+2+ivar)*vol end do end do #endif if(momentum_feedback>0)then do i=1,nleaf pp(i)=pstarold(ind_leaf(i)) end do endif ! Compute CFL time-step if(nleaf>0)then call cmpdt(uu,gg,pp,dx,dt_lev,nleaf) dt_loc=min(dt_loc,dt_lev) end if end do ! End loop over cells end do ! End loop over grids ! Compute global quantities #ifndef WITHOUTMPI comm_buffin(1)=mass_loc comm_buffin(2)=ekin_loc comm_buffin(3)=eint_loc call MPI_ALLREDUCE(comm_buffin,comm_buffout,3,MPI_DOUBLE_PRECISION,MPI_SUM,& &MPI_COMM_WORLD,info) call MPI_ALLREDUCE(dt_loc,dt_all,1,MPI_DOUBLE_PRECISION,MPI_MIN,& &MPI_COMM_WORLD,info) mass_all=comm_buffout(1) ekin_all=comm_buffout(2) eint_all=comm_buffout(3) #endif #ifdef WITHOUTMPI mass_all=mass_loc ekin_all=ekin_loc eint_all=eint_loc dt_all=dt_loc #endif mass_tot=mass_tot+mass_all ekin_tot=ekin_tot+ekin_all eint_tot=eint_tot+eint_all dtnew(ilevel)=MIN(dtnew(ilevel),dt_all) 111 format(' Entering courant_fine for level ',I2) end subroutine courant_fine subroutine check_cons(ilevel) use amr_commons use hydro_commons use poisson_commons use mpi_mod implicit none #ifndef WITHOUTMPI integer::info real(kind=8),dimension(3)::comm_buffin,comm_buffout #endif integer::ilevel !---------------------------------------------------------------------- ! Check mass and energy conservation !---------------------------------------------------------------------- integer::i,ivar,ind,ncache,igrid,iskip integer::nleaf,ngrid integer,dimension(1:nvector),save::ind_grid,ind_cell,ind_leaf real(dp)::dx,vol real(kind=8)::mass_loc,ekin_loc,eint_loc real(kind=8)::mass_all,ekin_all,eint_all real(dp),dimension(1:nvector,1:nvar),save::uu if(numbtot(1,ilevel)==0)return if(verbose)write(*,111)ilevel mass_all=0.0d0; mass_loc=0.0d0 ekin_all=0.0d0; ekin_loc=0.0d0 eint_all=0.0d0; eint_loc=0.0d0 ! Mesh spacing at that level dx=0.5D0**ilevel*boxlen vol=dx**ndim ! Loop over active grids by vector sweeps ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Loop over cells do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=ind_grid(i)+iskip end do ! Gather leaf cells nleaf=0 do i=1,ngrid if(son(ind_cell(i))==0)then nleaf=nleaf+1 ind_leaf(nleaf)=ind_cell(i) end if end do ! Gather hydro variables do ivar=1,nvar do i=1,nleaf uu(i,ivar)=uold(ind_leaf(i),ivar) end do end do ! Compute total mass do i=1,nleaf mass_loc=mass_loc+uu(i,1)*vol end do ! Compute total energy do i=1,nleaf ekin_loc=ekin_loc+uu(i,ndim+2)*vol end do ! Compute total internal energy do i=1,nleaf eint_loc=eint_loc+uu(i,ndim+2)*vol end do do ivar=1,ndim do i=1,nleaf eint_loc=eint_loc-0.5d0*uu(i,1+ivar)**2/max(uu(i,1),smallr)*vol end do end do #if NENER>0 do ivar=1,nener do i=1,nleaf eint_loc=eint_loc-uu(i,ndim+2+ivar)*vol end do end do #endif end do ! End loop over cells end do ! End loop over grids ! Compute global quantities #ifndef WITHOUTMPI comm_buffin(1)=mass_loc comm_buffin(2)=ekin_loc comm_buffin(3)=eint_loc call MPI_ALLREDUCE(comm_buffin,comm_buffout,3,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) mass_all=comm_buffout(1) ekin_all=comm_buffout(2) eint_all=comm_buffout(3) #endif #ifdef WITHOUTMPI mass_all=mass_loc ekin_all=ekin_loc eint_all=eint_loc #endif mass_tot=mass_tot+mass_all ekin_tot=ekin_tot+ekin_all eint_tot=eint_tot+eint_all 111 format(' Entering check_cons for level ',I2) end subroutine check_cons /home/hemv_g/srf_a/RAMSES_CODE_CORI/ramses06062019_grackle311Cosmo/patch/mom2/feedback.f90 !################################################################# !################################################################ !################################################################ !################################################################ #if NDIM==3 subroutine thermal_feedback(ilevel) use pm_commons use amr_commons use hydro_commons use mpi_mod implicit none #ifndef WITHOUTMPI integer::info2,dummy_io #endif integer::ilevel !------------------------------------------------------------------------ ! This routine computes the thermal energy, the kinetic energy and ! the metal mass dumped in the gas by stars (SNII, SNIa, winds). ! This routine is called every fine time step. !------------------------------------------------------------------------ real(dp)::scale_nH,scale_T2,scale_l,scale_d,scale_t,scale_v real(dp)::t_sn_cont,current_time integer::igrid,jgrid,ipart,jpart,next_part,ivar integer::ig,ip,npart1,npart2,icpu,ilun,idim integer,dimension(1:nvector),save::ind_grid,ind_part,ind_grid_part character(LEN=80)::filename,filedir,fileloc,filedirini character(LEN=5)::nchar,ncharcpu logical::file_exist integer,parameter::tag=1120 if(sf_log_properties) then call title(ifout-1,nchar) if(IOGROUPSIZEREP>0) then call title(((myid-1)/IOGROUPSIZEREP)+1,ncharcpu) filedirini='output_'//TRIM(nchar)//'/' filedir='output_'//TRIM(nchar)//'/group_'//TRIM(ncharcpu)//'/' else filedir='output_'//TRIM(nchar)//'/' endif filename=TRIM(filedir)//'stars_'//TRIM(nchar)//'.out' ilun=myid+103 call title(myid,nchar) fileloc=TRIM(filename)//TRIM(nchar) ! Wait for the token #ifndef WITHOUTMPI if(IOGROUPSIZE>0) then if (mod(myid-1,IOGROUPSIZE)/=0) then call MPI_RECV(dummy_io,1,MPI_INTEGER,myid-1-1,tag,& & MPI_COMM_WORLD,MPI_STATUS_IGNORE,info2) end if endif #endif inquire(file=fileloc,exist=file_exist) if(.not.file_exist) then open(ilun, file=fileloc, form='formatted') write(ilun,'(A24)',advance='no') '# event id ilevel mp ' do idim=1,ndim write(ilun,'(A2,I1,A2)',advance='no') 'xp',idim,' ' enddo do idim=1,ndim write(ilun,'(A2,I1,A2)',advance='no') 'vp',idim,' ' enddo do ivar=1,nvar if(ivar.ge.10) then write(ilun,'(A1,I2,A2)',advance='no') 'u',ivar,' ' else write(ilun,'(A1,I1,A2)',advance='no') 'u',ivar,' ' endif enddo write(ilun,'(A5)',advance='no') 'tag ' write(ilun,'(A1)') ' ' else open(ilun, file=fileloc, status='old', position='append', action='write', form='formatted') endif endif if(numbtot(1,ilevel)==0)return if(verbose)write(*,111)ilevel call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) ! Massive star lifetime from Myr to code units if(use_proper_time)then t_sn_cont=20.*1d6*(365.*24.*3600.)/(scale_t/aexp**2) current_time=texp else t_sn_cont=20.*1d6*(365.*24.*3600.)/scale_t current_time=t endif ! Gather star particles only. ! Loop over cpus do icpu=1,ncpu igrid=headl(icpu,ilevel) ig=0 ip=0 ! Loop over grids do jgrid=1,numbl(icpu,ilevel) npart1=numbp(igrid) ! Number of particles in the grid npart2=0 ! Count star particles if(npart1>0)then ipart=headp(igrid) ! Loop over particles do jpart=1,npart1 ! Save next particle <--- Very important !!! next_part=nextp(ipart) if ( is_star(typep(ipart)) .and. tp(ipart).ge.(current_time-t_sn_cont)) then npart2=npart2+1 endif ipart=next_part ! Go to next particle end do endif ! Gather star particles if(npart2>0)then ig=ig+1 ind_grid(ig)=igrid ipart=headp(igrid) ! Loop over particles do jpart=1,npart1 ! Save next particle <--- Very important !!! next_part=nextp(ipart) ! Select only star particles if ( is_star(typep(ipart)) .and. tp(ipart).ge.(current_time-t_sn_cont)) then if(ig==0)then ig=1 ind_grid(ig)=igrid end if ip=ip+1 ind_part(ip)=ipart ind_grid_part(ip)=ig endif if(ip==nvector)then call feedbk(ind_grid,ind_part,ind_grid_part,ig,ip,ilevel) ip=0 ig=0 end if ipart=next_part ! Go to next particle end do ! End loop over particles end if igrid=next(igrid) ! Go to next grid end do ! End loop over grids if(ip>0)call feedbk(ind_grid,ind_part,ind_grid_part,ig,ip,ilevel) end do ! End loop over cpus if(sf_log_properties) close(ilun) 111 format(' Entering thermal_feedback for level ',I2) end subroutine thermal_feedback #endif !################################################################ !################################################################ !################################################################ !################################################################ #if NDIM==3 subroutine feedbk(ind_grid,ind_part,ind_grid_part,ng,np,ilevel) use amr_commons use pm_commons use hydro_commons use random use constants, only: M_sun, Myr2sec, pc2cm implicit none integer::ng,np,ilevel integer,dimension(1:nvector)::ind_grid integer,dimension(1:nvector)::ind_grid_part,ind_part !----------------------------------------------------------------------- ! This routine is called by subroutine feedback. Each stellar particle ! dumps mass, momentum and energy in the nearest grid cell using array ! unew. !----------------------------------------------------------------------- integer::i,j,idim,nx_loc,ivar,ilun real(kind=8)::RandNum real(dp)::mstar,dx_min,vol_min real(dp)::t0,ESN,mejecta,zloss,e,uvar real(dp)::msne_min,mstar_max,FRAC_NT real(dp)::pressure,gas_density,metallicity real(dp)::p_SN,p_SN_z_exp,p_SN_n_exp,r_c,r_c_z_exp,r_c_n_exp real(dp)::M_SINGLE_SN,mpart_ini,cs_H2_2,p_boost,r_cool real(dp)::dx,dx_loc,scale,birth_time,current_time,t_sn_cont,avg_n,n_dot real(dp)::scale_nH,scale_T2,scale_l,scale_d,scale_t,scale_v ! Grid based arrays real(dp),dimension(1:nvector,1:ndim),save::x0 integer ,dimension(1:nvector),save::ind_cell integer ,dimension(1:nvector,1:threetondim),save::nbors_father_cells integer ,dimension(1:nvector,1:twotondim),save::nbors_father_grids ! Particle based arrays logical,dimension(1:nvector),save::ok real(dp),dimension(1:nvector),save::mloss,mzloss,ethermal,ekinetic,dteff real(dp),dimension(1:nvector),save::vol_loc real(dp),dimension(1:nvector,1:ndim),save::x integer ,dimension(1:nvector,1:ndim),save::id,igd,icd integer ,dimension(1:nvector),save::igrid,indp,n_SN integer ,dimension(1:nvector,1:twotondim),save::icell,kg real(dp),dimension(1:3)::skip_loc integer ,dimension(1:nvector),save::hra real(dp),dimension(1:nvector,1:ndim),save::dd,dg integer ,dimension(1:nvector,1:ndim),save::ig,igg,icg #if NENER>0 integer::irad #endif if(sf_log_properties) ilun=myid+103 ! Conversion factor from user units to cgs units call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) ! Mesh spacing in that level dx=0.5D0**ilevel nx_loc=(icoarse_max-icoarse_min+1) skip_loc=(/0.0d0,0.0d0,0.0d0/) if(ndim>0)skip_loc(1)=dble(icoarse_min) if(ndim>1)skip_loc(2)=dble(jcoarse_min) if(ndim>2)skip_loc(3)=dble(kcoarse_min) scale=boxlen/dble(nx_loc) dx_loc=dx*scale vol_loc(1:nvector)=dx_loc**ndim dx_min=(0.5D0**nlevelmax)*scale vol_min=dx_min**ndim ! Minimum star particle mass if(m_star < 0d0)then mstar=n_star/(scale_nH*aexp**3)*vol_min else mstar=m_star*mass_sph endif msne_min=mass_sne_min*M_sun/(scale_d*scale_l**3) mstar_max=mass_star_max*M_sun/(scale_d*scale_l**3) ! Massive star lifetime from Myr to code units if(use_proper_time)then t0=t_sne*1d6*(365.*24.*3600.)/(scale_t/aexp**2) t_sn_cont=20.*1d6*(365.*24.*3600.)/(scale_t/aexp**2) current_time=texp else t0=t_sne*1d6*(365.*24.*3600.)/scale_t t_sn_cont=20.*1d6*(365.*24.*3600.)/scale_t current_time=t endif ! Type II supernova specific energy from cgs to code units ESN=2.703d51/(10.*M_sun)/scale_v**2 ! Type II supernova average mass from cgs to code units M_SINGLE_SN=(10.*M_sun)/(scale_d*scale_l**3) if(momentum_feedback>0)then SELECT CASE (momentum_feedback) CASE (1) ! inhomogeneous medium (weak) ! momentum p_SN = 1.11*1d5*1d5*M_sun/(scale_v*scale_d*scale_l**3) p_SN_z_exp = -0.114 p_SN_n_exp = -0.190 ! cooling radius r_c = 6.3 * pc2cm / scale_l r_c_z_exp = -0.05 r_c_n_exp = -0.42 CASE (2) ! homogeneous medium (strong) ! momentum p_SN = 1.42*1d5*1d5*M_sun/(scale_v*scale_d*scale_l**3) p_SN_z_exp = -0.137 p_SN_n_exp = -0.160 ! cooling radius r_c = 3.0 * pc2cm / scale_l r_c_z_exp = -0.082 r_c_n_exp = -0.42 END SELECT endif ! Photoionization momentum injection from cgs to code units cs_H2_2=(22.0*1d5/scale_v)**2 ! 22 km/s ! Fraction of the SN energy into non-thermal component FRAC_NT=0.0 ! Lower left corner of 3x3x3 grid-cube do idim=1,ndim do i=1,ng x0(i,idim)=xg(ind_grid(i),idim)-3.0D0*dx end do end do ! Gather 27 neighboring father cells (should be present anytime !) do i=1,ng ind_cell(i)=father(ind_grid(i)) end do call get3cubefather(ind_cell,nbors_father_cells,nbors_father_grids,ng,ilevel) ! Rescale position at level ilevel do idim=1,ndim do j=1,np x(j,idim)=xp(ind_part(j),idim)/scale+skip_loc(idim) end do end do do idim=1,ndim do j=1,np x(j,idim)=x(j,idim)-x0(ind_grid_part(j),idim) end do end do do idim=1,ndim do j=1,np x(j,idim)=x(j,idim)/dx end do end do if(momentum_feedback>0)then ! CIC at level ilevel (dd: right cloud boundary; dg: left cloud boundary) do idim=1,ndim do j=1,np dd(j,idim)=x(j,idim)+0.5D0 id(j,idim)=int(dd(j,idim)) dd(j,idim)=dd(j,idim)-id(j,idim) dg(j,idim)=1.0D0-dd(j,idim) ig(j,idim)=id(j,idim)-1 end do end do ! Compute parent grids do idim=1,ndim do j=1,np igg(j,idim)=ig(j,idim)/2 igd(j,idim)=id(j,idim)/2 end do end do do j=1,np kg(j,1)=1+igg(j,1)+3*igg(j,2)+9*igg(j,3) kg(j,2)=1+igd(j,1)+3*igg(j,2)+9*igg(j,3) kg(j,3)=1+igg(j,1)+3*igd(j,2)+9*igg(j,3) kg(j,4)=1+igd(j,1)+3*igd(j,2)+9*igg(j,3) kg(j,5)=1+igg(j,1)+3*igg(j,2)+9*igd(j,3) kg(j,6)=1+igd(j,1)+3*igg(j,2)+9*igd(j,3) kg(j,7)=1+igg(j,1)+3*igd(j,2)+9*igd(j,3) kg(j,8)=1+igd(j,1)+3*igd(j,2)+9*igd(j,3) end do do j=1,np call ranf(localseed,RandNum) hra(j) = int(RandNum*8)+1 enddo do j=1,np igrid(j)=son(nbors_father_cells(ind_grid_part(j),kg(j,hra(j)))) end do ! Check if particles are entirely in level ilevel ok(1:np)=.true. do j=1,np ok(j)=ok(j).and.igrid(j)>0 end do ! Compute parent cell position do idim=1,ndim do j=1,np if(ok(j))then icg(j,idim)=ig(j,idim)-2*igg(j,idim) icd(j,idim)=id(j,idim)-2*igd(j,idim) endif end do end do do j=1,np if(ok(j))then icell(j,1)=1+icg(j,1)+2*icg(j,2)+4*icg(j,3) icell(j,2)=1+icd(j,1)+2*icg(j,2)+4*icg(j,3) icell(j,3)=1+icg(j,1)+2*icd(j,2)+4*icg(j,3) icell(j,4)=1+icd(j,1)+2*icd(j,2)+4*icg(j,3) icell(j,5)=1+icg(j,1)+2*icg(j,2)+4*icd(j,3) icell(j,6)=1+icd(j,1)+2*icg(j,2)+4*icd(j,3) icell(j,7)=1+icg(j,1)+2*icd(j,2)+4*icd(j,3) icell(j,8)=1+icd(j,1)+2*icd(j,2)+4*icd(j,3) endif end do ! Compute parent cell adresses do j=1,np if(ok(j))then indp(j)=ncoarse+(icell(j,hra(j))-1)*ngridmax+igrid(j) else indp(j) = nbors_father_cells(ind_grid_part(j),kg(j,hra(j))) vol_loc(j)=vol_loc(j)*2**ndim ! ilevel-1 cell volume end if end do else ! NGP at level ilevel do idim=1,ndim do j=1,np id(j,idim)=int(x(j,idim)) end do end do ! Compute parent grids do idim=1,ndim do j=1,np igd(j,idim)=id(j,idim)/2 end do end do do j=1,np kg(j,1)=1+igd(j,1)+3*igd(j,2)+9*igd(j,3) end do do j=1,np igrid(j)=son(nbors_father_cells(ind_grid_part(j),kg(j,1))) end do ! Check if particles are entirely in level ilevel ok(1:np)=.true. do j=1,np ok(j)=ok(j).and.igrid(j)>0 end do ! Compute parent cell position do idim=1,ndim do j=1,np if(ok(j))then icd(j,idim)=id(j,idim)-2*igd(j,idim) end if end do end do do j=1,np if(ok(j))then icell(j,1)=1+icd(j,1)+2*icd(j,2)+4*icd(j,3) end if end do ! Compute parent cell adresses do j=1,np if(ok(j))then indp(j)=ncoarse+(icell(j,1)-1)*ngridmax+igrid(j) else indp(j) = nbors_father_cells(ind_grid_part(j),kg(j,1)) vol_loc(j)=vol_loc(j)*2**ndim ! ilevel-1 cell volume end if end do endif ! Compute individual time steps do j=1,np dteff(j)=dtnew(levelp(ind_part(j))) end do if(use_proper_time)then do j=1,np dteff(j)=dteff(j)*aexp**2 end do endif ! Reset ejected mass, metallicity, thermal energy do j=1,np mloss(j)=0d0 mzloss(j)=0d0 ethermal(j)=0d0 end do ! Compute stellar mass loss and thermal feedback due to supernovae if(f_w==0)then do j=1,np n_SN(j)=0 birth_time=tp(ind_part(j)) if(birth_time.lt.(current_time-t0))then ! Guesstimate the initial star particle mass mpart_ini=mp(ind_part(j))/(1.0-eta_sn*(current_time-t0-birth_time)/(t_sn_cont-t0)) ! Compute the constant SN rate n_dot=eta_sn*mpart_ini/M_SINGLE_SN/(t_sn_cont-t0) ! Compute the mean SN count avg_n=n_dot*dteff(j) ! Draw a Poisson process for this mean call poissdev(localseed,avg_n,n_SN(j)) ! Stellar mass loss mejecta=n_SN(j)*M_SINGLE_SN mloss(j)=mloss(j)+mejecta/vol_loc(j) ! Thermal energy ethermal(j)=ethermal(j)+mejecta*ESN/vol_loc(j) ! Metallicity if(metal)then zloss=yield+(1d0-yield)*zp(ind_part(j)) mzloss(j)=mzloss(j)+mejecta*zloss/vol_loc(j) endif ! Reduce star particle mass mp(ind_part(j))=mp(ind_part(j))-mejecta if(sf_log_properties) then write(ilun,'(I10)',advance='no') 1 write(ilun,'(2I10,E24.12)',advance='no') idp(ind_part(j)),ilevel,mp(ind_part(j)) do idim=1,ndim write(ilun,'(E24.12)',advance='no') xp(ind_part(j),idim) enddo do idim=1,ndim write(ilun,'(E24.12)',advance='no') vp(ind_part(j),idim) enddo write(ilun,'(E24.12)',advance='no') unew(indp(j),1) do ivar=2,nvar if(ivar.eq.ndim+2)then e=0.0d0 do idim=1,ndim e=e+0.5d0*unew(ind_cell(i),idim+1)**2/max(unew(ind_cell(i),1),smallr) enddo #if NENER>0 do irad=0,nener-1 e=e+unew(ind_cell(i),inener+irad) enddo #endif #ifdef SOLVERmhd do idim=1,ndim e=e+0.125d0*(unew(ind_cell(i),idim+ndim+2)+unew(ind_cell(i),idim+nvar))**2 enddo #endif ! Temperature uvar=(gamma-1.0d0)*(unew(ind_cell(i),ndim+2)-e)*scale_T2 else uvar=unew(indp(j),ivar) endif write(ilun,'(E24.12)',advance='no') uvar/unew(indp(j),1) enddo write(ilun,'(I10)',advance='no') typep(ind_part(i))%tag write(ilun,'(A1)') ' ' endif endif end do ! Photo-ionization thermal feedback do j=1,np pressure=max(uold(indp(j),1),smallr)*cs_H2_2 ethermal(j)=ethermal(j)+pressure end do ! Use stellar momentum feedback if(momentum_feedback>0)then ! Momentum feedback from supernovae do j=1,np birth_time=tp(ind_part(j)) gas_density=max(uold(indp(j),1),smallr) ! Compute metallicity for cooling if(metal)then metallicity=max(uold(indp(j),imetal),smallr)/gas_density/0.02 else metallicity=z_ave endif metallicity=max(metallicity,0.01) p_boost=metallicity**(p_SN_z_exp)*(gas_density*scale_nH/100.0)**(p_SN_n_exp) r_cool=r_c*metallicity**(r_c_z_exp)*(gas_density*scale_nH/100.0)**(r_c_n_exp) if(birth_time.lt.(current_time-t0).and.r_cool.lt.(4.0*dx_min/aexp))then pstarnew(indp(j))=pstarnew(indp(j))+p_SN*n_SN(j)*p_boost*min(1.0,(dx_min/r_cool/aexp)**(3.0/2.0))/dx_loc**3 endif end do endif endif ! Update hydro variables due to feedback ! Loop on particles do j=1,np ! Specific kinetic energy of the star ekinetic(j)=0.5d0*(vp(ind_part(j),1)**2 & & +vp(ind_part(j),2)**2 & & +vp(ind_part(j),3)**2) ! Update hydro variable in NGP cell unew(indp(j),1)=unew(indp(j),1)+mloss(j) unew(indp(j),2)=unew(indp(j),2)+mloss(j)*vp(ind_part(j),1) unew(indp(j),3)=unew(indp(j),3)+mloss(j)*vp(ind_part(j),2) unew(indp(j),4)=unew(indp(j),4)+mloss(j)*vp(ind_part(j),3) unew(indp(j),5)=unew(indp(j),5)+mloss(j)*ekinetic(j)+ethermal(j) ! turbulent energy ! unew(indp(j),ivirial1)=unew(indp(j),ivirial1)+ethermal(j)*0.1 ! Update internal energy if(pressure_fix)then enew(indp(j))=enew(indp(j))+(1.0-FRAC_NT)*ethermal(j) endif ! Add a fraction of the SN energy to the non-thermal energy if(nener>0)then unew(indp(j),ndim+3)=unew(indp(j),ndim+3)+FRAC_NT*ethermal(j) endif end do ! Add metals if(metal)then do j=1,np unew(indp(j),imetal)=unew(indp(j),imetal)+mzloss(j) end do endif ! Add delayed cooling switch variable if(delayed_cooling)then do j=1,np unew(indp(j),idelay)=unew(indp(j),idelay)+mloss(j) end do endif end subroutine feedbk #endif !################################################################ !################################################################ !################################################################ !################################################################ subroutine kinetic_feedback use amr_commons use pm_commons use hydro_commons use constants, only:Myr2sec use mpi_mod implicit none #ifndef WITHOUTMPI integer::info integer,dimension(1:ncpu)::nSN_icpu_all real(dp),dimension(:),allocatable::mSN_all,sSN_all,ZSN_all real(dp),dimension(:,:),allocatable::xSN_all,vSN_all #endif !---------------------------------------------------------------------- ! This subroutine compute the kinetic feedback due to SNII and ! imolement this using exploding GMC particles. ! This routine is called only at coarse time step. ! Yohan Dubois !---------------------------------------------------------------------- ! local constants integer::ip,icpu,igrid,jgrid,npart1,npart2,ipart,jpart,next_part integer::nSN,nSN_loc,nSN_tot,iSN,ilevel,ivar integer,dimension(1:ncpu)::nSN_icpu real(dp)::scale_nH,scale_T2,scale_l,scale_d,scale_t,scale_v,t0 real(dp)::current_time real(dp)::scale,dx_min,vol_min,mstar integer::nx_loc integer,dimension(:),allocatable::ind_part,ind_grid logical,dimension(:),allocatable::ok_free integer,dimension(:),allocatable::indSN real(dp),dimension(:),allocatable::mSN,sSN,ZSN,m_gas,vol_gas,ekBlast real(dp),dimension(:,:),allocatable::xSN,vSN,u_gas,dq if(.not. hydro)return if(ndim.ne.3)return if(verbose)write(*,*)'Entering make_sn' ! Conversion factor from user units to cgs units call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) ! Mesh spacing in that level nx_loc=(icoarse_max-icoarse_min+1) scale=boxlen/dble(nx_loc) dx_min=(0.5D0**nlevelmax)*scale vol_min=dx_min**ndim ! Minimum star particle mass if(m_star < 0d0)then mstar=n_star/(scale_nH*aexp**3)*vol_min else mstar=m_star*mass_sph endif ! Lifetime of Giant Molecular Clouds from Myr to code units ! Massive star lifetime from Myr to code units if(use_proper_time)then t0=t_sne*Myr2sec/(scale_t/aexp**2) current_time=texp else t0=t_sne*Myr2sec/scale_t current_time=t endif !------------------------------------------------------ ! Gather GMC particles eligible for disruption !------------------------------------------------------ nSN_loc=0 ! Loop over levels do icpu=1,ncpu ! Loop over cpus igrid=headl(icpu,levelmin) ! Loop over grids do jgrid=1,numbl(icpu,levelmin) npart1=numbp(igrid) ! Number of particles in the grid npart2=0 ! Count old enough GMC particles if(npart1>0)then ipart=headp(igrid) ! Loop over particles do jpart=1,npart1 ! Save next particle <--- Very important !!! next_part=nextp(ipart) if ( is_debris(typep(ipart)) .and. tp(ipart).lt.(current_time-t0) ) then npart2=npart2+1 endif ipart=next_part ! Go to next particle end do endif nSN_loc=nSN_loc+npart2 ! Add SNe to the total igrid=next(igrid) ! Go to next grid end do end do ! End loop over levels nSN_icpu=0 nSN_icpu(myid)=nSN_loc #ifndef WITHOUTMPI ! Give an array of number of SN on each cpu available to all cpus call MPI_ALLREDUCE(nSN_icpu,nSN_icpu_all,ncpu,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,info) nSN_icpu=nSN_icpu_all #endif nSN_tot=sum(nSN_icpu(1:ncpu)) if (nSN_tot .eq. 0) return if(myid==1)then write(*,*)'-----------------------------------------------' write(*,*)'Number of GMC to explode=',nSN_tot write(*,*)'-----------------------------------------------' endif ! Allocate arrays for the position and the mass of the SN allocate(xSN(1:nSN_tot,1:3),vSN(1:nSN_tot,1:3)) allocate(mSN(1:nSN_tot),sSN(1:nSN_tot),ZSN(1:nSN_tot)) xSN=0; vSN=0; mSN=0; sSN=0; ZSN=0 ! Allocate arrays for particles index and parent grid if(nSN_loc>0)then allocate(ind_part(1:nSN_loc),ind_grid(1:nSN_loc),ok_free(1:nSN_loc)) endif !------------------------------------------------------ ! Store position and mass of the GMC into the SN array !------------------------------------------------------ if(myid==1)then iSN=0 else iSN=sum(nSN_icpu(1:myid-1)) endif ! Loop over levels ip=0 do icpu=1,ncpu igrid=headl(icpu,levelmin) ! Loop over grids do jgrid=1,numbl(icpu,levelmin) npart1=numbp(igrid) ! Number of particles in the grid ! Count old enough star particles that have not exploded if(npart1>0)then ipart=headp(igrid) ! Loop over particles do jpart=1,npart1 ! Save next particle <--- Very important !!! next_part=nextp(ipart) if ( is_debris(typep(ipart)) .and. tp(ipart).lt.(current_time-t0) ) then iSN=iSN+1 xSN(iSN,1)=xp(ipart,1) xSN(iSN,2)=xp(ipart,2) xSN(iSN,3)=xp(ipart,3) vSN(iSN,1)=vp(ipart,1) vSN(iSN,2)=vp(ipart,2) vSN(iSN,3)=vp(ipart,3) mSN(iSN)=mp(ipart) sSN(iSN)=dble(-idp(ipart))*mstar if(metal)ZSN(iSN)=zp(ipart) ip=ip+1 ind_grid(ip)=igrid ind_part(ip)=ipart endif ipart=next_part ! Go to next particle end do endif igrid=next(igrid) ! Go to next grid end do end do ! End loop over levels ! Remove GMC particle if(nSN_loc>0)then ok_free=.true. call remove_list(ind_part,ind_grid,ok_free,nSN_loc) call add_free_cond(ind_part,ok_free,nSN_loc) deallocate(ind_part,ind_grid,ok_free) endif #ifndef WITHOUTMPI allocate(xSN_all(1:nSN_tot,1:3),vSN_all(1:nSN_tot,1:3),mSN_all(1:nSN_tot),sSN_all(1:nSN_tot),ZSN_all(1:nSN_tot)) call MPI_ALLREDUCE(xSN,xSN_all,nSN_tot*3,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) call MPI_ALLREDUCE(vSN,vSN_all,nSN_tot*3,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) call MPI_ALLREDUCE(mSN,mSN_all,nSN_tot ,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) call MPI_ALLREDUCE(sSN,sSN_all,nSN_tot ,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) call MPI_ALLREDUCE(ZSN,ZSN_all,nSN_tot ,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) xSN=xSN_all vSN=vSN_all mSN=mSN_all sSN=sSN_all ZSN=ZSN_all deallocate(xSN_all,vSN_all,mSN_all,sSN_all,ZSN_all) #endif nSN=nSN_tot allocate(m_gas(1:nSN),u_gas(1:nSN,1:3),vol_gas(1:nSN),dq(1:nSN,1:3),ekBlast(1:nSN)) allocate(indSN(1:nSN)) ! Compute the grid discretization effects call average_SN(xSN,vol_gas,dq,ekBlast,indSN,nSN) ! Modify hydro quantities to account for a Sedov blast wave call Sedov_blast(xSN,vSN,mSN,sSN,ZSN,indSN,vol_gas,dq,ekBlast,nSN) deallocate(xSN,vSN,mSN,sSN,ZSN,indSN,m_gas,u_gas,vol_gas,dq,ekBlast) ! Update hydro quantities for split cells do ilevel=nlevelmax,levelmin,-1 call upload_fine(ilevel) do ivar=1,nvar call make_virtual_fine_dp(uold(1,ivar),ilevel) enddo enddo end subroutine kinetic_feedback !################################################################ !################################################################ !################################################################ !################################################################ subroutine average_SN(xSN,vol_gas,dq,ekBlast,ind_blast,nSN) use pm_commons use amr_commons use hydro_commons use constants, only: pc2cm use mpi_mod implicit none #ifndef WITHOUTMPI integer::info #endif !------------------------------------------------------------------------ ! This routine average the hydro quantities inside the SN bubble !------------------------------------------------------------------------ integer::ilevel,ncache,nSN,iSN,ind,ix,iy,iz,ngrid,iskip integer::i,nx_loc,igrid integer,dimension(1:nvector),save::ind_grid,ind_cell real(dp)::x,y,z,dr_SN,u,v,w,u2,v2,w2,dr_cell real(dp)::scale,dx,dxx,dyy,dzz,dx_min,dx_loc,vol_loc,rmax2,rmax real(dp)::scale_nH,scale_T2,scale_l,scale_d,scale_t,scale_v real(dp),dimension(1:3)::skip_loc real(dp),dimension(1:twotondim,1:3)::xc integer ,dimension(1:nSN)::ind_blast real(dp),dimension(1:nSN)::vol_gas,ekBlast real(dp),dimension(1:nSN,1:3)::xSN,dq,u2Blast #ifndef WITHOUTMPI real(dp),dimension(1:nSN)::vol_gas_all,ekBlast_all real(dp),dimension(1:nSN,1:3)::dq_all,u2Blast_all #endif logical ,dimension(1:nvector),save::ok if(nSN==0)return if(verbose)write(*,*)'Entering average_SN' ! Mesh spacing in that level nx_loc=(icoarse_max-icoarse_min+1) skip_loc=(/0.0d0,0.0d0,0.0d0/) skip_loc(1)=dble(icoarse_min) skip_loc(2)=dble(jcoarse_min) skip_loc(3)=dble(kcoarse_min) scale=boxlen/dble(nx_loc) dx_min=scale*0.5D0**nlevelmax ! Conversion factor from user units to cgs units call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) ! Maximum radius of the ejecta rmax=MAX(2.0d0*dx_min*scale_l/aexp,rbubble*pc2cm) rmax=rmax/scale_l rmax2=rmax*rmax ! Initialize the averaged variables vol_gas=0; dq=0; u2Blast=0; ekBlast=0; ind_blast=-1 ! Loop over levels do ilevel=levelmin,nlevelmax ! Computing local volume (important for averaging hydro quantities) dx=0.5D0**ilevel dx_loc=dx*scale vol_loc=dx_loc**ndim ! Cells center position relative to grid center position do ind=1,twotondim iz=(ind-1)/4 iy=(ind-1-4*iz)/2 ix=(ind-1-2*iy-4*iz) xc(ind,1)=(dble(ix)-0.5D0)*dx xc(ind,2)=(dble(iy)-0.5D0)*dx xc(ind,3)=(dble(iz)-0.5D0)*dx end do ! Loop over grids ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Loop over cells do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do ! Flag leaf cells do i=1,ngrid ok(i)=son(ind_cell(i))==0 end do do i=1,ngrid if(ok(i))then ! Get gas cell position x=(xg(ind_grid(i),1)+xc(ind,1)-skip_loc(1))*scale y=(xg(ind_grid(i),2)+xc(ind,2)-skip_loc(2))*scale z=(xg(ind_grid(i),3)+xc(ind,3)-skip_loc(3))*scale do iSN=1,nSN ! Check if the cell lies within the SN radius dxx=x-xSN(iSN,1) dyy=y-xSN(iSN,2) dzz=z-xSN(iSN,3) dr_SN=dxx**2+dyy**2+dzz**2 dr_cell=MAX(ABS(dxx),ABS(dyy),ABS(dzz)) if(dr_SN.lt.rmax2)then vol_gas(iSN)=vol_gas(iSN)+vol_loc ! Take account for grid effects on the conservation of the ! normalized linear momentum u=dxx/rmax v=dyy/rmax w=dzz/rmax ! Add the local normalized linear momentum to the total linear ! momentum of the blast wave (should be zero with no grid effect) dq(iSN,1)=dq(iSN,1)+u*vol_loc dq(iSN,2)=dq(iSN,2)+v*vol_loc dq(iSN,3)=dq(iSN,3)+w*vol_loc u2Blast(iSN,1)=u2Blast(iSN,1)+u*u*vol_loc u2Blast(iSN,2)=u2Blast(iSN,2)+v*v*vol_loc u2Blast(iSN,3)=u2Blast(iSN,3)+w*w*vol_loc endif if(dr_cell.le.dx_loc/2.0)then ind_blast(iSN)=ind_cell(i) ekBlast (iSN)=vol_loc endif end do endif end do end do ! End loop over cells end do ! End loop over grids end do ! End loop over levels #ifndef WITHOUTMPI call MPI_ALLREDUCE(vol_gas,vol_gas_all,nSN ,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) call MPI_ALLREDUCE(dq ,dq_all ,nSN*3,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) call MPI_ALLREDUCE(u2Blast,u2Blast_all,nSN*3,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) call MPI_ALLREDUCE(ekBlast,ekBlast_all,nSN ,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) vol_gas=vol_gas_all dq =dq_all u2Blast=u2Blast_all ekBlast=ekBlast_all #endif do iSN=1,nSN if(vol_gas(iSN)>0d0)then dq(iSN,1)=dq(iSN,1)/vol_gas(iSN) dq(iSN,2)=dq(iSN,2)/vol_gas(iSN) dq(iSN,3)=dq(iSN,3)/vol_gas(iSN) u2Blast(iSN,1)=u2Blast(iSN,1)/vol_gas(iSN) u2Blast(iSN,2)=u2Blast(iSN,2)/vol_gas(iSN) u2Blast(iSN,3)=u2Blast(iSN,3)/vol_gas(iSN) u2=u2Blast(iSN,1)-dq(iSN,1)**2 v2=u2Blast(iSN,2)-dq(iSN,2)**2 w2=u2Blast(iSN,3)-dq(iSN,3)**2 ekBlast(iSN)=max(0.5d0*(u2+v2+w2),0.0d0) endif end do if(verbose)write(*,*)'Exiting average_SN' end subroutine average_SN !################################################################ !################################################################ !################################################################ !################################################################ subroutine Sedov_blast(xSN,vSN,mSN,sSN,ZSN,indSN,vol_gas,dq,ekBlast,nSN) use pm_commons use amr_commons use hydro_commons use constants, only: M_sun, pc2cm use mpi_mod implicit none !------------------------------------------------------------------------ ! This routine merges SN using the FOF algorithm. !------------------------------------------------------------------------ integer::ilevel,iSN,nSN,ind,ix,iy,iz,ngrid,iskip integer::i,nx_loc,igrid,ncache integer,dimension(1:nvector),save::ind_grid,ind_cell real(dp)::x,y,z,dx,dxx,dyy,dzz,dr_SN,u,v,w,ESN,mstar,eta_sn2,msne_min,mstar_max real(dp)::scale,dx_min,dx_loc,vol_loc,rmax2,rmax,vol_min real(dp)::scale_nH,scale_T2,scale_l,scale_d,scale_t,scale_v real(dp),dimension(1:3)::skip_loc real(dp),dimension(1:twotondim,1:3)::xc real(dp),dimension(1:nSN)::mSN,sSN,ZSN,p_gas,d_gas,d_metal,vol_gas,uSedov,ekBlast real(dp),dimension(1:nSN,1:3)::xSN,vSN,dq integer ,dimension(1:nSN)::indSN logical ,dimension(1:nvector),save::ok if(nSN==0)return if(verbose)write(*,*)'Entering Sedov_blast' ! Mesh spacing in that level nx_loc=(icoarse_max-icoarse_min+1) skip_loc=(/0.0d0,0.0d0,0.0d0/) skip_loc(1)=dble(icoarse_min) skip_loc(2)=dble(jcoarse_min) skip_loc(3)=dble(kcoarse_min) scale=boxlen/dble(nx_loc) dx_min=scale*0.5D0**nlevelmax vol_min=dx_min**ndim ! Conversion factor from user units to cgs units call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) ! Maximum radius of the ejecta rmax=MAX(2.0d0*dx_min*scale_l/aexp,rbubble*pc2cm) rmax=rmax/scale_l rmax2=rmax*rmax ! Minimum star particle mass if(m_star < 0d0)then mstar=n_star/(scale_nH*aexp**3)*vol_min else mstar=m_star*mass_sph endif msne_min=mass_sne_min*M_sun/(scale_d*scale_l**3) mstar_max=mass_star_max*M_sun/(scale_d*scale_l**3) ! Supernova specific energy from cgs to code units ESN=(1d51/(10d0*M_sun))/scale_v**2 do iSN=1,nSN eta_sn2 = eta_sn if(sf_imf)then if(mSN(iSN).le.mstar_max)then if(mSN(iSN).ge.msne_min) eta_sn2 = eta_ssn if(mSN(iSN).lt.msne_min) eta_sn2 = 0 endif endif if(vol_gas(iSN)>0d0)then d_gas(iSN)=mSN(iSN)/vol_gas(iSN) if(metal)d_metal(iSN)=ZSN(iSN)*mSN(iSN)/vol_gas(iSN) if(ekBlast(iSN)==0d0)then p_gas(iSN)=eta_sn2*sSN(iSN)*ESN/vol_gas(iSN) uSedov(iSN)=0d0 else p_gas(iSN)=(1d0-f_ek)*eta_sn2*sSN(iSN)*ESN/vol_gas(iSN) uSedov(iSN)=sqrt(f_ek*eta_sn2*sSN(iSN)*ESN/mSN(iSN)/ekBlast(iSN)) endif else d_gas(iSN)=mSN(iSN)/ekBlast(iSN) p_gas(iSN)=eta_sn2*sSN(iSN)*ESN/ekBlast(iSN) if(metal)d_metal(iSN)=ZSN(iSN)*mSN(iSN)/ekBlast(iSN) endif end do ! Loop over levels do ilevel=levelmin,nlevelmax ! Computing local volume (important for averaging hydro quantities) dx=0.5D0**ilevel dx_loc=dx*scale vol_loc=dx_loc**ndim ! Cells center position relative to grid center position do ind=1,twotondim iz=(ind-1)/4 iy=(ind-1-4*iz)/2 ix=(ind-1-2*iy-4*iz) xc(ind,1)=(dble(ix)-0.5D0)*dx xc(ind,2)=(dble(iy)-0.5D0)*dx xc(ind,3)=(dble(iz)-0.5D0)*dx end do ! Loop over grids ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Loop over cells do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do ! Flag leaf cells do i=1,ngrid ok(i)=son(ind_cell(i))==0 end do do i=1,ngrid if(ok(i))then ! Get gas cell position x=(xg(ind_grid(i),1)+xc(ind,1)-skip_loc(1))*scale y=(xg(ind_grid(i),2)+xc(ind,2)-skip_loc(2))*scale z=(xg(ind_grid(i),3)+xc(ind,3)-skip_loc(3))*scale do iSN=1,nSN ! Check if the cell lies within the SN radius dxx=x-xSN(iSN,1) dyy=y-xSN(iSN,2) dzz=z-xSN(iSN,3) dr_SN=dxx**2+dyy**2+dzz**2 if(dr_SN.lt.rmax2)then ! Compute the mass density in the cell uold(ind_cell(i),1)=uold(ind_cell(i),1)+d_gas(iSN) ! Compute the metal density in the cell if(metal)uold(ind_cell(i),imetal)=uold(ind_cell(i),imetal)+d_metal(iSN) ! Velocity at a given dr_SN linearly interpolated between zero and uSedov u=uSedov(iSN)*(dxx/rmax-dq(iSN,1))+vSN(iSN,1) v=uSedov(iSN)*(dyy/rmax-dq(iSN,2))+vSN(iSN,2) w=uSedov(iSN)*(dzz/rmax-dq(iSN,3))+vSN(iSN,3) ! Add each momentum component of the blast wave to the gas uold(ind_cell(i),2)=uold(ind_cell(i),2)+d_gas(iSN)*u uold(ind_cell(i),3)=uold(ind_cell(i),3)+d_gas(iSN)*v uold(ind_cell(i),4)=uold(ind_cell(i),4)+d_gas(iSN)*w ! Finally update the total energy of the gas uold(ind_cell(i),5)=uold(ind_cell(i),5)+0.5d0*d_gas(iSN)*(u*u+v*v+w*w)+p_gas(iSN) endif end do endif end do end do ! End loop over cells end do ! End loop over grids end do ! End loop over levels do iSN=1,nSN if(vol_gas(iSN)==0d0)then u=vSN(iSN,1) v=vSN(iSN,2) w=vSN(iSN,3) if(indSN(iSN)>0)then uold(indSN(iSN),1)=uold(indSN(iSN),1)+d_gas(iSN) uold(indSN(iSN),2)=uold(indSN(iSN),2)+d_gas(iSN)*u uold(indSN(iSN),3)=uold(indSN(iSN),3)+d_gas(iSN)*v uold(indSN(iSN),4)=uold(indSN(iSN),4)+d_gas(iSN)*w uold(indSN(iSN),5)=uold(indSN(iSN),5)+d_gas(iSN)*0.5d0*(u*u+v*v+w*w)+p_gas(iSN) if(metal)uold(indSN(iSN),imetal)=uold(indSN(iSN),imetal)+d_metal(iSN) endif endif end do if(verbose)write(*,*)'Exiting Sedov_blast' end subroutine Sedov_blast !########################################################### !########################################################### !########################################################### !########################################################### /home/hemv_g/srf_a/RAMSES_CODE_CORI/ramses06062019_grackle311Cosmo/patch/mom2/godunov_fine.f90 !########################################################### !########################################################### !########################################################### !########################################################### subroutine godunov_fine(ilevel) use amr_commons use hydro_commons implicit none integer::ilevel !-------------------------------------------------------------------------- ! This routine is a wrapper to the second order Godunov solver. ! Small grids (2x2x2) are gathered from level ilevel and sent to the ! hydro solver. On entry, hydro variables are gathered from array uold. ! On exit, unew has been updated. !-------------------------------------------------------------------------- integer::i,igrid,ncache,ngrid integer,dimension(1:nvector),save::ind_grid if(numbtot(1,ilevel)==0)return if(static)return if(verbose)write(*,111)ilevel ! Loop over active grids by vector sweeps ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do call godfine1(ind_grid,ngrid,ilevel) end do 111 format(' Entering godunov_fine for level ',i2) end subroutine godunov_fine !########################################################### !########################################################### !########################################################### !########################################################### subroutine set_unew(ilevel) use amr_commons use hydro_commons implicit none integer::ilevel !-------------------------------------------------------------------------- ! This routine sets array unew to its initial value uold before calling ! the hydro scheme. unew is set to zero in virtual boundaries. !-------------------------------------------------------------------------- integer::i,ivar,ind,icpu,iskip real(dp)::d,u,v,w,e #if NENER>0 integer::irad #endif if(numbtot(1,ilevel)==0)return if(verbose)write(*,111)ilevel ! Set unew to uold for myid cells do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do ivar=1,nvar do i=1,active(ilevel)%ngrid unew(active(ilevel)%igrid(i)+iskip,ivar) = uold(active(ilevel)%igrid(i)+iskip,ivar) end do end do if(momentum_feedback>0)then do i=1,active(ilevel)%ngrid pstarnew(active(ilevel)%igrid(i)+iskip) = 0 end do endif if(pressure_fix)then do i=1,active(ilevel)%ngrid divu(active(ilevel)%igrid(i)+iskip) = 0 end do do i=1,active(ilevel)%ngrid d=max(uold(active(ilevel)%igrid(i)+iskip,1),smallr) u=0; v=0; w=0 if(ndim>0)u=uold(active(ilevel)%igrid(i)+iskip,2)/d if(ndim>1)v=uold(active(ilevel)%igrid(i)+iskip,3)/d if(ndim>2)w=uold(active(ilevel)%igrid(i)+iskip,4)/d e=uold(active(ilevel)%igrid(i)+iskip,ndim+2)-0.5d0*d*(u**2+v**2+w**2) #if NENER>0 do irad=1,nener e=e-uold(active(ilevel)%igrid(i)+iskip,ndim+2+irad) end do #endif enew(active(ilevel)%igrid(i)+iskip)=e end do end if end do ! Set unew to 0 for virtual boundary cells do icpu=1,ncpu do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do ivar=1,nvar do i=1,reception(icpu,ilevel)%ngrid unew(reception(icpu,ilevel)%igrid(i)+iskip,ivar)=0 end do end do if(momentum_feedback>0)then do i=1,reception(icpu,ilevel)%ngrid pstarnew(reception(icpu,ilevel)%igrid(i)+iskip) = 0 end do endif if(pressure_fix)then do i=1,reception(icpu,ilevel)%ngrid divu(reception(icpu,ilevel)%igrid(i)+iskip) = 0 enew(reception(icpu,ilevel)%igrid(i)+iskip) = 0 end do end if end do end do 111 format(' Entering set_unew for level ',i2) end subroutine set_unew !########################################################### !########################################################### !########################################################### !########################################################### subroutine set_uold(ilevel) use amr_commons use hydro_commons use poisson_commons implicit none integer::ilevel !--------------------------------------------------------- ! This routine sets array uold to its new value unew ! after the hydro step. !--------------------------------------------------------- integer::i,ivar,ind,iskip,nx_loc,ind_cell real(dp)::scale,d,u,v,w real(dp)::e_kin,e_cons,e_prim,e_trunc,div,dx #if NENER>0 integer::irad #endif if(numbtot(1,ilevel)==0)return if(verbose)write(*,111)ilevel nx_loc=icoarse_max-icoarse_min+1 scale=boxlen/dble(nx_loc) dx=0.5d0**ilevel*scale ! Add gravity source terms to unew if(poisson)then call add_gravity_source_terms(ilevel) end if ! Add non conservative pdV terms to unew ! for thermal and/or non-thermal energies if(pressure_fix.OR.nener>0)then call add_pdv_source_terms(ilevel) endif ! Add turbulent energy source and sink terms ! to the corresponding passive scalar ! Used for the sugrid scale turbulence model if(sf_model > 0)then call add_viscosity_source_terms(ilevel) endif ! Set uold to unew for myid cells do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do ivar=1,nvar do i=1,active(ilevel)%ngrid uold(active(ilevel)%igrid(i)+iskip,ivar) = unew(active(ilevel)%igrid(i)+iskip,ivar) end do end do if(momentum_feedback>0)then do i=1,active(ilevel)%ngrid pstarold(active(ilevel)%igrid(i)+iskip) = pstarnew(active(ilevel)%igrid(i)+iskip) end do endif if(pressure_fix)then ! Correct total energy if internal energy is too small do i=1,active(ilevel)%ngrid ind_cell=active(ilevel)%igrid(i)+iskip d=max(uold(ind_cell,1),smallr) u=0; v=0; w=0 if(ndim>0)u=uold(ind_cell,2)/d if(ndim>1)v=uold(ind_cell,3)/d if(ndim>2)w=uold(ind_cell,4)/d e_kin=0.5d0*d*(u**2+v**2+w**2) #if NENER>0 do irad=1,nener e_kin=e_kin+uold(ind_cell,ndim+2+irad) end do #endif e_cons=uold(ind_cell,ndim+2)-e_kin e_prim=enew(ind_cell) ! Note: here divu=-div.u*dt div=abs(divu(ind_cell))*dx/dtnew(ilevel) e_trunc=beta_fix*d*max(div,3.0d0*hexp*dx)**2 if(e_cons0)u=unew(ind_cell,2)/d if(ndim>1)v=unew(ind_cell,3)/d if(ndim>2)w=unew(ind_cell,4)/d e_kin=0.5d0*d*(u**2+v**2+w**2) e_prim=unew(ind_cell,ndim+2)-e_kin d_old=max(uold(ind_cell,1),smallr) fact=d_old/d*0.5d0*dtnew(ilevel) if(ndim>0)then u=u+f(ind_cell,1)*fact unew(ind_cell,2)=d*u endif if(ndim>1)then v=v+f(ind_cell,2)*fact unew(ind_cell,3)=d*v end if if(ndim>2)then w=w+f(ind_cell,3)*fact unew(ind_cell,4)=d*w endif e_kin=0.5d0*d*(u**2+v**2+w**2) unew(ind_cell,ndim+2)=e_prim+e_kin end do end do 111 format(' Entering add_gravity_source_terms for level ',i2) end subroutine add_gravity_source_terms !########################################################### !########################################################### !########################################################### !########################################################### subroutine add_pdv_source_terms(ilevel) use amr_commons use hydro_commons implicit none integer::ilevel !--------------------------------------------------------- ! This routine adds the pdV source term to the internal ! energy equation and to the non-thermal energy equations. !--------------------------------------------------------- integer::i,ind,iskip,nx_loc,ind_cell1 integer::ncache,igrid,ngrid,idim,id1,ig1,ih1,id2,ig2,ih2 integer,dimension(1:3,1:2,1:8)::iii,jjj real(dp)::scale,dx,dx_loc,d,u,v,w,eold integer ,dimension(1:nvector),save::ind_grid,ind_cell integer ,dimension(1:nvector,0:twondim),save::igridn integer ,dimension(1:nvector,1:ndim),save::ind_left,ind_right real(dp),dimension(1:nvector,1:ndim,1:ndim),save::velg,veld real(dp),dimension(1:nvector,1:ndim),save::dx_g,dx_d real(dp),dimension(1:nvector),save::divu_loc #if NENER>0 integer::irad #endif if(numbtot(1,ilevel)==0)return if(verbose)write(*,111)ilevel nx_loc=icoarse_max-icoarse_min+1 scale=boxlen/dble(nx_loc) dx=0.5d0**ilevel dx_loc=dx*scale iii(1,1,1:8)=(/1,0,1,0,1,0,1,0/); jjj(1,1,1:8)=(/2,1,4,3,6,5,8,7/) iii(1,2,1:8)=(/0,2,0,2,0,2,0,2/); jjj(1,2,1:8)=(/2,1,4,3,6,5,8,7/) iii(2,1,1:8)=(/3,3,0,0,3,3,0,0/); jjj(2,1,1:8)=(/3,4,1,2,7,8,5,6/) iii(2,2,1:8)=(/0,0,4,4,0,0,4,4/); jjj(2,2,1:8)=(/3,4,1,2,7,8,5,6/) iii(3,1,1:8)=(/5,5,5,5,0,0,0,0/); jjj(3,1,1:8)=(/5,6,7,8,1,2,3,4/) iii(3,2,1:8)=(/0,0,0,0,6,6,6,6/); jjj(3,2,1:8)=(/5,6,7,8,1,2,3,4/) ! Loop over myid grids by vector sweeps ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ! Gather nvector grids ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Gather neighboring grids do i=1,ngrid igridn(i,0)=ind_grid(i) end do do idim=1,ndim do i=1,ngrid ind_left (i,idim)=nbor(ind_grid(i),2*idim-1) ind_right(i,idim)=nbor(ind_grid(i),2*idim ) igridn(i,2*idim-1)=son(ind_left (i,idim)) igridn(i,2*idim )=son(ind_right(i,idim)) end do end do ! Loop over cells do ind=1,twotondim ! Compute central cell index iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do ! Gather all neighboring velocities do idim=1,ndim id1=jjj(idim,1,ind); ig1=iii(idim,1,ind) ih1=ncoarse+(id1-1)*ngridmax do i=1,ngrid if(igridn(i,ig1)>0)then velg(i,idim,1:ndim) = uold(igridn(i,ig1)+ih1,2:ndim+1)/max(uold(igridn(i,ig1)+ih1,1),smallr) dx_g(i,idim) = dx_loc else velg(i,idim,1:ndim) = uold(ind_left(i,idim),2:ndim+1)/max(uold(ind_left(i,idim),1),smallr) dx_g(i,idim) = dx_loc*1.5_dp end if enddo id2=jjj(idim,2,ind); ig2=iii(idim,2,ind) ih2=ncoarse+(id2-1)*ngridmax do i=1,ngrid if(igridn(i,ig2)>0)then veld(i,idim,1:ndim)= uold(igridn(i,ig2)+ih2,2:ndim+1)/max(uold(igridn(i,ig2)+ih2,1),smallr) dx_d(i,idim)=dx_loc else veld(i,idim,1:ndim)= uold(ind_right(i,idim),2:ndim+1)/max(uold(ind_right(i,idim),1),smallr) dx_d(i,idim)=dx_loc*1.5_dp end if enddo end do ! End loop over dimensions ! Compute divu = Trace G divu_loc(1:ngrid)=0.0d0 do i=1,ngrid do idim=1,ndim divu_loc(i) = divu_loc(i) + (veld(i,idim,idim)-velg(i,idim,idim)) & & / (dx_g(i,idim) +dx_d(i,idim)) enddo end do ! Update thermal internal energy if(pressure_fix)then do i=1,ngrid ! Compute old thermal energy d=max(uold(ind_cell(i),1),smallr) u=0; v=0; w=0 if(ndim>0)u=uold(ind_cell(i),2)/d if(ndim>1)v=uold(ind_cell(i),3)/d if(ndim>2)w=uold(ind_cell(i),4)/d eold=uold(ind_cell(i),ndim+2)-0.5d0*d*(u**2+v**2+w**2) #if NENER>0 do irad=1,nener eold=eold-uold(ind_cell(i),ndim+2+irad) end do #endif ! Add -pdV term enew(ind_cell(i))=enew(ind_cell(i)) & & -(gamma-1.0d0)*eold*divu_loc(i)*dtnew(ilevel) end do end if #if NENER>0 do irad=1,nener do i=1,ngrid ! Add -pdV term unew(ind_cell(i),ndim+2+irad)=unew(ind_cell(i),ndim+2+irad) & & -(gamma_rad(irad)-1.0d0)*uold(ind_cell(i),ndim+2+irad)*divu_loc(i)*dtnew(ilevel) end do end do #endif if(momentum_feedback>0)then ! Add +pdV term do i=1,ngrid unew(ind_cell(i),ndim+2)=unew(ind_cell(i),ndim+2) & & +pstarold(ind_cell(i))*divu_loc(i)*dx_loc/6.0 end do endif enddo ! End loop over cells end do ! End loop over grids return ! This is the old technique based on the 'pressure fix' option. ! Update thermal internal energy if(pressure_fix)then do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,active(ilevel)%ngrid ind_cell1=active(ilevel)%igrid(i)+iskip ! Compute old thermal energy d=max(uold(ind_cell1,1),smallr) u=0; v=0; w=0 if(ndim>0)u=uold(ind_cell1,2)/d if(ndim>1)v=uold(ind_cell1,3)/d if(ndim>2)w=uold(ind_cell1,4)/d eold=uold(ind_cell1,ndim+2)-0.5d0*d*(u**2+v**2+w**2) #if NENER>0 do irad=1,nener eold=eold-uold(ind_cell1,ndim+2+irad) end do #endif ! Add pdV term enew(ind_cell1)=enew(ind_cell1) & & +(gamma-1.0d0)*eold*divu(ind_cell1) ! Note: here divu=-div.u*dt end do end do end if #if NENER>0 do irad=1,nener do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,active(ilevel)%ngrid ind_cell1=active(ilevel)%igrid(i)+iskip unew(ind_cell1,ndim+2+irad)=unew(ind_cell1,ndim+2+irad) & & +(gamma_rad(irad)-1.0d0)*uold(ind_cell1,ndim+2+irad)*divu(ind_cell1) ! Note: here divu=-div.u*dt end do end do end do #endif 111 format(' Entering add_pdv_source_terms for level ',i2) end subroutine add_pdv_source_terms !########################################################### !########################################################### !########################################################### !########################################################### subroutine add_viscosity_source_terms(ilevel) use amr_commons use hydro_commons use poisson_commons use pm_commons implicit none integer::ilevel,levelmax !-------------------------------------------------------------------------- ! This routine adds to unew the viscosity terms ! with only half a time step. Only the momentum and the ! total energy are modified in array unew. !-------------------------------------------------------------------------- integer::i,ind,iskip,nx_loc integer::ncache,igrid,ngrid,idim,id1,ig1,ih1,id2,ig2,ih2,jdim integer,dimension(1:3,1:2,1:8)::iii,jjj real(dp)::scale,dx,dx_loc,dx_min real(dp)::Kturb,sigma,d_old real(dp)::scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2 integer ,dimension(1:nvector),save::ind_grid,ind_cell integer ,dimension(1:nvector,0:twondim),save::igridn integer ,dimension(1:nvector,1:ndim),save::ind_left,ind_right real(dp),dimension(1:nvector,1:ndim,1:ndim),save::velg,veld real(dp),dimension(1:nvector,1:ndim),save::dx_g,dx_d real(dp),dimension(1:nvector),save::divu_loc,phi_diss real(dp),dimension(1:nvector,1:ndim,1:ndim),save::gradu_loc,E_loc #if NENER>0 integer::irad #endif if(numbtot(1,ilevel)==0)return if(verbose)write(*,111)ilevel nx_loc=icoarse_max-icoarse_min+1 scale=boxlen/dble(nx_loc) dx=0.5d0**ilevel dx_loc=dx*scale dx_min=(0.5**levelmax)*scale call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) iii(1,1,1:8)=(/1,0,1,0,1,0,1,0/); jjj(1,1,1:8)=(/2,1,4,3,6,5,8,7/) iii(1,2,1:8)=(/0,2,0,2,0,2,0,2/); jjj(1,2,1:8)=(/2,1,4,3,6,5,8,7/) iii(2,1,1:8)=(/3,3,0,0,3,3,0,0/); jjj(2,1,1:8)=(/3,4,1,2,7,8,5,6/) iii(2,2,1:8)=(/0,0,4,4,0,0,4,4/); jjj(2,2,1:8)=(/3,4,1,2,7,8,5,6/) iii(3,1,1:8)=(/5,5,5,5,0,0,0,0/); jjj(3,1,1:8)=(/5,6,7,8,1,2,3,4/) iii(3,2,1:8)=(/0,0,0,0,6,6,6,6/); jjj(3,2,1:8)=(/5,6,7,8,1,2,3,4/) ! Loop over myid grids by vector sweeps ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ! Gather nvector grids ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Gather neighboring grids do i=1,ngrid igridn(i,0)=ind_grid(i) end do do idim=1,ndim do i=1,ngrid ind_left (i,idim)=nbor(ind_grid(i),2*idim-1) ind_right(i,idim)=nbor(ind_grid(i),2*idim ) igridn(i,2*idim-1)=son(ind_left (i,idim)) igridn(i,2*idim )=son(ind_right(i,idim)) end do end do ! Loop over cells do ind=1,twotondim ! Compute central cell index iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do ! Gather all neighboring velocities do idim=1,ndim id1=jjj(idim,1,ind); ig1=iii(idim,1,ind) ih1=ncoarse+(id1-1)*ngridmax do i=1,ngrid if(igridn(i,ig1)>0)then velg(i,idim,1:ndim) = uold(igridn(i,ig1)+ih1,2:ndim+1)/max(uold(igridn(i,ig1)+ih1,1),smallr) dx_g(i,idim) = dx_loc else velg(i,idim,1:ndim) = uold(ind_left(i,idim),2:ndim+1)/max(uold(ind_left(i,idim),1),smallr) dx_g(i,idim) = dx_loc*1.5_dp end if enddo id2=jjj(idim,2,ind); ig2=iii(idim,2,ind) ih2=ncoarse+(id2-1)*ngridmax do i=1,ngrid if(igridn(i,ig2)>0)then veld(i,idim,1:ndim)= uold(igridn(i,ig2)+ih2,2:ndim+1)/max(uold(igridn(i,ig2)+ih2,1),smallr) dx_d(i,idim)=dx_loc else veld(i,idim,1:ndim)= uold(ind_right(i,idim),2:ndim+1)/max(uold(ind_right(i,idim),1),smallr) dx_d(i,idim)=dx_loc*1.5_dp end if enddo end do ! End loop over dimensions ! Compute divu = Trace G divu_loc(1:ngrid)=0.0d0 do idim=1,ndim do i=1,ngrid divu_loc(i) = divu_loc(i) + (veld(i,idim,idim)-velg(i,idim,idim)) & & / (dx_g(i,idim) +dx_d(i,idim)) enddo end do ! Compute gradu = G gradu_loc(1:ngrid,1:ndim,1:ndim)=0.0d0 do idim=1,ndim do jdim=1,ndim do i=1,ngrid gradu_loc(i,idim,jdim) = (veld(i,idim,jdim)-velg(i,idim,jdim)) & & / (dx_g(i,idim) +dx_d(i,idim)) enddo enddo end do ! Compute 0.5*(gradu+gradu^T) = E E_loc(1:ngrid,1:ndim,1:ndim)=0.0d0 do idim=1,ndim do jdim=1,ndim do i=1,ngrid E_loc(i,idim,jdim) = 0.5*(gradu_loc(i,idim,jdim)+gradu_loc(i,jdim,idim)) enddo enddo end do ! Compute turbulent viscosity dissipation function phi_diss(1:ngrid)=0 do idim=1,ndim do jdim=1,ndim do i=1,ngrid phi_diss(i) = phi_diss(i)+2.0*E_loc(i,idim,jdim)**2 enddo enddo enddo do i=1,ngrid phi_diss(i)=phi_diss(i)-2.0/3.0*divu_loc(i)**2 end do ! Add viscosity term at time t with half time step do i=1,ngrid ! Compute old turbulent state d_old=max(uold(ind_cell(i),1),smallr) Kturb=uold(ind_cell(i),ivirial1) sigma=sqrt(max(2.0*Kturb/d_old,smallc**2)) ! Implicit solution unew(ind_cell(i),ivirial1)=(unew(ind_cell(i),ivirial1) & & +d_old*dx_loc*sigma*phi_diss(i)*dtnew(ilevel)) & & /(1.0+sigma/dx_loc*dtnew(ilevel)) ! Turbulence from SN ! unew(ind_cell(i),ivirial1) = unew(ind_cell(i),ivirial1)/(1.0+sigma/dx_loc*dtnew(ilevel)) ! Stationary solution ! unew(ind_cell(i),ivirial1)=d_old*dx_loc**2*phi_diss(i) end do ! End loop over grids end do ! End loop over cells end do ! End loop over sweeps 111 format(' Entering add_viscosity_terms for level ',i2) end subroutine add_viscosity_source_terms !########################################################### !########################################################### !########################################################### !########################################################### subroutine godfine1(ind_grid,ncache,ilevel) use amr_commons use hydro_commons use poisson_commons implicit none integer::ilevel,ncache integer,dimension(1:nvector)::ind_grid !------------------------------------------------------------------- ! This routine gathers first hydro variables from neighboring grids ! to set initial conditions in a 6x6x6 grid. It interpolate from ! coarser level missing grid variables. It then calls the ! Godunov solver that computes fluxes. These fluxes are zeroed at ! coarse-fine boundaries, since contribution from finer levels has ! already been taken into account. Conservative variables are updated ! and stored in array unew(:), both at the current level and at the ! coarser level if necessary. !------------------------------------------------------------------- integer ,dimension(1:nvector,1:threetondim ),save::nbors_father_cells integer ,dimension(1:nvector,1:twotondim ),save::nbors_father_grids integer ,dimension(1:nvector,0:twondim ),save::ibuffer_father real(dp),dimension(1:nvector,0:twondim ,1:nvar),save::u1 real(dp),dimension(1:nvector,1:twotondim,1:nvar),save::u2 real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar),save::uloc real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:ndim),save::gloc=0.0d0 real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2),save::ploc=0.0d0 real(dp),dimension(1:nvector,if1:if2,jf1:jf2,kf1:kf2,1:nvar,1:ndim),save::flux real(dp),dimension(1:nvector,if1:if2,jf1:jf2,kf1:kf2,1:2,1:ndim),save::tmp logical ,dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2),save::ok integer,dimension(1:nvector),save::igrid_nbor,ind_cell,ind_buffer,ind_exist,ind_nexist integer::i,j,ivar,idim,ind_son,ind_father,iskip,nbuffer integer::i0,j0,k0,i1,j1,k1,i2,j2,k2,i3,j3,k3,nx_loc,nb_noneigh,nexist integer::i1min,i1max,j1min,j1max,k1min,k1max integer::i2min,i2max,j2min,j2max,k2min,k2max integer::i3min,i3max,j3min,j3max,k3min,k3max real(dp)::dx,scale,oneontwotondim oneontwotondim = 1d0/dble(twotondim) ! Mesh spacing in that level nx_loc=icoarse_max-icoarse_min+1 scale=boxlen/dble(nx_loc) dx=0.5d0**ilevel*scale ! Integer constants i1min=0; i1max=0; i2min=0; i2max=0; i3min=1; i3max=1 j1min=0; j1max=0; j2min=0; j2max=0; j3min=1; j3max=1 k1min=0; k1max=0; k2min=0; k2max=0; k3min=1; k3max=1 if(ndim>0)then i1max=2; i2max=1; i3max=2 end if if(ndim>1)then j1max=2; j2max=1; j3max=2 end if if(ndim>2)then k1max=2; k2max=1; k3max=2 end if !------------------------------------------ ! Gather 3^ndim neighboring father cells !------------------------------------------ do i=1,ncache ind_cell(i)=father(ind_grid(i)) end do call get3cubefather(ind_cell,nbors_father_cells,nbors_father_grids,ncache,ilevel) !--------------------------- ! Gather 6x6x6 cells stencil !--------------------------- ! Loop over 3x3x3 neighboring father cells do k1=k1min,k1max do j1=j1min,j1max do i1=i1min,i1max ! Check if neighboring grid exists nbuffer=0 nexist=0 ind_father=1+i1+3*j1+9*k1 do i=1,ncache igrid_nbor(i)=son(nbors_father_cells(i,ind_father)) if(igrid_nbor(i)>0) then nexist=nexist+1 ind_exist(nexist)=i else nbuffer=nbuffer+1 ind_nexist(nbuffer)=i ind_buffer(nbuffer)=nbors_father_cells(i,ind_father) end if end do ! If not, interpolate hydro variables from parent cells if(nbuffer>0)then call getnborfather(ind_buffer,ibuffer_father,nbuffer,ilevel) do j=0,twondim do ivar=1,nvar do i=1,nbuffer u1(i,j,ivar)=uold(ibuffer_father(i,j),ivar) end do end do end do call interpol_hydro(u1,u2,nbuffer) endif ! Loop over 2x2x2 cells do k2=k2min,k2max do j2=j2min,j2max do i2=i2min,i2max ind_son=1+i2+2*j2+4*k2 iskip=ncoarse+(ind_son-1)*ngridmax do i=1,nexist ind_cell(i)=iskip+igrid_nbor(ind_exist(i)) end do i3=1; j3=1; k3=1 if(ndim>0)i3=1+2*(i1-1)+i2 if(ndim>1)j3=1+2*(j1-1)+j2 if(ndim>2)k3=1+2*(k1-1)+k2 ! Gather hydro variables do ivar=1,nvar do i=1,nexist uloc(ind_exist(i),i3,j3,k3,ivar)=uold(ind_cell(i),ivar) end do do i=1,nbuffer uloc(ind_nexist(i),i3,j3,k3,ivar)=u2(i,ind_son,ivar) end do end do ! Gather gravitational acceleration if(poisson)then do idim=1,ndim do i=1,nexist gloc(ind_exist(i),i3,j3,k3,idim)=f(ind_cell(i),idim) end do ! Use straight injection for buffer cells do i=1,nbuffer gloc(ind_nexist(i),i3,j3,k3,idim)=f(ibuffer_father(i,0),idim) end do end do end if ! Gather stellar momentum if(momentum_feedback>0)then do i=1,nexist ploc(ind_exist(i),i3,j3,k3)=pstarold(ind_cell(i))*dx/dtnew(ilevel)/6.0 end do ! Use straight injection for buffer cells do i=1,nbuffer ploc(ind_nexist(i),i3,j3,k3)=pstarold(ibuffer_father(i,0))*dx/dtnew(ilevel)/6.0 end do end if ! Gather refinement flag do i=1,nexist ok(ind_exist(i),i3,j3,k3)=son(ind_cell(i))>0 end do do i=1,nbuffer ok(ind_nexist(i),i3,j3,k3)=.false. end do end do end do end do ! End loop over cells end do end do end do ! End loop over neighboring grids !----------------------------------------------- ! Compute flux using second-order Godunov method !----------------------------------------------- call unsplit(uloc,ploc,gloc,flux,tmp,dx,dx,dx,dtnew(ilevel),ncache) !------------------------------------------------ ! Reset flux along direction at refined interface !------------------------------------------------ do idim=1,ndim i0=0; j0=0; k0=0 if(idim==1)i0=1 if(idim==2)j0=1 if(idim==3)k0=1 do k3=k3min,k3max+k0 do j3=j3min,j3max+j0 do i3=i3min,i3max+i0 do ivar=1,nvar do i=1,ncache if(ok(i,i3-i0,j3-j0,k3-k0) .or. ok(i,i3,j3,k3))then flux(i,i3,j3,k3,ivar,idim)=0.0d0 end if end do end do if(pressure_fix)then do ivar=1,2 do i=1,ncache if(ok(i,i3-i0,j3-j0,k3-k0) .or. ok(i,i3,j3,k3))then tmp (i,i3,j3,k3,ivar,idim)=0.0d0 end if end do end do end if end do end do end do end do !-------------------------------------- ! Conservative update at level ilevel !-------------------------------------- do idim=1,ndim i0=0; j0=0; k0=0 if(idim==1)i0=1 if(idim==2)j0=1 if(idim==3)k0=1 do k2=k2min,k2max do j2=j2min,j2max do i2=i2min,i2max ind_son=1+i2+2*j2+4*k2 iskip=ncoarse+(ind_son-1)*ngridmax do i=1,ncache ind_cell(i)=iskip+ind_grid(i) end do i3=1+i2 j3=1+j2 k3=1+k2 ! Update conservative variables new state vector do ivar=1,nvar do i=1,ncache unew(ind_cell(i),ivar)=unew(ind_cell(i),ivar)+ & & (flux(i,i3 ,j3 ,k3 ,ivar,idim) & & -flux(i,i3+i0,j3+j0,k3+k0,ivar,idim)) end do end do if(pressure_fix)then ! Update velocity divergence do i=1,ncache divu(ind_cell(i))=divu(ind_cell(i))+ & & (tmp(i,i3 ,j3 ,k3 ,1,idim) & & -tmp(i,i3+i0,j3+j0,k3+k0,1,idim)) end do ! Update internal energy do i=1,ncache enew(ind_cell(i))=enew(ind_cell(i))+ & & (tmp(i,i3 ,j3 ,k3 ,2,idim) & & -tmp(i,i3+i0,j3+j0,k3+k0,2,idim)) end do end if end do end do end do end do !-------------------------------------- ! Conservative update at level ilevel-1 !-------------------------------------- ! Loop over dimensions do idim=1,ndim i0=0; j0=0; k0=0 if(idim==1)i0=1 if(idim==2)j0=1 if(idim==3)k0=1 !---------------------- ! Left flux at boundary !---------------------- ! Check if grids sits near left boundary ! and gather neighbor father cells index nb_noneigh=0 do i=1,ncache if (son(nbor(ind_grid(i),2*idim-1))==0) then nb_noneigh = nb_noneigh + 1 ind_buffer(nb_noneigh) = nbor(ind_grid(i),2*idim-1) ind_cell(nb_noneigh) = i end if end do ! Conservative update of new state variables do ivar=1,nvar ! Loop over boundary cells do k3=k3min,k3max-k0 do j3=j3min,j3max-j0 do i3=i3min,i3max-i0 do i=1,nb_noneigh unew(ind_buffer(i),ivar)=unew(ind_buffer(i),ivar) & & -flux(ind_cell(i),i3,j3,k3,ivar,idim)*oneontwotondim end do end do end do end do end do if(pressure_fix)then ! Update velocity divergence do k3=k3min,k3max-k0 do j3=j3min,j3max-j0 do i3=i3min,i3max-i0 do i=1,nb_noneigh divu(ind_buffer(i))=divu(ind_buffer(i)) & & -tmp(ind_cell(i),i3,j3,k3,1,idim)*oneontwotondim end do end do end do end do ! Update internal energy do k3=k3min,k3max-k0 do j3=j3min,j3max-j0 do i3=i3min,i3max-i0 do i=1,nb_noneigh enew(ind_buffer(i))=enew(ind_buffer(i)) & & -tmp(ind_cell(i),i3,j3,k3,2,idim)*oneontwotondim end do end do end do end do end if !----------------------- ! Right flux at boundary !----------------------- ! Check if grids sits near right boundary ! and gather neighbor father cells index nb_noneigh=0 do i=1,ncache if (son(nbor(ind_grid(i),2*idim))==0) then nb_noneigh = nb_noneigh + 1 ind_buffer(nb_noneigh) = nbor(ind_grid(i),2*idim) ind_cell(nb_noneigh) = i end if end do ! Conservative update of new state variables do ivar=1,nvar ! Loop over boundary cells do k3=k3min+k0,k3max do j3=j3min+j0,j3max do i3=i3min+i0,i3max do i=1,nb_noneigh unew(ind_buffer(i),ivar)=unew(ind_buffer(i),ivar) & & +flux(ind_cell(i),i3+i0,j3+j0,k3+k0,ivar,idim)*oneontwotondim end do end do end do end do end do if(pressure_fix)then ! Update velocity divergence do k3=k3min+k0,k3max do j3=j3min+j0,j3max do i3=i3min+i0,i3max do i=1,nb_noneigh divu(ind_buffer(i))=divu(ind_buffer(i)) & & +tmp(ind_cell(i),i3+i0,j3+j0,k3+k0,1,idim)*oneontwotondim end do end do end do end do ! Update internal energy do k3=k3min+k0,k3max do j3=j3min+j0,j3max do i3=i3min+i0,i3max do i=1,nb_noneigh enew(ind_buffer(i))=enew(ind_buffer(i)) & & +tmp(ind_cell(i),i3+i0,j3+j0,k3+k0,2,idim)*oneontwotondim end do end do end do end do end if end do ! End loop over dimensions end subroutine godfine1 /home/hemv_g/srf_a/RAMSES_CODE_CORI/ramses06062019_grackle311Cosmo/patch/mom2/godunov_utils.f90 !########################################################### !########################################################### !########################################################### !########################################################### subroutine cmpdt(uu,gg,pp,dx,dt,ncell) use amr_parameters use hydro_parameters use const implicit none integer::ncell real(dp)::dx,dt real(dp),dimension(1:nvector,1:nvar)::uu real(dp),dimension(1:nvector,1:ndim)::gg real(dp),dimension(1:nvector)::pp real(dp)::dtcell,smallp integer::k,idim #if NENER>0 integer::irad #endif smallp = smallc**2/gamma ! Convert to primitive variables do k = 1,ncell uu(k,1)=max(uu(k,1),smallr) end do ! Velocity do idim = 1,ndim do k = 1, ncell uu(k,idim+1) = uu(k,idim+1)/uu(k,1) end do end do ! Internal energy do idim = 1,ndim do k = 1, ncell uu(k,ndim+2) = uu(k,ndim+2)-half*uu(k,1)*uu(k,idim+1)**2 end do end do #if NENER>0 do irad = 1,nener do k = 1, ncell uu(k,ndim+2) = uu(k,ndim+2)-uu(k,ndim+2+irad) end do end do #endif ! Debug if(debug)then do k = 1, ncell if(uu(k,ndim+2).le.0.or.uu(k,1).le.smallr)then write(*,*)'stop in cmpdt' write(*,*)'dx =',dx write(*,*)'ncell=',ncell write(*,*)'rho =',uu(k,1) write(*,*)'P =',uu(k,ndim+2) write(*,*)'vel =',uu(k,2:ndim+1) stop end if end do end if ! Compute pressure do k = 1, ncell uu(k,ndim+2) = max((gamma-one)*uu(k,ndim+2),uu(k,1)*smallp) end do #if NENER>0 do irad = 1,nener do k = 1, ncell uu(k,ndim+2+irad) = (gamma_rad(irad)-one)*uu(k,ndim+2+irad) end do end do #endif ! Compute sound speed do k = 1, ncell uu(k,ndim+2) = gamma*uu(k,ndim+2) end do #if NENER>0 do irad = 1,nener do k = 1, ncell uu(k,ndim+2) = uu(k,ndim+2) + gamma_rad(irad)*uu(k,ndim+2+irad) end do end do #endif do k = 1, ncell uu(k,ndim+2)=sqrt(uu(k,ndim+2)/uu(k,1)) end do ! Compute wave speed do k = 1, ncell uu(k,ndim+2)=dble(ndim)*uu(k,ndim+2) end do do idim = 1,ndim do k = 1, ncell uu(k,ndim+2)=uu(k,ndim+2)+abs(uu(k,idim+1)) end do end do if(momentum_feedback>0)then do k = 1, ncell uu(k,ndim+2) = uu(k,ndim+2) + abs(pp(k)/uu(k,1)) end do endif ! Compute gravity strength ratio do k = 1, ncell uu(k,1)=zero end do do idim = 1,ndim do k = 1, ncell uu(k,1)=uu(k,1)+abs(gg(k,idim)) end do end do do k = 1, ncell uu(k,1)=uu(k,1)*dx/uu(k,ndim+2)**2 uu(k,1)=MAX(uu(k,1),0.0001_dp) end do ! Compute maximum time step for each authorized cell dt = courant_factor*dx/smallc do k = 1,ncell dtcell = dx/uu(k,ndim+2)*(sqrt(one+two*courant_factor*uu(k,1))-one)/uu(k,1) dt = min(dt,dtcell) end do end subroutine cmpdt !########################################################### !########################################################### !########################################################### !########################################################### subroutine hydro_refine(ug,um,ud,ok,nn) use amr_parameters use hydro_parameters use const #ifdef RT use rt_parameters #endif implicit none ! dummy arguments integer nn real(dp)::ug(1:nvector,1:nvar) real(dp)::um(1:nvector,1:nvar) real(dp)::ud(1:nvector,1:nvar) logical ::ok(1:nvector) integer::k,idim real(dp),dimension(1:nvector),save::eking,ekinm,ekind real(dp)::dg,dm,dd,pg,pm,pd,vg,vm,vd,cg,cm,cd,error #if NENER>0 integer::irad #endif ! Convert to primitive variables do k = 1,nn ug(k,1) = max(ug(k,1),smallr) um(k,1) = max(um(k,1),smallr) ud(k,1) = max(ud(k,1),smallr) end do ! Velocity do idim = 1,ndim do k = 1,nn ug(k,idim+1) = ug(k,idim+1)/ug(k,1) um(k,idim+1) = um(k,idim+1)/um(k,1) ud(k,idim+1) = ud(k,idim+1)/ud(k,1) end do end do ! Pressure do k = 1,nn eking(k) = zero ekinm(k) = zero ekind(k) = zero end do do idim = 1,ndim do k = 1,nn eking(k) = eking(k) + half*ug(k,1)*ug(k,idim+1)**2 ekinm(k) = ekinm(k) + half*um(k,1)*um(k,idim+1)**2 ekind(k) = ekind(k) + half*ud(k,1)*ud(k,idim+1)**2 end do end do #if NENER>0 do irad = 1,nener do k = 1, nn eking(k) = eking(k) + ug(k,ndim+2+irad) ekinm(k) = ekinm(k) + um(k,ndim+2+irad) ekind(k) = ekind(k) + ud(k,ndim+2+irad) end do end do #endif do k = 1,nn ug(k,ndim+2) = (gamma-one)*(ug(k,ndim+2)-eking(k)) um(k,ndim+2) = (gamma-one)*(um(k,ndim+2)-ekinm(k)) ud(k,ndim+2) = (gamma-one)*(ud(k,ndim+2)-ekind(k)) end do ! Passive scalars #if NVAR > NDIM + 2 do idim = ndim+3,nvar do k = 1,nn ug(k,idim) = ug(k,idim)/ug(k,1) um(k,idim) = um(k,idim)/um(k,1) ud(k,idim) = ud(k,idim)/ud(k,1) end do end do #endif ! Compute errors if(err_grad_d >= 0.)then do k=1,nn dg=ug(k,1); dm=um(k,1); dd=ud(k,1) error=2.0d0*MAX( & & ABS((dd-dm)/(dd+dm+floor_d)) , & & ABS((dm-dg)/(dm+dg+floor_d)) ) ok(k) = ok(k) .or. error > err_grad_d end do end if if(err_grad_p >= 0.)then do k=1,nn pg=ug(k,ndim+2); pm=um(k,ndim+2); pd=ud(k,ndim+2) error=2.0d0*MAX( & & ABS((pd-pm)/(pd+pm+floor_p)), & & ABS((pm-pg)/(pm+pg+floor_p)) ) ok(k) = ok(k) .or. error > err_grad_p end do end if if(err_grad_u >= 0.)then do idim = 1,ndim do k=1,nn vg=ug(k,idim+1); vm=um(k,idim+1); vd=ud(k,idim+1) cg=sqrt(max(gamma*ug(k,ndim+2)/ug(k,1),floor_u**2)) cm=sqrt(max(gamma*um(k,ndim+2)/um(k,1),floor_u**2)) cd=sqrt(max(gamma*ud(k,ndim+2)/ud(k,1),floor_u**2)) error=2.0d0*MAX( & & ABS((vd-vm)/(cd+cm+ABS(vd)+ABS(vm)+floor_u)) , & & ABS((vm-vg)/(cm+cg+ABS(vm)+ABS(vg)+floor_u)) ) ok(k) = ok(k) .or. error > err_grad_u end do end do end if #ifdef RT ! Ionization state (only Hydrogen) if(rt_err_grad_xHII >= 0.) then !--------------------------------------- do k=1,nn dg=min(1d0,max(0d0,ug(k,iIons))) dm=min(1d0,max(0d0,um(k,iIons))) dd=min(1d0,max(0d0,ud(k,iIons))) error=2.0d0*MAX( & & ABS((dd-dm)/(dd+dm+rt_floor_xHII)) , & & ABS((dm-dg)/(dm+dg+rt_floor_xHII)) ) ok(k) = ok(k) .or. error > rt_err_grad_xHII end do end if ! Neutral state (only Hydrogen) if(rt_err_grad_xHI >= 0.) then !--------------------------------------- do k=1,nn dg=min(1d0,max(0d0,1d0 - ug(k,iIons))) dm=min(1d0,max(0d0,1d0 - um(k,iIons))) dd=min(1d0,max(0d0,1d0 - ud(k,iIons))) error=2.0d0*MAX( & & ABS((dd-dm)/(dd+dm+rt_floor_xHI)) , & & ABS((dm-dg)/(dm+dg+rt_floor_xHI)) ) ok(k) = ok(k) .or. error > rt_err_grad_xHI end do end if #endif end subroutine hydro_refine !########################################################### !########################################################### !########################################################### !########################################################### subroutine riemann_approx(qleft,qright,fgdnv,ngrid) use amr_parameters use hydro_parameters use const implicit none ! dummy arguments integer::ngrid real(dp),dimension(1:nvector,1:nvar)::qleft,qright real(dp),dimension(1:nvector,1:nvar+1)::fgdnv ! local arrays real(dp),dimension(1:nvector,1:nvar+1),save::qgdnv real(dp),dimension(1:nvector),save::rl ,ul ,pl ,cl real(dp),dimension(1:nvector),save::rr ,ur ,pr ,cr real(dp),dimension(1:nvector),save::ro ,uo ,po ,co real(dp),dimension(1:nvector),save::rstar,ustar,pstar,cstar real(dp),dimension(1:nvector),save::wl ,wr ,wo real(dp),dimension(1:nvector),save::sgnm ,spin ,spout,ushock real(dp),dimension(1:nvector),save::frac ,delp ,pold integer ,dimension(1:nvector),save::ind ,ind2 ! local variables real(dp)::smallp, gamma6, ql, qr, usr, usl, wwl, wwr, smallpp, entho, etot integer ::i, j, n, iter, n_new ! constants smallp = smallc**2/gamma smallpp = smallr*smallp gamma6 = (gamma+one)/(two*gamma) entho = one/(gamma-one) ! Pressure, density and velocity do i=1,ngrid rl(i)=MAX(qleft (i,1),smallr) ul(i)= qleft (i,2) pl(i)=MAX(qleft (i,3),rl(i)*smallp) rr(i)=MAX(qright(i,1),smallr) ur(i)= qright(i,2) pr(i)=MAX(qright(i,3),rr(i)*smallp) end do ! Lagrangian sound speed do i=1,ngrid cl(i) = gamma*pl(i)*rl(i) cr(i) = gamma*pr(i)*rr(i) end do ! First guess do i=1,ngrid wl(i)= sqrt(cl(i)); wr(i) = sqrt(cr(i)) pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i)))/(wl(i)+wr(i)) pstar(i) = MAX(pstar(i),0.0_dp) pold (i)= pstar(i) end do n = ngrid do i = 1,n ind(i)=i end do ! Newton-Raphson iterations to find pstar at the required accuracy ! for a two-shock Riemann problem do iter = 1,niter_riemann do i=1,n wwl=sqrt(cl(ind(i))*(one+gamma6*(pold(i)-pl(ind(i)))/pl(ind(i)))) wwr=sqrt(cr(ind(i))*(one+gamma6*(pold(i)-pr(ind(i)))/pr(ind(i)))) ql=two*wwl**3/(wwl**2+cl(ind(i))) qr=two*wwr**3/(wwr**2+cr(ind(i))) usl=ul(ind(i))-(pold(i)-pl(ind(i)))/wwl usr=ur(ind(i))+(pold(i)-pr(ind(i)))/wwr delp(i)=MAX(qr*ql/(qr+ql)*(usl-usr),-pold(i)) end do do i=1,n pold(i)=pold(i)+delp(i) end do ! Convergence indicator do i=1,n uo(i)=ABS(delp(i)/(pold(i)+smallpp)) end do n_new=0 do i=1,n if(uo(i)>1d-06)then n_new=n_new+1 ind2(n_new)=ind (i) po (n_new)=pold(i) end if end do j=n_new do i=1,n if(uo(i)<=1d-06)then n_new=n_new+1 ind2(n_new)=ind (i) po (n_new)=pold(i) end if end do ind (1:n)=ind2(1:n) pold(1:n)=po (1:n) n=j end do ! Star region pressure ! for a two-shock Riemann problem do i=1,ngrid pstar(ind(i))=pold(i) end do do i=1,ngrid wl(i) = sqrt(cl(i)*(one+gamma6*(pstar(i)-pl(i))/pl(i))) wr(i) = sqrt(cr(i)*(one+gamma6*(pstar(i)-pr(i))/pr(i))) end do ! Star region velocity ! for a two shock Riemann problem do i=1,ngrid ustar(i) = half*(ul(i) + (pl(i)-pstar(i))/wl(i) + & & ur(i) - (pr(i)-pstar(i))/wr(i) ) end do ! Left going or right going contact wave do i=1,ngrid sgnm(i) = sign(one,ustar(i)) end do ! Left or right unperturbed state do i=1,ngrid if(sgnm(i)==one)then ro(i) = rl(i) uo(i) = ul(i) po(i) = pl(i) wo(i) = wl(i) else ro(i) = rr(i) uo(i) = ur(i) po(i) = pr(i) wo(i) = wr(i) end if end do do i=1,ngrid co(i) = max(smallc,sqrt(abs(gamma*po(i)/ro(i)))) end do ! Star region density do i=1,ngrid if(pstar(i)>= po(i))then ! Shock rstar(i) = ro(i)/(one+ro(i)*(po(i)-pstar(i))/wo(i)**2) else ! Rarefaction rstar(i) = ro(i)*(pstar(i)/po(i))**(one/gamma) end if end do do i=1,ngrid ! Prevent vacuum formation in star region rstar(i) = max(rstar(i),smallr) ! Star region sound speed cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i))) cstar(i) = max(cstar(i),smallc) ! Compute rarefaction head and tail speed spout(i) = co (i)-sgnm(i)*uo (i) spin (i) = cstar(i)-sgnm(i)*ustar(i) ! Compute shock speed ushock(i) = wo(i)/ro(i)-sgnm(i)*uo(i) end do do i=1,ngrid if(pstar(i)>=po(i))then spout(i)=ushock(i) spin (i)=spout (i) end if end do ! Sample the solution at x/t=0 do i=1,ngrid if(spout(i)<=zero)then qgdnv(i,1) = ro(i) qgdnv(i,2) = uo(i) qgdnv(i,3) = po(i) else if(spin(i)>=zero)then qgdnv(i,1) = rstar(i) qgdnv(i,2) = ustar(i) qgdnv(i,3) = pstar(i) else frac(i)=spout(i)/(spout(i)-spin(i)) qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i) qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i) qgdnv(i,1) = ro(i)*(qgdnv(i,3)/po(i))**(one/gamma) end if end do ! Passive scalars #if NVAR > 3 do n = 4,nvar do i=1,ngrid if(sgnm(i)==one)then qgdnv(i,n) = qleft(i,n) else qgdnv(i,n) = qright(i,n) end if end do end do #endif ! Specific internal energy do i=1,ngrid qgdnv(i,nvar+1) = po(i)/ro(i)*entho end do ! Compute fluxes do i = 1, ngrid fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2) ! Mass density fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2 ! Normal momentum etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2 #if NDIM>1 etot = etot + half*qgdnv(i,1)*qgdnv(i,4)**2 #endif #if NDIM>2 etot = etot + half*qgdnv(i,1)*qgdnv(i,5)**2 #endif fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3)) ! Total energy end do ! Other advected quantities do n = 4, nvar+1 do i = 1, ngrid fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n) end do end do end subroutine riemann_approx !########################################################### !########################################################### !########################################################### !########################################################### subroutine riemann_acoustic(qleft,qright,fgdnv,ngrid) use amr_parameters use hydro_parameters use const implicit none ! dummy arguments integer::ngrid real(dp),dimension(1:nvector,1:nvar)::qleft,qright real(dp),dimension(1:nvector,1:nvar+1)::fgdnv ! local variables integer::i,n real(dp)::smallp, entho, etot ! local arrays real(dp),dimension(1:nvector,1:nvar+1),save::qgdnv real(dp),dimension(1:nvector),save::rl ,ul ,pl ,cl real(dp),dimension(1:nvector),save::rr ,ur ,pr ,cr real(dp),dimension(1:nvector),save::ro ,uo ,po ,co real(dp),dimension(1:nvector),save::rstar,ustar,pstar,cstar real(dp),dimension(1:nvector),save::wl ,wr ,wo real(dp),dimension(1:nvector),save::sgnm ,spin ,spout,ushock real(dp),dimension(1:nvector),save::frac ! constants smallp = smallc**2/gamma entho = one/(gamma-one) ! Initial states pressure, density and velocity do i=1,ngrid rl(i)=max(qleft (i,1),smallr) ul(i)= qleft (i,2) pl(i)=max(qleft (i,3),rl(i)*smallp) rr(i)=max(qright(i,1),smallr) ur(i)= qright(i,2) pr(i)=max(qright(i,3),rr(i)*smallp) end do ! Acoustic star state do i=1,ngrid cl(i) = sqrt(gamma*pl(i)/rl(i)) cr(i) = sqrt(gamma*pr(i)/rr(i)) wl(i) = cl(i)*rl(i) wr(i) = cr(i)*rr(i) pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i))) & & / (wl(i)+wr(i)) ustar(i) = ((wr(i)*ur(i)+wl(i)*ul(i))+(pl(i)-pr(i))) & & / (wl(i)+wr(i)) !! pstar(i) = MAX(pstar(i),zero) end do ! Left going or right going contact wave do i=1,ngrid sgnm(i) = sign(one,ustar(i)) end do ! Left or right unperturbed state do i=1,ngrid if(sgnm(i)==one)then ro(i) = rl(i) uo(i) = ul(i) po(i) = pl(i) wo(i) = wl(i) co(i) = cl(i) else ro(i) = rr(i) uo(i) = ur(i) po(i) = pr(i) wo(i) = wr(i) co(i) = cr(i) end if end do ! Star region density and sound speed do i=1,ngrid rstar(i) = ro(i)+(pstar(i)-po(i))/co(i)**2 rstar(i) = max(rstar(i),smallr) cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i))) cstar(i) = max(cstar(i),smallc) end do ! Head and tail speed of rarefaction do i=1,ngrid spout(i) = co (i)-sgnm(i)*uo (i) spin (i) = cstar(i)-sgnm(i)*ustar(i) end do ! Shock speed do i=1,ngrid ushock(i) = half*(spin(i)+spout(i)) ushock(i) = max(ushock(i),-sgnm(i)*ustar(i)) end do do i=1,ngrid if(pstar(i)>=po(i))then spout(i)=ushock(i) spin (i)=spout (i) end if end do ! Sample the solution at x/t=0 do i=1,ngrid if(spout(i)=zero)then ! Star region qgdnv(i,1) = rstar(i) qgdnv(i,2) = ustar(i) qgdnv(i,3) = pstar(i) else ! Rarefaction frac(i) = spout(i)/(spout(i)-spin(i)) qgdnv(i,1) = frac(i)*rstar(i) + (one - frac(i))*ro(i) qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i) qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i) end if end do ! Passive scalars #if NVAR > 3 do n = 4,nvar do i=1,ngrid if(sgnm(i)==one)then qgdnv(i,n) = qleft (i,n) else qgdnv(i,n) = qright(i,n) end if end do end do #endif ! Specific internal energy do i=1,ngrid qgdnv(i,nvar+1) = po(i)/ro(i)*entho end do ! Compute fluxes do i = 1, ngrid fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2) ! Mass density fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2 ! Normal momentum etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2 #if NDIM>1 etot = etot + half*qgdnv(i,1)*qgdnv(i,4)**2 #endif #if NDIM>2 etot = etot + half*qgdnv(i,1)*qgdnv(i,5)**2 #endif fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3)) ! Total energy end do ! Other advected quantities do n = 4, nvar+1 do i = 1, ngrid fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n) end do end do end subroutine riemann_acoustic !########################################################### !########################################################### !########################################################### !########################################################### subroutine riemann_llf(qleft,qright,fgdnv,ngrid) use amr_parameters use hydro_parameters use const implicit none ! dummy arguments integer::ngrid real(dp),dimension(1:nvector,1:nvar)::qleft,qright real(dp),dimension(1:nvector,1:nvar+1)::fgdnv ! local arrays real(dp),dimension(1:nvector,1:nvar+1),save::fleft,fright real(dp),dimension(1:nvector,1:nvar+1),save::uleft,uright real(dp),dimension(1:nvector),save::cmax ! local variables integer::i,n real(dp)::smallp, entho real(dp)::rl ,ul ,pl ,cl real(dp)::rr ,ur ,pr ,cr ! Constants smallp = smallc**2/gamma entho = one/(gamma-one) !=========================== ! Compute maximum wave speed !=========================== do i = 1,ngrid ! Left states rl = max(qleft (i,1),smallr) ul = qleft (i,2) pl = max(qleft (i,3),rl*smallp) cl = gamma*pl #if NENER>0 do n = 1,nener cl = cl + gamma_rad(n)*qleft(i,ndim+2+n) end do #endif cl = sqrt(cl/rl) ! Right states rr = max(qright(i,1),smallr) ur = qright(i,2) pr = max(qright(i,3),rr*smallp) cr = gamma*pr #if NENER>0 do n = 1,nener cr = cr + gamma_rad(n)*qright(i,ndim+2+n) end do #endif cr = sqrt(cr/rr) ! Local max. wave speed cmax(i) = max(abs(ul)+cl,abs(ur)+cr) end do !=============================== ! Compute conservative variables !=============================== do i = 1, ngrid ! Mass density uleft (i,1) = qleft (i,1) uright(i,1) = qright(i,1) ! Normal momentum uleft (i,2) = qleft (i,1)*qleft (i,2) uright(i,2) = qright(i,1)*qright(i,2) ! Total energy uleft (i,3) = qleft (i,3)*entho + half*qleft (i,1)*qleft (i,2)**2 uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2 #if NDIM>1 uleft (i,3) = uleft (i,3) + half*qleft (i,1)*qleft (i,4)**2 uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,4)**2 #endif #if NDIM>2 uleft (i,3) = uleft (i,3) + half*qleft (i,1)*qleft (i,5)**2 uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,5)**2 #endif #if NENER>0 do n = 1,nener uleft (i,3) = uleft (i,3) + qleft (i,ndim+2+n)/(gamma_rad(n)-one) uright(i,3) = uright(i,3) + qright(i,ndim+2+n)/(gamma_rad(n)-one) end do #endif end do ! Transverse velocities #if NDIM > 1 do n = 4, ndim+2 do i = 1, ngrid uleft (i,n) = qleft (i,1)*qleft (i,n) uright(i,n) = qright(i,1)*qright(i,n) end do end do #endif ! Non-thermal energies #if NENER>0 do n = 1, nener do i = 1, ngrid uleft (i,ndim+2+n) = qleft (i,ndim+2+n)/(gamma_rad(n)-one) uright(i,ndim+2+n) = qright(i,ndim+2+n)/(gamma_rad(n)-one) end do end do #endif ! Other passively advected quantities #if NVAR > 2+NDIM+NENER do n = 3+ndim+nener, nvar do i = 1, ngrid uleft (i,n) = qleft (i,1)*qleft (i,n) uright(i,n) = qright(i,1)*qright(i,n) end do end do #endif ! Thermal energy do i = 1, ngrid uleft (i,nvar+1) = qleft (i,3)*entho uright(i,nvar+1) = qright(i,3)*entho end do !============================== ! Compute left and right fluxes !============================== do i = 1, ngrid ! Mass density fleft (i,1) = qleft (i,2)*uleft (i,1) fright(i,1) = qright(i,2)*uright(i,1) ! Normal momentum fleft (i,2) = qleft (i,2)*uleft (i,2) + qleft (i,3) fright(i,2) = qright(i,2)*uright(i,2) + qright(i,3) #if NENER>0 do n = 1,nener fleft (i,2) = fleft (i,2) + qleft (i,ndim+2+n) fright(i,2) = fright(i,2) + qright(i,ndim+2+n) end do #endif ! Total energy fleft (i,3) = qleft (i,2)*(uleft (i,3)+qleft (i,3)) fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3)) #if NENER>0 do n = 1,nener fleft (i,3) = fleft (i,3) + qleft (i,2)*qleft (i,ndim+2+n) fright(i,3) = fright(i,3) + qright(i,2)*qright(i,ndim+2+n) end do #endif end do ! Other passively advected quantities do n = 4, nvar+1 do i = 1, ngrid fleft (i,n) = qleft (i,2)*uleft (i,n) fright(i,n) = qright(i,2)*uright(i,n) end do end do !============================= ! Compute Lax-Friedrich fluxes !============================= do n = 1, nvar+1 do i = 1, ngrid fgdnv(i,n) = half*(fleft(i,n)+fright(i,n)-cmax(i)*(uright(i,n)-uleft(i,n))) end do end do end subroutine riemann_llf !########################################################### !########################################################### !########################################################### !########################################################### subroutine riemann_hll(qleft,qright,fgdnv,ngrid) USE amr_parameters USE const USE hydro_parameters ! 1D HLL Riemann solver IMPLICIT NONE integer::ngrid real(dp),dimension(1:nvector,1:nvar)::qleft,qright real(dp),dimension(1:nvector,1:nvar+1)::fgdnv real(dp),dimension(1:nvector,1:nvar+1),save::fleft,fright real(dp),dimension(1:nvector,1:nvar+1),save::uleft,uright real(dp),dimension(1:nvector),save::SL,SR integer::i,n real(dp)::smallp, entho real(dp)::rl ,ul ,pl ,cl real(dp)::rr ,ur ,pr ,cr ! Constants smallp = smallc**2/gamma entho = one/(gamma-one) !=========================== ! Compute maximum wave speed !=========================== do i=1,ngrid ! Left states rl = max(qleft (i,1),smallr) ul = qleft (i,2) pl = max(qleft (i,3),rl*smallp) cl = gamma*pl #if NENER>0 do n = 1,nener cl = cl + gamma_rad(n)*qleft(i,ndim+2+n) end do #endif cl = sqrt(cl/rl) ! Right states rr = max(qright(i,1),smallr) ur = qright(i,2) pr = max(qright(i,3),rr*smallp) cr = gamma*pr #if NENER>0 do n = 1,nener cr = cr + gamma_rad(n)*qright(i,ndim+2+n) end do #endif cr = sqrt(cr/rr) ! Left and right max. wave speed SL(i)=min(min(ul,ur)-max(cl,cr),zero) SR(i)=max(max(ul,ur)+max(cl,cr),zero) end do !=============================== ! Compute conservative variables !=============================== do i = 1, ngrid ! Mass density uleft (i,1) = qleft (i,1) uright(i,1) = qright(i,1) ! Normal momentum uleft (i,2) = qleft (i,1)*qleft (i,2) uright(i,2) = qright(i,1)*qright(i,2) ! Total energy uleft (i,3) = qleft (i,3)*entho + half*qleft (i,1)*qleft (i,2)**2 uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2 #if NDIM>1 uleft (i,3) = uleft (i,3) + half*qleft (i,1)*qleft (i,4)**2 uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,4)**2 #endif #if NDIM>2 uleft (i,3) = uleft (i,3) + half*qleft (i,1)*qleft (i,5)**2 uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,5)**2 #endif #if NENER>0 do n = 1,nener uleft (i,3) = uleft (i,3) + qleft (i,ndim+2+n)/(gamma_rad(n)-one) uright(i,3) = uright(i,3) + qright(i,ndim+2+n)/(gamma_rad(n)-one) end do #endif end do ! Transverse velocities #if NDIM > 1 do n = 4, ndim+2 do i = 1, ngrid uleft (i,n) = qleft (i,1)*qleft (i,n) uright(i,n) = qright(i,1)*qright(i,n) end do end do #endif ! Non-thermal energies #if NENER>0 do n = 1, nener do i = 1, ngrid uleft (i,ndim+2+n) = qleft (i,ndim+2+n)/(gamma_rad(n)-one) uright(i,ndim+2+n) = qright(i,ndim+2+n)/(gamma_rad(n)-one) end do end do #endif ! Other passively advected quantities #if NVAR > 2+NDIM+NENER do n = 3+ndim+nener, nvar do i = 1, ngrid uleft (i,n) = qleft (i,1)*qleft (i,n) uright(i,n) = qright(i,1)*qright(i,n) end do end do #endif ! Thermal energy do i = 1, ngrid uleft (i,nvar+1) = qleft (i,3)*entho uright(i,nvar+1) = qright(i,3)*entho end do !============================== ! Compute left and right fluxes !============================== do i = 1, ngrid ! Mass density fleft (i,1) = uleft (i,2) fright(i,1) = uright(i,2) ! Normal momentum fleft (i,2) = qleft (i,3)+uleft (i,2)*qleft (i,2) fright(i,2) = qright(i,3)+uright(i,2)*qright(i,2) #if NENER>0 do n = 1,nener fleft (i,2) = fleft (i,2) + qleft (i,ndim+2+n) fright(i,2) = fright(i,2) + qright(i,ndim+2+n) end do #endif ! Total energy fleft (i,3) = qleft (i,2)*(uleft (i,3)+qleft (i,3)) fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3)) #if NENER>0 do n = 1,nener fleft (i,3) = fleft (i,3) + qleft (i,2)*qleft (i,ndim+2+n) fright(i,3) = fright(i,3) + qright(i,2)*qright(i,ndim+2+n) end do #endif end do ! Other advected quantities do n = 4, nvar+1 do i = 1, ngrid fleft (i,n) = qleft (i,2)*uleft (i,n) fright(i,n) = qright(i,2)*uright(i,n) end do end do !=================== ! Compute HLL fluxes !=================== do n = 1, nvar+1 do i = 1, ngrid fgdnv(i,n) = (SR(i)*fleft(i,n)-SL(i)*fright(i,n) & & + SR(i)*SL(i)*(uright(i,n)-uleft(i,n)))/(SR(i)-SL(i)) end do end do end subroutine riemann_hll !########################################################### !########################################################### !########################################################### !########################################################### subroutine riemann_hllc(qleft,qright,snleft,snright,fgdnv,ngrid) use amr_parameters use hydro_parameters use const implicit none ! HLLC Riemann solver (Toro) integer::ngrid real(dp),dimension(1:nvector)::snleft,snright real(dp),dimension(1:nvector,1:nvar)::qleft,qright real(dp),dimension(1:nvector,1:nvar+1)::fgdnv REAL(dp)::SL,SR REAL(dp)::entho REAL(dp)::rl,pl,ul,ecinl,etotl,el,ptotl REAL(dp)::rr,pr,ur,ecinr,etotr,er,ptotr REAL(dp)::cfastl,rcl,rstarl,estarl REAL(dp)::cfastr,rcr,rstarr,estarr REAL(dp)::etotstarl,etotstarr REAL(dp)::ustar,ptotstar REAL(dp)::ro,uo,ptoto,etoto,eo REAL(dp)::smallp INTEGER::ivar,i #if NENER>0 REAL(dp),dimension(1:nener)::eradl,eradr,erado REAL(dp),dimension(1:nener)::eradstarl,eradstarr INTEGER::irad #endif ! constants smallp = smallc**2/gamma entho = one/(gamma-one) do i=1,ngrid ! Left variables rl=max(qleft (i,1),smallr) Pl=max(qleft (i,3),rl*smallp) ul= qleft (i,2) el=Pl*entho ecinl = half*rl*ul*ul #if NDIM>1 ecinl=ecinl+half*rl*qleft(i,4)**2 #endif #if NDIM>2 ecinl=ecinl+half*rl*qleft(i,5)**2 #endif etotl = el+ecinl #if NENER>0 do irad=1,nener eradl(irad)=qleft(i,2+ndim+irad)/(gamma_rad(irad)-one) etotl=etotl+eradl(irad) end do #endif Ptotl = Pl #if NENER>0 do irad=1,nener Ptotl=Ptotl+qleft(i,2+ndim+irad) end do #endif Ptotl=Ptotl+snleft(i) ! Right variables rr=max(qright(i,1),smallr) Pr=max(qright(i,3),rr*smallp) ur= qright(i,2) er=Pr*entho ecinr = half*rr*ur*ur #if NDIM>1 ecinr=ecinr+half*rr*qright(i,4)**2 #endif #if NDIM>2 ecinr=ecinr+half*rr*qright(i,5)**2 #endif etotr = er+ecinr #if NENER>0 do irad=1,nener eradr(irad)=qright(i,2+ndim+irad)/(gamma_rad(irad)-one) etotr=etotr+eradr(irad) end do #endif Ptotr = Pr #if NENER>0 do irad=1,nener Ptotr=Ptotr+qright(i,2+ndim+irad) end do #endif Ptotr=Ptotr+snright(i) ! Find the largest eigenvalues in the normal direction to the interface cfastl=gamma*Pl #if NENER>0 do irad = 1,nener cfastl = cfastl + gamma_rad(irad)*qleft(i,ndim+2+irad) end do #endif cfastl = cfastl + snleft(i) cfastl=sqrt(max(cfastl/rl,smallc**2)) cfastr=gamma*Pr #if NENER>0 do irad = 1,nener cfastr = cfastr + gamma_rad(irad)*qright(i,ndim+2+irad) end do #endif cfastr = cfastr + snright(i) cfastr=sqrt(max(cfastr/rr,smallc**2)) ! Compute HLL wave speed SL=min(ul,ur)-max(cfastl,cfastr) SR=max(ul,ur)+max(cfastl,cfastr) ! Compute lagrangian sound speed rcl=rl*(ul-SL) rcr=rr*(SR-ur) ! Compute acoustic star state ustar =(rcr*ur +rcl*ul + (Ptotl-Ptotr))/(rcr+rcl) Ptotstar=(rcr*Ptotl+rcl*Ptotr+rcl*rcr*(ul-ur))/(rcr+rcl) ! Left star region variables rstarl=rl*(SL-ul)/(SL-ustar) etotstarl=((SL-ul)*etotl-Ptotl*ul+Ptotstar*ustar)/(SL-ustar) estarl=el*(SL-ul)/(SL-ustar) #if NENER>0 do irad = 1,nener eradstarl(irad)=eradl(irad)*(SL-ul)/(SL-ustar) end do #endif ! Right star region variables rstarr=rr*(SR-ur)/(SR-ustar) etotstarr=((SR-ur)*etotr-Ptotr*ur+Ptotstar*ustar)/(SR-ustar) estarr=er*(SR-ur)/(SR-ustar) #if NENER>0 do irad = 1,nener eradstarr(irad)=eradr(irad)*(SR-ur)/(SR-ustar) end do #endif ! Sample the solution at x/t=0 if(SL>0d0)then ro=rl uo=ul Ptoto=Ptotl etoto=etotl eo=el #if NENER>0 do irad = 1,nener erado(irad)=eradl(irad) end do #endif else if(ustar>0d0)then ro=rstarl uo=ustar Ptoto=Ptotstar etoto=etotstarl eo=estarl #if NENER>0 do irad = 1,nener erado(irad)=eradstarl(irad) end do #endif else if (SR>0d0)then ro=rstarr uo=ustar Ptoto=Ptotstar etoto=etotstarr eo=estarr #if NENER>0 do irad = 1,nener erado(irad)=eradstarr(irad) end do #endif else ro=rr uo=ur Ptoto=Ptotr etoto=etotr eo=er #if NENER>0 do irad = 1,nener erado(irad)=eradr(irad) end do #endif end if !========================= ! Compute the Godunov flux !========================= fgdnv(i,1) = ro*uo fgdnv(i,2) = ro*uo*uo+Ptoto fgdnv(i,3) = (etoto+Ptoto)*uo ! Transverse velocities #if NDIM > 1 do ivar = 4,ndim+2 if(ustar>0)then fgdnv(i,ivar) = ro*uo*qleft(i,ivar) else fgdnv(i,ivar) = ro*uo*qright(i,ivar) endif end do #endif ! Non-thermal energies #if NENER>0 do irad = 1,nener fgdnv(i,ndim+2+irad) = uo*erado(irad) end do #endif ! Other passively advected quantities #if NVAR > 2+NDIM+NENER do ivar = 3+ndim+nener,nvar if(ustar>0)then fgdnv(i,ivar) = ro*uo*qleft(i,ivar) else fgdnv(i,ivar) = ro*uo*qright(i,ivar) endif end do #endif ! Thermal energy fgdnv(i,nvar+1) = uo*eo end do end subroutine riemann_hllc !########################################################### !########################################################### !########################################################### !###########################################################/home/hemv_g/srf_a/RAMSES_CODE_CORI/ramses06062019_grackle311Cosmo/patch/mom2/star_formation.f90 !################################################################ !################################################################ !################################################################ !################################################################ subroutine star_formation(ilevel) use amr_commons use pm_commons use hydro_commons use poisson_commons use cooling_module, ONLY: XH=>X, rhoc, mH , twopi use random use mpi_mod implicit none #ifndef WITHOUTMPI integer::info,info2,dummy_io integer,parameter::tag=1120 #endif integer::ilevel !---------------------------------------------------------------------- ! Description: This subroutine spawns star-particle of constant mass ! using a Poisson probability law if some gas condition are fulfilled. ! It modifies hydrodynamic variables according to mass conservation ! and assumes an isothermal transformation... ! On exit, the gas velocity and sound speed are unchanged. ! New star particles are synchronized with other collisionless particles. ! Array flag2 is used as temporary work space. ! Yann Rasera 10/2002-01/2003 !---------------------------------------------------------------------- ! local constants real(dp)::t0,d0,d00,mgas,mcell real(dp)::scale_nH,scale_T2,scale_l,scale_d,scale_t,scale_v real(dp),dimension(1:twotondim,1:3)::xc ! other variables integer ::ncache,nnew,ivar,ngrid,icpu,index_star,ndebris_tot,ilun integer ::igrid,ix,iy,iz,ind,i,n,iskip,nx_loc,idim integer ::ntot,ntot_all,nstar_corrected logical ::ok_free real(dp)::d,x,y,z,u,v,w,e,tg,zg real(dp)::mstar,dstar,tstar,nISM,nCOM,phi_t,phi_x,theta,sigs,scrit,b_turb real(dp)::T2,nH,T_poly,cs2,cs2_poly,trel,t_dyn,t_ff,uvar real(dp)::alpha0,npoly real(dp)::sigma2 real(dp)::birth_epoch,factG,M2 real(kind=8)::mlost,mtot,mlost_all,mtot_all real(kind=8)::PoissMean real(dp),parameter::pi=0.5*twopi real(dp),dimension(1:3)::skip_loc real(dp)::dx,dx_loc,scale,vol_loc,dx_min,vol_min real(dp),dimension(1:nvector)::sfr_ff integer ,dimension(1:ncpu,1:IRandNumSize)::allseed integer ,dimension(1:nvector),save::ind_grid,ind_cell,nstar integer ,dimension(1:nvector),save::ind_grid_new,ind_cell_new,ind_part logical ,dimension(1:nvector),save::ok,ok_new=.true. integer ,dimension(1:ncpu)::ntot_star_cpu,ntot_star_all character(LEN=80)::filename,filedir,fileloc,filedirini character(LEN=5)::nchar,ncharcpu logical::file_exist #ifdef SOLVERmhd real(dp)::bx1,bx2,by1,by2,bz1,bz2,A,B,C,emag,beta,fbeta #endif #if NENER>0 integer::irad #endif if(numbtot(1,ilevel)==0) return if(.not. hydro)return if(ndim.ne.3)return if(static)return if(verbose)write(*,*)' Entering star_formation' if(sf_log_properties.and.ifout.gt.1) then call title(ifout-1,nchar) if(IOGROUPSIZEREP>0) then filedirini='output_'//TRIM(nchar)//'/' filedir='output_'//TRIM(nchar)//'/group_'//TRIM(ncharcpu)//'/' else filedir='output_'//TRIM(nchar)//'/' endif filename=TRIM(filedir)//'stars_'//TRIM(nchar)//'.out' ilun=myid+10 call title(myid,nchar) fileloc=TRIM(filename)//TRIM(nchar) ! Wait for the token #ifndef WITHOUTMPI if(IOGROUPSIZE>0) then if (mod(myid-1,IOGROUPSIZE)/=0) then call MPI_RECV(dummy_io,1,MPI_INTEGER,myid-1-1,tag,& & MPI_COMM_WORLD,MPI_STATUS_IGNORE,info2) end if endif #endif inquire(file=fileloc,exist=file_exist) if((.not.file_exist).or.(abs(t-trestart).lt.dtnew(ilevel))) then open(ilun, file=fileloc, form='formatted') write(ilun,'(A24)',advance='no') '# event id ilevel mp ' do idim=1,ndim write(ilun,'(A2,I1,A2)',advance='no') 'xp',idim,' ' enddo do idim=1,ndim write(ilun,'(A2,I1,A2)',advance='no') 'vp',idim,' ' enddo do ivar=1,nvar if(ivar.ge.10) then write(ilun,'(A1,I2,A2)',advance='no') 'u',ivar,' ' else write(ilun,'(A1,I1,A2)',advance='no') 'u',ivar,' ' endif enddo write(ilun,'(A5)',advance='no') 'tag ' write(ilun,'(A1)') ' ' else open(ilun, file=fileloc, status='old', position='append', action='write', form='formatted') endif endif ! Conversion factor from user units to cgs units call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2) ! Mesh spacing in that level dx=0.5D0**ilevel nx_loc=(icoarse_max-icoarse_min+1) skip_loc=(/0.0d0,0.0d0,0.0d0/) if(ndim>0)skip_loc(1)=dble(icoarse_min) if(ndim>1)skip_loc(2)=dble(jcoarse_min) if(ndim>2)skip_loc(3)=dble(kcoarse_min) scale=boxlen/dble(nx_loc) dx_loc=dx*scale vol_loc=dx_loc**ndim dx_min=(0.5D0**nlevelmax)*scale vol_min=dx_min**ndim ! Star formation time scale from Gyr to code units ! SFR apply here for long lived stars only t0=t_star*(1d9*365.*24.*3600.)/scale_t trel=sf_trelax*1d6*(365.*24.*3600.)/scale_t ! ISM density threshold from H/cc to code units nISM = n_star if(cosmo)then nCOM = del_star*omega_b*rhoc*(h0/100.)**2/aexp**3*XH/mH nISM = MAX(nCOM,nISM) endif d0 = nISM/scale_nH d00 = n_star/scale_nH npoly=n_poly_fact*n_star ! Initial star particle mass if(m_star < 0d0)then mstar=n_star/(scale_nH*aexp**3)*vol_min else mstar=m_star*mass_sph endif dstar=mstar/vol_loc factG = 1d0 if(cosmo) factG = 3d0/4d0/twopi*omega_m*aexp ! Birth epoch as proper time if(use_proper_time)then birth_epoch=texp else birth_epoch=t endif ! Cells center position relative to grid center position do ind=1,twotondim iz=(ind-1)/4 iy=(ind-1-4*iz)/2 ix=(ind-1-2*iy-4*iz) xc(ind,1)=(dble(ix)-0.5D0)*dx xc(ind,2)=(dble(iy)-0.5D0)*dx xc(ind,3)=(dble(iz)-0.5D0)*dx end do ! If necessary, initialize random number generator if(localseed(1)==-1)then call rans(ncpu,iseed,allseed) localseed=allseed(myid,1:IRandNumSize) end if #if NDIM==3 !------------------------------------------------ ! Convert hydro variables to primitive variables !------------------------------------------------ ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do do i=1,ngrid d=uold(ind_cell(i),1) u=uold(ind_cell(i),2)/d v=uold(ind_cell(i),3)/d w=uold(ind_cell(i),4)/d e=uold(ind_cell(i),5) #ifdef SOLVERmhd bx1=uold(ind_cell(i),6) by1=uold(ind_cell(i),7) bz1=uold(ind_cell(i),8) bx2=uold(ind_cell(i),nvar+1) by2=uold(ind_cell(i),nvar+2) bz2=uold(ind_cell(i),nvar+3) e=e-0.125d0*((bx1+bx2)**2+(by1+by2)**2+(bz1+bz2)**2) #endif e=e-0.5d0*d*(u**2+v**2+w**2) #if NENER>0 do irad=0,nener-1 e=e-uold(ind_cell(i),inener+irad) end do #endif uold(ind_cell(i),1)=d uold(ind_cell(i),2)=u uold(ind_cell(i),3)=v uold(ind_cell(i),4)=w uold(ind_cell(i),5)=e/d end do do ivar=imetal,nvar do i=1,ngrid d=uold(ind_cell(i),1) w=uold(ind_cell(i),ivar)/d uold(ind_cell(i),ivar)=w end do end do end do end do ! get values of uold for density and velocities in virtual boundaries #ifndef WITHOUTMPI do ivar=1,4 call make_virtual_fine_dp(uold(1,ivar),ilevel) end do #endif !------------------------------------------------ ! Compute number of new stars in each cell !------------------------------------------------ ntot=0 ndebris_tot=0 ! Loop over grids ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Star formation criterion ---> logical array ok(i) do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do ! Flag leaf cells do i=1,ngrid ok(i)=son(ind_cell(i))==0 end do if(sf_model > 0)then do i=1,ngrid ! if cell is a leaf cell if (ok(i)) then ! Gas density d = uold(ind_cell(i),1) ! Compute temperature in K/mu T2 = (gamma-1.0)*uold(ind_cell(i),5)*scale_T2 ! Correct from polytrope T_poly = T2_star*(uold(ind_cell(i),1)*scale_nH/npoly)**(g_star-1.0) T2 = T2-T_poly ! Compute sound speed squared cs2 = (gamma-1.0)*uold(ind_cell(i),5) ! prevent numerical crash due to negative temperature cs2 = max(cs2,smallc**2) ! Correct from polytrope cs2_poly = (T2_star/scale_T2)*(uold(ind_cell(i),1)*scale_nH/npoly)**(g_star-1.0) cs2 = cs2-cs2_poly ! Turbulence 1D velocity dispersion sigma2 = uold(ind_cell(i),ivirial1)*2.0/3.0 ! Density criterion if(d<=d0) ok(i)=.false. if(ok(i)) then SELECT CASE (sf_model) ! Multi-ff KM model CASE (1) ! Virial parameter alpha0 = (5.0*(sigma2+cs2))/(pi*factG*d*dx_loc**2) M2 = max(sigma2/cs2,1.0) ! Turbulent forcing parameter (Federrath 2008 & 2010) b_turb = 0.4 ! Fudge for alpha dependence (KM 2005). phi_x = 1.12 ! The prefered value for eps_star = 1.0, ! which represents the theoretical maximum efficiency. sigs = log(1.0+(b_turb**2)*(M2)) scrit = log(((pi**2)/5)*(phi_x**2)*alpha0*(M2)) sfr_ff(i) = (eps_star/2.0)*exp(3.0/8.0*sigs)*(2.0-erfc((sigs-scrit)/sqrt(2.0*sigs))) ! Multi-ff PN model CASE (2) ! Virial parameter alpha0 = (5.0*(sigma2+cs2))/(pi*factG*d*dx_loc**2) b_turb = 0.4 ! Best fit values to the Multi-ff PN model (Hydro) phi_t = 1.0/0.49 theta = 0.97 sigs = log(1.0+(b_turb**2)*(sigma2/cs2)) scrit = log(0.067/(theta**2)*alpha0*(sigma2/cs2)) sfr_ff(i) = (eps_star/(2.0*phi_t))*exp(3.0/8.0*sigs)*(2.0-erfc((sigs-scrit)/sqrt(2.0*sigs))) ! Virial criterion threshold a la Hopkins CASE (3) alpha0 = (5.0*(sigma2+cs2))/(pi*factG*d*dx_loc**2) if(alpha0<1.0) then sfr_ff(i) = eps_star else sfr_ff(i) = 0.0 ok(i) = .false. endif ! Padoan 2012 a la Semenov CASE (4) ! Feedback efficiency t_dyn = dx_loc/(2.0*sqrt(sigma2+cs2)) t_ff = 0.5427*sqrt(1.0/(factG*max(d,smallr))) sfr_ff(i) = eps_star*exp(-1.6*t_ff/t_dyn) CASE (5) ! Virial parameter alpha0 = (5.0*(sigma2+cs2))/(pi*factG*d*dx_loc**2) M2 = max(sigma2/cs2,smallr) ! Turbulent forcing parameter (Federrath 2008 & 2010) b_turb = 0.4 ! Fudge for alpha dependence (KM 2005). ! phi_x = 1.12 ! The prefered value for eps_star = 1.0, ! which represents the theoretical maximum efficiency. sigs = log(1.0+(b_turb**2)*(M2)) ! scrit = log(alpha0*(1.0+(M2*(pi**2)/5.0)*(phi_x**2))) scrit = log(alpha0*(1.0+(2*(M2**2)/(1.0+M2)))) sfr_ff(i) = (eps_star/2.0)*exp(3.0/8.0*sigs)*(2.0-erfc((sigs-scrit)/sqrt(2.0*sigs))) END SELECT endif endif end do else ! Density criterion do i=1,ngrid sfr_ff(i) = eps_star d=uold(ind_cell(i),1) if(d<=d0)ok(i)=.false. end do ! Temperature criterion do i=1,ngrid T2=uold(ind_cell(i),5)*scale_T2*(gamma-1.0) nH=max(uold(ind_cell(i),1),smallr)*scale_nH T_poly=T2_star*(nH/npoly)**(g_star-1.0) T2=T2-T_poly if(T2>2e4)ok(i)=.false. end do endif ! Geometrical criterion if(ivar_refine>0)then do i=1,ngrid d=uold(ind_cell(i),ivar_refine) if(d<=var_cut_refine)ok(i)=.false. end do endif ! Calculate number of new stars in each cell using Poisson statistics do i=1,ngrid nstar(i)=0 if(ok(i))then ! Compute mean number of events d=uold(ind_cell(i),1) mcell=d*vol_loc ! Free fall time of an homogeneous sphere tstar= .5427*sqrt(1.0/(factG*max(d,smallr))) ! Gas mass to be converted into stars mgas=dtnew(ilevel)*(sfr_ff(i)/tstar)*mcell ! Poisson mean PoissMean=mgas/mstar if((trel>0.).and.(.not.cosmo)) PoissMean = PoissMean*min((t/trel),1.0) ! Compute Poisson realisation call poissdev(localseed,PoissMean,nstar(i)) ! Compute depleted gas mass mgas=nstar(i)*mstar ! Security to prevent more than 90% of gas depletion if (mgas > 0.9*mcell) then nstar_corrected=int(0.9*mcell/mstar) mstar_lost=mstar_lost+(nstar(i)-nstar_corrected)*mstar nstar(i)=nstar_corrected endif ! Compute new stars local statistics mstar_tot=mstar_tot+nstar(i)*mstar if(nstar(i)>0)then ntot=ntot+1 endif endif enddo ! Store nstar in array flag2 do i=1,ngrid flag2(ind_cell(i))=nstar(i) end do end do end do !--------------------------------- ! Check for free particle memory !--------------------------------- ok_free=(numbp_free-ntot-ndebris_tot)>=0 #ifndef WITHOUTMPI call MPI_ALLREDUCE(numbp_free,numbp_free_tot,1,MPI_INTEGER,MPI_MIN,MPI_COMM_WORLD,info) #endif #ifdef WITHOUTMPI numbp_free_tot=numbp_free #endif if(.not. ok_free)then write(*,*)'No more free memory for particles' write(*,*)'Increase npartmax' #ifndef WITHOUTMPI call MPI_ABORT(MPI_COMM_WORLD,1,info) #else stop #endif end if !--------------------------------- ! Compute global stars statistics !--------------------------------- #ifndef WITHOUTMPI mlost=mstar_lost; mtot=mstar_tot call MPI_ALLREDUCE(ntot,ntot_all,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,info) call MPI_ALLREDUCE(mtot,mtot_all,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) call MPI_ALLREDUCE(mlost,mlost_all,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info) #endif #ifdef WITHOUTMPI ntot_all=ntot mtot_all=mstar_tot mlost_all=mstar_lost #endif ntot_star_cpu=0; ntot_star_all=0 ntot_star_cpu(myid)=ntot #ifndef WITHOUTMPI call MPI_ALLREDUCE(ntot_star_cpu,ntot_star_all,ncpu,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,info) ntot_star_cpu(1)=ntot_star_all(1) #endif do icpu=2,ncpu ntot_star_cpu(icpu)=ntot_star_cpu(icpu-1)+ntot_star_all(icpu) end do nstar_tot=nstar_tot+ntot_all if(myid==1)then if(ntot_all.gt.0)then write(*,'(' Level=',I6,' New star=',I6,' Tot=',I10,' Mass=',1PE10.3,' Lost=',0PF5.1,'%')')& & ilevel,ntot_all,nstar_tot,mtot_all,mlost_all/(mlost_all+mtot_all)*100. endif end if !------------------------------ ! Create new star particles !------------------------------ ! Starting identity number if(myid==1)then index_star=nstar_tot-ntot_all else index_star=nstar_tot-ntot_all+ntot_star_cpu(myid-1) end if ! Loop over grids ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do ! Loop over cells do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do ! Flag cells with at least one new star do i=1,ngrid ok(i)=flag2(ind_cell(i))>0 end do ! Gather new star arrays nnew=0 do i=1,ngrid if (ok(i))then nnew=nnew+1 ind_grid_new(nnew)=ind_grid(i) ind_cell_new(nnew)=ind_cell(i) end if end do ! Update linked list for stars call remove_free(ind_part,nnew) call add_list(ind_part,ind_grid_new,ok_new,nnew) ! Calculate new star particle and modify gas density do i=1,nnew index_star=index_star+1 ! Get gas variables n=flag2(ind_cell_new(i)) d=uold(ind_cell_new(i),1) u=uold(ind_cell_new(i),2) v=uold(ind_cell_new(i),3) w=uold(ind_cell_new(i),4) x=(xg(ind_grid_new(i),1)+xc(ind,1)-skip_loc(1))*scale y=(xg(ind_grid_new(i),2)+xc(ind,2)-skip_loc(2))*scale z=(xg(ind_grid_new(i),3)+xc(ind,3)-skip_loc(3))*scale tg=uold(ind_cell_new(i),5)*(gamma-1)*scale_T2 if(metal)zg=uold(ind_cell_new(i),imetal) ! Set star particle variables tp(ind_part(i))=birth_epoch ! Birth epoch mp(ind_part(i))=n*mstar ! Mass levelp(ind_part(i))=ilevel ! Level idp(ind_part(i))=index_star ! Star identity typep(ind_part(i))%family=FAM_STAR typep(ind_part(i))%tag=0 xp(ind_part(i),1)=x xp(ind_part(i),2)=y xp(ind_part(i),3)=z vp(ind_part(i),1)=u vp(ind_part(i),2)=v vp(ind_part(i),3)=w if(metal)zp(ind_part(i))=zg ! Initial star metallicity if(sf_log_properties) then write(ilun,'(I10)',advance='no') 0 write(ilun,'(2I10,E24.12)',advance='no') idp(ind_part(i)),ilevel,mp(ind_part(i)) do idim=1,ndim write(ilun,'(E24.12)',advance='no') xp(ind_part(i),idim) enddo do idim=1,ndim write(ilun,'(E24.12)',advance='no') vp(ind_part(i),idim) enddo write(ilun,'(E24.12)',advance='no') uold(ind_cell_new(i),1) do ivar=2,nvar if(ivar.eq.ndim+2)then ! Temperature uvar=(gamma-1.0)*(uold(ind_cell_new(i),ndim+2))*scale_T2 else uvar=uold(ind_cell_new(i),ivar) endif write(ilun,'(E24.12)',advance='no') uvar enddo write(ilun,'(I10)',advance='no') typep(ind_part(i))%tag write(ilun,'(A1)') ' ' endif end do ! End loop over new star particles ! Modify gas density according to mass depletion do i=1,ngrid if(flag2(ind_cell(i))>0)then n=flag2(ind_cell(i)) d=uold(ind_cell(i),1) uold(ind_cell(i),1)=max(d-n*dstar,0.5*d) endif end do end do ! End loop over cells end do ! End loop over grids !--------------------------------------------------------- ! Convert hydro variables back to conservative variables !--------------------------------------------------------- ncache=active(ilevel)%ngrid do igrid=1,ncache,nvector ngrid=MIN(nvector,ncache-igrid+1) do i=1,ngrid ind_grid(i)=active(ilevel)%igrid(igrid+i-1) end do do ind=1,twotondim iskip=ncoarse+(ind-1)*ngridmax do i=1,ngrid ind_cell(i)=iskip+ind_grid(i) end do do i=1,ngrid d=uold(ind_cell(i),1) u=uold(ind_cell(i),2) v=uold(ind_cell(i),3) w=uold(ind_cell(i),4) e=uold(ind_cell(i),5)*d #ifdef SOLVERmhd bx1=uold(ind_cell(i),6) by1=uold(ind_cell(i),7) bz1=uold(ind_cell(i),8) bx2=uold(ind_cell(i),nvar+1) by2=uold(ind_cell(i),nvar+2) bz2=uold(ind_cell(i),nvar+3) e=e+0.125d0*((bx1+bx2)**2+(by1+by2)**2+(bz1+bz2)**2) #endif e=e+0.5d0*d*(u**2+v**2+w**2) #if NENER>0 do irad=0,nener-1 e=e+uold(ind_cell(i),inener+irad) end do #endif uold(ind_cell(i),1)=d uold(ind_cell(i),2)=d*u uold(ind_cell(i),3)=d*v uold(ind_cell(i),4)=d*w uold(ind_cell(i),5)=e end do do ivar=imetal,nvar do i=1,ngrid d=uold(ind_cell(i),1) w=uold(ind_cell(i),ivar) uold(ind_cell(i),ivar)=d*w end do end do end do end do #endif if(sf_log_properties) close(ilun) end subroutine star_formation !################################################################ !################################################################ !################################################################ !################################################################ subroutine getnbor(ind_cell,ind_father,ncell,ilevel) use amr_commons implicit none integer::ncell,ilevel integer,dimension(1:nvector)::ind_cell integer,dimension(1:nvector,0:twondim)::ind_father !----------------------------------------------------------------- ! This subroutine determines the 2*ndim neighboring cells ! cells of the input cell (ind_cell). ! If for some reasons they don't exist, the routine returns ! the input cell. !----------------------------------------------------------------- integer::i,j,iok,ind integer,dimension(1:nvector),save::ind_grid_father,pos integer,dimension(1:nvector,0:twondim),save::igridn,igridn_ok integer,dimension(1:nvector,1:twondim),save::icelln_ok if(ilevel==1)then write(*,*) 'Warning: attempting to form stars on level 1 --> this is not allowed ...' return endif ! Get father cell do i=1,ncell ind_father(i,0)=ind_cell(i) end do ! Get father cell position in the grid do i=1,ncell pos(i)=(ind_father(i,0)-ncoarse-1)/ngridmax+1 end do ! Get father grid do i=1,ncell ind_grid_father(i)=ind_father(i,0)-ncoarse-(pos(i)-1)*ngridmax end do ! Get neighboring father grids call getnborgrids(ind_grid_father,igridn,ncell) ! Loop over position do ind=1,twotondim ! Select father cells that sit at position ind do j=0,twondim iok=0 do i=1,ncell if(pos(i)==ind)then iok=iok+1 igridn_ok(iok,j)=igridn(i,j) end if end do end do ! Get neighboring cells for selected cells if(iok>0)call getnborcells(igridn_ok,ind,icelln_ok,iok) ! Update neighboring father cells for selected cells do j=1,twondim iok=0 do i=1,ncell if(pos(i)==ind)then iok=iok+1 if(icelln_ok(iok,j)>0)then ind_father(i,j)=icelln_ok(iok,j) else ind_father(i,j)=ind_cell(i) end if end if end do end do end do end subroutine getnbor !############################################################## !############################################################## !############################################################## !############################################################## function erfc(x) ! complementary error function use amr_commons, ONLY: dp implicit none real(dp) erfc real(dp) x, y real(kind=8) pv, ph real(kind=8) q0, q1, q2, q3, q4, q5, q6, q7 real(kind=8) p0, p1, p2, p3, p4, p5, p6, p7 parameter(pv= 1.26974899965115684d+01, ph= 6.10399733098688199d+00) parameter(p0= 2.96316885199227378d-01, p1= 1.81581125134637070d-01) parameter(p2= 6.81866451424939493d-02, p3= 1.56907543161966709d-02) parameter(p4= 2.21290116681517573d-03, p5= 1.91395813098742864d-04) parameter(p6= 9.71013284010551623d-06, p7= 1.66642447174307753d-07) parameter(q0= 6.12158644495538758d-02, q1= 5.50942780056002085d-01) parameter(q2= 1.53039662058770397d+00, q3= 2.99957952311300634d+00) parameter(q4= 4.95867777128246701d+00, q5= 7.41471251099335407d+00) parameter(q6= 1.04765104356545238d+01, q7= 1.48455557345597957d+01) y = x*x y = exp(-y)*x*(p7/(y+q7)+p6/(y+q6) + p5/(y+q5)+p4/(y+q4)+p3/(y+q3) & & + p2/(y+q2)+p1/(y+q1)+p0/(y+q0)) if (x < ph) y = y+2d0/(exp(pv*x)+1.0) erfc = y return end function erfc /home/hemv_g/srf_a/RAMSES_CODE_CORI/ramses06062019_grackle311Cosmo/patch/mom2/umuscl.f90 ! --------------------------------------------------------------- ! UNSPLIT Unsplit second order Godunov integrator for ! polytropic gas dynamics using either ! MUSCL-HANCOCK scheme or Collela's PLMDE scheme ! with various slope limiters. ! ! inputs/outputs ! uin => (const) input state ! gravin => (const) input gravitational acceleration ! iu1,iu2 => (const) first and last index of input array, ! ju1,ju2 => (const) cell centered, ! ku1,ku2 => (const) including buffer cells. ! flux <= (modify) return fluxes in the 3 coord directions ! if1,if2 => (const) first and last index of output array, ! jf1,jf2 => (const) edge centered, ! kf1,kf2 => (const) for active cells only. ! dx,dy,dz => (const) (dx,dy,dz) ! dt => (const) time step ! ngrid => (const) number of sub-grids ! ndim => (const) number of dimensions ! ---------------------------------------------------------------- subroutine unsplit(uin,pin,gravin,flux,tmp,dx,dy,dz,dt,ngrid) use amr_parameters use const use hydro_parameters implicit none integer ::ngrid real(dp)::dx,dy,dz,dt ! Input states real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2)::pin real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar)::uin real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:ndim)::gravin ! Output fluxes real(dp),dimension(1:nvector,if1:if2,jf1:jf2,kf1:kf2,1:nvar,1:ndim)::flux real(dp),dimension(1:nvector,if1:if2,jf1:jf2,kf1:kf2,1:2 ,1:ndim)::tmp ! Primitive variables real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar),save::qin real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2 ),save::cin ! Slopes real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim),save::dq ! Left and right state arrays real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim),save::qm real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim),save::qp ! Intermediate fluxes real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar),save::fx real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:2 ),save::tx ! Velocity divergence real(dp),dimension(1:nvector,if1:if2,jf1:jf2,kf1:kf2)::divu ! Local scalar variables integer::i,j,k,l,ivar integer::ilo,ihi,jlo,jhi,klo,khi ilo=MIN(1,iu1+2); ihi=MAX(1,iu2-2) jlo=MIN(1,ju1+2); jhi=MAX(1,ju2-2) klo=MIN(1,ku1+2); khi=MAX(1,ku2-2) ! Translate to primative variables, compute sound speeds call ctoprim(uin,qin,cin,gravin,dt,ngrid) ! Compute TVD slopes call uslope(qin,dq,dx,dt,ngrid) ! Compute 3D traced-states in all three directions if(scheme=='muscl')then #if NDIM==1 call trace1d(qin,dq,qm,qp,dx ,dt,ngrid) #endif #if NDIM==2 call trace2d(qin,dq,qm,qp,dx,dy ,dt,ngrid) #endif #if NDIM==3 call trace3d(qin,dq,qm,qp,dx,dy,dz,dt,ngrid) #endif endif if(scheme=='plmde')then #if NDIM==1 call tracex (qin,dq,cin,qm,qp,dx ,dt,ngrid) #endif #if NDIM==2 call tracexy (qin,dq,cin,qm,qp,dx,dy ,dt,ngrid) #endif #if NDIM==3 call tracexyz(qin,dq,cin,qm,qp,dx,dy,dz,dt,ngrid) #endif endif ! Solve for 1D flux in X direction call cmpflxm(qm,iu1+1,iu2+1,ju1 ,ju2 ,ku1 ,ku2 , & & qp,iu1 ,iu2 ,ju1 ,ju2 ,ku1 ,ku2 , pin, & & if1 ,if2 ,jlo ,jhi ,klo ,khi , 2,3,4,fx,tx,ngrid) ! Save flux in output array do i=if1,if2 do j=jlo,jhi do k=klo,khi do ivar=1,nvar do l=1,ngrid flux(l,i,j,k,ivar,1)=fx(l,i,j,k,ivar)*dt/dx end do end do do ivar=1,2 do l=1,ngrid tmp (l,i,j,k,ivar,1)=tx(l,i,j,k,ivar)*dt/dx end do end do end do end do end do ! Solve for 1D flux in Y direction #if NDIM>1 call cmpflxm(qm,iu1 ,iu2 ,ju1+1,ju2+1,ku1 ,ku2 , & & qp,iu1 ,iu2 ,ju1 ,ju2 ,ku1 ,ku2 , pin, & & ilo ,ihi ,jf1 ,jf2 ,klo ,khi , 3,2,4,fx,tx,ngrid) ! Save flux in output array do i=ilo,ihi do j=jf1,jf2 do k=klo,khi do ivar=1,nvar do l=1,ngrid flux(l,i,j,k,ivar,2)=fx(l,i,j,k,ivar)*dt/dy end do end do do ivar=1,2 do l=1,ngrid tmp (l,i,j,k,ivar,2)=tx(l,i,j,k,ivar)*dt/dy end do end do end do end do end do #endif ! Solve for 1D flux in Z direction #if NDIM>2 call cmpflxm(qm,iu1 ,iu2 ,ju1 ,ju2 ,ku1+1,ku2+1, & & qp,iu1 ,iu2 ,ju1 ,ju2 ,ku1 ,ku2 , pin, & & ilo ,ihi ,jlo ,jhi ,kf1 ,kf2 , 4,2,3,fx,tx,ngrid) ! Save flux in output array do i=ilo,ihi do j=jlo,jhi do k=kf1,kf2 do ivar=1,nvar do l=1,ngrid flux(l,i,j,k,ivar,3)=fx(l,i,j,k,ivar)*dt/dz end do end do do ivar=1,2 do l=1,ngrid tmp (l,i,j,k,ivar,3)=tx(l,i,j,k,ivar)*dt/dz end do end do end do end do end do #endif if(difmag>0.0)then call cmpdivu(qin,divu,dx,dy,dz,ngrid) call consup(uin,flux,divu,dt,ngrid) endif end subroutine unsplit !########################################################### !########################################################### !########################################################### !########################################################### subroutine trace1d(q,dq,qm,qp,dx,dt,ngrid) use amr_parameters use hydro_parameters use const implicit none integer ::ngrid real(dp)::dx, dt real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar)::q real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::dq real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::qm real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::qp ! Local variables integer ::i, j, k, l integer ::ilo,ihi,jlo,jhi,klo,khi integer ::ir, iu, ip real(dp)::dtdx real(dp)::r, u, p real(dp)::drx, dux, dpx real(dp)::sr0, su0, sp0 #if NENER>0 integer::irad real(dp),dimension(1:nener)::e, dex, se0 #endif #if NVAR > NDIM + 2 + NENER integer::n real(dp)::a, dax, sa0 #endif dtdx = dt/dx ilo=MIN(1,iu1+1); ihi=MAX(1,iu2-1) jlo=MIN(1,ju1+1); jhi=MAX(1,ju2-1) klo=MIN(1,ku1+1); khi=MAX(1,ku2-1) ir=1; iu=2; ip=3 do k = klo, khi do j = jlo, jhi do i = ilo, ihi do l = 1, ngrid ! Cell centered values r = q(l,i,j,k,ir) u = q(l,i,j,k,iu) p = q(l,i,j,k,ip) #if NENER>0 do irad=1,nener e(irad) = q(l,i,j,k,ip+irad) end do #endif ! TVD slopes in X direction drx = dq(l,i,j,k,ir,1) dux = dq(l,i,j,k,iu,1) dpx = dq(l,i,j,k,ip,1) #if NENER>0 do irad=1,nener dex(irad) = dq(l,i,j,k,ip+irad,1) end do #endif ! Source terms (including transverse derivatives) sr0 = -u*drx - (dux)*r sp0 = -u*dpx - (dux)*gamma*p su0 = -u*dux - (dpx)/r #if NENER>0 do irad=1,nener su0 = su0 - (dex(irad))/r se0(irad) = -u*dex(irad) & & - (dux)*gamma_rad(irad)*e(irad) end do #endif ! Right state qp(l,i,j,k,ir,1) = r - half*drx + sr0*dtdx*half qp(l,i,j,k,iu,1) = u - half*dux + su0*dtdx*half qp(l,i,j,k,ip,1) = p - half*dpx + sp0*dtdx*half ! qp(l,i,j,k,ir,1) = max(smallr, qp(l,i,j,k,ir,1)) if(qp(l,i,j,k,ir,1)0 do irad=1,nener qp(l,i,j,k,ip+irad,1) = e(irad) - half*dex(irad) + se0(irad)*dtdx*half end do #endif ! Left state qm(l,i,j,k,ir,1) = r + half*drx + sr0*dtdx*half qm(l,i,j,k,iu,1) = u + half*dux + su0*dtdx*half qm(l,i,j,k,ip,1) = p + half*dpx + sp0*dtdx*half ! qm(l,i,j,k,ir,1) = max(smallr, qm(l,i,j,k,ir,1)) if(qm(l,i,j,k,ir,1)0 do irad=1,nener qm(l,i,j,k,ip+irad,1) = e(irad) + half*dex(irad) + se0(irad)*dtdx*half end do #endif end do end do end do end do #if NVAR > NDIM + 2 + NENER ! Passive scalars do n = ndim+nener+3, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi do l = 1, ngrid a = q(l,i,j,k,n) ! Cell centered values u = q(l,i,j,k,iu) dax = dq(l,i,j,k,n,1) ! TVD slopes sa0 = -u*dax ! Source terms qp(l,i,j,k,n,1) = a - half*dax + sa0*dtdx*half ! Right state qm(l,i,j,k,n,1) = a + half*dax + sa0*dtdx*half ! Left state end do end do end do end do end do #endif end subroutine trace1d !########################################################### !########################################################### !########################################################### !########################################################### #if NDIM>1 subroutine trace2d(q,dq,qm,qp,dx,dy,dt,ngrid) use amr_parameters use hydro_parameters use const implicit none integer ::ngrid real(dp)::dx, dy, dt real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar)::q real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::dq real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::qm real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::qp ! declare local variables integer ::i, j, k, l integer ::ilo,ihi,jlo,jhi,klo,khi integer ::ir, iu, iv, ip real(dp)::dtdx, dtdy real(dp)::r, u, v, p real(dp)::drx, dux, dvx, dpx real(dp)::dry, duy, dvy, dpy real(dp)::sr0, su0, sv0, sp0 #if NENER>0 integer ::irad real(dp),dimension(1:nener)::e, dex, dey, se0 #endif #if NVAR > NDIM + 2 + NENER integer ::n real(dp)::a, dax, day, sa0 #endif dtdx = dt/dx dtdy = dt/dy ilo=MIN(1,iu1+1); ihi=MAX(1,iu2-1) jlo=MIN(1,ju1+1); jhi=MAX(1,ju2-1) klo=MIN(1,ku1+1); khi=MAX(1,ku2-1) ir=1; iu=2; iv=3; ip=4 do k = klo, khi do j = jlo, jhi do i = ilo, ihi do l = 1, ngrid ! Cell centered values r = q(l,i,j,k,ir) u = q(l,i,j,k,iu) v = q(l,i,j,k,iv) p = q(l,i,j,k,ip) #if NENER>0 do irad=1,nener e(irad) = q(l,i,j,k,ip+irad) end do #endif ! TVD slopes in all directions drx = dq(l,i,j,k,ir,1) dux = dq(l,i,j,k,iu,1) dvx = dq(l,i,j,k,iv,1) dpx = dq(l,i,j,k,ip,1) #if NENER>0 do irad=1,nener dex(irad) = dq(l,i,j,k,ip+irad,1) end do #endif dry = dq(l,i,j,k,ir,2) duy = dq(l,i,j,k,iu,2) dvy = dq(l,i,j,k,iv,2) dpy = dq(l,i,j,k,ip,2) #if NENER>0 do irad=1,nener dey(irad) = dq(l,i,j,k,ip+irad,2) end do #endif ! source terms (with transverse derivatives) sr0 = -u*drx-v*dry - (dux+dvy)*r sp0 = -u*dpx-v*dpy - (dux+dvy)*gamma*p su0 = -u*dux-v*duy - (dpx )/r sv0 = -u*dvx-v*dvy - (dpy )/r #if NENER>0 do irad=1,nener su0 = su0 - (dex(irad))/r sv0 = sv0 - (dey(irad))/r se0(irad) = -u*dex(irad)-v*dey(irad) & & - (dux+dvy)*gamma_rad(irad)*e(irad) end do #endif ! Right state at left interface qp(l,i,j,k,ir,1) = r - half*drx + sr0*dtdx*half qp(l,i,j,k,iu,1) = u - half*dux + su0*dtdx*half qp(l,i,j,k,iv,1) = v - half*dvx + sv0*dtdx*half qp(l,i,j,k,ip,1) = p - half*dpx + sp0*dtdx*half ! qp(l,i,j,k,ir,1) = max(smallr, qp(l,i,j,k,ir,1)) if(qp(l,i,j,k,ir,1)0 do irad=1,nener qp(l,i,j,k,ip+irad,1) = e(irad) - half*dex(irad) + se0(irad)*dtdx*half end do #endif ! Left state at right interface qm(l,i,j,k,ir,1) = r + half*drx + sr0*dtdx*half qm(l,i,j,k,iu,1) = u + half*dux + su0*dtdx*half qm(l,i,j,k,iv,1) = v + half*dvx + sv0*dtdx*half qm(l,i,j,k,ip,1) = p + half*dpx + sp0*dtdx*half ! qm(l,i,j,k,ir,1) = max(smallr, qm(l,i,j,k,ir,1)) if(qm(l,i,j,k,ir,1)0 do irad=1,nener qm(l,i,j,k,ip+irad,1) = e(irad) + half*dex(irad) + se0(irad)*dtdx*half end do #endif ! Top state at bottom interface qp(l,i,j,k,ir,2) = r - half*dry + sr0*dtdy*half qp(l,i,j,k,iu,2) = u - half*duy + su0*dtdy*half qp(l,i,j,k,iv,2) = v - half*dvy + sv0*dtdy*half qp(l,i,j,k,ip,2) = p - half*dpy + sp0*dtdy*half ! qp(l,i,j,k,ir,2) = max(smallr, qp(l,i,j,k,ir,2)) if(qp(l,i,j,k,ir,2)0 do irad=1,nener qp(l,i,j,k,ip+irad,2) = e(irad) - half*dey(irad) + se0(irad)*dtdy*half end do #endif ! Bottom state at top interface qm(l,i,j,k,ir,2) = r + half*dry + sr0*dtdy*half qm(l,i,j,k,iu,2) = u + half*duy + su0*dtdy*half qm(l,i,j,k,iv,2) = v + half*dvy + sv0*dtdy*half qm(l,i,j,k,ip,2) = p + half*dpy + sp0*dtdy*half ! qm(l,i,j,k,ir,2) = max(smallr, qm(l,i,j,k,ir,2)) if(qm(l,i,j,k,ir,2)0 do irad=1,nener qm(l,i,j,k,ip+irad,2) = e(irad) + half*dey(irad) + se0(irad)*dtdy*half end do #endif end do end do end do end do #if NVAR > NDIM + 2 + NENER ! passive scalars do n = ndim+nener+3, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi do l = 1, ngrid a = q(l,i,j,k,n) ! Cell centered values u = q(l,i,j,k,iu) v = q(l,i,j,k,iv) dax = dq(l,i,j,k,n,1) ! TVD slopes day = dq(l,i,j,k,n,2) sa0 = -u*dax-v*day ! Source terms qp(l,i,j,k,n,1) = a - half*dax + sa0*dtdx*half ! Right state qm(l,i,j,k,n,1) = a + half*dax + sa0*dtdx*half ! Left state qp(l,i,j,k,n,2) = a - half*day + sa0*dtdy*half ! Top state qm(l,i,j,k,n,2) = a + half*day + sa0*dtdy*half ! Bottom state end do end do end do end do end do #endif end subroutine trace2d #endif !########################################################### !########################################################### !########################################################### !########################################################### #if NDIM>2 subroutine trace3d(q,dq,qm,qp,dx,dy,dz,dt,ngrid) use amr_parameters use hydro_parameters use const implicit none integer ::ngrid real(dp)::dx, dy, dz, dt real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar)::q real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::dq real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::qm real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::qp ! declare local variables integer ::i, j, k, l integer ::ilo,ihi,jlo,jhi,klo,khi integer ::ir, iu, iv, iw, ip real(dp)::dtdx, dtdy, dtdz real(dp)::r, u, v, w, p real(dp)::drx, dux, dvx, dwx, dpx real(dp)::dry, duy, dvy, dwy, dpy real(dp)::drz, duz, dvz, dwz, dpz real(dp)::sr0, su0, sv0, sw0, sp0 #if NENER>0 integer ::irad real(dp),dimension(1:nener)::e, dex, dey, dez, se0 #endif #if NVAR > NDIM + 2 + NENER integer ::n real(dp)::a, dax, day, daz, sa0 #endif dtdx = dt/dx dtdy = dt/dy dtdz = dt/dz ilo=MIN(1,iu1+1); ihi=MAX(1,iu2-1) jlo=MIN(1,ju1+1); jhi=MAX(1,ju2-1) klo=MIN(1,ku1+1); khi=MAX(1,ku2-1) ir=1; iu=2; iv=3; iw=4; ip=5 do k = klo, khi do j = jlo, jhi do i = ilo, ihi do l = 1, ngrid ! Cell centered values r = q(l,i,j,k,ir) u = q(l,i,j,k,iu) v = q(l,i,j,k,iv) w = q(l,i,j,k,iw) p = q(l,i,j,k,ip) #if NENER>0 do irad=1,nener e(irad) = q(l,i,j,k,ip+irad) end do #endif ! TVD slopes in all 3 directions drx = dq(l,i,j,k,ir,1) dpx = dq(l,i,j,k,ip,1) dux = dq(l,i,j,k,iu,1) dvx = dq(l,i,j,k,iv,1) dwx = dq(l,i,j,k,iw,1) #if NENER>0 do irad=1,nener dex(irad) = dq(l,i,j,k,ip+irad,1) end do #endif dry = dq(l,i,j,k,ir,2) dpy = dq(l,i,j,k,ip,2) duy = dq(l,i,j,k,iu,2) dvy = dq(l,i,j,k,iv,2) dwy = dq(l,i,j,k,iw,2) #if NENER>0 do irad=1,nener dey(irad) = dq(l,i,j,k,ip+irad,2) end do #endif drz = dq(l,i,j,k,ir,3) dpz = dq(l,i,j,k,ip,3) duz = dq(l,i,j,k,iu,3) dvz = dq(l,i,j,k,iv,3) dwz = dq(l,i,j,k,iw,3) #if NENER>0 do irad=1,nener dez(irad) = dq(l,i,j,k,ip+irad,3) end do #endif ! Source terms (including transverse derivatives) sr0 = -u*drx-v*dry-w*drz - (dux+dvy+dwz)*r sp0 = -u*dpx-v*dpy-w*dpz - (dux+dvy+dwz)*gamma*p su0 = -u*dux-v*duy-w*duz - (dpx )/r sv0 = -u*dvx-v*dvy-w*dvz - (dpy )/r sw0 = -u*dwx-v*dwy-w*dwz - (dpz )/r #if NENER>0 do irad=1,nener su0 = su0 - (dex(irad))/r sv0 = sv0 - (dey(irad))/r sw0 = sw0 - (dez(irad))/r se0(irad) = -u*dex(irad)-v*dey(irad)-w*dez(irad) & & - (dux+dvy+dwz)*gamma_rad(irad)*e(irad) end do #endif ! Right state at left interface qp(l,i,j,k,ir,1) = r - half*drx + sr0*dtdx*half qp(l,i,j,k,ip,1) = p - half*dpx + sp0*dtdx*half qp(l,i,j,k,iu,1) = u - half*dux + su0*dtdx*half qp(l,i,j,k,iv,1) = v - half*dvx + sv0*dtdx*half qp(l,i,j,k,iw,1) = w - half*dwx + sw0*dtdx*half ! qp(l,i,j,k,ir,1) = max(smallr, qp(l,i,j,k,ir,1)) if(qp(l,i,j,k,ir,1)0 do irad=1,nener qp(l,i,j,k,ip+irad,1) = e(irad) - half*dex(irad) + se0(irad)*dtdx*half end do #endif ! Left state at left interface qm(l,i,j,k,ir,1) = r + half*drx + sr0*dtdx*half qm(l,i,j,k,ip,1) = p + half*dpx + sp0*dtdx*half qm(l,i,j,k,iu,1) = u + half*dux + su0*dtdx*half qm(l,i,j,k,iv,1) = v + half*dvx + sv0*dtdx*half qm(l,i,j,k,iw,1) = w + half*dwx + sw0*dtdx*half ! qm(l,i,j,k,ir,1) = max(smallr, qm(l,i,j,k,ir,1)) if(qm(l,i,j,k,ir,1)0 do irad=1,nener qm(l,i,j,k,ip+irad,1) = e(irad) + half*dex(irad) + se0(irad)*dtdx*half end do #endif ! Top state at bottom interface qp(l,i,j,k,ir,2) = r - half*dry + sr0*dtdy*half qp(l,i,j,k,ip,2) = p - half*dpy + sp0*dtdy*half qp(l,i,j,k,iu,2) = u - half*duy + su0*dtdy*half qp(l,i,j,k,iv,2) = v - half*dvy + sv0*dtdy*half qp(l,i,j,k,iw,2) = w - half*dwy + sw0*dtdy*half ! qp(l,i,j,k,ir,2) = max(smallr, qp(l,i,j,k,ir,2)) if(qp(l,i,j,k,ir,2)0 do irad=1,nener qp(l,i,j,k,ip+irad,2) = e(irad) - half*dey(irad) + se0(irad)*dtdy*half end do #endif ! Bottom state at top interface qm(l,i,j,k,ir,2) = r + half*dry + sr0*dtdy*half qm(l,i,j,k,ip,2) = p + half*dpy + sp0*dtdy*half qm(l,i,j,k,iu,2) = u + half*duy + su0*dtdy*half qm(l,i,j,k,iv,2) = v + half*dvy + sv0*dtdy*half qm(l,i,j,k,iw,2) = w + half*dwy + sw0*dtdy*half ! qm(l,i,j,k,ir,2) = max(smallr, qm(l,i,j,k,ir,2)) if(qm(l,i,j,k,ir,2)0 do irad=1,nener qm(l,i,j,k,ip+irad,2) = e(irad) + half*dey(irad) + se0(irad)*dtdy*half end do #endif ! Back state at front interface qp(l,i,j,k,ir,3) = r - half*drz + sr0*dtdz*half qp(l,i,j,k,ip,3) = p - half*dpz + sp0*dtdz*half qp(l,i,j,k,iu,3) = u - half*duz + su0*dtdz*half qp(l,i,j,k,iv,3) = v - half*dvz + sv0*dtdz*half qp(l,i,j,k,iw,3) = w - half*dwz + sw0*dtdz*half ! qp(l,i,j,k,ir,3) = max(smallr, qp(l,i,j,k,ir,3)) if(qp(l,i,j,k,ir,3)0 do irad=1,nener qp(l,i,j,k,ip+irad,3) = e(irad) - half*dez(irad) + se0(irad)*dtdz*half end do #endif ! Front state at back interface qm(l,i,j,k,ir,3) = r + half*drz + sr0*dtdz*half qm(l,i,j,k,ip,3) = p + half*dpz + sp0*dtdz*half qm(l,i,j,k,iu,3) = u + half*duz + su0*dtdz*half qm(l,i,j,k,iv,3) = v + half*dvz + sv0*dtdz*half qm(l,i,j,k,iw,3) = w + half*dwz + sw0*dtdz*half ! qm(l,i,j,k,ir,3) = max(smallr, qm(l,i,j,k,ir,3)) if(qm(l,i,j,k,ir,3)0 do irad=1,nener qm(l,i,j,k,ip+irad,3) = e(irad) + half*dez(irad) + se0(irad)*dtdz*half end do #endif end do end do end do end do #if NVAR > NDIM + 2 + NENER ! Passive scalars do n = ndim+nener+3, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi do l = 1, ngrid a = q(l,i,j,k,n) ! Cell centered values u = q(l,i,j,k,iu) v = q(l,i,j,k,iv) w = q(l,i,j,k,iw) dax = dq(l,i,j,k,n,1) ! TVD slopes day = dq(l,i,j,k,n,2) daz = dq(l,i,j,k,n,3) sa0 = -u*dax-v*day-w*daz ! Source terms qp(l,i,j,k,n,1) = a - half*dax + sa0*dtdx*half ! Right state qm(l,i,j,k,n,1) = a + half*dax + sa0*dtdx*half ! Left state qp(l,i,j,k,n,2) = a - half*day + sa0*dtdy*half ! Bottom state qm(l,i,j,k,n,2) = a + half*day + sa0*dtdy*half ! Upper state qp(l,i,j,k,n,3) = a - half*daz + sa0*dtdz*half ! Front state qm(l,i,j,k,n,3) = a + half*daz + sa0*dtdz*half ! Back state end do end do end do end do end do #endif end subroutine trace3d #endif !########################################################### !########################################################### !########################################################### !########################################################### subroutine cmpflxm(qm,im1,im2,jm1,jm2,km1,km2, & & qp,ip1,ip2,jp1,jp2,kp1,kp2, pin, & & ilo,ihi,jlo,jhi,klo,khi, ln,lt1,lt2, & & flx,tmp,ngrid) use amr_parameters use hydro_parameters use const implicit none integer ::ngrid integer ::ln,lt1,lt2 integer ::im1,im2,jm1,jm2,km1,km2 integer ::ip1,ip2,jp1,jp2,kp1,kp2 integer ::ilo,ihi,jlo,jhi,klo,khi real(dp),dimension(1:nvector,im1:im2,jm1:jm2,km1:km2,1:nvar,1:ndim)::qm real(dp),dimension(1:nvector,ip1:ip2,jp1:jp2,kp1:kp2,1:nvar,1:ndim)::qp real(dp),dimension(1:nvector,ip1:ip2,jp1:jp2,kp1:kp2)::pin real(dp),dimension(1:nvector,ip1:ip2,jp1:jp2,kp1:kp2,1:nvar)::flx real(dp),dimension(1:nvector,ip1:ip2,jp1:jp2,kp1:kp2,1:2)::tmp ! local variables integer ::i, j, k, l, xdim real(dp)::entho real(dp),dimension(1:nvector),save::snleft,snright real(dp),dimension(1:nvector,1:nvar),save::qleft,qright real(dp),dimension(1:nvector,1:nvar+1),save::fgdnv #if NVAR > NDIM + 2 integer ::n #endif entho=one/(gamma-one) xdim=ln-1 do k = klo, khi do j = jlo, jhi do i = ilo, ihi ! Mass density do l = 1, ngrid qleft (l,1) = qm(l,i,j,k,1,xdim) qright(l,1) = qp(l,i,j,k,1,xdim) end do ! Normal velocity do l = 1, ngrid qleft (l,2) = qm(l,i,j,k,ln,xdim) qright(l,2) = qp(l,i,j,k,ln,xdim) end do ! Pressure do l = 1, ngrid qleft (l,3) = qm(l,i,j,k,ndim+2,xdim) qright(l,3) = qp(l,i,j,k,ndim+2,xdim) end do ! Supernovae if(ln==2)then do l = 1, ngrid snleft (l) = pin(l,i-1,j,k) snright(l) = pin(l,i,j,k) end do endif if(ln==3)then do l = 1, ngrid snleft (l) = pin(l,i,j-1,k) snright(l) = pin(l,i,j,k) end do endif if(ln==4)then do l = 1, ngrid snleft (l) = pin(l,i,j,k-1) snright(l) = pin(l,i,j,k) end do endif ! Tangential velocity 1 #if NDIM>1 do l = 1, ngrid qleft (l,4) = qm(l,i,j,k,lt1,xdim) qright(l,4) = qp(l,i,j,k,lt1,xdim) end do #endif ! Tangential velocity 2 #if NDIM>2 do l = 1, ngrid qleft (l,5) = qm(l,i,j,k,lt2,xdim) qright(l,5) = qp(l,i,j,k,lt2,xdim) end do #endif #if NVAR > NDIM + 2 ! Other advected quantities do n = ndim+3, nvar do l = 1, ngrid qleft (l,n) = qm(l,i,j,k,n,xdim) qright(l,n) = qp(l,i,j,k,n,xdim) end do end do #endif ! Solve Riemann problem if(riemann.eq.'acoustic')then call riemann_acoustic(qleft,qright,fgdnv,ngrid) else if (riemann.eq.'exact')then call riemann_approx (qleft,qright,fgdnv,ngrid) else if (riemann.eq.'llf')then call riemann_llf (qleft,qright,fgdnv,ngrid) else if (riemann.eq.'hllc')then call riemann_hllc (qleft,qright,snleft,snright,fgdnv,ngrid) else if (riemann.eq.'hll')then call riemann_hll (qleft,qright,fgdnv,ngrid) else write(*,*)'unknown Riemann solver' stop end if ! Compute fluxes ! Mass density do l = 1, ngrid flx(l,i,j,k,1) = fgdnv(l,1) end do ! Normal momentum do l = 1, ngrid flx(l,i,j,k,ln) = fgdnv(l,2) end do ! Transverse momentum 1 #if NDIM>1 do l = 1, ngrid flx(l,i,j,k,lt1) = fgdnv(l,4) end do #endif ! Transverse momentum 2 #if NDIM>2 do l = 1, ngrid flx(l,i,j,k,lt2) = fgdnv(l,5) end do #endif ! Total energy do l = 1, ngrid flx(l,i,j,k,ndim+2) = fgdnv(l,3) end do #if NVAR > NDIM + 2 ! Other advected quantities do n = ndim+3, nvar do l = 1, ngrid flx(l,i,j,k,n) = fgdnv(l,n) end do end do #endif ! Normal velocity do l = 1, ngrid tmp(l,i,j,k,1) = half*(qleft(l,2)+qright(l,2)) end do ! Internal energy flux do l = 1,ngrid tmp(l,i,j,k,2) = fgdnv(l,nvar+1) end do end do end do end do end subroutine cmpflxm !########################################################### !########################################################### !########################################################### !########################################################### subroutine ctoprim(uin,q,c,gravin,dt,ngrid) use amr_parameters use hydro_parameters use const implicit none integer ::ngrid real(dp)::dt real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar)::uin real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:ndim)::gravin real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar)::q real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2)::c integer ::i, j, k, l real(dp)::eint, smalle, dtxhalf, oneoverrho real(dp)::eken, erad #if NVAR > NDIM + 2 + NENER integer ::n #endif #if NENER>0 integer ::irad #endif smalle = smallc**2/gamma/(gamma-one) dtxhalf = dt*half ! Convert to primitive variable do k = ku1, ku2 do j = ju1, ju2 do i = iu1, iu2 do l = 1, ngrid ! Compute density q(l,i,j,k,1) = max(uin(l,i,j,k,1),smallr) ! Compute velocities oneoverrho = one/q(l,i,j,k,1) q(l,i,j,k,2) = uin(l,i,j,k,2)*oneoverrho #if NDIM>1 q(l,i,j,k,3) = uin(l,i,j,k,3)*oneoverrho #endif #if NDIM>2 q(l,i,j,k,4) = uin(l,i,j,k,4)*oneoverrho #endif ! Compute specific kinetic energy eken = half*q(l,i,j,k,2)*q(l,i,j,k,2) #if NDIM>1 eken = eken + half*q(l,i,j,k,3)*q(l,i,j,k,3) #endif #if NDIM>2 eken = eken + half*q(l,i,j,k,4)*q(l,i,j,k,4) #endif ! Compute non-thermal pressure erad = zero #if NENER>0 do irad = 1,nener q(l,i,j,k,ndim+2+irad) = (gamma_rad(irad)-one)*uin(l,i,j,k,ndim+2+irad) erad = erad+uin(l,i,j,k,ndim+2+irad)*oneoverrho enddo #endif ! Compute thermal pressure eint = MAX(uin(l,i,j,k,ndim+2)*oneoverrho-eken-erad,smalle) q(l,i,j,k,ndim+2) = (gamma-one)*q(l,i,j,k,1)*eint ! Compute sound speed c(l,i,j,k)=gamma*q(l,i,j,k,ndim+2) #if NENER>0 do irad=1,nener c(l,i,j,k)=c(l,i,j,k)+gamma_rad(irad)*q(l,i,j,k,ndim+2+irad) enddo #endif c(l,i,j,k)=sqrt(c(l,i,j,k)*oneoverrho) ! Gravity predictor step q(l,i,j,k,2) = q(l,i,j,k,2) + gravin(l,i,j,k,1)*dtxhalf #if NDIM>1 q(l,i,j,k,3) = q(l,i,j,k,3) + gravin(l,i,j,k,2)*dtxhalf #endif #if NDIM>2 q(l,i,j,k,4) = q(l,i,j,k,4) + gravin(l,i,j,k,3)*dtxhalf #endif end do end do end do end do #if NVAR > NDIM + 2 + NENER ! Passive scalar do n = ndim+nener+3, nvar do k = ku1, ku2 do j = ju1, ju2 do i = iu1, iu2 do l = 1, ngrid oneoverrho = one/q(l,i,j,k,1) q(l,i,j,k,n) = uin(l,i,j,k,n)*oneoverrho end do end do end do end do end do #endif end subroutine ctoprim !########################################################### !########################################################### !########################################################### !########################################################### subroutine uslope(q,dq,dx,dt,ngrid) use amr_parameters use hydro_parameters use const implicit none integer::ngrid real(dp)::dx,dt real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar)::q real(dp),dimension(1:nvector,iu1:iu2,ju1:ju2,ku1:ku2,1:nvar,1:ndim)::dq ! local arrays integer::i, j, k, l, n real(dp)::dsgn, dlim, dcen, dlft, drgt, slop #if NDIM==2 real(dp)::dfll,dflm,dflr,dfml,dfmm,dfmr,dfrl,dfrm,dfrr #endif #if NDIM==3 real(dp)::dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl real(dp)::dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm real(dp)::dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr real(dp)::dfz #endif #if NDIM>1 real(dp)::vmin,vmax,dfx,dfy,dff #endif integer::ilo,ihi,jlo,jhi,klo,khi ilo=MIN(1,iu1+1); ihi=MAX(1,iu2-1) jlo=MIN(1,ju1+1); jhi=MAX(1,ju2-1) klo=MIN(1,ku1+1); khi=MAX(1,ku2-1) if(slope_type==0)then dq=zero return end if #if NDIM==1 do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi if(slope_type==1.or.slope_type==2.or.slope_type==3)then ! minmod or average do l = 1, ngrid dlft = MIN(slope_type,2)*(q(l,i ,j,k,n) - q(l,i-1,j,k,n)) drgt = MIN(slope_type,2)*(q(l,i+1,j,k,n) - q(l,i ,j,k,n)) dcen = half*(dlft+drgt)/MIN(slope_type,2) dsgn = sign(one, dcen) slop = min(abs(dlft),abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,1) = dsgn*min(dlim,abs(dcen)) end do else if(slope_type==4)then ! superbee do l = 1, ngrid dcen = q(l,i,j,k,2)*dt/dx dlft = two/(one+dcen)*(q(l,i,j,k,n)-q(l,i-1,j,k,n)) drgt = two/(one-dcen)*(q(l,i+1,j,k,n)-q(l,i,j,k,n)) dcen = half*(q(l,i+1,j,k,n)-q(l,i-1,j,k,n)) dsgn = sign(one, dlft) slop = min(abs(dlft),abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,1) = dsgn*dlim !min(dlim,abs(dcen)) end do else if(slope_type==5)then ! ultrabee if(n==1)then do l = 1, ngrid dcen = q(l,i,j,k,2)*dt/dx if(dcen>=0)then dlft = two/(zero+dcen+1d-10)*(q(l,i,j,k,n)-q(l,i-1,j,k,n)) drgt = two/(one -dcen )*(q(l,i+1,j,k,n)-q(l,i,j,k,n)) else dlft = two/(one +dcen )*(q(l,i,j,k,n)-q(l,i-1,j,k,n)) drgt = two/(zero-dcen+1d-10)*(q(l,i+1,j,k,n)-q(l,i,j,k,n)) endif dsgn = sign(one, dlft) slop = min(abs(dlft),abs(drgt)) dlim = slop dcen = half*(q(l,i+1,j,k,n)-q(l,i-1,j,k,n)) if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,1) = dsgn*dlim !min(dlim,abs(dcen)) end do else do l = 1, ngrid dq(l,i,j,k,n,1) = 0 end do end if else if(slope_type==6)then ! unstable if(n==1)then do l = 1, ngrid dlft = (q(l,i,j,k,n)-q(l,i-1,j,k,n)) drgt = (q(l,i+1,j,k,n)-q(l,i,j,k,n)) slop = 0.5d0*(dlft+drgt) dlim = slop dq(l,i,j,k,n,1) = dlim end do else do l = 1, ngrid dq(l,i,j,k,n,1) = 0 end do end if else if(slope_type==7)then ! van Leer do l = 1, ngrid dlft = (q(l,i ,j,k,n) - q(l,i-1,j,k,n)) drgt = (q(l,i+1,j,k,n) - q(l,i ,j,k,n)) if((dlft*drgt)<=zero) then dq(l,i,j,k,n,1)=zero else dq(l,i,j,k,n,1)=(2*dlft*drgt/(dlft+drgt)) end if end do else if(slope_type==8)then ! generalized moncen/minmod parameterisation (van Leer 1979) do l = 1, ngrid dlft = (q(l,i ,j,k,n) - q(l,i-1,j,k,n)) drgt = (q(l,i+1,j,k,n) - q(l,i ,j,k,n)) dcen = half*(dlft+drgt) dsgn = sign(one, dcen) slop = min(slope_theta*abs(dlft),slope_theta*abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,1) = dsgn*min(dlim,abs(dcen)) end do else write(*,*)'Unknown slope type' stop end if end do end do end do end do #endif #if NDIM==2 if(slope_type==1.or.slope_type==2)then ! minmod or average do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi ! slopes in first coordinate direction do l = 1, ngrid dlft = slope_type*(q(l,i ,j,k,n) - q(l,i-1,j,k,n)) drgt = slope_type*(q(l,i+1,j,k,n) - q(l,i ,j,k,n)) dcen = half*(dlft+drgt)/slope_type dsgn = sign(one, dcen) slop = min(abs(dlft),abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,1) = dsgn*min(dlim,abs(dcen)) end do ! slopes in second coordinate direction do l = 1, ngrid dlft = slope_type*(q(l,i,j ,k,n) - q(l,i,j-1,k,n)) drgt = slope_type*(q(l,i,j+1,k,n) - q(l,i,j ,k,n)) dcen = half*(dlft+drgt)/slope_type dsgn = sign(one,dcen) slop = min(abs(dlft),abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,2) = dsgn*min(dlim,abs(dcen)) end do end do end do end do end do else if(slope_type==3)then ! positivity preserving 2d unsplit slope do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi do l = 1, ngrid dfll = q(l,i-1,j-1,k,n)-q(l,i,j,k,n) dflm = q(l,i-1,j ,k,n)-q(l,i,j,k,n) dflr = q(l,i-1,j+1,k,n)-q(l,i,j,k,n) dfml = q(l,i ,j-1,k,n)-q(l,i,j,k,n) dfmm = q(l,i ,j ,k,n)-q(l,i,j,k,n) dfmr = q(l,i ,j+1,k,n)-q(l,i,j,k,n) dfrl = q(l,i+1,j-1,k,n)-q(l,i,j,k,n) dfrm = q(l,i+1,j ,k,n)-q(l,i,j,k,n) dfrr = q(l,i+1,j+1,k,n)-q(l,i,j,k,n) vmin = min(dfll,dflm,dflr,dfml,dfmm,dfmr,dfrl,dfrm,dfrr) vmax = max(dfll,dflm,dflr,dfml,dfmm,dfmr,dfrl,dfrm,dfrr) dfx = half*(q(l,i+1,j,k,n)-q(l,i-1,j,k,n)) dfy = half*(q(l,i,j+1,k,n)-q(l,i,j-1,k,n)) dff = half*(abs(dfx)+abs(dfy)) if(dff>zero)then slop = min(one,min(abs(vmin),abs(vmax))/dff) else slop = one endif dlim = slop dq(l,i,j,k,n,1) = dlim*dfx dq(l,i,j,k,n,2) = dlim*dfy end do end do end do end do end do else if(slope_type==7)then ! van Leer do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi ! slopes in first coordinate direction do l = 1, ngrid dlft = (q(l,i ,j,k,n) - q(l,i-1,j,k,n)) drgt = (q(l,i+1,j,k,n) - q(l,i ,j,k,n)) if((dlft*drgt)<=zero) then dq(l,i,j,k,n,1)=zero else dq(l,i,j,k,n,1)=(2*dlft*drgt/(dlft+drgt)) end if end do ! slopes in second coordinate direction do l = 1, ngrid dlft = (q(l,i,j ,k,n) - q(l,i,j-1,k,n)) drgt = (q(l,i,j+1,k,n) - q(l,i,j ,k,n)) if((dlft*drgt)<=zero) then dq(l,i,j,k,n,2)=zero else dq(l,i,j,k,n,2)=(2*dlft*drgt/(dlft+drgt)) end if end do end do end do end do end do else if(slope_type==8)then ! generalized moncen/minmod parameterisation (van Leer 1979) do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi ! slopes in first coordinate direction do l = 1, ngrid dlft = (q(l,i ,j,k,n) - q(l,i-1,j,k,n)) drgt = (q(l,i+1,j,k,n) - q(l,i ,j,k,n)) dcen = half*(dlft+drgt) dsgn = sign(one, dcen) slop = min(slope_theta*abs(dlft),slope_theta*abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,1) = dsgn*min(dlim,abs(dcen)) end do ! slopes in second coordinate direction do l = 1, ngrid dlft = (q(l,i,j ,k,n) - q(l,i,j-1,k,n)) drgt = (q(l,i,j+1,k,n) - q(l,i,j ,k,n)) dcen = half*(dlft+drgt) dsgn = sign(one,dcen) slop = min(slope_theta*abs(dlft),slope_theta*abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,2) = dsgn*min(dlim,abs(dcen)) end do end do end do end do end do else write(*,*)'Unknown slope type' stop endif #endif #if NDIM==3 if(slope_type==1)then ! minmod do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi ! slopes in first coordinate direction do l = 1, ngrid dlft = q(l,i ,j,k,n) - q(l,i-1,j,k,n) drgt = q(l,i+1,j,k,n) - q(l,i ,j,k,n) if((dlft*drgt)<=zero) then dq(l,i,j,k,n,1) = zero else if(dlft>0) then dq(l,i,j,k,n,1) = min(dlft,drgt) else dq(l,i,j,k,n,1) = max(dlft,drgt) end if end do ! slopes in second coordinate direction do l = 1, ngrid dlft = q(l,i,j ,k,n) - q(l,i,j-1,k,n) drgt = q(l,i,j+1,k,n) - q(l,i,j ,k,n) if((dlft*drgt)<=zero) then dq(l,i,j,k,n,2) = zero else if(dlft>0) then dq(l,i,j,k,n,2) = min(dlft,drgt) else dq(l,i,j,k,n,2) = max(dlft,drgt) end if end do ! slopes in third coordinate direction do l = 1, ngrid dlft = q(l,i,j,k ,n) - q(l,i,j,k-1,n) drgt = q(l,i,j,k+1,n) - q(l,i,j,k ,n) if((dlft*drgt)<=zero) then dq(l,i,j,k,n,3) = zero else if(dlft>0) then dq(l,i,j,k,n,3) = min(dlft,drgt) else dq(l,i,j,k,n,3) = max(dlft,drgt) end if end do end do end do end do end do else if(slope_type==2)then ! moncen do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi ! slopes in first coordinate direction do l = 1, ngrid dlft = slope_type*(q(l,i ,j,k,n) - q(l,i-1,j,k,n)) drgt = slope_type*(q(l,i+1,j,k,n) - q(l,i ,j,k,n)) dcen = half*(dlft+drgt)/slope_type dsgn = sign(one, dcen) slop = min(abs(dlft),abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,1) = dsgn*min(dlim,abs(dcen)) end do ! slopes in second coordinate direction do l = 1, ngrid dlft = slope_type*(q(l,i,j ,k,n) - q(l,i,j-1,k,n)) drgt = slope_type*(q(l,i,j+1,k,n) - q(l,i,j ,k,n)) dcen = half*(dlft+drgt)/slope_type dsgn = sign(one,dcen) slop = min(abs(dlft),abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,2) = dsgn*min(dlim,abs(dcen)) end do ! slopes in third coordinate direction do l = 1, ngrid dlft = slope_type*(q(l,i,j,k ,n) - q(l,i,j,k-1,n)) drgt = slope_type*(q(l,i,j,k+1,n) - q(l,i,j,k ,n)) dcen = half*(dlft+drgt)/slope_type dsgn = sign(one,dcen) slop = min(abs(dlft),abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,3) = dsgn*min(dlim,abs(dcen)) end do end do end do end do end do else if(slope_type==3)then ! positivity preserving 3d unsplit slope do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi do l = 1, ngrid dflll = q(l,i-1,j-1,k-1,n)-q(l,i,j,k,n) dflml = q(l,i-1,j ,k-1,n)-q(l,i,j,k,n) dflrl = q(l,i-1,j+1,k-1,n)-q(l,i,j,k,n) dfmll = q(l,i ,j-1,k-1,n)-q(l,i,j,k,n) dfmml = q(l,i ,j ,k-1,n)-q(l,i,j,k,n) dfmrl = q(l,i ,j+1,k-1,n)-q(l,i,j,k,n) dfrll = q(l,i+1,j-1,k-1,n)-q(l,i,j,k,n) dfrml = q(l,i+1,j ,k-1,n)-q(l,i,j,k,n) dfrrl = q(l,i+1,j+1,k-1,n)-q(l,i,j,k,n) dfllm = q(l,i-1,j-1,k ,n)-q(l,i,j,k,n) dflmm = q(l,i-1,j ,k ,n)-q(l,i,j,k,n) dflrm = q(l,i-1,j+1,k ,n)-q(l,i,j,k,n) dfmlm = q(l,i ,j-1,k ,n)-q(l,i,j,k,n) dfmmm = q(l,i ,j ,k ,n)-q(l,i,j,k,n) dfmrm = q(l,i ,j+1,k ,n)-q(l,i,j,k,n) dfrlm = q(l,i+1,j-1,k ,n)-q(l,i,j,k,n) dfrmm = q(l,i+1,j ,k ,n)-q(l,i,j,k,n) dfrrm = q(l,i+1,j+1,k ,n)-q(l,i,j,k,n) dfllr = q(l,i-1,j-1,k+1,n)-q(l,i,j,k,n) dflmr = q(l,i-1,j ,k+1,n)-q(l,i,j,k,n) dflrr = q(l,i-1,j+1,k+1,n)-q(l,i,j,k,n) dfmlr = q(l,i ,j-1,k+1,n)-q(l,i,j,k,n) dfmmr = q(l,i ,j ,k+1,n)-q(l,i,j,k,n) dfmrr = q(l,i ,j+1,k+1,n)-q(l,i,j,k,n) dfrlr = q(l,i+1,j-1,k+1,n)-q(l,i,j,k,n) dfrmr = q(l,i+1,j ,k+1,n)-q(l,i,j,k,n) dfrrr = q(l,i+1,j+1,k+1,n)-q(l,i,j,k,n) vmin = min(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, & & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, & & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr) vmax = max(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, & & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, & & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr) dfx = half*(q(l,i+1,j,k,n)-q(l,i-1,j,k,n)) dfy = half*(q(l,i,j+1,k,n)-q(l,i,j-1,k,n)) dfz = half*(q(l,i,j,k+1,n)-q(l,i,j,k-1,n)) dff = half*(abs(dfx)+abs(dfy)+abs(dfz)) if(dff>zero)then slop = min(one,min(abs(vmin),abs(vmax))/dff) else slop = one endif dlim = slop dq(l,i,j,k,n,1) = dlim*dfx dq(l,i,j,k,n,2) = dlim*dfy dq(l,i,j,k,n,3) = dlim*dfz end do end do end do end do end do else if(slope_type==7)then ! van Leer do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi ! slopes in first coordinate direction do l = 1, ngrid dlft = (q(l,i ,j,k,n) - q(l,i-1,j,k,n)) drgt = (q(l,i+1,j,k,n) - q(l,i ,j,k,n)) if((dlft*drgt)<=zero) then dq(l,i,j,k,n,1)=zero else dq(l,i,j,k,n,1)=(2*dlft*drgt/(dlft+drgt)) end if end do ! slopes in second coordinate direction do l = 1, ngrid dlft = (q(l,i,j ,k,n) - q(l,i,j-1,k,n)) drgt = (q(l,i,j+1,k,n) - q(l,i,j ,k,n)) if((dlft*drgt)<=zero) then dq(l,i,j,k,n,2)=zero else dq(l,i,j,k,n,2)=(2*dlft*drgt/(dlft+drgt)) end if end do ! slopes in third coordinate direction do l = 1, ngrid dlft = (q(l,i,j,k ,n) - q(l,i,j,k-1,n)) drgt = (q(l,i,j,k+1,n) - q(l,i,j,k ,n)) if((dlft*drgt)<=zero) then dq(l,i,j,k,n,3)=zero else dq(l,i,j,k,n,3)=(2*dlft*drgt/(dlft+drgt)) end if end do end do end do end do end do else if(slope_type==8)then ! generalized moncen/minmod parameterisation (van Leer 1979) do n = 1, nvar do k = klo, khi do j = jlo, jhi do i = ilo, ihi ! slopes in first coordinate direction do l = 1, ngrid dlft = (q(l,i ,j,k,n) - q(l,i-1,j,k,n)) drgt = (q(l,i+1,j,k,n) - q(l,i ,j,k,n)) dcen = half*(dlft+drgt) dsgn = sign(one, dcen) slop = min(slope_theta*abs(dlft),slope_theta*abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,1) = dsgn*min(dlim,abs(dcen)) end do ! slopes in second coordinate direction do l = 1, ngrid dlft = (q(l,i,j ,k,n) - q(l,i,j-1,k,n)) drgt = (q(l,i,j+1,k,n) - q(l,i,j ,k,n)) dcen = half*(dlft+drgt) dsgn = sign(one,dcen) slop = min(slope_theta*abs(dlft),slope_theta*abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,2) = dsgn*min(dlim,abs(dcen)) end do ! slopes in third coordinate direction do l = 1, ngrid dlft = (q(l,i,j,k ,n) - q(l,i,j,k-1,n)) drgt = (q(l,i,j,k+1,n) - q(l,i,j,k ,n)) dcen = half*(dlft+drgt) dsgn = sign(one,dcen) slop = min(slope_theta*abs(dlft),slope_theta*abs(drgt)) dlim = slop if((dlft*drgt)<=zero)dlim=zero dq(l,i,j,k,n,3) = dsgn*min(dlim,abs(dcen)) end do end do end do end do end do else write(*,*)'Unknown slope type' stop endif #endif end subroutine uslope