! $Id$ ! ! 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 ! ! Note that Qrad is the heating rate (Q=I-S) and *not* the cooling rate used ! in Heinemann et al. 2006. ! Furthermore, Frad is the radiative flux multiplied by kappa*rho, ! not actually the radiative flux itself. ! ! TODO: Calculate weights properly ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! CPARAM logical, parameter :: lradiation = .true. ! ! MVAR CONTRIBUTION 0 ! MAUX CONTRIBUTION 2 ! !*************************************************************** module Radiation ! use Cparam use Cdata 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 ! integer, parameter :: mnu=2 integer, parameter :: maxdir=26 ! real, dimension (mx,my,mz) :: Srad, tau=0, Qrad=0, Qrad0=0 real, dimension (nx,ny,nz) :: Srad_noghost, kapparho_noghost real, dimension (mx,my) :: Irad_refl_xy real, target, dimension (:,:,:), allocatable :: Jrad_xy real, target, dimension (:,:,:), allocatable :: Jrad_xy2 real, target, dimension (:,:,:), allocatable :: Jrad_xy3 real, target, dimension (:,:,:), allocatable :: Jrad_xy4 real, target, dimension (:,:,:), allocatable :: Jrad_xz real, target, dimension (:,:,:), allocatable :: Jrad_yz real, target, dimension (:,:,:), allocatable :: Jrad_xz2 real, target, dimension (:,:,:,:,:,:), allocatable :: Jrad_r type (radslice), dimension (maxdir), target :: Isurf real, dimension (maxdir,3) :: unit_vec real, dimension (maxdir) :: weight, weightn, mu real, dimension (mnu) :: scalefactor_Srad=1.0, scalefactor_kappa=1.0 real, dimension (mnu) :: kappa_cst=1.0, kappa20_cst=0.0 real, dimension (:), allocatable :: lnTT_table real, dimension (:,:), allocatable :: lnSS_table real :: arad real :: dtau_thresh_min, dtau_thresh_max real :: tau_top=0.0, TT_top=0.0 real :: tau_bot=0.0, TT_bot=0.0 real :: kapparho_cst=1.0, kappa_Kconst=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=0.1, cdtrad_thin=1.0, cdtrad_thick=0.25, cdtrad_cgam=0.25 real :: scalefactor_cooling=1.0, scalefactor_radpressure=1.0 real :: expo_rho_opa=0.0, expo_temp_opa=0.0, expo_temp_opa_buff=0.0 real :: expo2_rho_opa=0.0, expo2_temp_opa=0.0 real :: ref_rho_opa=1.0, ref_temp_opa=1.0 real :: knee_temp_opa=0.0, width_temp_opa=1.0 real :: ampl_Isurf=0.0, radius_Isurf=0.0 real :: lnTT_table0=0.0, dlnTT_table=0.0, kapparho_floor=0.0 real :: TT_bump=0.0, sigma_bump=1.0, ampl_bump=0.0 real :: z_cutoff=impossible,cool_wid=impossible, qrad_max=0.0, & zclip_dwn=-max_real,zclip_up=max_real,kappa_ceiling=max_real ! integer :: radx=0, rady=0, radz=1, rad2max=1, nnu=1 integer, dimension (maxdir,3) :: dir integer, dimension (3) :: single_ray=0 integer :: lrad, mrad, nrad, rad2 integer :: idir, ndir integer :: l integer :: llstart, llstop, ll1, ll2, lsign integer :: mmstart, mmstop, mm1, mm2, msign integer :: nnstart, nnstop, nn1, nn2, nsign integer :: ipzstart, ipzstop, ipystart, ipystop, ipxstart, ipxstop integer :: nIsurf=1 integer :: nlnTT_table=1 ! logical :: lperiodic_ray, lperiodic_ray_x, lperiodic_ray_y logical :: lfix_radweight_1d=.true. logical :: lcooling=.true., lrad_debug=.false. logical :: lno_rad_heating=.false. logical :: lintrinsic=.true., lcommunicate=.true., lrevision=.true. logical :: lradpressure=.false., lradflux=.false., lsingle_ray=.false. logical :: lrad_cool_diffus=.false., lrad_pres_diffus=.false. logical :: lcheck_tau_division=.false., lread_source_function=.false. logical :: lcutoff_opticallythin=.false.,lcutoff_zconst=.false. logical :: lcdtrad_old=.true. ! character (len=2*bclen+1), dimension(3) :: bc_rad=(/'0:0','0:0','S:0'/) character (len=bclen), dimension(3) :: bc_rad1, bc_rad2 character (len=bclen) :: bc_ray_x, bc_ray_y, bc_ray_z character (len=labellen) :: source_function_type='LTE', opacity_type='Hminus' character (len=labellen) :: angle_weight='corrected' character :: lrad_str, mrad_str, nrad_str character (len=3) :: raydir_str ! 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 ! integer :: idiag_Qradrms=0, idiag_Qradmax=0 integer :: idiag_Fradzm=0, idiag_kapparhom=0, idiag_Sradm=0 integer :: idiag_Fradzmz=0, idiag_kapparhomz=0, idiag_taumz=0 integer :: idiag_dtchi=0, idiag_dtrad=0 integer :: ivid_Jrad=0, ivid_Isurf=0 ! namelist /radiation_init_pars/ & radx, rady, radz, rad2max, bc_rad, lrad_debug, kapparho_cst, & kappa_cst, kappa20_cst, & TT_top, TT_bot, tau_top, tau_bot, source_function_type, opacity_type, & nnu, lsingle_ray, single_ray, Srad_const, amplSrad, radius_Srad, nIsurf, & kappa_Kconst, kapparho_const, amplkapparho, radius_kapparho, lintrinsic, & lcommunicate, lrevision, lradflux, Frad_boundary_ref, lrad_cool_diffus, & lrad_pres_diffus, scalefactor_Srad, scalefactor_kappa, & angle_weight, lcheck_tau_division, & lfix_radweight_1d, expo_rho_opa, expo_temp_opa, expo_temp_opa_buff, & expo2_rho_opa, expo2_temp_opa, & ref_rho_opa, ref_temp_opa, knee_temp_opa, width_temp_opa, & lread_source_function, kapparho_floor,lcutoff_opticallythin, & lcutoff_zconst,z_cutoff,cool_wid, TT_bump, sigma_bump, ampl_bump, & kappa_ceiling ! namelist /radiation_run_pars/ & radx, rady, radz, rad2max, bc_rad, lrad_debug, kapparho_cst, & kappa_cst, kappa20_cst, & TT_top, TT_bot, tau_top, tau_bot, source_function_type, opacity_type, & nnu, lsingle_ray, single_ray, Srad_const, amplSrad, radius_Srad, nIsurf, & kx_Srad, ky_Srad, kz_Srad, kx_kapparho, ky_kapparho, kz_kapparho, & kappa_Kconst, kapparho_const, amplkapparho, radius_kapparho, lintrinsic, & lcommunicate, lrevision, lcooling, lradflux, lradpressure, & Frad_boundary_ref, lrad_cool_diffus, lrad_pres_diffus, lcdtrad_old, & cdtrad, cdtrad_thin, cdtrad_thick, cdtrad_cgam, & scalefactor_Srad, scalefactor_kappa, & angle_weight, & lcheck_tau_division, lfix_radweight_1d, expo_rho_opa, expo_temp_opa, & expo2_rho_opa, expo2_temp_opa, & ref_rho_opa, expo_temp_opa_buff, ref_temp_opa, knee_temp_opa, & width_temp_opa, ampl_Isurf, radius_Isurf, & scalefactor_cooling, scalefactor_radpressure, & lread_source_function, kapparho_floor, lcutoff_opticallythin, & lcutoff_zconst,z_cutoff,cool_wid,lno_rad_heating,qrad_max,zclip_dwn, & zclip_up, TT_bump, sigma_bump, ampl_bump, kappa_ceiling ! contains !*********************************************************************** subroutine register_radiation ! ! Initialize radiation flags. ! ! 24-mar-03/axel+tobi: coded ! use FArrayManager ! lradiation_ray=.true. ! ! Set indices for auxiliary variables. ! call farray_register_auxiliary('Qrad',iQrad) call farray_register_auxiliary('kapparho',ikapparho) ! ! Allocated auxiliary arrays for radiative flux only if lradflux=T ! Remember putting "! MAUX CONTRIBUTION 3" (or adding 3) in cparam.local! ! if (lradflux) then call farray_register_auxiliary('KR_Frad',iKR_Frad,vector=3) iKR_Fradx = iKR_Frad iKR_Frady = iKR_Frad+1 iKR_Fradz = iKR_Frad+2 endif ! ! Identify version number (generated automatically by SVN). ! 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)=',kapparho $' aux_count=aux_count+1 ! ! Frad is used only if lradflux=.true. ! if (lradflux) then if (naux < maux) aux_var(aux_count)=',KR_Frad $' if (naux == maux) aux_var(aux_count)=',KR_Frad' aux_count=aux_count+3 endif ! ! Output. ! if (lroot) then write(15,*) 'Qrad = fltarr(mx,my,mz)*one' write(15,*) 'kapparho = fltarr(mx,my,mz)*one' if (lradflux) write(15,*) 'KR_Frad = 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. ! ! 16-jun-03/axel+tobi: coded ! 03-jul-03/tobi: position array added ! 19-feb-14/axel: read tabulated source function ! use Sub, only: parse_bc_rad use Slices_methods, only: alloc_slice_buffers ! integer :: itable real :: radlength,arad_normal logical :: periodic_xy_plane,bad_ray,ray_good character (len=labellen) :: header ! ! Check that the number of rays does not exceed maximum. ! if (radx>1) call fatal_error('initialize_radiation', & 'radx currently must not be greater than 1') if (rady>1) call fatal_error('initialize_radiation', & 'rady currently must not be greater than 1') if (radz>1) call fatal_error('initialize_radiation', & 'radz currently must not be greater than 1') ! ! Empirically we have found that cdtrad>0.1 is unsafe. ! But in Perri & Brandenburg, for example, cdtrad=0.5 was still ok. ! It may be useful to define 2 different cdtrad for the cases below. ! !if (ldt.and.cdtrad>0.1) then ! call fatal_error('initialize_radiation', & ! 'cdtrad is larger than 0.1 - do you really want this?') !endif ! ! Check boundary conditions. ! if (lroot.and.ip<14) print*,'initialize_radiation: bc_rad =',bc_rad ! call parse_bc_rad(bc_rad,bc_rad1,bc_rad2) ! ! Count. ! idir=1 ! 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 ! ! Check whether the horizontal plane is fully periodic. ! periodic_xy_plane=all(bc_rad1(1:2)=='p').and.all(bc_rad2(1:2)=='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) ! ! For single rays: single_ray must match (lrad,mrad,nrad). ! Otherwise, check for length of rays. ! if (lsingle_ray) then ray_good=(lrad==single_ray(1) & .and.mrad==single_ray(2) & .and.nrad==single_ray(3) ) else ray_good=(rad2>0.and.rad2<=rad2max).and.(.not.bad_ray) endif ! ! Proceed with good rays. ! if (ray_good) 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 ! ! Print all unit vectors (if ray_good=T). ! if (lroot.and.ip<14) write(*,'(1x,a,i4,3f6.2)') & 'initialize_radiation: idir, unit_vec=', & idir, unit_vec(idir,:) ! idir=idir+1 ! endif ! enddo enddo enddo ! ! Total number of directions; correct for latest idir+1 operation. ! 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 ! if (source_function_type=='LTE') then arad=sigmaSB/pi arad_normal=4*sigmaSB/c_light else arad=0.0 arad_normal=0.0 endif ! ! Debug output. ! NOTE: arad is only used when S=(c/4pi)*aT^4=(sigmaSB/pi)*T^4 ! if (lroot.and.ip<9) then print*, 'initialize_radiation: arad=', arad print*, 'initialize_radiation: arad_normal=', arad_normal print*, 'initialize_radiation: sigmaSB=', sigmaSB endif ! ! Write constants to disk. ! if (lroot) then open (1,file=trim(datadir)//'/pc_constants.pro',position="append") write (1,*) 'arad_normal=',arad_normal write (1,*) 'arad=',arad write (1,*) 'c_light=',c_light write (1,*) 'sigmaSB=',sigmaSB write (1,*) 'kappa_es=',kappa_es close (1) endif ! ! Calculate weights. ! call calc_angle_weights ! if (lroot.and.ip<=14) print*, 'initialize_radiation: ndir =', ndir ! ! Read tabulated source function. ! if (lread_source_function) then open (1,file='nnu2.dat') read (1,*) nlnTT_table read (1,*) header allocate(lnTT_table(nlnTT_table),lnSS_table(nlnTT_table,nnu)) do itable=1,nlnTT_table read(1,*) lnTT_table(itable),lnSS_table(itable,:) enddo close (1) lnTT_table0=lnTT_table(1) dlnTT_table=(nlnTT_table-1)/(lnTT_table(nlnTT_table)-lnTT_table(1)) endif ! if (ivid_Jrad/=0) & call alloc_slice_buffers(Jrad_xy,Jrad_xz,Jrad_yz,Jrad_xy2,Jrad_xy3,Jrad_xy4,Jrad_xz2,ncomp=nnu) endsubroutine initialize_radiation !*********************************************************************** subroutine calc_angle_weights ! ! This subroutine calculates the weights needed for the discrete version ! of the angular integration of the radiative transfer equation. By now ! the code was using only constant weights, which just works for very special ! cases. Using spherical harmonics just works for some special cases as well, ! but in a slightly broader range. ! ! 18-may-07/wlad: coded ! real :: xyplane,yzplane,xzplane,room,xaxis real :: yaxis,zaxis,aspect_ratio,mu2,correction_factor=1. ! ! Calculate weights for weighed integrals involving one unit vector nhat ! select case (angle_weight) ! ! Weight factors corrected for dimensionality less than 3. ! In 1-D, we need to correct by a factor 1/3, while in 2-D by 2/3. ! Thus, we multiply by the number of dimensions (=radx+rady+radz) ! and divide by 3. This has been verified in 1-D with 2 rays, ! in 2-D with 4 (xy periodic) or 8 (xz semi-periodic) rays, ! as well as in 3-D with 14 or 22 rays; for details, see appendix ! of Barekat & Brandenburg (2014, Astron. Astrophys. 571, A68). ! case ('corrected') correction_factor=(radx+rady+radz)/3. if (ndir>0) weight=(4*pi/ndir)*correction_factor if (radx>1.or.rady>1.or.radz>1) call fatal_error( & 'calc_angle_weights','radx, rady, and radz cannot be larger than 1') if (lroot.and.ip<=14) print*, & 'calc_angle_weights: correction_factor =', correction_factor weightn=weight ! ! Constant weight factors, ignore dimensionality, which is wrong, ! but we keep it for now for compatibility reasons. ! case ('constant') if (ndir>0) weight=4*pi/ndir weightn=weight ! if (lfix_radweight_1d.and.ndir==2) weightn=weightn/3. ! ! Spherical harmonics calculation of weight factors; see the manual. ! 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 fatal_error('initialize_radiation', & 'weights not calculated for dx/dy != 1') endif aspect_ratio=dx/dz if (aspect_ratio<0.69.or.aspect_ratio>sqrt(3.0)) & call fatal_error('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 fatal_error('calc_angle_weights', & '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. ! ! 16-jun-03/axel+tobi: coded ! 5-dec-13/axel: alterations to allow non-gray opacities ! real, dimension(mx,my,mz,mfarray) :: f ! integer :: j,k,inu ! ! Identifier. ! if (ldebug.and.headt) print*, 'radtransfer' ! ! Continue only if we either have more than a single ray, or, if we do have a ! single ray, when also lvideo.and.lfirst are true so that the result is used ! for visualization. ! if ((.not.lsingle_ray) .or. (lsingle_ray.and.lvideo.and.lfirst)) then ! ! Do loop over all frequency bins. ! do inu=1,nnu ! ! Calculate source function and opacity. ! call source_function(f,inu) call opacity(f,inu) ! ! Do the rest only if we do not do diffusion approximation. ! If *either* lrad_cool_diffus *or* lrad_pres_diffus, no rays are computed. ! if (lrad_cool_diffus.or.lrad_pres_diffus) then if (headt) print*, 'radtransfer: do diffusion approximation, no rays' else ! ! Initialize heating rate and radiative flux only at first frequency point. ! if (inu==1) then f(:,:,:,iQrad)=0.0 if (lradflux) f(:,:,:,iKR_Fradx:iKR_Fradz)=0.0 endif ! ! Loop over rays. ! do idir=1,ndir call raydirection ! ! Do the 3 steps: first intrinsic (compute intensive), ! then communication (not compute intensive), ! and finally revision (again compute intensive). ! if (lintrinsic) call Qintrinsic(f) ! 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. ! Turn this now into heating rate by multiplying with opacity. ! This allows the opacity to be frequency-dependent. ! f(:,:,:,iQrad)=f(:,:,:,iQrad)+weight(idir)*Qrad*f(:,:,:,ikapparho) ! ! Calculate radiative flux. Multiply it here by opacity to have the correct ! frequency-dependent contributions from all frequencies. ! if (lradflux) then do j=1,3 k=iKR_Frad+(j-1) f(:,:,:,k)=f(:,:,:,k)+weightn(idir)*unit_vec(idir,j) & *(Qrad+Srad)*f(:,:,:,ikapparho) enddo endif ! ! Store outgoing intensity in case of lower reflective boundary condition. ! We need to set Iup=Idown+I0 at the lower boundary. We must first integrate ! Idown into the lower boundary following: ! ! Idown(n1-1) = Idown(n1) - dtau*Qdown(n1) ! ! [using simply Idown(n1) as the outgoing intensity is not consistent] ! ! This corresponds to ! ! Idown(n1-1) = S(n1) + (1-dtau)*Qdown(n1) ! ! where dtau=kappa*rho(n1)*dz. ! if ((bc_rad1(3)=='R').or.(bc_rad1(3)=='R+F')) then if ((ndir==2).and.(idir==1).and.(ipz==ipzstop)) & Irad_refl_xy=Srad(:,:,nnstop)+Qrad(:,:,nnstop)* & (1.0-f(:,:,nnstop,ikapparho)/dz_1(nnstop)) endif ! ! enddo from idir. ! enddo ! ! End of no-diffusion approximation query. ! endif ! ! Calculate slices of J=S+Q/(4pi). ! if (lvideo.and.lfirst.and.ivid_Jrad/=0) then if (lwrite_slice_yz) & Jrad_yz(:,:,inu)= Qrad(ix_loc,m1:m2,n1:n2) +Srad(ix_loc,m1:m2,n1:n2) if (lwrite_slice_xz) & Jrad_xz(:,:,inu)= Qrad(l1:l2,iy_loc,n1:n2) +Srad(l1:l2,iy_loc,n1:n2) if (lwrite_slice_xz2) & Jrad_xz2(:,:,inu)=Qrad(l1:l2,iy2_loc,n1:n2)+Srad(l1:l2,iy2_loc,n1:n2) if (lwrite_slice_xy) & Jrad_xy(:,:,inu)= Qrad(l1:l2,m1:m2,iz_loc) +Srad(l1:l2,m1:m2,iz_loc) if (lwrite_slice_xy2) & Jrad_xy2(:,:,inu)=Qrad(l1:l2,m1:m2,iz2_loc)+Srad(l1:l2,m1:m2,iz2_loc) if (lwrite_slice_xy3) & Jrad_xy3(:,:,inu)=Qrad(l1:l2,m1:m2,iz3_loc)+Srad(l1:l2,m1:m2,iz3_loc) if (lwrite_slice_xy4) & Jrad_xy4(:,:,inu)=Qrad(l1:l2,m1:m2,iz4_loc)+Srad(l1:l2,m1:m2,iz4_loc) endif ! ! End of frequency loop (inu). ! enddo ! ! endif from single ray check ! endif ! endsubroutine radtransfer !*********************************************************************** subroutine raydirection ! ! Determine certain variables depending on the ray direction. ! ! 10-nov-03/tobi: coded ! 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=bc_rad1(1) if (lrad<0) bc_ray_x=bc_rad2(1) if (mrad>0) bc_ray_y=bc_rad1(2) if (mrad<0) bc_ray_y=bc_rad2(2) if (nrad>0) bc_ray_z=bc_rad1(3) if (nrad<0) bc_ray_z=bc_rad2(3) ! ! 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 (lrad>0) then; ipxstart=0; ipxstop=nprocx-1; endif if (lrad<0) then; ipxstart=nprocx-1; ipxstop=0; endif 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(f) ! ! Integration radiation transfer equation along rays. ! ! This routine is called before the communication part. ! All rays start with zero intensity. ! ! 16-jun-03/axel+tobi: coded ! 3-aug-03/axel: added max(dtau,dtaumin) construct ! use Debug_IO, only: output ! real, dimension(mx,my,mz,mfarray), intent(in) :: f real,dimension(mz) :: dlength real :: Srad1st,Srad2nd,emdtau1,emdtau2,emdtau real :: dtau_m,dtau_p,dSdtau_m,dSdtau_p ! ! Identifier. ! if (ldebug.and.headt) print*,'Qintrinsic' ! ! Line elements. ! if (nzgrid/=1) then dlength=sqrt((dx*lrad)**2+(dy*mrad)**2+(nrad/dz_1)**2) else dlength=sqrt((dx*lrad)**2+(dy*mrad)**2+(nrad/dz)**2) endif ! ! Set optical depth and intensity initially to zero. ! tau=0.0 Qrad=0.0 ! ! Loop over all meshpoints. ! do n=nnstart,nnstop,nsign do m=mmstart,mmstop,msign do l=llstart,llstop,lsign ! dtau_m=sqrt(f(l-lrad,m-mrad,n-nrad,ikapparho)* & f(l,m,n,ikapparho))*0.5*(dlength(n-nrad)+dlength(n)) dtau_p=sqrt(f(l,m,n,ikapparho)* & f(l+lrad,m+mrad,n+nrad,ikapparho))* & 0.5*(dlength(n)+dlength(n+nrad)) ! ! Avoid divisions by zero when the optical depth is such. ! dtau_m = max(dtau_m,epsi) dtau_p = max(dtau_p,epsi) ! 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 ! 3-oct-20/axel: initializing Qrecv, Qrad, and Qsend to zero ! 15-apr-22/felipe: fixed the communication on the x and y directions ! use Mpicomm, only: radboundary_xy_recv,radboundary_xy_send use Mpicomm, only: radboundary_yz_sendrecv,radboundary_zx_sendrecv use Mpicomm, only: radboundary_yz_send, radboundary_yz_recv use Mpicomm, only: radboundary_zx_send, radboundary_zx_recv ! real, dimension (my,mz) :: emtau_yz=0, Qrad_yz=0 real, dimension (mx,mz) :: emtau_zx=0, Qrad_zx=0 real, dimension (mx,my) :: emtau_xy=0, Qrad_xy=0 real, dimension (my,mz) :: Qrecv_yz=0, Qsend_yz=0 real, dimension (mx,mz) :: Qrecv_zx=0, Qsend_zx=0 real, dimension (mx,my) :: Qrecv_xy=0, Qsend_xy=0 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 (ipx==ipxstart) then if (bc_ray_x/='p') then call radboundary_yz_set(Qrecv_yz) else call radboundary_yz_recv(lrad,idir,Qrecv_yz) endif ! Qbc_yz(mm1:mm2,nnstart-nrad)%val = Qrecv_yz(llstart-lrad,mm1:mm2) ! Qbc_yz(mm1:mm2,nnstart-nrad)%set = .true. ! else call radboundary_yz_recv(lrad,idir,Qrecv_yz) 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 ! ! Copy the above heating rates to the yz-target arrays which are then set. ! 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 all_yz=.true. endif ! if (mrad/=0) then if (ipy==ipystart) then if (bc_ray_y/='p') then call radboundary_zx_set(Qrecv_zx) else call radboundary_zx_recv(mrad,idir,Qrecv_zx) endif ! Qbc_yz(mm1:mm2,nnstart-nrad)%val = Qrecv_yz(llstart-lrad,mm1:mm2) ! Qbc_yz(mm1:mm2,nnstart-nrad)%set = .true. ! else call radboundary_zx_recv(mrad,idir,Qrecv_zx) 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 ! ! Copy the above heating rates to the zx-target arrays which are then set. ! 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 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 ! ! x-direction. ! 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) Qsend_yz(m,n) = Qpt_yz(m,n)%val*emtau_yz(m,n)+Qrad_yz(m,n) endforall ! if (nprocx>1) then call radboundary_yz_sendrecv(lrad,idir,Qsend_yz,Qrecv_yz) else Qrecv_yz=Qsend_yz endif ! forall (m=mm1:mm2,n=nn1:nn2,Qpt_yz(m,n)%set.and..not.Qbc_yz(m,n)%set) Qbc_yz(m,n)%val = Qrecv_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 ! ! y-direction. ! 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 if (lrad/=0.and.ipx/=ipxstop) then forall (m=mm1:mm2,n=nn1:nn2) emtau_yz(m,n) = exp(-tau(llstop,m,n)) Qrad_yz(m,n) = Qrad(llstop,m,n) Qsend_yz(m,n) = Qpt_yz(m,n)%val*emtau_yz(m,n)+Qrad_yz(m,n) endforall ! call radboundary_yz_send(lrad,idir,Qsend_yz) endif if (mrad/=0.and.ipy/=ipystop) then forall (l=ll1:ll2,n=nn1:nn2) emtau_zx(l,n) = exp(-tau(l,mmstop,n)) Qrad_zx(l,n) = Qrad(l,mmstop,n) Qsend_zx(l,n) = Qpt_zx(l,n)%val*emtau_zx(l,n)+Qrad_zx(l,n) endforall ! call radboundary_zx_send(mrad,idir,Qsend_zx) endif endsubroutine Qcommunicate !*********************************************************************** subroutine Qperiodic ! ! This routine deals with communication for horizontal rays. ! use Mpicomm, only: radboundary_yz_periodic_ray,radboundary_zx_periodic_ray use Debug_IO, only: output ! real, dimension(ny,nz) :: Qrad_yz,tau_yz real, dimension(nx,nz) :: Qrad_zx,tau_zx real, dimension(ny,nz) :: Qrad_tot_yz,tau_tot_yz,emtau1_tot_yz real, dimension(nx,nz) :: Qrad_tot_zx,tau_tot_zx,emtau1_tot_zx real, dimension(ny,nz,0:nprocx-1) :: Qrad_yz_all,tau_yz_all real, dimension(nx,nz,0:nprocy-1) :: Qrad_zx_all,tau_zx_all integer :: ipxstart,ipxstop,ipl integer :: ipystart,ipystop,ipm ! ! x-direction ! if (lrad/=0) then ! ! Intrinsic heating rate and optical depth at the downstream boundary of ! each processor. ! Qrad_yz=Qrad(llstop,m1:m2,n1:n2) tau_yz=tau(llstop,m1:m2,n1:n2) ! ! Gather intrinsic heating rates and optical depths from all processors ! into one rank-3 array available on each processor. ! call radboundary_yz_periodic_ray(Qrad_yz,tau_yz,Qrad_yz_all,tau_yz_all) ! ! Find out in which direction we want to loop over processors. ! if (lrad>0) then; ipxstart=0; ipxstop=nprocx-1; endif if (lrad<0) then; ipxstart=nprocx-1; ipxstop=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_yz=0.0 tau_tot_yz=0.0 ! ! Do the sum from this processor to the last one in the downstream direction. ! do ipl=ipx,ipxstop,lsign Qrad_tot_yz=Qrad_tot_yz*exp(-tau_yz_all(:,:,ipl))+Qrad_yz_all(:,:,ipl) tau_tot_yz=tau_tot_yz+tau_yz_all(:,:,ipl) enddo ! ! Do the sum from the first processor in the upstream direction to the one ! before this one. ! do ipl=ipxstart,ipx-lsign,lsign Qrad_tot_yz=Qrad_tot_yz*exp(-tau_yz_all(:,:,ipl))+Qrad_yz_all(:,:,ipl) tau_tot_yz=tau_tot_yz+tau_yz_all(:,:,ipl) 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_yz>dtau_thresh_max) emtau1_tot_yz=1.0 elsewhere (tau_tot_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 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) & Isurf(idir)%xy2=Qrad(l1:l2,m1:m2,nnstop)+Srad(l1:l2,m1:m2,nnstop) endif ! if (lrad_debug) & call output(trim(directory_dist)//'/Qrev-'//raydir_str//'.dat',Qrad,1) ! ! xy-averages ! if (l1davgfirst) then do n=n1,n2 do m=m1,m2 call xysum_mn_name_z(tau(l1:l2,m,n),idiag_taumz) enddo enddo 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 ! real, dimension(my,mz), intent(out) :: Qrad0_yz ! ! No incoming intensity. ! if (bc_ray_x=='0') then Qrad0_yz=-Srad(llstart-lrad,:,:) endif ! ! Set intensity equal to unity (as a numerical experiment for now). ! if (bc_ray_x=='1') then Qrad0_yz=1.0-Srad(llstart-lrad,:,:) endif ! ! Set intensity equal to source function. ! if (bc_ray_x=='S') then Qrad0_yz=0 endif ! ! Set intensity equal to S-F. ! if (bc_ray_x=='S-F') then Qrad0_yz=-Frad_boundary_ref/(2*weightn(idir)) endif ! ! Set intensity equal to S+F. ! if (bc_ray_x=='S+F') then Qrad0_yz=+Frad_boundary_ref/(2*weightn(idir)) endif ! ! Set intensity equal to a pre-defined incoming intensity. ! if (bc_ray_x=='F') then Qrad0_yz=-Srad(llstart-lrad,:,:)+Frad_boundary_ref/(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 ! real, dimension(mx,mz), intent(out) :: Qrad0_zx ! ! No incoming intensity. ! if (bc_ray_y=='0') then Qrad0_zx=-Srad(:,mmstart-mrad,:) endif ! ! Set intensity equal to unity (as a numerical experiment for now). ! if (bc_ray_y=='1') then Qrad0_zx=1.0-Srad(:,mmstart-mrad,:) endif ! ! Set intensity equal to source function. ! if (bc_ray_y=='S') then Qrad0_zx=0 endif ! ! Set intensity equal to S-F. ! if (bc_ray_y=='S-F') then Qrad0_zx=-Frad_boundary_ref/(2*weightn(idir)) endif ! ! Set intensity equal to S+F. ! if (bc_ray_y=='S+F') then Qrad0_zx=+Frad_boundary_ref/(2*weightn(idir)) endif ! ! Set intensity equal to a pre-defined incoming intensity. ! if (bc_ray_y=='F') then Qrad0_zx=-Srad(:,mmstart-mrad,:)+Frad_boundary_ref/(2*weightn(idir)) 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 ! real, dimension(mx,my), intent(out) :: Qrad0_xy real :: Irad_xy, fact ! ! No incoming intensity. ! if (bc_ray_z=='0') then Qrad0_xy=-Srad(:,:,nnstart-nrad) endif ! ! Set intensity equal to unity (as a numerical experiment for now). ! if (bc_ray_z=='1') then Qrad0_xy=1.0-Srad(:,:,nnstart-nrad) endif ! ! Set intensity equal to a gaussian (as a numerical experiment). ! if (bc_ray_z=='blo') then fact=.5/radius_Isurf**2 Qrad0_xy=ampl_Isurf*exp(-fact*(spread(x**2,2,my)+spread(y**2,1,mx))) endif ! ! Incoming intensity from a layer of constant temperature TT_top. ! if (bc_ray_z=='c') then if (nrad<0) then Irad_xy=arad*TT_top**4*(1-exp(tau_top/mu(idir))) elseif (nrad>0) then Irad_xy=arad*TT_bot**4*(1-exp(tau_bot/mu(idir))) else call fatal_error('radboundary_xy_set','nrad=0 should not happen') Irad_xy=impossible endif 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 outgoing intensity. ! if (bc_ray_z=='R') then Qrad0_xy=-Srad(:,:,nnstart-nrad)+Irad_refl_xy endif ! ! Set intensity equal to outgoing intensity plus pre-defined flux. ! if (bc_ray_z=='R+F') then Qrad0_xy=-Srad(:,:,nnstart-nrad)+ & Frad_boundary_ref/(2*weightn(idir))+Irad_refl_xy endif ! endsubroutine radboundary_xy_set !*********************************************************************** subroutine radiative_cooling(f,df,p) ! ! Add radiative cooling to entropy/temperature equation. ! ! 25-mar-03/axel+tobi: coded ! 5-dec-13/axel: removed opacity multiplicaton to allow non-gray opacities ! use Diagnostics use Sub ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! real, dimension (nx) :: cooling, kappa, dt1_rad, Qrad2 real, dimension (nx) :: cgam, ell, chi, dtrad_thick, dtrad_thin real, dimension (nx) :: dt1rad_cgam integer :: l ! ! Add radiative cooling, either from the intensity or in the diffusion ! approximation (if either lrad_cool_diffus=F or lrad_pres_diffus=F). ! if (lrad_cool_diffus.or.lrad_pres_diffus) call calc_rad_diffusion(f,p) if (lno_rad_heating .and. (qrad_max > 0)) then ! ! Upper limit radiative heating by qrad_max ! do l=l1-radx, l2+radx if (f(l,m,n,iqrad) .gt. qrad_max) & f(l,m,n,iQrad)=qrad_max enddo endif cooling=f(l1:l2,m,n,iQrad) ! ! Possibility of rescaling the radiative cooling term. ! if (scalefactor_cooling/=1.) then cooling=cooling*scalefactor_cooling 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 ! ! Time-step contribution from cooling. ! !if (lfirst.and.ldt) then ! ! Choose less stringent time-scale of optically thin or thick cooling. ! This is currently not correct in the non-gray case! ! Instead of a factor 4, one should have a factor of 16; see BB14, Eq.(A.2). ! kapparho^2 > dxyz_2 means short mean-free path, so optically thick. ! Then, dt1_min = 4*kappa*sigmaSB*T^3/(cv*dx^2*kapparho^2*cdtrad), so ! dt1_min = cgam*ell/(4*dx^2); otherwise, dt_min = cgam*ell/(4*ell). ! In the optically thin regime: no constraint for z > z_cutoff. ! kappa=f(l1:l2,m,n,ikapparho)*p%rho1 if (lcdtrad_old) then do l=1,nx if (f(l1-1+l,m,n,ikapparho)**2>dxyz_2(l)) then dt1_rad(l)=4*kappa(l)*sigmaSB*p%TT(l)**3*p%cv1(l)* & dxyz_2(l)/f(l1-1+l,m,n,ikapparho)**2/cdtrad if (z_cutoff/=impossible .and. cool_wid/=impossible) & dt1_rad(l)=0.5*dt1_rad(l)*(1.-tanh((z(n)-z_cutoff)/cool_wid)) else dt1_rad(l)=4*kappa(l)*sigmaSB*p%TT(l)**3*p%cv1(l)/cdtrad if (z_cutoff/=impossible .and. cool_wid/=impossible) & dt1_rad(l)=0.5*dt1_rad(l)*(1.-tanh((z(n)-z_cutoff)/cool_wid)) endif enddo else cgam=16.*sigmaSB*p%TT**3*p%rho1*p%cp1 ell=1./f(l1:l2,m,n,ikapparho) chi=cgam*ell/3. dtrad_thick=cdtrad_thick/(dxyz_2*chi*dimensionality) !dtrad_cgam=cdtrad_cgam/(sqrt(dxyz_2)*cgam) dt1rad_cgam=sqrt(dxyz_2)*cgam/cdtrad_cgam dtrad_thin=cdtrad_thin*ell/cgam dt1_rad=1./(dtrad_thick+dtrad_thin) endif dt1_max=max(dt1_max,dt1_rad) !dt1_max=max(dt1_max,dt1_rad,dt1rad_cgam) !endif endif if (lfirst.and.ldt) then if (idiag_dtrad/=0) & call max_mn_name(dt1_rad,idiag_dtrad,l_dt=.true.) endif ! ! Diagnostics. ! if (ldiagnos) then Qrad2=f(l1:l2,m,n,iQrad)**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 Diagnostics use Sub ! 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 ! ! Radiative pressure force = kappa*Frad/c = (kappa*rho)*rho1*Frad/c. ! Moved multiplication with opacity to earlier Frad calculation to ! allow opacities to be non-gray. This implies that Frad is not actually ! the radiative flux, but the radiative flux multiplied by kappa*rho. ! The auxiliary array is therefore now called KR_Frad. ! if (lradpressure) then do j=1,3 k=iKR_Frad+(j-1) radpressure(:,j)=scalefactor_radpressure*p%rho1*f(l1:l2,m,n,k)/c_light enddo df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)+radpressure endif ! ! Diagnostics. Here we divide again by kapparho, so we get actual Frad. ! if (ldiagnos) then if (lradflux) & call sum_mn_name(f(l1:l2,m,n,iKR_Fradz)/f(l1:l2,m,n,ikapparho),idiag_Fradzm) call sum_mn_name(f(l1:l2,m,n,ikapparho),idiag_kapparhom) endif ! if (l1davgfirst) then if (lradflux) & call xysum_mn_name_z(f(l1:l2,m,n,iKR_Fradz)/f(l1:l2,m,n,ikapparho),idiag_Fradzmz) call xysum_mn_name_z(f(l1:l2,m,n,ikapparho),idiag_kapparhomz) endif ! endsubroutine radiative_pressure !*********************************************************************** subroutine source_function(f,inu) ! ! Calculates source function. ! ! This module is currently ignored if diffusion approximation is used. ! ! 3-apr-04/tobi: coded ! 8-feb-09/axel: added B2 for visualisation purposes ! 5-dec-13/axel: alterations to allow non-gray opacities ! use EquationOfState, only: eoscalc use Debug_IO, only: output use SharedVariables, only: put_shared_variable use Mpicomm, only: stop_it ! real, dimension(mx,my,mz,mfarray), intent(in) :: f logical, save :: lfirst=.true. integer, dimension(mx) :: ilnTT_table real, dimension(mx,my) :: z_cutoff1 real, dimension(mx) :: lnTT integer :: lun_input = 1 integer :: inu integer :: ierr ! select case (source_function_type) ! ! usual Planck function, S=sigmaSB*TT^4 ! case ('LTE') if (lcutoff_opticallythin) then ! ! This works for stratification in the z-direction ! if (z_cutoff==impossible .or. cool_wid==impossible) & call fatal_error("source_function:","z_cutoff or cool_wid is not set") call put_shared_variable('z_cutoff',z_cutoff,ierr) if (ierr/=0) call stop_it("source_function: "//& "there was a problem when putting z_cutoff") call put_shared_variable('cool_wid',cool_wid,ierr) if (ierr/=0) call stop_it("source_function: "//& "there was a problem when putting cool_wid") z_cutoff1=z_cutoff if (.not. lcutoff_zconst) then do l=l1-radx,l2+radx do m=m1-rady,m2+rady do n=n1-radz,n2+radz ! ! Put Srad smoothly to zero for z above which the ! photon mean free path kappa*rho > 1/(1000*dz) ! if (abs(f(l,m,n,ikapparho)-1.0e-3*dz_1(n)) .lt. epsi) then z_cutoff1(l,m)=min(max(z(n),zclip_dwn),zclip_up) exit endif enddo enddo enddo endif do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnTT=lnTT) do l=l1-radx,l2+radx Srad(l,m,n)=arad*exp(4*lnTT(l))*scalefactor_Srad(inu)* & 0.5*(1.-tanh((z(n)-z_cutoff1(l,m))/cool_wid)) enddo enddo enddo else 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(inu) enddo enddo endif ! ! read tabulated 2-colored source function for inu=1 and 2. ! case ('two-colored') do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnTT=lnTT) ilnTT_table=max(min(int((lnTT-lnTT_table0)*dlnTT_table),nlnTT_table-1),1) Srad(:,m,n)=exp(lnSS_table(ilnTT_table,inu) & +(lnSS_table(ilnTT_table+1,inu)-lnSS_table(ilnTT_table,inu)) & *(lnTT-lnTT_table(ilnTT_table))/(lnTT_table(ilnTT_table+1)-lnTT_table(ilnTT_table)))/unit_flux enddo enddo ! ! for testing purposes, a blob-like spatial profile ! 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 ! ! for testing purposes, a cosine-like spatial profile. ! 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 ! ! Source function proportional to magnetic energy density ! (used for visualization purposes). ! Needs to be done at every time step (not just when lfirst=.true.). ! case ('B2') call calc_Srad_B2(f,Srad) ! ! Combination of magnetic energy density and enstrophy. ! case ('B2+W2') if (inu==1) then call calc_Srad_B2(f,Srad) elseif (inu==2) then call calc_Srad_W2(f,Srad) endif ! ! Read from file ! case ('read_file') open (lun_input, file=trim(directory_prestart)//'/Srad.dat', form='unformatted') read (lun_input) Srad_noghost close (lun_input) Srad(l1:l2,m1:m2,n1:n2)=Srad_noghost Srad(l1:l2,m1:m2,n1-1)=impossible Srad(l1:l2,m1:m2,n2+1)=impossible ! ! Nothing. ! case ('nothing') Srad=0.0 ! case default call fatal_error('source_function', & 'no such source function type: '//trim(source_function_type)) ! endselect ! if (lrad_debug) then call output(trim(directory_dist)//'/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. ! ! 3-apr-04/tobi: coded ! 8-feb-09/axel: added B2 for visualisation purposes ! 5-dec-13/axel: alterations to allow non-gray opacities ! use Debug_IO, only: output use EquationOfState, only: eoscalc use Sub, only: cubic_step ! real, dimension(mx,my,mz,mfarray), intent(inout) :: f real, dimension(mx) :: tmp,lnrho,lnTT,yH,rho,TT,profile=0 real, dimension(mx) :: kappa1,kappa2,kappae real, dimension(mx) :: kappa_rad,kappa_cond,kappa_tot real :: kappa0, kappa0_cgs,k1,k2 logical, save :: lfirst=.true. integer :: lun_input = 1 integer :: i,inu ! select case (opacity_type) ! case ('Hminus') do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,kapparho=tmp) f(:,m,n,ikapparho)=kapparho_floor+tmp*scalefactor_kappa(inu) enddo enddo ! case ('total_Rosseland_mean') ! ! All coefficients valid for cgs units ONLY ! We use solar abundances X=0.7381, Y=0.2485, metallicity Z=0.0134 ! kappa1 is Kramer's opacity, kappa2 is Hminus opacity, kappa_cond is conductive opacity ! total kappa is calculated by taking the harmonic mean of radiative nd conductive opacities ! kappa_rad=1/(1/kappa2+1/kappa1) ! kappa=1/(1/kappa_rad+1/kappa3) ! do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnrho=lnrho) call eoscalc(f,mx,lntt=lntt) kappa1=4.0d25*1.7381*0.0135*unit_density**2*unit_length* & exp(lnrho)*(exp(lnTT)*unit_temperature)**(-3.5) kappa2=1.25d-29*0.0134*unit_density**1.5*unit_length*& unit_temperature**9*exp(0.5*lnrho)*exp(9.0*lnTT) kappae=0.2*1.7381*(1.+2.7d11*exp(lnrho-2*lnTT)*unit_density/unit_temperature**2)**(-1.) kappa_cond=2.6d-7*unit_length*unit_temperature**2*exp(2*lnTT)*exp(-lnrho) kappa_rad=kapparho_floor+1./(1./(kappa1+kappae)+1./kappa2) kappa_tot=1./(1./kappa_rad+1./kappa_cond) if (lcutoff_opticallythin) & kappa_tot=0.5*(1.-tanh((z(n)-0.5*z_cutoff)/(2*cool_wid)))/ & (1./kappa_rad+1./kappa_cond) do i=1,mx kappa_tot(i)=min(kappa_tot(i),kappa_ceiling) enddo f(:,m,n,ikapparho)=exp(lnrho)*kappa_tot*scalefactor_kappa(inu) 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)=kapparho_floor+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)=kapparho_floor+kappa_cst(inu)*exp(lnrho) enddo enddo ! case ('kapparho_cst') do n=n1-radz,n2+radz do m=m1-rady,m2+rady f(:,m,n,ikapparho)=kapparho_floor+kapparho_cst enddo enddo ! ! To have a constant K, kappa=kappa0*T^3/rho, ! where kappa0=16./3.*sigmaSB/K ! case ('kappa_Kconst') kappa0=16./3.*sigmaSB/kappa_Kconst do n=n1-radz,n2+radz do m=m1-rady,m2+rady call eoscalc(f,mx,lnTT=lnTT) TT=exp(lnTT) f(:,m,n,ikapparho)=kapparho_floor+kappa0*TT**3 enddo enddo ! ! Rescaled (generalized) Kramers opacity ! case ('kappa_power_law') do n=n1-radz,n2+radz; do m=m1-rady,m2+rady call eoscalc(f,mx,lnrho=lnrho,lnTT=lnTT) rho=exp(lnrho) TT=exp(lnTT) if (knee_temp_opa==0.0) then f(:,m,n,ikapparho)=kapparho_floor+rho*kappa_cst(inu)* & (rho/ref_rho_opa)**expo_rho_opa* & (TT/ref_temp_opa)**expo_temp_opa else ! ! Use erf profile to connect smoothly to constant opacity beyond ``knee'' ! temperature. ! profile=1.0-cubic_step(TT,knee_temp_opa,width_temp_opa) f(:,m,n,ikapparho)=kapparho_floor+profile*rho*kappa_cst(inu)* & (rho/ref_rho_opa)**expo_rho_opa* & (TT/ref_temp_opa)**expo_temp_opa + & (1.0-profile)*rho*kappa_cst(inu)* & (knee_temp_opa/ref_temp_opa)**expo_temp_opa* & (TT/knee_temp_opa)**expo_temp_opa_buff endif enddo; enddo ! ! Rescaled (generalized) Kramers opacity ! kappa^{-1} = kappa1^{-1} + (kappa2+kappa20)^{-1} ! Here, kappa1 corresponds to exponents for Hminus ! and kappa2 corresponds to exponents for Kramers ! case ('kappa_double_power_law') do n=n1-radz,n2+radz; do m=m1-rady,m2+rady call eoscalc(f,mx,lnrho=lnrho,lnTT=lnTT) rho=exp(lnrho) TT=exp(lnTT) kappa1=kappa_cst(inu)* & (rho/ref_rho_opa)**expo_rho_opa* & (TT/ref_temp_opa)**expo_temp_opa kappa2=kappa20_cst(inu)+kappa_cst(inu)* & (rho/ref_rho_opa)**expo2_rho_opa* & (TT/ref_temp_opa)**expo2_temp_opa f(:,m,n,ikapparho)=kapparho_floor+rho/(1./kappa1+1./kappa2) & *(1.+ampl_bump*exp(-.5*((TT-TT_bump)/sigma_bump)**2)) enddo; enddo ! !AB: the following looks suspicious with *_cgs case ('Tsquare') !! Morfill et al. 1985 kappa0_cgs=2e-4 ! 2e-6 in D'Angelo 2003 (wrong!) 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)=kapparho_floor+exp(lnrho)*kappa0*((exp(lnTT))**2) enddo enddo ! case ('Kramers') !! as in Frank et al. 1992 (D'Angelo 2003) kappa0_cgs=6.6e22 !! (7.5e22 in Prialnik) 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)=kapparho_floor+kappa0*(exp(lnrho)**2)*(exp(lnTT))**(-3.5) enddo enddo ! case ('dust-infrared') 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))<=150) then tmp(i)=2e-4*exp(lnTT(i))**2 elseif (exp(lnTT(i))>=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 f(:,m,n,ikapparho)=kapparho_floor+exp(lnrho)*tmp enddo enddo ! case ('blob') if (lfirst) then f(:,:,:,ikapparho)=kapparho_floor+kapparho_const + amplkapparho & *spread(spread(exp(-(x/radius_kapparho)**2),2,my),3,mz) & *spread(spread(exp(-(y/radius_kapparho)**2),1,mx),3,mz) & *spread(spread(exp(-(z/radius_kapparho)**2),1,mx),2,my) lfirst=.false. endif ! case ('cos') if (lfirst) then f(:,:,:,ikapparho)=kapparho_floor+kapparho_const + amplkapparho & *spread(spread(cos(kx_kapparho*x),2,my),3,mz) & *spread(spread(cos(ky_kapparho*y),1,mx),3,mz) & *spread(spread(cos(kz_kapparho*z),1,mx),2,my) lfirst=.false. endif ! case ('rad_ionization') !! as in Miao et al. 2006 ! TODO: The code still needs to be able to calculate the right ! ionization fraction for this. This does not work right without it, ! but does no harm either to anyone else. - mvaisala do n= n1-radz, n2+radz do m= m1-rady, m2+rady call eoscalc(f,mx,lnrho=lnrho,yH=yH) f(:,m,n,ikapparho)=kapparho_floor+ sigmaH_*(1 - yH)*exp(2*lnrho+log(m_H)) enddo enddo ! ! Opacity proportional to magnetic energy density ! (used for visualization purposes) ! Needs to be done at every time step (not just when lfirst=.true.) ! case ('B2') !! magnetic field call calc_kapparho_B2(f) ! case ('B2+W2') !! magnetic field and vorticity call calc_kapparho_B2_W2(f) ! ! Read from file ! case ('read_file') open (lun_input, file=trim(directory_prestart)//'/kapparho.dat', form='unformatted') read (lun_input) kapparho_noghost close (lun_input) f(l1:l2,m1:m2,n1:n2,ikapparho)=kapparho_noghost ! case ('nothing') do n=n1-radz,n2+radz do m=m1-rady,m2+rady f(l1:l2,m,n,ikapparho)=0.0 enddo enddo ! case default call fatal_error('opacity','no such opacity type: '//trim(opacity_type)) ! endselect ! endsubroutine opacity !*********************************************************************** subroutine calc_Srad_B2(f,Srad) ! ! Calculate B2 for special prescriptions of source function. ! ! 21-feb-2009/axel: coded ! use Sub, only: gij, curl_mn, dot2_mn ! real, dimension(mx,my,mz,mfarray), intent(in) :: f real, dimension(mx,my,mz), intent(out) :: Srad real, dimension(nx,3,3) :: aij real, dimension(nx,3) :: aa,bb real, dimension(nx) :: b2 ! ! Check whether B-field is available, and then calculate (curlA)^2. ! if (iaa==0) then call fatal_error('calc_Srad_B2','no magnetic field available') else do n=n1,n2 do m=m1,m2 aa=f(l1:l2,m,n,iax:iaz) call gij(f,iaa,aij,1) call curl_mn(aij,bb,aa) call dot2_mn(bb,b2) Srad(l1:l2,m,n)=b2 enddo enddo endif ! endsubroutine calc_Srad_B2 !*********************************************************************** subroutine calc_Srad_W2(f,Srad) ! ! Calculate W2 for special prescriptions of source function. ! ! 21-feb-2009/axel: coded ! use Sub, only: gij, curl_mn, dot2_mn ! real, dimension(mx,my,mz,mfarray), intent(in) :: f real, dimension(mx,my,mz), intent(out) :: Srad real, dimension(nx,3,3) :: uij real, dimension(nx,3) :: uu,oo real, dimension(nx) :: o2 ! ! Check whether u-field is available, and then calculate (curl u)^2. ! if (iuu==0) then call fatal_error('calc_Srad_W2','no velocity field available') else do n=n1,n2 do m=m1,m2 uu=f(l1:l2,m,n,iux:iuz) call gij(f,iuu,uij,1) call curl_mn(uij,oo,uu) call dot2_mn(oo,o2) Srad(l1:l2,m,n)=o2 enddo enddo endif ! endsubroutine calc_Srad_W2 !*********************************************************************** subroutine calc_kapparho_B2(f) ! ! Calculate B2 for special prescriptions of opacity ! ! 21-feb-2009/axel: coded ! use Sub, only: gij, curl_mn, dot2_mn ! real, dimension(mx,my,mz,mfarray), intent(inout) :: f real, dimension(nx,3,3) :: aij real, dimension(nx,3) :: aa,bb real, dimension(nx) :: b2 integer :: ighost ! ! Check whether B-field is available, and then calculate (curlA)^2 ! if (iaa==0) then call fatal_error('calc_kapparho_B2','no magnetic field available') else do n=n1,n2 do m=m1,m2 aa=f(l1:l2,m,n,iax:iaz) call gij(f,iaa,aij,1) call curl_mn(aij,bb,aa) call dot2_mn(bb,b2) f(l1:l2,m,n,ikapparho)=kapparho_floor+b2 ! ! in the ghost zones, put the magnetic field equal to the value ! in the layer layer inside the domain. ! do ighost=1,nghost f(l1-ighost,:,:,ikapparho)=f(l1,:,:,ikapparho) f(l2+ighost,:,:,ikapparho)=f(l2,:,:,ikapparho) f(:,m1-ighost,:,ikapparho)=f(:,m1,:,ikapparho) f(:,m2+ighost,:,ikapparho)=f(:,m2,:,ikapparho) f(:,:,n1-ighost,ikapparho)=f(:,:,n1,ikapparho) f(:,:,n2+ighost,ikapparho)=f(:,:,n2,ikapparho) enddo enddo enddo endif ! endsubroutine calc_kapparho_B2 !*********************************************************************** subroutine calc_kapparho_B2_W2(f) ! ! Calculate B2+W2 for special prescriptions of opacity. ! Assume that both are opaque in all colors. ! ! 21-feb-2009/axel: coded ! use Sub, only: gij, curl_mn, dot2_mn ! real, dimension(mx,my,mz,mfarray), intent(inout) :: f real, dimension(nx,3,3) :: aij,uij real, dimension(nx,3) :: aa,bb,uu,oo real, dimension(nx) :: b2,o2 integer :: ighost ! ! Check whether B-field is available, and then calculate (curlA)^2. ! if (iaa==0) then call fatal_error('calc_kapparho_B2_W2','no magnetic field available') else do n=n1,n2 do m=m1,m2 aa=f(l1:l2,m,n,iax:iaz) uu=f(l1:l2,m,n,iux:iuz) call gij(f,iaa,aij,1) call gij(f,iuu,uij,1) call curl_mn(aij,bb,aa) call curl_mn(uij,oo,uu) call dot2_mn(bb,b2) call dot2_mn(oo,o2) f(l1:l2,m,n,ikapparho)=kapparho_floor+b2+o2 ! ! In the ghost zones, put the magnetic field equal to the value ! in the layer layer inside the domain. ! do ighost=1,nghost f(l1-ighost,:,:,ikapparho)=f(l1,:,:,ikapparho) f(l2+ighost,:,:,ikapparho)=f(l2,:,:,ikapparho) f(:,m1-ighost,:,ikapparho)=f(:,m1,:,ikapparho) f(:,m2+ighost,:,ikapparho)=f(:,m2,:,ikapparho) f(:,:,n1-ighost,ikapparho)=f(:,:,n1,ikapparho) f(:,:,n2+ighost,ikapparho)=f(:,:,n2,ikapparho) enddo enddo enddo endif ! endsubroutine calc_kapparho_B2_W2 !*********************************************************************** subroutine init_rad(f) ! ! Dummy routine for Flux Limited Diffusion routine ! initialise radiation; called from start.f90 ! ! 15-jul-2002/nils: dummy routine ! real, dimension (mx,my,mz,mfarray) :: f ! call keep_compiler_quiet(f) ! endsubroutine init_rad !*********************************************************************** subroutine pencil_criteria_radiation ! ! All pencils that the Radiation module depends on are specified here. ! ! 21-11-04/anders: coded ! if (lcooling) then if (ldt) then lpenc_requested(i_rho1)=.true. lpenc_requested(i_cp1)=.true. lpenc_requested(i_cv1)=.true. lpenc_requested(i_cp)=.true. lpenc_requested(i_cv)=.true. lpenc_requested(i_TT)=.true. endif if (lrad_cool_diffus.or.lrad_pres_diffus) then lpenc_requested(i_glnrho)=.true. lpenc_requested(i_TT)=.true. lpenc_requested(i_cp1)=.true. lpenc_requested(i_rho1)=.true. lpenc_requested(i_glnTT)=.true. lpenc_requested(i_del2lnTT)=.true. lpenc_requested(i_cv1)=.true. endif lpenc_requested(i_TT1)=.true. lpenc_requested(i_rho1)=.true. if (ltemperature) then lpenc_requested(i_cv1)=.true. endif endif ! endsubroutine pencil_criteria_radiation !*********************************************************************** subroutine pencil_interdep_radiation(lpencil_in) ! ! Interdependency among pencils provided by the Radiation module ! is specified here. ! ! 21-11-04/anders: coded ! logical, dimension (npencils) :: lpencil_in ! call keep_compiler_quiet(lpencil_in) ! endsubroutine pencil_interdep_radiation !*********************************************************************** subroutine calc_pencils_radiation(f,p) ! ! Calculate Radiation pencils. ! Most basic pencils should come first, as others may depend on them. ! ! 21-11-04/anders: coded ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! intent(in) :: f,p ! call keep_compiler_quiet(f) call keep_compiler_quiet(p) ! endsubroutine calc_pencils_radiation !*********************************************************************** subroutine de_dt(f,df,p) ! ! Dummy routine for Flux Limited Diffusion routine ! ! 15-jul-2002/nils: dummy routine ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! call keep_compiler_quiet(f,df) call keep_compiler_quiet(p) ! endsubroutine de_dt !*********************************************************************** subroutine read_radiation_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=radiation_init_pars, IOSTAT=iostat) ! endsubroutine read_radiation_init_pars !*********************************************************************** subroutine write_radiation_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=radiation_init_pars) ! endsubroutine write_radiation_init_pars !*********************************************************************** subroutine read_radiation_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=radiation_run_pars, IOSTAT=iostat) ! endsubroutine read_radiation_run_pars !*********************************************************************** subroutine write_radiation_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=radiation_run_pars) ! endsubroutine write_radiation_run_pars !*********************************************************************** subroutine rprint_radiation(lreset,lwrite) ! ! Dummy routine for Flux Limited Diffusion routine ! reads and registers print parameters relevant for radiative part ! ! 16-jul-02/nils: adapted from rprint_hydro ! use Diagnostics, only: parse_name use FArrayManager, only: farray_index_append ! integer :: iname,inamez logical :: lreset,lwr logical, optional :: lwrite ! lwr = .false. if (present(lwrite)) lwr=lwrite ! ! reset everything in case of RELOAD ! (this needs to be consistent with what is defined above!) ! if (lreset) then idiag_Qradrms=0; idiag_Qradmax=0; idiag_Fradzm=0; idiag_kapparhom=0; idiag_Sradm=0 idiag_Fradzmz=0; idiag_kapparhomz=0; idiag_taumz=0 idiag_dtchi=0; idiag_dtrad=0 ivid_Jrad=0; ivid_Isurf=0 endif ! ! check for those quantities that we want to evaluate online ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),'Qradrms',idiag_Qradrms) call parse_name(iname,cname(iname),cform(iname),'Qradmax',idiag_Qradmax) call parse_name(iname,cname(iname),cform(iname),'Fradzm',idiag_Fradzm) call parse_name(iname,cname(iname),cform(iname),'kapparhom',idiag_kapparhom) call parse_name(iname,cname(iname),cform(iname),'Sradm',idiag_Sradm) call parse_name(iname,cname(iname),cform(iname),'dtchi',idiag_dtchi) call parse_name(iname,cname(iname),cform(iname),'dtrad',idiag_dtrad) enddo ! ! check for those quantities for which we want xy-averages ! do inamez=1,nnamez call parse_name(inamez,cnamez(inamez),cformz(inamez),'Fradzmz',idiag_Fradzmz) call parse_name(inamez,cnamez(inamez),cformz(inamez),'taumz',idiag_taumz) call parse_name(inamez,cnamez(inamez),cformz(inamez),'kapparhomz',idiag_kapparhomz) enddo ! ! check for those quantities for which we want video slices ! if (lwrite_slices) then where(cnamev=='Qrad'.or.cnamev=='Frad'.or. & cnamev=='Srad'.or.cnamev=='kapparho') cformv='DEFINED' endif do iname=1,nnamev call parse_name(iname,cnamev(iname),cformv(iname),'Isurf',ivid_Isurf) call parse_name(iname,cnamev(iname),cformv(iname),'Jrad',ivid_Jrad) enddo ! ! write column where which radiative variable is stored ! if (lwr) then call farray_index_append('iQrad',iQrad) call farray_index_append('ikapparho',ikapparho) call farray_index_append('iKR_Frad',iKR_Frad) call farray_index_append('iKR_Fradx',iKR_Fradx) call farray_index_append('iKR_Frady',iKR_Frady) call farray_index_append('iKR_Fradz',iKR_Fradz) endif ! endsubroutine rprint_radiation !*********************************************************************** subroutine get_slices_radiation(f,slices) ! ! Write slices for animation of radiation variables. ! ! 26-jul-06/tony: coded ! use Slices_methods, only: assign_slices_scal, assign_slices_vec, & assign_slices_f_scal, nullify_slice_pointers ! real, dimension (mx,my,mz,mfarray) :: f type (slice_data) :: slices ! integer, save :: idir_last=0 ! ! Loop over slices ! select case (trim(slices%name)) ! ! Surface intensity. If nIsurf==1, then we only want the vertical ray. ! Otherwise, for bigger nIsurf, we get up to nIsurf upward pointing rays, ! but possibly not the vertical one, if nIsurf is not big enough. ! Count which one is which by putting ip<14 and look at output. ! case ('Isurf') call nullify_slice_pointers(slices) if (lwrite_slice_xy2) then if (slices%index>=nIsurf) then slices%ready=.false. else if (slices%index==0) idir_last=0 do idir=idir_last+1,ndir nrad=dir(idir,3) if ((nIsurf>1.and.nrad>0).or. & (nIsurf==1.and.lrad==0.and.mrad==0.and.nrad==1)) 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 endif ! ! Heating rate ! case ('Qrad'); call assign_slices_scal(slices,f,iQrad) ! ! Mean intensity ! case ('Jrad') !J = S + Q/(4pi) !slices%yz=.25*pi_1*f(ix_loc,m1:m2,n1:n2,iQrad)+Srad(ix_loc,m1:m2,n1:n2) !slices%xz=.25*pi_1*f(l1:l2,iy_loc,n1:n2,iQrad)+Srad(l1:l2,iy_loc,n1:n2) !slices%xy=.25*pi_1*f(l1:l2,m1:m2,iz_loc,iQrad)+Srad(l1:l2,m1:m2,iz_loc) !slices%xy2=.25*pi_1*f(l1:l2,m1:m2,iz2_loc,iQrad)+Srad(l1:l2,m1:m2,iz2_loc) ! call assign_slices_vec(slices,Jrad_xy,Jrad_xz,Jrad_yz,Jrad_xy2,Jrad_xy3,Jrad_xy4,Jrad_xz2,ncomp=nnu) ! ! Source function ! case ('Srad'); call assign_slices_f_scal(slices,Srad,1) ! ! Opacity ! case ('kapparho'); call assign_slices_scal(slices,f,ikapparho) ! ! Radiative Flux ! case ('Frad'); call assign_slices_vec(slices,f,iKR_Fradx) ! 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 ! 6-sep-16/MR: replaced dxmax by dxmax_pencil as suggested ! use Diagnostics use EquationOfState use Sub use General, only: notanumber ! real, dimension (mx,my,mz,mvar+maux) :: f type (pencil_case) :: p real, dimension (nx) :: Krad,chi_rad,g2, diffus_chi,advec_crad2 real, dimension (nx) :: local_optical_depth,opt_thin,opt_thick,dt1_rad real :: fact integer :: j,k ! intent(inout) :: f intent(in) :: p ! ! Calculate diffusion coefficient, Krad=16*sigmaSB*T^3/(3*kappa*rho). ! fact=16.*sigmaSB/3. Krad=fact*p%TT**3/f(l1:l2,m,n,ikapparho) ! ! 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=iKR_Frad+(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 advec2=advec2+advec_crad2 if (notanumber(advec_crad2)) print*, 'advec_crad2=',advec_crad2 ! ! Check maximum diffusion from thermal diffusion. ! With heat conduction, the second-order term for leading entropy term ! is gamma*chi_rad*del2ss. ! ! 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) ! diffus_chi=0. 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 ! local_optical_depth=dxmax_pencil*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) maxdiffus=max(maxdiffus,diffus_chi) ! !-- 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