! $Id: baroclinic_run.f90 19193 2012-06-30 12:55:46Z wdobler $ !** AUTOMATIC CPARAM.INC GENERATION **************************** ! Declare (for generation of special_dummies.inc) the number of f array ! variables and auxiliary variables added by this module ! ! CPARAM logical, parameter :: lspecial = .true. ! ! MVAR CONTRIBUTION 0 ! MAUX CONTRIBUTION 0 ! !*************************************************************** !------------------------------------------------------------------- ! ! HOW TO USE THIS FILE ! -------------------- ! ! The rest of this file may be used as a template for your own ! special module. Lines which are double commented are intended ! as examples of code. Simply fill out the prototypes for the ! features you want to use. ! ! Save the file with a meaningful name, eg. geo_kws.f90 and place ! it in the $PENCIL_HOME/src/special directory. This path has ! been created to allow users ot optionally check their contributions ! in to the Pencil-Code CVS repository. This may be useful if you ! are working on/using the additional physics with somebodyelse or ! may require some assistance from one of the main Pencil-Code team. ! ! To use your additional physics code edit the Makefile.local in ! the src directory under the run directory in which you wish to ! use your additional physics. Add a line with all the module ! selections to say something like: ! ! SPECIAL=special/nstar ! ! Where nstar it replaced by the filename of your new module ! upto and not including the .f90 ! !-------------------------------------------------------------------- ! ! Equation of state from Lyra & Kuchner (2013), due to photoelectric ! heating of the gas by the dust. A polytropic term is added for ! completeness and to provide support against gravity in stratified ! simulations. The pressure gradient itself is ! ! gradP = cs20/(gamma*rho0)*(p%rho*p%grhop(:,j) + p%rhop*p%grho(:,j)) ! ! 17-may-12/wlad: coded ! 20-dec-12/wlad: streamlined, and moved here to special ! module Special ! use Cparam use Cdata use General, only: keep_compiler_quiet use Messages ! implicit none ! include '../special.h' ! real :: mu=1.0, kappa=0.0, factor_localiso=0.0, factor_photoelectric=1.0 real, dimension(ndustspec) :: const_pr=0.0 logical :: ldust_pressureforce=.true. logical :: lradiation_PRdrag=.false. logical :: ldivrhop_by_vol=.false. ! namelist /special_init_pars/ mu, kappa, factor_localiso, factor_photoelectric, & ldust_pressureforce, lradiation_PRdrag, const_pr, ldivrhop_by_vol ! namelist /special_run_pars/ mu, kappa, factor_localiso, factor_photoelectric, & ldust_pressureforce, lradiation_PRdrag, const_pr, ldivrhop_by_vol ! ! integer, parameter :: nmode_max = 50 ! real, dimension(nmode_max) :: gauss_ampl, rcenter, phicenter ! type InternalPencils real, dimension(nx,3) :: fpres real, dimension(nx,3) :: fpres_photoelectric real, dimension(nx,3) :: fpres_polytropic real, dimension(nx,3) :: fpres_localisotropic endtype InternalPencils type (InternalPencils) :: q ! real :: const1,const2,const3 integer :: idiag_photom=0,idiag_photomax=0 integer :: idiag_photomin=0,idiag_polym=0 integer :: idiag_polymax=0,idiag_polymin=0 integer :: idiag_dtcp=0 ! logical, pointer :: lpressuregradient_gas ! contains !*********************************************************************** subroutine register_special() ! ! Configure pre-initialised (i.e. before parameter read) variables ! which should be know to be able to evaluate ! ! 14-jul-09/wlad: coded ! use SharedVariables integer :: ierr ! call get_shared_variable('lpressuregradient_gas',lpressuregradient_gas,ierr) if (ierr/=0) call fatal_error('register_special','lpressuregradient_gas') ! endsubroutine register_special !*********************************************************************** subroutine initialize_special(f) ! ! called by run.f90 after reading parameters, but before the time loop ! ! 14-jul-09/wlad: coded ! use Cdata use Mpicomm use EquationOfState, only: gamma1,cs20 use SharedVariables ! real, dimension (mx,my,mz,mfarray) :: f integer :: ierr ! ! Initializing the pressures. Seeting mu=1 makes the polytrope a linear barotrope. ! const1=kappa*mu const2=factor_photoelectric * cs20 * gamma1 const3=factor_localiso ! if (lrun) then call get_shared_variable('lpressuregradient_gas',lpressuregradient_gas,ierr) if (ierr/=0) call fatal_error('initialize_special','lpressuregradient_gas') if (lpressuregradient_gas) then call fatal_error('initialize_special','switch lpressuregradient_gas=F in hydro_run_pars') endif endif ! if (lentropy) & call fatal_error('initialize_special','This code should be used with noentropy.') ! call keep_compiler_quiet(f) ! endsubroutine initialize_special !*********************************************************************** subroutine pencil_criteria_special() ! use Mpicomm ! if (const1 /= 0.0) then lpenc_requested(i_rho)=.true. lpenc_requested(i_glnrho)=.true. endif ! if (const2 /= 0.0) then lpenc_requested(i_glnrho)=.true. if (lparticles) then lpenc_requested(i_rhop)=.true. lpenc_requested(i_grhop)=.true. elseif (ldustdensity) then lpenc_requested(i_rhodsum)=.true. lpenc_requested(i_glnrhodsum)=.true. else call fatal_error("pencil_criteria_special",& "the world is flat, and we never got here") endif endif ! if (const3 /= 0.0) then lpenc_requested(i_glnTT)=.true. lpenc_requested(i_glnrho)=.true. endif ! lpenc_requested(i_cs2)=.true. ! if (lradiation_PRdrag) then lpenc_requested(i_uud)=.true. endif ! endsubroutine pencil_criteria_special !*********************************************************************** subroutine calc_pencils_special(f,p) ! ! Calculate the pressure gradient. For the polytropic, P=S*rho**mu, ! so fpres = - 1./rho * grad(P) = - S*mu * rho**(mu-2) * grad(rho) ! ! For the photoelectric pressure, P=cs20/(gamma*/rho0) * rho*rhod ! !use Particles_sub, only: find_grid_volume real, dimension(mx,my,mz,mfarray) :: f real, dimension(nx) :: areas, rhop_tmp real, dimension(nx,3) :: grhop_tmp type (pencil_case) :: p integer :: j, l ! if (ldivrhop_by_vol) then areas = x(l1:l2)*y(m)/dy_1(m) grhop_tmp(:,1) = p%grhop(:,1)/areas grhop_tmp(:,2) = p%grhop(:,2)/areas grhop_tmp(:,3) = p%grhop(:,3)/areas rhop_tmp = p%rhop(:)/areas else grhop_tmp = p%grhop rhop_tmp = p%rhop endif ! do j=1,3 if (const1/=0.0) then q%fpres_polytropic(:,j) = -const1 * p%rho**(mu-1) * p%glnrho(:,j) else q%fpres_polytropic(:,j) = 0. endif ! if (const2/=0.0) then if (lparticles) then q%fpres_photoelectric(:,j) = -const2 * (grhop_tmp(:,j) + rhop_tmp*p%glnrho(:,j)) else q%fpres_photoelectric(:,j) = -const2 * p%rhodsum * (p%glnrhodsum(:,j) + p%glnrho(:,j)) endif else q%fpres_photoelectric(:,j) = 0. endif ! if (const3/=0.0) then q%fpres_localisotropic(:,j) = -const3*p%cs2*(p%glnrho(:,j)+p%glnTT(:,j)) endif enddo ! q%fpres = q%fpres_localisotropic + q%fpres_photoelectric + q%fpres_polytropic ! call keep_compiler_quiet(f) ! endsubroutine calc_pencils_special !*********************************************************************** subroutine special_before_boundary(f) ! ! This subroutine calculates the full potential due to the turbulence. ! ! 03-oct-12/wlad: coded ! real, dimension(mx,my,mz,mfarray), intent(inout) :: f ! call keep_compiler_quiet(f) ! endsubroutine special_before_boundary !*********************************************************************** subroutine read_special_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=special_init_pars, IOSTAT=iostat) ! endsubroutine read_special_init_pars !*********************************************************************** subroutine write_special_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=special_init_pars) ! endsubroutine write_special_init_pars !*********************************************************************** subroutine read_special_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=special_run_pars, IOSTAT=iostat) ! endsubroutine read_special_run_pars !*********************************************************************** subroutine write_special_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=special_run_pars) ! endsubroutine write_special_run_pars !*********************************************************************** subroutine rprint_special(lreset,lwrite) ! use Diagnostics ! ! reads and registers print parameters relevant to special ! ! 14-jul-09/wlad: coded ! use FArrayManager, only: farray_index_append ! integer :: iname logical :: lreset,lwr logical, optional :: lwrite ! lwr = .false. if (present(lwrite)) lwr=lwrite ! ! Write information to index.pro ! if (lreset) then idiag_photom=0;idiag_photomax=0;idiag_photomin=0 idiag_polym=0;idiag_polymax=0;idiag_polymin=0 idiag_dtcp=0 endif ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),'photom',idiag_photom) call parse_name(iname,cname(iname),cform(iname),'photomax',idiag_photomax) call parse_name(iname,cname(iname),cform(iname),'photomin',idiag_photomin) call parse_name(iname,cname(iname),cform(iname),'polym',idiag_polym) call parse_name(iname,cname(iname),cform(iname),'polymax',idiag_polymax) call parse_name(iname,cname(iname),cform(iname),'polymin',idiag_polymin) call parse_name(iname,cname(iname),cform(iname),'dtcp',idiag_dtcp) enddo ! if (lwr) then call farray_index_append('i_photom',idiag_photom) call farray_index_append('i_photomax',idiag_photomax) call farray_index_append('i_photomin',idiag_photomin) call farray_index_append('i_polym',idiag_polym) call farray_index_append('i_polymax',idiag_polymax) call farray_index_append('i_polymin',idiag_polymin) call farray_index_append('i_dtcp',idiag_dtcp) endif ! endsubroutine rprint_special !*********************************************************************** subroutine special_calc_hydro(f,df,p) ! ! calculate a additional 'special' term on the right hand side of the ! momentum equation. ! ! Some precalculated pencils of data are passed in for efficiency ! others may be calculated directly from the f array ! use Cdata use Diagnostics use EquationOfState, only: gamma1 ! real, dimension (mx,my,mz,mvar+maux), intent(in) :: f real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p integer :: j,ju,k ! if (lfirst.and.ldt) advec_cs2 = (const3*p%cs2 + const2*gamma1 + const1) * dxyz_2 if (headtt.or.ldebug) & print*, 'dss_dt: max(advec_cs2) =', maxval(advec_cs2) ! ! Modified momentum equation ! if (ldust_pressureforce) then do j=1,3 ju=j+iuu-1 df(l1:l2,m,n,ju) = df(l1:l2,m,n,ju) + q%fpres(:,j) enddo endif ! if (lradiation_PRdrag) then do k=1,ndustspec df(l1:l2,m,n,iudx(k)) = df(l1:l2,m,n,iudx(k)) - 2*const_pr(k)*p%uud(:,1,k) df(l1:l2,m,n,iudy(k)) = df(l1:l2,m,n,iudy(k)) - const_pr(k)*p%uud(:,2,k) enddo endif ! if (ldiagnos) then if (idiag_photom/=0) call sum_mn_name( q%fpres_photoelectric(:,1),idiag_photom) if (idiag_photomax/=0) call max_mn_name( q%fpres_photoelectric,idiag_photomax) if (idiag_photomin/=0) call max_mn_name(-q%fpres_photoelectric,idiag_photomin,lneg=.true.) if (idiag_polym/=0) call sum_mn_name( q%fpres_polytropic(:,1),idiag_polym) if (idiag_polymax/=0) call max_mn_name( q%fpres_polytropic,idiag_polymax) if (idiag_polymin/=0) call max_mn_name(-q%fpres_polytropic,idiag_polymin,lneg=.true.) if (idiag_dtcp/=0) call max_mn_name(sqrt(advec_cs2)/cdt,idiag_dtcp,l_dt=.true.) endif ! call keep_compiler_quiet(f) call keep_compiler_quiet(p) ! endsubroutine special_calc_hydro !*********************************************************************** ! !*********************************************************************** !*********************************************************************** !******************************************************************** ! !******************************************************************** !************ DO NOT DELETE THE FOLLOWING ************** !******************************************************************** !** This is an automatically generated include file that creates ** !** copies dummy routines from nospecial.f90 for any Special ** !** routines not implemented in this file ** !** ** include '../special_dummies.inc' !******************************************************************** endmodule Special