! *********************************************************************** ! ! Copyright (C) 2010 Bill Paxton ! ! this file is part of mesa. ! ! mesa is free software; you can redistribute it and/or modify ! it under the terms of the gnu general library public license as published ! by the free software foundation; either version 2 of the license, or ! (at your option) any later version. ! ! mesa is distributed in the hope that it will be useful, ! but without any warranty; without even the implied warranty of ! merchantability or fitness for a particular purpose. see the ! gnu library general public license for more details. ! ! you should have received a copy of the gnu library general public license ! along with this software; if not, write to the free software ! foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa ! ! *********************************************************************** module run_star_extras use star_lib use star_def use const_def use alert_lib use mlt_def implicit none ! these routines are called by the standard run_star check_model contains ! include 'standard_run_star_extras.inc' subroutine extras_controls(s, ierr) type (star_info), pointer :: s integer, intent(out) :: ierr ierr = 0 ! this is the place to set any procedure pointers you want to change ! e.g., other_wind, other_mixing, other_energy (see star_data.inc) s % other_mixing => pavel_mixing_micro end subroutine extras_controls integer function extras_startup(s, id, restart, ierr) type (star_info), pointer :: s integer, intent(in) :: id logical, intent(in) :: restart integer, intent(out) :: ierr ierr = 0 extras_startup = 0 if (.not. restart) then call alloc_extra_info(s) else ! it is a restart call unpack_extra_info(s) end if end function extras_startup ! returns either keep_going, retry, backup, or terminate. integer function extras_check_model(s, id, id_extra) type (star_info), pointer :: s integer, intent(in) :: id, id_extra extras_check_model = keep_going if (.false. .and. s% star_mass_h1 < 0.35d0) then ! stop when star hydrogen mass drops to specified level extras_check_model = terminate write(*, *) 'have reached desired hydrogen mass' return end if end function extras_check_model integer function how_many_extra_log_columns(s, id, id_extra) type (star_info), pointer :: s integer, intent(in) :: id, id_extra how_many_extra_log_columns = 0 end function how_many_extra_log_columns subroutine data_for_extra_log_columns(s, id, id_extra, n, names, vals, ierr) type (star_info), pointer :: s integer, intent(in) :: id, id_extra, n character (len=maxlen_log_column_name) :: names(n) real(dp) :: vals(n) integer, intent(out) :: ierr !note: do NOT add these names to log_columns.list ! the log_columns.list is only for the built-in log column options. ! it must not include the new column names you are adding here. ierr = 0 end subroutine data_for_extra_log_columns integer function how_many_extra_profile_columns(s, id, id_extra) type (star_info), pointer :: s integer, intent(in) :: id, id_extra how_many_extra_profile_columns = 0 end function how_many_extra_profile_columns subroutine data_for_extra_profile_columns(s, id, id_extra, n, nz, names, vals, ierr) type (star_info), pointer :: s integer, intent(in) :: id, id_extra, n, nz character (len=maxlen_profile_column_name) :: names(n) real(dp) :: vals(nz,n) integer, intent(out) :: ierr integer :: k ierr = 0 !note: do NOT add these names to profile_columns.list ! the profile_columns.list is only for the built-in profile column options. ! it must not include the new column names you are adding here. ! here is an example for adding a profile column !if (n /= 1) stop 'data_for_extra_profile_columns' !names(1) = 'beta' !do k = 1, nz ! vals(k,1) = s% Pgas(k)/s% P(k) !end do end subroutine data_for_extra_profile_columns ! returns either keep_going or terminate. ! note: cannot request retry or backup; extras_check_model can do that. integer function extras_finish_step(s, id, id_extra) type (star_info), pointer :: s integer, intent(in) :: id, id_extra integer :: ierr extras_finish_step = keep_going call store_extra_info(s) ! to save a profile, ! s% need_to_save_profiles_now = .true. ! to update the star log, ! s% need_to_update_logfile_now = .true. end function extras_finish_step subroutine extras_after_evolve(s, id, id_extra, ierr) type (star_info), pointer :: s integer, intent(in) :: id, id_extra integer, intent(out) :: ierr ierr = 0 end subroutine extras_after_evolve ! routines for saving and restoring extra data so can do restarts ! put these defs at the top and delete from the following routines !integer, parameter :: extra_info_alloc = 1 !integer, parameter :: extra_info_get = 2 !integer, parameter :: extra_info_put = 3 subroutine alloc_extra_info(s) integer, parameter :: extra_info_alloc = 1 type (star_info), pointer :: s call move_extra_info(s,extra_info_alloc) end subroutine alloc_extra_info subroutine unpack_extra_info(s) integer, parameter :: extra_info_get = 2 type (star_info), pointer :: s call move_extra_info(s,extra_info_get) end subroutine unpack_extra_info subroutine store_extra_info(s) integer, parameter :: extra_info_put = 3 type (star_info), pointer :: s call move_extra_info(s,extra_info_put) end subroutine store_extra_info subroutine move_extra_info(s,op) integer, parameter :: extra_info_alloc = 1 integer, parameter :: extra_info_get = 2 integer, parameter :: extra_info_put = 3 type (star_info), pointer :: s integer, intent(in) :: op integer :: i, j, num_ints, num_dbls, ierr i = 0 ! call move_int or move_flg num_ints = i i = 0 ! call move_dbl num_dbls = i if (op /= extra_info_alloc) return if (num_ints == 0 .and. num_dbls == 0) return ierr = 0 call star_alloc_extras(s% id, num_ints, num_dbls, ierr) if (ierr /= 0) then write(*,*) 'failed in star_alloc_extras' write(*,*) 'alloc_extras num_ints', num_ints write(*,*) 'alloc_extras num_dbls', num_dbls stop 1 end if contains subroutine move_dbl(dbl) real(dp) :: dbl i = i+1 select case (op) case (extra_info_get) dbl = s% extra_work(i) case (extra_info_put) s% extra_work(i) = dbl end select end subroutine move_dbl subroutine move_int(int) integer :: int i = i+1 select case (op) case (extra_info_get) int = s% extra_iwork(i) case (extra_info_put) s% extra_iwork(i) = int end select end subroutine move_int subroutine move_flg(flg) logical :: flg i = i+1 select case (op) case (extra_info_get) flg = (s% extra_iwork(i) /= 0) case (extra_info_put) if (flg) then s% extra_iwork(i) = 1 else s% extra_iwork(i) = 0 end if end select end subroutine move_flg end subroutine move_extra_info ! examples of other_mixing subroutines ! examples subroutine pavel_mixing_micro(id, ierr) use chem_def, only: ih1, ihe4, in14, ic12, io16 integer, intent(in) :: id integer, intent(out) :: ierr integer :: k, nz, h1, he4, n14, c12, o16, h2, he3, li7, be7, & b8, c13, n13, n15, o17, o18, f19, ne22, lnlam, indice double precision :: x, dq00, dqm1, dqsum, T, kap, rho, cp, & molecular_diffusivity, radiative_diffusivity, & micro_diffusivity, y, z, car, oxy, maxabu, mol_dif, amas, & zat, x2, y3, lit, ber, bor, car3, z3, z5, oxy7, oxy8, flu, & neo type (star_info), pointer :: s ierr = 0 call get_star_ptr(id, s, ierr) if (ierr /= 0) return nz = s% nz h1 = s% net_iso(ih1) he4 = s% net_iso(ihe4) n14 = s% net_iso(in14) c12 = s% net_iso(ic12) o16 = s% net_iso(io16) h2 = s% net_iso(ih2) he3 = s% net_iso(ihe3) lit7 = s% net_iso(ilit7) be7 = s% net_iso(ibe7) b8 = s% net_iso(ib8) c13 = s% net_iso(ic13) n13 = s% net_iso(in13) n15 = s% net_iso(in15) o17 = s% net_iso(io17) o18 = s% net_iso(io18) f19 = s% net_iso(if19) ne22 = s% net_iso(ine22) ! add microscopic diffusivity to D_mix ! clip_D_limit = 1 ! zero mixing diffusion coeffs that are smaller than this ! in inlist controls make clip_D_limt = 0.0001 do k = 2, nz ! start at 2 since no mixing across surface ! interpolate values at face x = s% xa(h1,k) y = s% xa(he4,k) z = s% xa(n14,k) car = s% xa(c12,k) oxy = s% xa(o16,k) x2 = s% xa(h2,k) y3 = s% xa(he3,k) lit = s% xa(li7,k) ber = s% xa(be7,k) bor = s% xa(b8,k) car3 = s% xa(c13,k) z3 = s% xa(n13,k) z5 = s% xa(n15,k) oxy7 = s% xa(o17,k) oxy8 = s% xa(o18,k) flu = s% xa(f19,k) neo = s% xa(ne22,k) amas=((x*1)+(y*4)+(z*14)+(car*12)+(oxy*16)+(x2*2)+(y3*3)+ & (lit*7)+(ber*7)+(bor*8)+(car3*13)+(z3*13)+(z5*15)+ & (oxy7*17)+(oxy8*18)+(flu*19)+(neo*22))/(x+y+z+car+oxy+ & zat+x2+y3+lit+ber+bor+car3+z3+z5+oxy7+oxy8+flu+neo) zat=((x*1)+(y*2)+(z*7)+(car*6)+(oxy*8)+(x2*1)+(y3*2)+ & (lit*3)+(ber*4)+(bor*5)+(car3*6)+(z3*7)+(z5*7)+(oxy7*8)+ & (oxy8*8)+(flu*9)+(neo*10))/(x+y+z+car+oxy+zat+x2+y3+lit+ & ber+bor+car3+z3+z5+oxy7+oxy8+flu+neo) maxabu=max(x,y,z,car,oxy, x2, y3, lit, ber, bor, car3, & z3, z5, oxy7, oxy8, flu,neo) else if (maxabu .eq. y ) then lnlam=40 else if (maxabu .eq. car ) then lnlam=38 else if (maxabu .eq. oxy ) then lnlam=35 else lnlam=15 endif dq00 = s% dq(k) dqm1 = s% dq(k-1) dqsum = dq00 + dqm1 kap = (dqm1*s% opacity(k) + dq00*s% opacity(k-1))/dqsum rho = (dqm1*s% rho(k) + dq00*s% rho(k-1))/dqsum cp = (dqm1*s% cp(k) + dq00*s% cp(k-1))/dqsum T = (dqm1*s% T(k) + dq00*s% T(k-1))/dqsum !!!radiative_diffusivity = 4*crad*T**4/(15*kap*clight*rho**2) !!!molecular_diffusivity = 1.84d-17*(1+7*x)*(T**2.5/rho) mol_dif=2.21d-15*(amas**0.5/((zat**4)*lnlam)*(T**2.5/rho) !!!micro_diffusivity = mol_dif s% D_mix(k) = s% D_mix(k) + mol_diff !!print*, 'PAVEL MIXING' !! For debugging end do end subroutine pavel_mixing_micro end module run_star_extras