! $Id$ !!! NOTE: this routine will perhaps be renamed to radiation_feautrier !!! or it may be combined with radiation_ray. !** AUTOMATIC CPARAM.INC GENERATION **************************** ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! MVAR CONTRIBUTION 0 ! MAUX CONTRIBUTION 10 ! !*************************************************************** module Radiation ! Radiation (solves transfer equation along rays) ! The direction of the ray is given by the vector (lrad,mrad,nrad), ! and the parameters radx0,rady0,radz0 gives the maximum number of ! steps of the direction vector in the corresponding direction. ! ! This module currently does not work for fully periodic domains. ! The z-direction has to be non-periodic ! ! This non-grey module is EXPERIMENTAL! ! By now, it handles only two wavelenghts, very far apart (visible ! and infrared), such that one wavelength radiates all of the energy ! and the other one contributes to heating only. It reproduces in a ! simplistic way a cold disk that is illuminated by a star ! in visible wavelengths and re-radiates all the energy in infrared. ! use Cparam use General, only: keep_compiler_quiet use Messages ! implicit none include 'radiation.h' ! type Qbound !Qbc real :: val logical :: set endtype Qbound type Qpoint !Qpt real, pointer :: val logical, pointer :: set endtype Qpoint type radslice real, dimension (nx,ny) :: xy2 endtype radslice ! ! useful constant - solar luminosity ! double precision, parameter :: solar_luminosity_cgs=3.827d+33 double precision, parameter :: solar_radius_cgs=6.960d+10 double precision, parameter :: AU_cgs=1.496d+13 real :: solar_luminosity=impossible,solar_radius=impossible real :: solar_flux,solar_flux_cgs ! real, dimension (mx,my,mz) :: Srad,tau,Qrad,Qrad0,divF real, dimension (mx,my,mz,3) :: Frad type (Qbound), dimension (my,mz), target :: Qbc_yz type (Qbound), dimension (mx,mz), target :: Qbc_zx type (Qbound), dimension (mx,my), target :: Qbc_xy type (Qpoint), dimension (my,mz) :: Qpt_yz type (Qpoint), dimension (mx,mz) :: Qpt_zx type (Qpoint), dimension (mx,my) :: Qpt_xy character (len=2*bclen+1), dimension(nnu) :: bcx_rad,bcy_rad,bcz_rad !character (len=2*bclen+1), dimension(3,nnu) :: bc_rad!=(/'0:0','0:0','S:0'/) character (len=bclen), dimension(nnu) :: bcx_rad1,bcx_rad2,bcy_rad1 character (len=bclen), dimension(nnu) :: bcy_rad2,bcz_rad1,bcz_rad2 character (len=bclen) :: bc_ray_x,bc_ray_y,bc_ray_z integer, parameter :: maxdir=26 integer, dimension (maxdir,3) :: dir real, dimension (maxdir,3) :: unit_vec real, dimension (maxdir) :: weight,weightn,mu type (radslice), dimension (maxdir), target :: Isurf real :: arad real :: dtau_thresh_min,dtau_thresh_max integer :: lrad,mrad,nrad,rad2 integer :: idir,ndir integer :: l,m,n integer :: llstart,llstop,ll1,ll2,lsign integer :: mmstart,mmstop,mm1,mm2,msign integer :: nnstart,nnstop,nn1,nn2,nsign integer :: ipzstart,ipzstop,ipystart,ipystop integer :: nIsurf logical :: lperiodic_ray,lperiodic_ray_x,lperiodic_ray_y,lperiodic_ray_z character (len=labellen),dimension(nnu) :: source_function_type='LTE',opacity_type='Planck' character (len=labellen) :: angle_weight='constant' real :: tau_top=0.0,TT_top=0.0 real :: tau_bot=0.0,TT_bot=0.0 real :: kappa_cst=1.0 real :: Srad_const=1.0,amplSrad=1.0,radius_Srad=1.0 real :: kx_Srad=0.0,ky_Srad=0.0,kz_Srad=0.0 real :: kapparho_const=1.0,amplkapparho=1.0,radius_kapparho=1.0 real :: kx_kapparho=0.0,ky_kapparho=0.0,kz_kapparho=0.0 real :: Frad_boundary_ref=0.0 real :: cdtrad_thin=1., cdtrad_thick=0.8 real :: scalefactor_Srad=1. ! ! Default values for one pair of vertical rays ! integer :: radx=0,rady=0,radz=1,rad2max=1 ! logical :: lcooling=.true.,lrad_debug=.false. logical :: lintrinsic=.true.,lcommunicate=.true.,lrevision=.true. logical :: lradpressure=.false.,lradflux=.false. logical :: lrad_cool_diffus=.false., lrad_pres_diffus=.false. character :: lrad_str,mrad_str,nrad_str character(len=3) :: raydir_str ! ! Definition of dummy variables for FLD routine ! real :: DFF_new=0. !(dum) integer :: idiag_frms=0,idiag_fmax=0,idiag_Erad_rms=0,idiag_Erad_max=0 integer :: idiag_Egas_rms=0,idiag_Egas_max=0,idiag_Qradrms=0,idiag_Qradmax=0 integer :: idiag_Fradzm=0,idiag_Sradm=0,idiag_xyFradzm=0 integer :: idiag_dtchi=0,idiag_dtrad=0 real, dimension(nnu) :: nufreq namelist /radiation_init_pars/ & radx,rady,radz,rad2max,bcx_rad,bcy_rad,bcz_rad,lrad_debug,kappa_cst, & TT_top,TT_bot,tau_top,tau_bot,source_function_type,opacity_type, & Srad_const,amplSrad,radius_Srad, & kapparho_const,amplkapparho,radius_kapparho, & lintrinsic,lcommunicate,lrevision,lradflux, & Frad_boundary_ref,lrad_cool_diffus, lrad_pres_diffus, & scalefactor_Srad,nufreq,angle_weight namelist /radiation_run_pars/ & radx,rady,radz,rad2max,bcx_rad,bcy_rad,bcz_rad,lrad_debug,kappa_cst, & TT_top,TT_bot,tau_top,tau_bot,source_function_type,opacity_type, & Srad_const,amplSrad,radius_Srad, & kx_Srad,ky_Srad,kz_Srad,kx_kapparho,ky_kapparho,kz_kapparho, & kapparho_const,amplkapparho,radius_kapparho, & lintrinsic,lcommunicate,lrevision,lcooling,lradflux,lradpressure, & Frad_boundary_ref,lrad_cool_diffus,lrad_pres_diffus, & cdtrad_thin,cdtrad_thick, & scalefactor_Srad,nufreq,angle_weight contains !*********************************************************************** subroutine register_radiation() ! ! Initialize radiation flags ! ! 18-apr-07/wlad+heidar: apadted from radiation_ray ! use Cdata use FArrayManager use Mpicomm, only: stop_it ! lradiation=.true. lradiation_ray=.true. ! ! Set indices for auxiliary variables ! call farray_register_auxiliary('Qrad',iQrad) call farray_register_auxiliary('Qrad2',iQrad2) ! call farray_register_auxiliary('Qkapparho',iQkapparho) call farray_register_auxiliary('Qkapparho2',iQkapparho2) ! call farray_register_auxiliary('Frad',iFrad,vector=3) iFradx = iFrad iFrady = iFrad+1 iFradz = iFrad+2 ! call farray_register_auxiliary('Frad2',iFrad,vector=3) iFradx2 = iFrad2 iFrady2 = iFrad2+1 iFradz2 = iFrad2+2 ! if ((ip<=8) .and. lroot) then print*, 'register_radiation: radiation naux = ', naux print*, 'iQrad = ', iQrad print*, 'iQrad2 = ', iQrad2 print*, 'ikapparho = ', ikapparho print*, 'ikapparho2 = ', ikapparho2 print*, 'iFrad = ', iFrad print*, 'iFrad2 = ', iFrad2 endif ! ! Identify version number (generated automatically by CVS) ! if (lroot) call svn_id( & "$Id$") ! ! Writing files for use with IDL ! aux_var(aux_count)=',Qrad $' aux_count=aux_count+1 aux_var(aux_count)=',Qrad2 $' aux_count=aux_count+1 aux_var(aux_count)=',kapparho $' aux_count=aux_count+1 aux_var(aux_count)=',kapparho2 $' aux_count=aux_count+1 ! if (naux < maux) aux_var(aux_count)=',Frad $' if (naux == maux) aux_var(aux_count)=',Frad' aux_count=aux_count+3 ! if (naux < maux) aux_var(aux_count)=',Frad2 $' if (naux == maux) aux_var(aux_count)=',Frad2' aux_count=aux_count+3 ! if (lroot) then write(15,*) 'Qrad = fltarr(mx,my,mz)*one' write(15,*) 'Qrad2 = fltarr(mx,my,mz)*one' write(15,*) 'kapparho = fltarr(mx,my,mz)*one' write(15,*) 'kapparho2 = fltarr(mx,my,mz)*one' write(15,*) 'Frad = fltarr(mx,my,mz,3)*one' write(15,*) 'Frad2 = fltarr(mx,my,mz,3)*one' endif ! endsubroutine register_radiation !*********************************************************************** subroutine initialize_radiation() ! ! Calculate number of directions of rays ! Do this in the beginning of each run ! ! 18-apr-07/wlad+heidar: adapted from radiation_ray ! use Cdata, only: lroot,sigmaSB,c_light,pi,datadir use Cdata, only: dx,dy,dz use Cdata, only: unit_system,unit_flux,unit_length use Sub, only: parse_bc_radg use Mpicomm, only: stop_it real :: radlength,arad_normal logical :: periodic_xy_plane,bad_ray integer :: inu ! ! Check that the number of rays does not exceed maximum ! if (radx>1) call stop_it("radx currently must not be greater than 1") if (rady>1) call stop_it("rady currently must not be greater than 1") if (radz>1) call stop_it("radz currently must not be greater than 1") ! ! Check boundary conditions ! if (lroot.and.ip<14) then print*,'initialize_radiation: bcx_rad =',bcx_rad print*,'initialize_radiation: bcy_rad =',bcy_rad print*,'initialize_radiation: bcz_rad =',bcz_rad endif ! do inu=1,nnu call parse_bc_radg(bcx_rad(inu),bcx_rad1(inu),bcx_rad2(inu)) call parse_bc_radg(bcy_rad(inu),bcy_rad1(inu),bcy_rad2(inu)) call parse_bc_radg(bcz_rad(inu),bcz_rad1(inu),bcz_rad2(inu)) enddo ! ! Count ! idir=1 nIsurf=0 ! do nrad=-radz,radz do mrad=-rady,rady do lrad=-radx,radx ! ! The value of rad2 determines whether a ray is along a coordinate axis (1), ! a face diagonal (2), or a room diagonal (3) ! rad2=lrad**2+mrad**2+nrad**2 !if (rad2.eq.2) then ! print*,lrad,mrad,nrad,rad2 !endif ! ! Check whether the horizontal plane is fully periodic ! The all here refers to frequencies, not directions as ! in radiation ray. ! This should be ok, since it doesn't make much sense to make ! some frequencies periodic while the others are not. ! periodic_xy_plane=all(bcx_rad1=='p').and.all(bcx_rad2=='p')& .and.all(bcy_rad1=='p').and.all(bcy_rad2=='p') ! ! If it is, we currently don't want to calculate rays along face diagonals in ! the horizontal plane because a ray can pass through the computational domain ! many, many times before `biting itself in its tail'. ! bad_ray=(rad2==2.and.nrad==0.and.periodic_xy_plane) if ((rad2>0.and.rad2<=rad2max).and.(.not.bad_ray)) then dir(idir,1)=lrad dir(idir,2)=mrad dir(idir,3)=nrad ! ! Ray length from one grid point to the next one ! radlength=sqrt((lrad*dx)**2+(mrad*dy)**2+(nrad*dz)**2) ! ! mu = cos(theta) ! mu(idir)=nrad*dz/radlength ! ! define unit vector ! unit_vec(idir,1)=lrad*dx/radlength unit_vec(idir,2)=mrad*dy/radlength unit_vec(idir,3)=nrad*dz/radlength ! idir=idir+1 ! endif ! enddo enddo enddo ! ! Total number of directions ! ndir=idir-1 ! ! Determine when terms like exp(-dtau)-1 are to be evaluated ! as a power series ! ! Experimentally determined optimum ! Relative errors for (emdtau1, emdtau2) will be ! (1e-6, 1.5e-4) for floats and (3e-13, 1e-8) for doubles ! dtau_thresh_min=1.6*epsilon(dtau_thresh_min)**0.25 dtau_thresh_max=-log(tiny(dtau_thresh_max)) ! ! Calculate arad for LTE source function ! Note that this arad is *not* the usual radiation-density constant. ! so S = arad*TT^4 = (c/4pi)*arad_normal*TT^4, so ! arad_normal = 4pi/c*arad ! arad=SigmaSB/pi if (lroot) then arad_normal=4.*SigmaSB/c_light print*,'initialize_radiation: arad=',arad print*,'initialize_radiation: arad_normal=',arad_normal print*,'initialize_radiation: sigmaSB=',sigmaSB endif ! ! Calculate solar flux ! solar_flux_cgs=solar_luminosity_cgs/(4*pi*solar_radius_cgs**2) if (unit_system=='cgs') then solar_flux=solar_flux_cgs/unit_flux solar_radius=solar_radius_cgs/unit_length endif ! ! Calculate weights ! call calc_angle_weights ! if (lroot.and.ip<14) print*,'initialize_radiation: ndir =',ndir ! endsubroutine initialize_radiation !*********************************************************************** subroutine calc_angle_weights ! ! Weights for angle integration. ! Began coding the weigths by spherical harmonics ! ! 18-may-07/wlad+heidar: coded ! use Cdata, only: dx,dy,dz,pi,lroot use Mpicomm, only: stop_it ! real :: xyplane,yzplane,xzplane,room,xaxis,yaxis,zaxis,aspect_ratio,mu2 ! select case (angle_weight) ! case ('constant') if (ndir>0) weight=4*pi/ndir weightn=weight if (ndir==2) weightn=weightn/3. ! case ('spherical-harmonics') ! ! Check if dx==dy and that dx/dz lies in the positive range ! if (dx/=dy) then print*,'dx,dy=',dx,dy call stop_it("initialize_radiation: weights not "//& "calculated for dx/dy != 1") endif aspect_ratio=dx/dz if (aspect_ratio.lt.0.69.or.aspect_ratio.gt.sqrt(3.)) & call stop_it("initialize_radiation: weights go "//& "negative for this dx/dz ratio") ! ! Calculate the weights ! mu2=dx**2/(dx**2+dz**2) xaxis=1/42.*(4.-1./mu2) ; yaxis=xaxis zaxis=(21.-54.*mu2+34.*mu2**2)/(210.*(mu2-1)**2) xyplane=2./105.*(4.-1./mu2) yzplane=(5.-6.*mu2)/(420.*mu2*(mu2-1)**2) ; xzplane=yzplane room=(2.-mu2)**3/(840.*mu2*(mu2-1)**2) ! ! Allocate the weights on the appropriate rays ! do idir=1,ndir !axes if (dir(idir,1)/=0.and.dir(idir,2)==0.and.dir(idir,3)==0) weight(idir)=xaxis if (dir(idir,1)==0.and.dir(idir,2)/=0.and.dir(idir,3)==0) weight(idir)=yaxis if (dir(idir,1)==0.and.dir(idir,2)==0.and.dir(idir,3)/=0) weight(idir)=zaxis !face diagonals if (dir(idir,1)==0.and.dir(idir,2)/=0.and.dir(idir,3)/=0) weight(idir)=yzplane if (dir(idir,1)/=0.and.dir(idir,2)==0.and.dir(idir,3)/=0) weight(idir)=xzplane if (dir(idir,1)/=0.and.dir(idir,2)/=0.and.dir(idir,3)==0) weight(idir)=xyplane !room diagonal if (dir(idir,1)/=0.and.dir(idir,2)/=0.and.dir(idir,3)/=0) weight(idir)=room if (lroot.and.ip<11) & print*,'initialize_radiation: dir(idir,1:3),weight(idir) =',& dir(idir,1:3),weight(idir) enddo weightn=weight ! case default call stop_it('no such angle-weighting: '//& trim(angle_weight)) endselect ! endsubroutine calc_angle_weights !*********************************************************************** subroutine radtransfer(f) ! ! Integration of the radiative transfer equation along rays ! ! This routine is called before the communication part ! (certainly needs to be given a better name) ! All rays start with zero intensity ! ! 18-apr-07/wlad+heidar: adpated from radiation_ray ! use Cdata, only: ldebug,headt,iQrad,ikapparho use Cdata, only: iFrad,iFradx,iFradz,iFrad2,iFradx2,iFradz2 ! real, dimension(mx,my,mz,mfarray) :: f integer :: j,k,nnu,inu real, dimension(2) :: nuvector integer :: iQ,ikr ! ! Identifier ! if (ldebug.and.headt) print*,'radtransfer' ! ! Calculate source function and opacity ! nnu=2 ! ! Initialize (bolometric) radiative flux and cooling function ! if (lradflux) then f(:,:,:,iFradx:iFradz)=0 f(:,:,:,iFradx2:iFradz2)=0 endif ! ! Initialize cooling function ! divF(:,:,:)=0. ! do inu=1,nnu ! call source_function(f,inu) call opacity(f,inu) ! ! do the rest only if we don't do diffusion approximation ! if (lrad_cool_diffus.or.lrad_pres_diffus) then if (headt) print*,'do diffusion approximation, no rays' else ! ! Initialize heating rates ! iQ=iQrad-1+inu ikr=ikapparho-1+inu ! f(:,:,:,iQ)=0 ! ! loop over rays ! do idir=1,ndir ! call raydirection(inu) ! ! Qintrinsic: no need to send the full f array if just the ! current opacity is needed ! if (lintrinsic) call Qintrinsic(f(:,:,:,ikr)) ! if (lcommunicate) then if (lperiodic_ray) then call Qperiodic else call Qpointers call Qcommunicate endif endif ! if (lrevision) call Qrevision ! ! calculate heating rate, so at the end of the loop ! f(:,:,:,iQrad) = \int_{4\pi} (I-S) d\Omega, not divided by 4pi. ! In the paper the directional Q is defined with the opposite sign. ! f(:,:,:,iQ)=f(:,:,:,iQ)+weight(idir)*Qrad ! ! calculate bolometric radiative flux ! if (lradflux) then do j=1,3 !which frequency? k=iFrad+3*(inu-1)+(j-1) f(:,:,:,k)=f(:,:,:,k)+weightn(idir)*unit_vec(idir,j)*(Qrad+Srad) enddo endif ! enddo ! divF=divF + f(:,:,:,ikr)*f(:,:,:,iQ) ! ! end of (ifnot) diffusion approximation ! endif ! ! end of frequency loop ! enddo ! endsubroutine radtransfer !*********************************************************************** subroutine raydirection(inu) ! ! Determine certain variables depending on the ray direction ! ! 18-apr-07/wlad+heidar: apapted from radiation_ray ! use Cdata, only: ldebug,headt integer :: inu ! ! Identifier ! if (ldebug.and.headt) print*,'raydirection' ! ! Get direction components ! lrad=dir(idir,1) mrad=dir(idir,2) nrad=dir(idir,3) ! ! Determine start and stop positions ! llstart=l1; llstop=l2; ll1=l1; ll2=l2; lsign=+1 mmstart=m1; mmstop=m2; mm1=m1; mm2=m2; msign=+1 nnstart=n1; nnstop=n2; nn1=n1; nn2=n2; nsign=+1 if (lrad>0) then; llstart=l1; llstop=l2; ll1=l1-lrad; lsign=+1; endif if (lrad<0) then; llstart=l2; llstop=l1; ll2=l2-lrad; lsign=-1; endif if (mrad>0) then; mmstart=m1; mmstop=m2; mm1=m1-mrad; msign=+1; endif if (mrad<0) then; mmstart=m2; mmstop=m1; mm2=m2-mrad; msign=-1; endif if (nrad>0) then; nnstart=n1; nnstop=n2; nn1=n1-nrad; nsign=+1; endif if (nrad<0) then; nnstart=n2; nnstop=n1; nn2=n2-nrad; nsign=-1; endif ! ! Determine boundary conditions ! if (lrad>0) bc_ray_x=bcx_rad1(inu) if (lrad<0) bc_ray_x=bcx_rad2(inu) if (mrad>0) bc_ray_y=bcy_rad1(inu) if (mrad<0) bc_ray_y=bcy_rad2(inu) if (nrad>0) bc_ray_z=bcz_rad1(inu) if (nrad<0) bc_ray_z=bcz_rad2(inu) ! ! Are we dealing with a periodic ray? ! lperiodic_ray_x=(lrad/=0.and.mrad==0.and.nrad==0.and.bc_ray_x=='p') lperiodic_ray_y=(lrad==0.and.mrad/=0.and.nrad==0.and.bc_ray_y=='p') lperiodic_ray=(lperiodic_ray_x.or.lperiodic_ray_y) ! ! Determine start and stop processors ! if (mrad>0) then; ipystart=0; ipystop=nprocy-1; endif if (mrad<0) then; ipystart=nprocy-1; ipystop=0; endif if (nrad>0) then; ipzstart=0; ipzstop=nprocz-1; endif if (nrad<0) then; ipzstart=nprocz-1; ipzstop=0; endif ! ! Label for debug output ! if (lrad_debug) then lrad_str='0'; mrad_str='0'; nrad_str='0' if (lrad>0) lrad_str='p' if (lrad<0) lrad_str='m' if (mrad>0) mrad_str='p' if (mrad<0) mrad_str='m' if (nrad>0) nrad_str='p' if (nrad<0) nrad_str='m' raydir_str=lrad_str//mrad_str//nrad_str endif endsubroutine raydirection !*********************************************************************** subroutine Qintrinsic(fikr) ! ! Integration radiation transfer equation along rays ! ! This routine is called before the communication part ! All rays start with zero intensity ! ! 18-apr-07/wlad+heidar: adapted from radiation_ray ! use Cdata, only: ldebug,headt,dx,dy,dz,directory_snap,ikapparho use IO, only: output integer :: ikr ! real, dimension(mx,my,mz), intent(in) :: fikr real :: Srad1st,Srad2nd,dlength,emdtau1,emdtau2,emdtau real :: dtau_m,dtau_p,dSdtau_m,dSdtau_p character(len=3) :: raydir ! ! identifier ! if (ldebug.and.headt) print*,'Qintrinsic' ! ! line elements ! dlength=sqrt((dx*lrad)**2+(dy*mrad)**2+(dz*nrad)**2) ! ! set optical depth and intensity initially to zero ! tau=0 Qrad=0 ! ! loop over all meshpoints ! do n=nnstart,nnstop,nsign do m=mmstart,mmstop,msign do l=llstart,llstop,lsign dtau_m=sqrt(fikr(l-lrad,m-mrad,n-nrad)* & fikr(l,m,n))*dlength dtau_p=sqrt(fikr(l,m,n)* & fikr(l+lrad,m+mrad,n+nrad))*dlength dSdtau_m=(Srad(l,m,n)-Srad(l-lrad,m-mrad,n-nrad))/dtau_m dSdtau_p=(Srad(l+lrad,m+mrad,n+nrad)-Srad(l,m,n))/dtau_p Srad1st=(dSdtau_p*dtau_m+dSdtau_m*dtau_p)/(dtau_m+dtau_p) Srad2nd=2*(dSdtau_p-dSdtau_m)/(dtau_m+dtau_p) if (dtau_m>dtau_thresh_max) then emdtau=0.0 emdtau1=1.0 emdtau2=-1.0 elseif (dtau_m set Qpt_yz(m,n)%val => val enddo enddo endif ! ! zx-plane ! if (mrad/=0) then m=mmstop msteps=(m+mrad-mmstart)/mrad nsteps=huge(nsteps) lsteps=huge(lsteps) do n=nn1,nn2 do l=ll1,ll2 if (nrad/=0) nsteps = (n+nrad-nnstart)/nrad if (lrad/=0) lsteps = (l+lrad-llstart)/lrad call assign_pointer(lsteps,msteps,nsteps,val,set) Qpt_zx(l,n)%set => set Qpt_zx(l,n)%val => val enddo enddo endif ! ! xy-plane ! if (nrad/=0) then n=nnstop nsteps=(n+nrad-nnstart)/nrad lsteps=huge(lsteps) msteps=huge(msteps) do l=ll1,ll2 do m=mm1,mm2 if (lrad/=0) lsteps = (l+lrad-llstart)/lrad if (mrad/=0) msteps = (m+mrad-mmstart)/mrad call assign_pointer(lsteps,msteps,nsteps,val,set) Qpt_xy(l,m)%set => set Qpt_xy(l,m)%val => val enddo enddo endif endsubroutine Qpointers !*********************************************************************** subroutine assign_pointer(lsteps,msteps,nsteps,val,set) integer, intent(in) :: lsteps,msteps,nsteps real, pointer :: val logical, pointer :: set integer :: steps steps=min(lsteps,msteps,nsteps) if (steps==lsteps) then val => Qbc_yz(m-mrad*steps,n-nrad*steps)%val set => Qbc_yz(m-mrad*steps,n-nrad*steps)%set endif if (steps==msteps) then val => Qbc_zx(l-lrad*steps,n-nrad*steps)%val set => Qbc_zx(l-lrad*steps,n-nrad*steps)%set endif if (steps==nsteps) then val => Qbc_xy(l-lrad*steps,m-mrad*steps)%val set => Qbc_xy(l-lrad*steps,m-mrad*steps)%set endif endsubroutine assign_pointer !*********************************************************************** subroutine Qcommunicate ! ! Determine the boundary heating rates at all upstream boundaries. ! ! First the boundary heating rates at the non-periodic xy-boundary ! are set either through the boundary condition for the entire ! computational domain (ipz==ipzstart) or through communication with ! the neighboring processor in the upstream z-direction (ipz/=ipzstart). ! ! The boundary heating rates at the periodic yz- and zx-boundaries ! are then obtained by repetitive communication along the y-direction ! until both boundaries are entirely set with the correct values. ! ! 30-jul-05/tobi: coded ! use Cdata, only: ipy,ipz use Mpicomm, only: radboundary_xy_recv,radboundary_xy_send use Mpicomm, only: radboundary_zx_recv,radboundary_zx_send use Mpicomm, only: radboundary_zx_sendrecv real, dimension (my,mz) :: emtau_yz,Qrad_yz real, dimension (mx,mz) :: emtau_zx,Qrad_zx real, dimension (mx,my) :: emtau_xy,Qrad_xy real, dimension (my,mz) :: Qsend_yz,Qrecv_yz real, dimension (mx,mz) :: Qsend_zx,Qrecv_zx real, dimension (mx,my) :: Qrecv_xy,Qsend_xy logical :: all_yz,all_zx ! ! Initially no boundaries are set ! Qbc_xy%set=.false. Qbc_yz%set=.false. Qbc_zx%set=.false. all_yz=.false. all_zx=.false. ! ! Either receive or set xy-boundary heating rate ! if (nrad/=0) then if (ipz==ipzstart) then call radboundary_xy_set(Qrecv_xy) else call radboundary_xy_recv(nrad,idir,Qrecv_xy) endif ! ! Copy the above heating rates to the xy-target arrays which are then set ! Qbc_xy(ll1:ll2,mm1:mm2)%val = Qrecv_xy(ll1:ll2,mm1:mm2) Qbc_xy(ll1:ll2,mm1:mm2)%set = .true. endif ! ! do the same for the yz- and zx-target arrays where those boundaries ! overlap with the xy-boundary and calculate exp(-tau) and Qrad at the ! downstream boundaries. ! if (lrad/=0) then if (bc_ray_x/='p') then call radboundary_yz_set(Qrecv_yz) Qbc_yz(mm1:mm2,nn1:nn2)%val = Qrecv_yz(mm1:mm2,nn1:nn2) Qbc_yz(mm1:mm2,nn1:nn2)%set = .true. all_yz=.true. else Qbc_yz(mm1:mm2,nnstart-nrad)%val = Qrecv_xy(llstart-lrad,mm1:mm2) Qbc_yz(mm1:mm2,nnstart-nrad)%set = .true. emtau_yz(mm1:mm2,nn1:nn2) = exp(-tau(llstop,mm1:mm2,nn1:nn2)) Qrad_yz(mm1:mm2,nn1:nn2) = Qrad(llstop,mm1:mm2,nn1:nn2) endif else all_yz=.true. endif if (mrad/=0) then if (bc_ray_y/='p') then call radboundary_zx_set(Qrecv_zx) Qbc_zx(ll1:ll2,nn1:nn2)%val = Qrecv_zx(ll1:ll2,nn1:nn2) Qbc_zx(ll1:ll2,nn1:nn2)%set = .true. all_zx=.true. else Qbc_zx(ll1:ll2,nnstart-nrad)%val = Qrecv_xy(ll1:ll2,mmstart-mrad) Qbc_zx(ll1:ll2,nnstart-nrad)%set = .true. emtau_zx(ll1:ll2,nn1:nn2) = exp(-tau(ll1:ll2,mmstop,nn1:nn2)) Qrad_zx(ll1:ll2,nn1:nn2) = Qrad(ll1:ll2,mmstop,nn1:nn2) endif else all_zx=.true. endif ! ! communicate along the y-direction until all upstream heating rates at ! the yz- and zx-boundaries are determined. ! if ((lrad/=0.and..not.all_yz).or.(mrad/=0.and..not.all_zx)) then; do if (lrad/=0.and..not.all_yz) then forall (m=mm1:mm2,n=nn1:nn2,Qpt_yz(m,n)%set.and..not.Qbc_yz(m,n)%set) Qbc_yz(m,n)%val = Qpt_yz(m,n)%val*emtau_yz(m,n)+Qrad_yz(m,n) Qbc_yz(m,n)%set = Qpt_yz(m,n)%set endforall all_yz=all(Qbc_yz(mm1:mm2,nn1:nn2)%set) if (all_yz.and.all_zx) exit endif if (mrad/=0.and..not.all_zx) then forall (l=ll1:ll2,n=nn1:nn2,Qpt_zx(l,n)%set.and..not.Qbc_zx(l,n)%set) Qsend_zx(l,n) = Qpt_zx(l,n)%val*emtau_zx(l,n)+Qrad_zx(l,n) endforall if (nprocy>1) then call radboundary_zx_sendrecv(mrad,idir,Qsend_zx,Qrecv_zx) else Qrecv_zx=Qsend_zx endif forall (l=ll1:ll2,n=nn1:nn2,Qpt_zx(l,n)%set.and..not.Qbc_zx(l,n)%set) Qbc_zx(l,n)%val = Qrecv_zx(l,n) Qbc_zx(l,n)%set = Qpt_zx(l,n)%set endforall all_zx=all(Qbc_zx(ll1:ll2,nn1:nn2)%set) if (all_yz.and.all_zx) exit endif enddo; endif ! ! copy all heating rates at the upstream boundaries to Qrad0 which is used in ! Qrevision below. ! if (lrad/=0) then Qrad0(llstart-lrad,mm1:mm2,nn1:nn2)=Qbc_yz(mm1:mm2,nn1:nn2)%val endif if (mrad/=0) then Qrad0(ll1:ll2,mmstart-mrad,nn1:nn2)=Qbc_zx(ll1:ll2,nn1:nn2)%val endif if (nrad/=0) then Qrad0(ll1:ll2,mm1:mm2,nnstart-nrad)=Qbc_xy(ll1:ll2,mm1:mm2)%val endif ! ! If this is not the last processor in ray direction (z-component) then ! calculate the downstream heating rates at the xy-boundary and send them ! to the next processor. ! if (nrad/=0.and.ipz/=ipzstop) then forall (l=ll1:ll2,m=mm1:mm2) emtau_xy(l,m) = exp(-tau(l,m,nnstop)) Qrad_xy(l,m) = Qrad(l,m,nnstop) Qsend_xy(l,m) = Qpt_xy(l,m)%val*emtau_xy(l,m)+Qrad_xy(l,m) endforall call radboundary_xy_send(nrad,idir,Qsend_xy) endif endsubroutine Qcommunicate !*********************************************************************** subroutine Qperiodic ! ! DOCUMENT ME! ! use Cdata, only: ipy,iproc use Mpicomm, only: radboundary_zx_periodic_ray use IO, only: output real, dimension(ny,nz) :: Qrad_yz,tau_yz,emtau1_yz real, dimension(nx,nz) :: Qrad_zx,tau_zx,emtau1_zx real, dimension(nx,nz) :: Qrad_tot_zx,tau_tot_zx,emtau1_tot_zx real, dimension(nx,nz,0:nprocy-1) :: Qrad_zx_all,tau_zx_all integer :: ipystart,ipystop,ipm ! ! x-direction ! if (lrad/=0) then ! ! Intrinsic heating rate and optical depth at the downstream boundary. ! Qrad_yz=Qrad(llstop,m1:m2,n1:n2) tau_yz=tau(llstop,m1:m2,n1:n2) ! ! Try to avoid time consuming exponentials and loss of precision. ! where (tau_yz>dtau_thresh_max) emtau1_yz=1.0 elsewhere (tau_yz0) then; ipystart=0; ipystop=nprocy-1; endif if (mrad<0) then; ipystart=nprocy-1; ipystop=0; endif ! ! We need the sum of all intrinsic optical depths and the attenuated sum of ! all intrinsic heating rates. The latter needs to be summed in the ! downstream direction starting at the current processor. Set both to zero ! initially. ! Qrad_tot_zx=0.0 tau_tot_zx=0.0 ! ! Do the sum from the this processor to the last one in the downstream ! direction. ! do ipm=ipy,ipystop,msign Qrad_tot_zx=Qrad_tot_zx*exp(-tau_zx_all(:,:,ipm))+Qrad_zx_all(:,:,ipm) tau_tot_zx=tau_tot_zx+tau_zx_all(:,:,ipm) enddo ! ! Do the sum from the first processor in the upstream direction to the one ! before this one. ! do ipm=ipystart,ipy-msign,msign Qrad_tot_zx=Qrad_tot_zx*exp(-tau_zx_all(:,:,ipm))+Qrad_zx_all(:,:,ipm) tau_tot_zx=tau_tot_zx+tau_zx_all(:,:,ipm) enddo ! ! To calculate the boundary heating rate we need to compute an exponential ! term involving the total optical depths across all processors. ! Try to avoid time consuming exponentials and loss of precision. ! where (tau_tot_zx>dtau_thresh_max) emtau1_tot_zx=1.0 elsewhere (tau_tot_zx0) then Isurf(idir)%xy2=Qrad(l1:l2,m1:m2,nnstop)+Srad(l1:l2,m1:m2,nnstop) endif ! if (lrad_debug) then call output(trim(directory_snap)//'/Qrev-'//raydir_str//'.dat',Qrad,1) endif ! endsubroutine Qrevision !*********************************************************************** subroutine radboundary_yz_set(Qrad0_yz) ! ! Sets the physical boundary condition on yz plane ! ! 6-jul-03/axel+tobi: coded ! 26-may-06/axel: added S+F and S-F to model interior more accurately ! use Mpicomm, only: stop_it use Cdata, only: x ! real, dimension(my,mz), intent(out) :: Qrad0_yz real :: Irad_yz real :: rad,Frad_bound ! ! No incoming intensity ! if (bc_ray_z=='0') then Qrad0_yz=-Srad(llstart-lrad,:,:) endif ! ! Set intensity equal to source function ! if (bc_ray_z=='S') then Qrad0_yz=0 endif ! ! Set intensity equal to S-F ! if (bc_ray_z=='S-F') then Qrad0_yz=-Frad_boundary_ref/(2.*weightn(idir)) endif ! ! Set intensity equal to S+F ! if (bc_ray_z=='S+F') then Qrad0_yz=+Frad_boundary_ref/(2.*weightn(idir)) endif ! ! Set intensity equal to a pre-defined incoming intensity ! if (bc_ray_z=='F') then Qrad0_yz=-Srad(llstart-lrad,:,:)+Frad_boundary_ref/(2.*weightn(idir)) endif ! ! Set intensity equal to a pre-defined radially variable intensity ! if (bc_ray_z=='Fr') then rad = x(llstart-lrad) Frad_bound=solar_flux*(solar_radius/rad)**2 Qrad0_yz=-Srad(llstart-lrad,:,:)+Frad_bound/(2.*weightn(idir)) endif ! endsubroutine radboundary_yz_set !*********************************************************************** subroutine radboundary_zx_set(Qrad0_zx) ! ! Sets the physical boundary condition on zx plane ! ! 6-jul-03/axel+tobi: coded ! 26-may-06/axel: added S+F and S-F to model interior more accurately ! use Cdata, only: x use Mpicomm, only: stop_it ! real, dimension(mx,mz), intent(out) :: Qrad0_zx real :: Irad_zx real, dimension(mx) :: rad,Frad_bound integer :: in ! ! No incoming intensity ! if (bc_ray_z=='0') then Qrad0_zx=-Srad(:,mmstart-mrad,:) endif ! ! Set intensity equal to source function ! if (bc_ray_z=='S') then Qrad0_zx=0 endif ! ! Set intensity equal to S-F ! if (bc_ray_z=='S-F') then Qrad0_zx=-Frad_boundary_ref/(2.*weightn(idir)) endif ! ! Set intensity equal to S+F ! if (bc_ray_z=='S+F') then Qrad0_zx=+Frad_boundary_ref/(2.*weightn(idir)) endif ! ! Set intensity equal to a pre-defined incoming intensity ! if (bc_ray_z=='F') then Qrad0_zx=-Srad(:,mmstart-mrad,:)+Frad_boundary_ref/(2.*weightn(idir)) endif ! ! Set intensity equal to a pre-defined radially variable incoming intensity ! if (bc_ray_z=='Fr') then rad = x Frad_bound=solar_flux*(solar_radius/rad)**2 do in=1,mz Qrad0_zx(:,in)=-Srad(:,mmstart-mrad,in)+Frad_bound/(2.*weightn(idir)) enddo endif ! endsubroutine radboundary_zx_set !*********************************************************************** subroutine radboundary_xy_set(Qrad0_xy) ! ! Sets the physical boundary condition on xy plane ! ! 6-jul-03/axel+tobi: coded ! 26-may-06/axel: added S+F and S-F to model interior more accurately ! use Cdata, only: x,unit_length use Mpicomm, only: stop_it ! real, dimension(mx,my), intent(out) :: Qrad0_xy real :: Irad_xy real, dimension(mx) :: rad,Frad_bound,angle_recipe,rau integer :: im ! ! No incoming intensity ! if (bc_ray_z=='0') then Qrad0_xy=-Srad(:,:,nnstart-nrad) endif ! ! incoming intensity from a layer of constant temperature TT_top ! if (bc_ray_z=='c') then if (nrad<0) Irad_xy=arad*TT_top**4*(1-exp(tau_top/mu(idir))) if (nrad>0) Irad_xy=arad*TT_bot**4*(1-exp(tau_top/mu(idir))) Qrad0_xy=Irad_xy-Srad(:,:,nnstart-nrad) endif ! ! Set intensity equal to source function ! if (bc_ray_z=='S') then Qrad0_xy=0 endif ! ! Set intensity equal to S-F ! if (bc_ray_z=='S-F') then Qrad0_xy=-Frad_boundary_ref/(2.*weightn(idir)) endif ! ! Set intensity equal to S+F ! if (bc_ray_z=='S+F') then Qrad0_xy=+Frad_boundary_ref/(2.*weightn(idir)) endif ! ! Set intensity equal to a pre-defined incoming intensity ! if (bc_ray_z=='F') then Qrad0_xy=-Srad(:,:,nnstart-nrad)+Frad_boundary_ref/(2.*weightn(idir)) endif ! ! Set intensity equal to a pre-defined radially variable incoming intensity ! The stellar flux falls with the r^2 law and in addition ! the angle recipe accounts for the effect that the flaring of the disk ! has on the intercepted stellar radiation, based on Chiang & Goldreich '97 ! The factor is 0.5*0.05*r[AU]^(2/7), valid for 0.4 to 84 AU, ! where 0.5 comes in because only half the star is seen ! if (bc_ray_z=='Fr') then rad = x rau=rad*unit_length/AU_cgs angle_recipe=0.025*rau**0.286 Frad_bound=solar_flux*(solar_radius/rad)**2*angle_recipe*.25 do im=1,my Qrad0_xy(:,im)=-Srad(:,im,nnstart-nrad)+Frad_bound/(2.*weightn(idir)) enddo endif ! endsubroutine radboundary_xy_set !*********************************************************************** subroutine radiative_cooling(f,df,p) ! ! Add radiative cooling to entropy/temperature equation ! ! 18-apr-07/wlad+heidar: adapted from radiation_ray ! use Cdata use Diagnostics use Mpicomm, only: stop_it ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p real, dimension (nx) :: cooling,Qrad2 integer :: i ! ! add radiative cooling, either from the intensity ! or in the diffusion approximation ! AB: Tobi, would it be better/possible to redefine Qrad so as ! AB: to include the kapparho factor. Then we'd have Qrad=-divFrad. ! if (lrad_cool_diffus.or.lrad_pres_diffus) then call stop_it("radiative_cooling: diffusion not implemented for non-grey RT") !call calc_rad_diffusion(f,p) !cooling=f(l1:l2,m,n,iQrad) else if (any(opacity_type=='kappa_es')) then call stop_it("radiative_cooling: scattering not implemented for non-grey RT") call calc_rad_diffusion(f,p) endif cooling=divF(l1:l2,m,n) !cooling=f(l1:l2,m,n,iQrad)*f(l1:l2,m,n,ikapparho) !just the 1st frequency (visible) endif ! ! Add radiative cooling ! if (lcooling) then if (lentropy) then df(l1:l2,m,n,iss)=df(l1:l2,m,n,iss)+p%rho1*p%TT1*cooling endif if (ltemperature) then df(l1:l2,m,n,ilnTT)=df(l1:l2,m,n,ilnTT)+p%rho1*p%cv1*p%TT1*cooling endif endif ! ! diagnostics ! if (ldiagnos) then Qrad2=f(l1:l2,m,n,iQrad )**2 !Qrad2=f(l1:l2,m,n,iQrad2)**2 if (idiag_Sradm/=0) call sum_mn_name(Srad(l1:l2,m,n),idiag_Sradm) if (idiag_Qradrms/=0) then call sum_mn_name(Qrad2,idiag_Qradrms,lsqrt=.true.) endif if (idiag_Qradmax/=0) then call max_mn_name(Qrad2,idiag_Qradmax,lsqrt=.true.) endif endif ! endsubroutine radiative_cooling !*********************************************************************** subroutine radiative_pressure(f,df,p) ! ! Add radiative pressure to equation of motion ! ! 17-may-06/axel: coded ! use Cdata use Diagnostics use Mpicomm, only: stop_it real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p real, dimension (nx,3) :: radpressure integer :: j,k,iflux,ik,inu ! ! radiative pressure force = kappa*Frad/c = (kappa*rho)*rho1*Frad/c ! if (lradpressure) then call stop_it("radiative_pressure: not implemented for non-grey RT") do inu=1,nnu do j=1,3 iflux=iFrad+(j-1)+3*(inu-1) ik=ikapparho+inu-1 radpressure(:,j)=p%rho1*f(l1:l2,m,n,ik)*f(l1:l2,m,n,iflux)/c_light enddo df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)+radpressure enddo endif ! ! diagnostics ! if (ldiagnos) then if (idiag_Fradzm/=0) then call sum_mn_name(f(l1:l2,m,n,iFradz),idiag_Fradzm) endif endif ! if (l1davgfirst) then if (idiag_xyFradzm/=0) then call xysum_mn_name_z(f(l1:l2,m,n,iFradz),idiag_xyFradzm) endif endif ! endsubroutine radiative_pressure !*********************************************************************** subroutine source_function(f,inu) ! ! calculates source function ! (This module is currently ignored if diffusion approximation is used) ! ! 18-apr-07/wlad+heidar: adapted from radiation_ray ! use Cdata, only: m,n,x,y,z,Lx,Ly,Lz,dx,dy,dz,pi,directory_snap,hbar,k_B,c_light use Mpicomm, only: stop_it use EquationOfState, only: eoscalc use IO, only: output real, dimension(mx,my,mz,mfarray), intent(in) :: f logical, save :: lfirst=.true. real, dimension(mx) :: lnTT,rr real :: nu,width,rrp integer :: inu,i !if (source_function_type(1)/='zero') & ! call stop_it("source_function: thou shall not have cold gas emitting in visible wavelengths") select case (source_function_type(inu)) case ('Planck') nu=nufreq(inu) do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnTT=lnTT) Srad(:,m,n)=4*pi*hbar*nu**3/& (c_light**2*(exp((2*pi*hbar*nu)/(k_B*exp(lnTT)))-1)) enddo enddo case ('LTE') do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnTT=lnTT) Srad(:,m,n)=arad*exp(4*lnTT)*scalefactor_Srad enddo enddo case ('blob') if (lfirst) then Srad=Srad_const & +amplSrad*spread(spread(exp(-(x/radius_Srad)**2),2,my),3,mz) & *spread(spread(exp(-(y/radius_Srad)**2),1,mx),3,mz) & *spread(spread(exp(-(z/radius_Srad)**2),1,mx),2,my) lfirst=.false. endif case ('cos') if (lfirst) then Srad=Srad_const & +amplSrad*spread(spread(cos(kx_Srad*x),2,my),3,mz) & *spread(spread(cos(ky_Srad*y),1,mx),3,mz) & *spread(spread(cos(kz_Srad*z),1,mx),2,my) lfirst=.false. endif case ('zero') Srad(:,:,:)=0. case default call stop_it('no such source function type: '//& trim(source_function_type(inu))) end select if (lrad_debug) then call output(trim(directory_snap)//'/Srad.dat',Srad,1) endif endsubroutine source_function !*********************************************************************** subroutine opacity(f,inu) ! ! calculates opacity ! ! Note that currently the diffusion approximation does not take ! into account any gradient terms of kappa. ! ! 18-apr-07/wlad+heidar: adapted from radiation_ray ! use Cdata, only: ilnrho,x,y,z,m,n,Lx,Ly,Lz,dx,dy,dz,pi,directory_snap use Cdata, only: kappa_es,ikapparho,unit_density,unit_length,unit_temperature use EquationOfState, only: eoscalc use Mpicomm, only: stop_it use IO, only: output real, dimension(mx,my,mz,mfarray), intent(inout) :: f real, dimension(mx) :: tmp,lnrho,lnTT,qq,kapparho,lnkappa real :: kappa0, kappa0_cgs,k1,k2,rhoref=1e-9 real :: TT7,TT8,TT9,lnTT7,lnTT8,lnTT9,qqedge,a,b !! logical, save :: lfirst=.true. integer :: i,inu,ikr select case (opacity_type(inu)) case ('Hminus') do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,kapparho=tmp) f(:,m,n,ikapparho)=tmp enddo enddo case ('kappa_es') do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnrho=lnrho) f(:,m,n,ikapparho)=kappa_es*exp(lnrho) enddo enddo case ('kappa_cst') do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnrho=lnrho) f(:,m,n,ikapparho)=kappa_cst*exp(lnrho) enddo enddo case ('Tsquare') ! ! HT: Case of Morfill et al. 1985 ! The coefficient is wrongly stated as 2e-6 in D'Angelo et al 2003 ! kappa0_cgs=2e-4 kappa0=kappa0_cgs*unit_density*unit_length do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnrho=lnrho,lnTT=lnTT) f(:,m,n,ikapparho)=exp(lnrho)*kappa0*((exp(lnTT))**2) enddo enddo case ('Kramers') ! ! HT: Case of Frank et al. 1992, also used in D'Angelo et al 2003 ! The coefficient is 7.5e22 in Prialnik ! kappa0_cgs=6.6e22 kappa0=kappa0_cgs*unit_density*unit_length do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnrho=lnrho,lnTT=lnTT) f(:,m,n,ikapparho)=kappa0*(exp(lnrho)**2)*(exp(lnTT))**(-3.5) enddo enddo case ('dust-visible', 'dust-infrared') ! ! WL: What's the interval of validity? ! if (unit_temperature/=1) & call stop_it("opacity: dust opacity just works for unit_temperature=1") ikr=ikapparho-1+inu do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnrho=lnrho,lnTT=lnTT) do i=1,mx if (exp(lnTT(i)).le.150) then tmp(i)=2e-4*exp(lnTT(i))**2 elseif (exp(lnTT(i)).ge.200) then k1=0.861353*lnTT(i)-4.56372 tmp(i)=exp(k1) else k2=-5.22826*lnTT(i)+27.7010 tmp(i)=exp(k2) endif enddo kapparho=exp(lnrho)*tmp if (opacity_type(inu)=='dust-infrared') & f(:,m,n,ikr)=kapparho if (opacity_type(inu)=='dust-visible') then ! ! HT: q-recipe of Calvet et al 1991 for the opacity ! increase on visible wavelengths ! valid for 100K < T < 200K, assumed constant for T < 100K ! for 200-1400K, qq=exp(-8.93e-4*exp(lnTT)+3.42) ! qq=min(exp(-1.2e-2*exp(lnTT)+5.8),1e2) f(:,m,n,ikr)=kapparho*qq endif enddo enddo case ('dust-visible2', 'dust-infrared2') ! ! HT: extended to higher temperatures, following table in Bell et al. 97 ! if (unit_temperature/=1) & call stop_it("opacity: dust opacity just works for unit_temperature=1") ikr=ikapparho-1+inu do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnrho=lnrho,lnTT=lnTT) do i=1,mx lnkappa(i)=1. ! skip? if (exp(lnTT(i)).le.132) then tmp(i)=1e-4*exp(lnTT(i))**(2.1) elseif (exp(lnTT(i)).le.170) then tmp(i)=3.0*exp(lnTT(i))**(-0.01) elseif (exp(lnTT(i)).le.375) then tmp(i)=1e-2*exp(lnTT(i))**(1.1) elseif (exp(lnTT(i)).le.390) then tmp(i)=5e4*exp(lnTT(i))**(-1.5) elseif (exp(lnTT(i)).le.580) then tmp(i)=1e-1*exp(lnTT(i))**(0.7) elseif (exp(lnTT(i)).le.680) then lnkappa(i)=23.94-5.2*lnTT(i) tmp(i)=exp(lnkappa(i)) ! ! density dependence comes in... ! elseif (exp(lnTT(i)).le.960) then tmp(i)=2e-2*exp(lnTT(i))**(0.8) elseif (exp(lnTT(i)).le.1570) then lnkappa(i)=187.2+log(rhoref)-24.*lnTT(i) tmp(i)=exp(lnkappa(i)) elseif (exp(lnTT(i)).le.3730) then tmp(i)=1e-8*rhoref**(2./3.)*exp(lnTT(i))**(3.0) else lnkappa(i)=-82.9+(log(rhoref)/3.)+(10.*lnTT(i)) tmp(i)=exp(lnkappa(i)) endif enddo kapparho=exp(lnrho)*tmp if (opacity_type(inu)=='dust-infrared2') & f(:,m,n,ikr)=kapparho if (opacity_type(inu)=='dust-visible2') then ! ! q-recipe of Calvet et al 1991 for the opacity ! increase on visible wavelengths ! really only valid for 100K= nIsurf) then slices%ready=.false. else if (slices%index == 0) idir_last=0 nullify(slices%yz) nullify(slices%xz) nullify(slices%xy) do idir=idir_last+1,ndir nrad=dir(idir,3) if (nrad>0) then slices%xy2=>Isurf(idir)%xy2 slices%index=slices%index+1 if (slices%index<=nIsurf) slices%ready=.true. idir_last=idir exit endif enddo endif ! ! Heating rate ! case ('Qrad') slices%yz=f(ix_loc,m1:m2,n1:n2,iQrad) slices%xz=f(l1:l2,iy_loc,n1:n2,iQrad) slices%xy=f(l1:l2,m1:m2,iz_loc,iQrad) slices%xy2=f(l1:l2,m1:m2,iz2_loc,iQrad) slices%ready = .true. ! ! Opacity ! case ('kapparho') slices%yz=f(ix_loc,m1:m2,n1:n2,ikapparho) slices%xz=f(l1:l2,iy_loc,n1:n2,ikapparho) slices%xy=f(l1:l2,m1:m2,iz_loc,ikapparho) slices%xy2=f(l1:l2,m1:m2,iz2_loc,ikapparho) slices%ready = .true. ! ! Radiative Flux ! case ('Frad') if (slices%index>=3) then slices%ready=.false. else slices%index=slices%index+1 slices%yz =abs(f(ix_loc,m1:m2 ,n1:n2 ,iFradx-1+slices%index)) slices%xz =abs(f(l1:l2 ,iy_loc,n1:n2 ,iFradx-1+slices%index)) slices%xy =abs(f(l1:l2 ,m1:m2 ,iz_loc ,iFradx-1+slices%index)) slices%xy2=abs(f(l1:l2 ,m1:m2 ,iz2_loc,iFradx-1+slices%index)) if (slices%index<=3) slices%ready = .true. endif ! endselect ! endsubroutine get_slices_radiation !*********************************************************************** subroutine calc_rad_diffusion(f,p) ! ! radiation in the diffusion approximation ! ! 12-apr-06/natalia: adapted from Wolfgang's more complex version ! 3-nov-06/axel: included gradient of conductivity, gradK.gradT ! use Sub, only: max_mn_name,dot use Cdata use Cparam use EquationOfState ! real, dimension (mx,my,mz,mvar+maux) :: f type (pencil_case) :: p real, dimension (mx,my,mz,mvar) :: df real, dimension (nx,3) :: glnThcond,duu real, dimension (nx) :: Krad,chi_rad,g2,diffus_chi1 real, dimension (nx) :: local_optical_depth,opt_thin,opt_thick,dt1_rad real :: fact, rho_max=0. !, dl_max=0. integer :: j,k, zone_size intent(inout) :: f,df intent(in) :: p ! ! calculate diffusion coefficient, Krad=16*sigmaSB*T^3/(3*kappa*rho) ! fact=16.*sigmaSB/(3.*kappa_es) Krad=fact*p%TT**3*p%rho1 ! ! calculate Qrad = div(K*gradT) = KT*[del2lnTT + (4*glnTT-glnrho).glnTT] ! Note: this is only correct for constant kappa (and here kappa=kappa_es) ! if (lrad_cool_diffus.and.lcooling) then call dot(4*p%glnTT-p%glnrho,p%glnTT,g2) f(l1:l2,m,n,iQrad)=Krad*p%TT*(p%del2lnTT+g2) endif ! ! radiative flux, Frad = -K*gradT; note that -div(Frad)=Qrad ! if (lrad_pres_diffus.and.lradflux) then do j=1,3 k=iFrad+(j-1) f(l1:l2,m,n,k)=-Krad*p%TT*p%glnTT(:,j) enddo endif ! ! include constraint from radiative time step ! (has to do with radiation pressure waves) ! if (lfirst.and.ldt) then advec_crad2=(16./3.)*p%rho1*(sigmaSB/c_light)*p%TT**4 endif ! ! check maximum diffusion from thermal diffusion ! With heat conduction, the second-order term for leading entropy term ! is gamma*chi_rad*del2ss ! if (lfirst.and.ldt) then ! ! Calculate timestep limitation. In the diffusion approximation the ! time step is just the diffusive time step, but with full radiation ! transfer there is a similar constraint resulting from the finite ! propagation speed of the radiation field, Vrad=Frad/Egas, which limits ! the time step to dt_rad = Vrad/lphoton, where lphoton=1/(kappa*rho), ! Frad=sigmaSB*T^4, and Egas=rho*cv*T. (At the moment we use cp.) ! cdtrad=0.8 is an empirical coefficient (harcoded for the time being) ! chi_rad=Krad*p%rho1*p%cp1 if (lrad_cool_diffus .and. lrad_pres_diffus) then diffus_chi=max(diffus_chi,gamma*chi_rad*dxyz_2) else ! ! calculate switches for optically thin/thick regions ! (should really make dxmax a pencil) ! local_optical_depth=dxmax*f(l1:l2,m,n,ikapparho) opt_thick=sign(.5,local_optical_depth-1.)+.5 opt_thin=1.-opt_thick ! !-- rho_max=maxval(p%rho) !-- dl_max=max(dx,dy,dz) !-- if (dxmax .GT. sqrt(gamma)/(rho_max*kappa_es)) then !-- diffus_chi=max(diffus_chi,gamma*chi_rad*dxyz_2) !-- else !-- diffus_chi1=min(gamma*chi_rad*dxyz_2, & !-- real(sigmaSB*kappa_es*p%TT**3*p%cv1/cdtrad)) !-- diffus_chi1=real(sigmaSB*kappa_es*p%TT**3*p%cv1/cdtrad) !-- diffus_chi=max(diffus_chi,diffus_chi1) dt1_rad=opt_thin*sigmaSB*kappa_es*p%TT**3*p%cv1/cdtrad_thin dt1_max=max(dt1_max,dt1_rad) diffus_chi=max(diffus_chi,opt_thick*gamma*chi_rad*dxyz_2/cdtrad_thick) !-- endif endif endif ! ! check radiative time step ! if (lfirst.and.ldt) then if (ldiagnos) then if (idiag_dtchi/=0) then call max_mn_name(diffus_chi/cdtv,idiag_dtchi,l_dt=.true.) endif if (idiag_dtrad/=0) then call max_mn_name(dt1_rad,idiag_dtrad,l_dt=.true.) endif endif endif ! endsubroutine calc_rad_diffusion !*********************************************************************** endmodule Radiation