1c1 < ! $Id: radiation_ray.f90,v 1.86 2006-05-17 17:13:16 theine Exp $ --- > ! $Id: radiation_ray.f90,v 1.90 2006-05-22 23:22:57 brandenb Exp $ 20a21,25 > ! > ! This module currently does not work for fully periodic domains. > ! The z-direction has to be non-periodic > ! > ! TODO: Calculate weights properly 29a35,53 > type Qbound !Qbc > real :: val > logical :: set > endtype Qbound > > type Qpoint !Qpt > real, pointer :: val > logical, pointer :: set > endtype Qpoint > > real, dimension (mx,my,mz) :: Srad,kapparho,tau,Qrad,Qrad0 > 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 > 34,36d57 < real, dimension (mx,my,mz) :: Srad,kapparho,tau,Qrad,Qrad0 < real, dimension (mx,my,mz,3) :: Frad < real, dimension (maxdir,3) :: unit_vec 37a59 > real, dimension (maxdir,3) :: unit_vec 47c69,70 < integer :: ipystart,ipystop,ipzstart,ipzstop --- > integer :: ipzstart,ipzstop,ipystart,ipystop > logical :: lperiodic_ray,lperiodic_ray_x,lperiodic_ray_y,lperiodic_ray_z 62d84 < logical :: lradpressure=.false.,lradflux=.false. 63a86 > logical :: lradpressure=.false.,lradflux=.false. 80c103 < lintrinsic,lcommunicate,lrevision --- > lintrinsic,lcommunicate,lrevision,lradflux 88c111 < lintrinsic,lcommunicate,lrevision,lcooling,lradpressure,lradflux --- > lintrinsic,lcommunicate,lrevision,lcooling,lradflux,lradpressure 139c162 < "$Id: radiation_ray.f90,v 1.86 2006-05-17 17:13:16 theine Exp $") --- > "$Id: radiation_ray.f90,v 1.90 2006-05-22 23:22:57 brandenb Exp $") 175a199 > logical :: periodic_xy_plane,bad_ray 182a207,212 > ! 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) > ! 190c220,223 < --- > ! > ! The value of rad2 determines whether a ray is along a coordinate axis (1), > ! a face diagonal (2), or a room diagonal (3) > ! 191a225,234 > ! > ! 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) 193c236 < if (rad2>0.and.rad2<=rad2max) then --- > if ((rad2>0.and.rad2<=rad2max).and.(.not.bad_ray)) then 244c287 < if (ndir>0) weight=1.0/ndir --- > if (ndir>0) weight=4*pi/ndir 248,253d290 < ! 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) < ! 296c333,340 < if (lcommunicate) call Qcommunicate --- > if (lcommunicate) then > if (lperiodic_ray) then > call Qperiodic > else > call Qpointers > call Qcommunicate > endif > endif 355a400,405 > ! 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) > ! 450,459c500 < subroutine Qcommunicate < ! < ! set boundary intensities or receive from neighboring processors < ! < ! 11-jul-03/tobi: coded < ! < use Cparam, only: nprocy,nprocz < use Mpicomm, only: ipy,ipz < use Mpicomm, only: radboundary_zx_recv,radboundary_zx_send < use Mpicomm, only: radboundary_xy_recv,radboundary_xy_send --- > subroutine Qpointers 461,465c502,516 < real, dimension(my,mz) :: Qrad0_yz < real, dimension(mx,mz) :: Qrad0_zx < real, dimension(mx,my) :: Qrad0_xy < integer :: raysteps < integer :: l,m,n --- > ! For each gridpoint at the downstream boundaries, set up a > ! pointer (Qpt_{yz,zx,xy}) that points to a unique location > ! at the upstream boundaries (Qbc_{yz,zx,xy}). Both > ! Qpt_{yz,zx,xy} and Qbc_{yz,zx,xy} are derived types > ! containing at each grid point the value of the heating rate > ! (...%val) and whether the heating rate at that point has > ! been already set or not (...%set). > ! > ! 30-jul-05/tobi: coded > ! > integer :: steps > integer :: minsteps > integer :: lsteps,msteps,nsteps > real, pointer :: val > logical, pointer :: set 467,486c518 < ! set boundary values < ! < if (lrad/=0) then < call radboundary_yz_set(Qrad0_yz) < Qrad0(llstart-lrad,mm1:mm2,nn1:nn2)=Qrad0_yz(mm1:mm2,nn1:nn2) < endif < ! < if (mrad/=0) then < if (ipy==ipystart) call radboundary_zx_set(Qrad0_zx) < if (ipy/=ipystart) call radboundary_zx_recv(mrad,idir,Qrad0_zx) < Qrad0(ll1:ll2,mmstart-mrad,nn1:nn2)=Qrad0_zx(ll1:ll2,nn1:nn2) < endif < ! < if (nrad/=0) then < if (ipz==ipzstart) call radboundary_xy_set(Qrad0_xy) < if (ipz/=ipzstart) call radboundary_xy_recv(nrad,idir,Qrad0_xy) < Qrad0(ll1:ll2,mm1:mm2,nnstart-nrad)=Qrad0_xy(ll1:ll2,mm1:mm2) < endif < ! < ! propagate boundary values --- > ! yz-plane 488a521,527 > > l=llstop > lsteps=(l+lrad-llstart)/lrad > > msteps=huge(msteps) > nsteps=huge(nsteps) > 491,495c530,538 < raysteps=(llstop-llstart)/lrad < if (mrad/=0) raysteps=min(raysteps,(mmstop-m)/mrad) < if (nrad/=0) raysteps=min(raysteps,(nnstop-n)/nrad) < Qrad0(llstart+lrad*raysteps,m+mrad*raysteps,n+nrad*raysteps) & < =Qrad0(llstart-lrad, m-mrad, n-nrad ) --- > > if (mrad/=0) msteps = (m+mrad-mmstart)/mrad > if (nrad/=0) nsteps = (n+nrad-nnstart)/nrad > > call assign_pointer(lsteps,msteps,nsteps,val,set) > > Qpt_yz(m,n)%set => set > Qpt_yz(m,n)%val => val > 497a541 > 499a544,545 > ! zx-plane > ! 500a547,553 > > m=mmstop > msteps=(m+mrad-mmstart)/mrad > > nsteps=huge(nsteps) > lsteps=huge(lsteps) > 503,507c556,564 < raysteps=(mmstop-mmstart)/mrad < if (nrad/=0) raysteps=min(raysteps,(nnstop-n)/nrad) < if (lrad/=0) raysteps=min(raysteps,(llstop-l)/lrad) < Qrad0(l+lrad*raysteps,mmstart+mrad*raysteps,n+nrad*raysteps) & < =Qrad0(l-lrad, mmstart-mrad, n-nrad ) --- > > 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 > 509a567 > 511a570,571 > ! xy-plane > ! 512a573,579 > > n=nnstop > nsteps=(n+nrad-nnstart)/nrad > > lsteps=huge(lsteps) > msteps=huge(msteps) > 515,519c582,590 < raysteps=(nnstop-nnstart)/nrad < if (lrad/=0) raysteps=min(raysteps,(llstop-l)/lrad) < if (mrad/=0) raysteps=min(raysteps,(mmstop-m)/mrad) < Qrad0(l+lrad*raysteps,m+mrad*raysteps,nnstart+nrad*raysteps) & < =Qrad0(l-lrad, m-mrad, nnstart-nrad ) --- > > 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 > 521a593 > 522a595,624 > > 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 524c626 < ! send boundary values --- > ! Determine the boundary heating rates at all upstream boundaries. 526,530c628,676 < if (mrad/=0.and.ipy/=ipystop) then < Qrad0_zx(ll1:ll2,nn1:nn2)=Qrad0(ll1:ll2,mmstop,nn1:nn2) & < *exp(-tau(ll1:ll2,mmstop,nn1:nn2)) & < +Qrad(ll1:ll2,mmstop,nn1:nn2) < call radboundary_zx_send(mrad,idir,Qrad0_zx) --- > ! 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. > 533,537c679,794 < if (nrad/=0.and.ipz/=ipzstop) then < Qrad0_xy(ll1:ll2,mm1:mm2)=Qrad0(ll1:ll2,mm1:mm2,nnstop) & < *exp(-tau(ll1:ll2,mm1:mm2,nnstop)) & < +Qrad(ll1:ll2,mm1:mm2,nnstop) < call radboundary_xy_send(nrad,idir,Qrad0_xy) --- > ! 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 538a796,803 > > 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. 539a805,832 > 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 > > 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 > 541a835,939 > 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_yz emtau1_yz=tau_yz*(1-0.5*tau_yz*(1-tau_yz/3)) > elsewhere > emtau1_yz=1-exp(-tau_yz) > endwhere > ! > ! The requirement of periodicity gives the following heating rate at the > ! upstream boundary. > ! > Qrad0(llstart-lrad,m1:m2,n1:n2)=Qrad_yz/emtau1_yz > > endif > ! > ! y-direction > ! > if (mrad/=0) then > ! > ! Intrinsic heating rate and optical depth at the downstream boundary of > ! each processor. > ! > Qrad_zx=Qrad(l1:l2,mmstop,n1:n2) > tau_zx=tau(l1:l2,mmstop,n1:n2) > ! > ! Gather intrinsic heating rates and optical depths from all processors > ! into one rank-3 array available on each processor. > ! > call radboundary_zx_periodic_ray(Qrad_zx,tau_zx,Qrad_zx_all,tau_zx_all) > ! > ! Find out in which direction we want to loop over processors. > ! > if (mrad>0) 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_zx emtau1_tot_zx=tau_tot_zx*(1-0.5*tau_tot_zx*(1-tau_tot_zx/3)) > elsewhere > emtau1_tot_zx=1-exp(-tau_tot_zx) > endwhere > ! > ! The requirement of periodicity gives the following heating rate at the > ! upstream boundary of this processor. > ! > Qrad0(l1:l2,mmstart-mrad,n1:n2)=Qrad_tot_zx/emtau1_tot_zx > > endif > > endsubroutine Qperiodic > !*********************************************************************** 558,561d955 < ! avoid underflows < ! < !tau=min(tau,-log((1+10*epsilon(tau))*tiny(tau))) < ! 568,570c962 < if (tau(l,m,n) < -log(tiny(dtau_thresh_max))) then < Qrad(l,m,n)=Qrad(l,m,n)+Qrad0(l,m,n)*exp(-tau(l,m,n)) < endif --- > Qrad(l,m,n)=Qrad(l,m,n)+Qrad0(l,m,n)*exp(-tau(l,m,n)) 589,591c981 < ! sets the physical boundary condition on yz plane < ! < ! 6-jul-03/axel: coded --- > ! Sets the physical boundary condition on yz plane 593c983 < real, dimension(my,mz) :: Qrad0_yz --- > ! 6-jul-03/axel+tobi: coded 595c985 < ! no incoming intensity --- > use Mpicomm, only: stop_it 597,599c987,988 < if (bc_ray_x=='0') then < Qrad0_yz=-Srad(llstart-lrad,:,:) < endif --- > real, dimension(my,mz), intent(out) :: Qrad0_yz > real :: Irad_yz 601c990 < ! periodic boundary consition --- > ! No incoming intensity 603,604c992,993 < if (bc_ray_x=='p') then < Qrad0_yz=Qrad(llstop-lrad,:,:) --- > if (bc_ray_z=='0') then > Qrad0_yz=-Srad(llstart-lrad,:,:) 607c996 < ! set intensity equal to source function --- > ! Set intensity equal to source function 609c998 < if (bc_ray_x=='S') then --- > if (bc_ray_z=='S') then 617c1006 < ! sets the physical boundary condition on zx plane --- > ! Sets the physical boundary condition on zx plane 619c1008 < ! 6-jul-03/axel: coded --- > ! 6-jul-03/axel+tobi: coded 623c1012,1013 < real, dimension(mx,mz) :: Qrad0_zx --- > real, dimension(mx,mz), intent(out) :: Qrad0_zx > real :: Irad_zx 625c1015 < ! no incoming intensity --- > ! No incoming intensity 627c1017 < if (bc_ray_y=='0') then --- > if (bc_ray_z=='0') then 631,641c1021 < ! periodic boundary consition (currently not implemented for < ! multiple processors in the y-direction) < ! < if (bc_ray_y=='p') then < if (nprocy>1) then < call stop_it("radboundary_zx_set: periodic bc not implemented for nprocy>1") < endif < Qrad0_zx=Qrad(:,mmstop-mrad,:) < endif < ! < ! set intensity equal to source function --- > ! Set intensity equal to source function 643c1023 < if (bc_ray_y=='S') then --- > if (bc_ray_z=='S') then 674,684c1054 < ! periodic boundary consition (currently not implemented for < ! multiple processors in the z-direction) < ! < if (bc_ray_z=='p') then < if (nprocz>1) then < call stop_it("radboundary_xy_set: periodic bc not implemented for nprocz>1") < endif < Qrad0_xy=Qrad(:,:,nnstop-nrad) < endif < ! < ! set intensity equal to source function --- > ! Set intensity equal to source function 706c1076 < cooling=4*pi*kapparho(l1:l2,m,n)*f(l1:l2,m,n,iQrad) --- > cooling=kapparho(l1:l2,m,n)*f(l1:l2,m,n,iQrad) 741c1111 < ! --- > 768c1138 < ! --- > 832a1203 > use Cdata, only: kappa_es 850a1222,1229 > case ('kappa_es') > do n=n1-radz,n2+radz > do m=m1-rady,m2+rady > call eoscalc(f,mx,lnrho=lnrho) > kapparho(:,m,n)=kappa_es*exp(lnrho) > enddo > enddo > 1101d1479 < ------------------------------------------------------------------------ r5674 | axel.brandenburg | 2006-05-23 01:22:57 +0200 (Tue, 23 May 2006) | 3 lines added kappa_es ------------------------------------------------------------------------ r5650 | tobias.heinemann | 2006-05-19 19:03:30 +0200 (Fri, 19 May 2006) | 6 lines Merged radiation_ray and radiation_ray_periodic, the latter being obsolete now and soon to be removed I hope. Please do make harsh complaints if I don't document this module properly in the near future. ------------------------------------------------------------------------ r5649 | tobias.heinemann | 2006-05-19 15:05:23 +0200 (Fri, 19 May 2006) | 3 lines Added lradflux to radiation_init_pars ------------------------------------------------------------------------ r5648 | tobias.heinemann | 2006-05-19 13:56:17 +0200 (Fri, 19 May 2006) | 3 lines Changed weights such that sum(wights)=4*pi. Axel: Please check whether we calculate the radiative flux correctly now. ------------------------------------------------------------------------ r5638 | tobias.heinemann | 2006-05-17 19:13:16 +0200 (Wed, 17 May 2006) | 3 lines Added radiative flux to radiation_ray_periodic. Also renamed lradpressure_calc to lradflux.