Commit c8840bca authored by Lena Noack's avatar Lena Noack
Browse files

Update el cond

parent 2cac3544
......@@ -885,7 +885,8 @@ dr_r1=dr_r(nt);
%%path_write = '/home/noacklen/ownCloud/Literature_Projects/Eigene_Projekte/Project_Ortenzi_Dorn/OctaveResultsNew3/';
%%path_write = '/scratch-2/PublishedProjects/Project_CaroDorn/Results_New/';
%path_write = '/scratch-2/PublishedProjects/Project_CaroDorn/Results_Revision_New/';
path_write = '/scratch-2/PublishedProjects/Project_CaroDorn/Results_Revision_New_Fast/';
%path_write = '/scratch-2/PublishedProjects/Project_CaroDorn/Results_Revision_New_Fast/';
path_write = '/scratch-2/PublishedProjects/Project_CaroDorn/Results_Revision_New_Fast_SolNew/';
%path_write = 'C:\Arbeit\NextCloud_FUB\Literature_Projects\Eigene_Projekte\Project_Ortenzi_Dorn\_Revision\TestCases\';
%%path_write = sprintf('/home/noacklen/ownCloud/Literature_Projects/Eigene_Projekte/Project_Ortenzi_Dorn/OctaveResultsNew/%s/',redox);
%%path_write = sprintf('/home/noacklen/Arbeit/Ortenzi_Dorn_Paper/%s/',redox); % to avoid update conflicts with NextCloud
......@@ -1079,7 +1080,8 @@ fac=1;
while ((abs(H2O_dissolved-H2O_dissolved_old)>tol || abs(CO2_dissolved-CO2_dissolved_old)>tol) && (i<i_max))
H2O_dissolved_old = H2O_dissolved;
CO2_dissolved_old = CO2_dissolved;
[H2O_dissolved,CO2_dissolved] = solubility(p_bar,T_lid,H2O_fl_mole,X_H2O,hydrous);
[H2O_dissolved,CO2_dissolved] = solubility(p_bar,T_lid,H2O_fl_mole,H2O_dissolved,hydrous);
%[H2O_dissolved,CO2_dissolved] = solubility(p_bar,T_lid,H2O_fl_mole,X_H2O,hydrous);
dX_H2O = max(0,X_H2O-H2O_dissolved); % what would go into gas bubble in wt-%
dX_CO2 = max(0,X_CO2-CO2_dissolved); % what would go into gas bubble in wt-ppm
if (H2O_fl_mole==1), dX_CO2=0;, CO2_dissolved=CO2_dissolved_old;, end
......
......@@ -346,6 +346,7 @@ contains
Tm_surf = T_K*exp(-pX%alpha_melt*pF%alpha/(pX%density_melt*pF%rho*pX%Cp_melt*pF%Cp)*(P_bar-Psurf)*10.0**5)
if (pX%solubility.eq.1) then
write(*,*) "Pre-sol:",m,part%X_CO2(m),part%X_H2O(m)
call solubility_calc(mc%P_N2+mc%P_CO2+mc%P_H2O+mc%P_CO2+mc%P_H2O,& ! psurf total
&mc%P_CO2,mc%P_H2O,Tm_surf,part%X_H2O(m)/part%df(m),part%X_CO2(m)/part%df(m),&
&H2O_diss_melt,CO2_diss_melt,H2O_gas,CO2_gas)
......@@ -354,6 +355,7 @@ contains
part%X_H2O(m) = H2O_gas*part%df(m) ! only volatiles that are released into atmosphere; calculated back to cell volume fraction
part%X_CO2(m) = CO2_gas*part%df(m)
write(*,*) "Post-sol:",m,part%X_CO2(m),part%X_H2O(m)
endif
end if
......
......@@ -901,7 +901,9 @@ module initialize
else if (pH%use_magn_H.ge.2.0_dp) then
! 2: time-dependent version, values in Qr_file in erg g^-1 s^-1
! 3: time-dependent and angle-dependent version, dQ nicht in erg g^-1 s^-1 sondern in erg/(s cm^3)! % ToDo: max wert am äquator hier, noch winkelabhängig machen mit sin so dass 0 am Pol
! 3: time-dependent, use angle-dep version with max value, divide by 2/rho, should be then the same as version 2 (only for testing)
! 4: time-dependent and angle-dependent version, dQ nicht in erg g^-1 s^-1 sondern in erg/(s cm^3)!
! max value given at equator, vec_magn_H(k) is multiplied with angle in thermal.f90
n_t = 0 !prof_nm ! profiles calculated for my input files, therefore same number of entries for mantle!
if (pH%use_magn_H.eq.2.0_dp) then
......@@ -947,14 +949,22 @@ module initialize
enddo
close(70)
pH%vals_H(1,:) = pH%vals_H(1,:)/(100.0*pF%D) ! radius from cm -> m -> nondim
pH%vals_H(2:n_z+1,:) = pH%vals_H(2:n_z+1,:) * 0.1_dp * pF%D**2/(pF%k*pF%DeltaT) ! values given in erg/(s cm³) = 1e-7 W/cm³ = 1e-1 W/m³ / rho -> 1e-1 W/kg (divide here by rho, later multiply again with rho, hence no pF%rho)
if (pH%use_magn_H.lt.5.0_dp) then
! units not in SI
pH%vals_H(1,:) = pH%vals_H(1,:)/(100.0_dp*pF%D) ! radius from cm -> m -> nondim
pH%vals_H(2:n_z+1,:) = pH%vals_H(2:n_z+1,:) * 0.1_dp * pF%D**2/(pF%k*pF%DeltaT) ! values given in erg/(s cm³) = 1e-7 W/cm³ = 1e-1 W/m³ / rho -> 1e-1 W/kg (divide here by rho, later multiply again with rho, hence no pF%rho)
else
! was supposed to be in m, by mistake given in dm -> divide by 10
pH%vals_H(1,:) = pH%vals_H(1,:)/(10.0_dp*pF%D)
! already in W/kg
pH%vals_H(2:n_z+1,:) = pH%vals_H(2:n_z+1,:) * pF%D**2/(pF%k*pF%DeltaT)
endif
endif
Q = 0.0_dp
! interpolate values, here time is zero, hence second row (first heating row)
do k=1,nr
do kp=1,n_t-1
do k=1,nr ! nr is number of elements in radial direction in mesh
do kp=1,n_t-1 ! n_t is number of elements in radial direction in heating file
if ((pH%vals_H(1,kp).le.rc(k)).and.(pH%vals_H(1,kp+1).ge.rc(k))) then
pH%vec_magn_H(k)=pH%vals_H(2,kp+1)-(pH%vals_H(2,kp+1)-pH%vals_H(2,kp))&
& *(pH%vals_H(1,kp+1)-rc(k))/(pH%vals_H(1,kp+1)-pH%vals_H(1,kp))
......@@ -970,12 +980,20 @@ module initialize
endif
endif
if (pH%use_magn_H.eq.3.0_dp) then
pH%vec_magn_H(k)=pH%vec_magn_H(k)/cell%ref_r(1,1,k) ! in version 2 division by rho was already done by Kristina
! no angle-dependence, for average value divide by 2
pH%vec_magn_H(k)=0.5_dp*pH%vec_magn_H(k)/cell%ref_r(1,1,k) ! in V2, division by rho was already done by Kristina
elseif (pH%use_magn_H.ge.4.0_dp) then
! this is max value -> in thermal.f90 multiply with angle
pH%vec_magn_H(k)=pH%vec_magn_H(k)/cell%ref_r(1,1,k)
endif
Q = Q + pH%vec_magn_H(k) * cell%ref_r(1,1,k) * 4.0_dp*acos(-1.0_dp)/3.0_dp * &
& ((0.5_dp*(rc(k)+rc(k+1)))**3-(0.5_dp*(rc(k)+rc(k-1)))**3)
! Q here multiplied with dimensionless ref_r and radii, to account for variations
! with depth -> in main.f90 Q_eddy is then calculated by multiplying with D, rho_ref etc.
! only question: should it be time rho or not? Want to have pW/kg???
enddo
if (pH%use_magn_H.ge.4.0_dp) Q = Q/2.0_dp ! actually is angle-dependent, but divide here by 2 to get correct mantle-averaged Q
pH%Q_eddy = Q
else if (pH%use_magn_H.lt.0.0_dp) then
......@@ -1083,13 +1101,24 @@ module initialize
endif
endif
if (pH%use_magn_H.eq.3.0_dp) then
! no angle-dependence, for average value divide by 2
pH%vec_magn_H(k)=0.5_dp*pH%vec_magn_H(k)/cell%ref_r(1,1,k) ! in V2, division by rho was already done by Kristina
elseif (pH%use_magn_H.ge.4.0_dp) then
! this is max value -> in thermal.f90 multiply with angle
pH%vec_magn_H(k)=pH%vec_magn_H(k)/cell%ref_r(1,1,k)
endif
Q = Q + pH%vec_magn_H(k) * cell%ref_r(1,1,k) * 4.0_dp*acos(-1.0_dp)/3.0_dp * &
& ((0.5_dp*(rc(k)+rc(k+1)))**3-(0.5_dp*(rc(k)+rc(k-1)))**3)
enddo
!write(*,*) time/(1.0e+6_dp*pF%time_yr), Q*1.0e+7_dp*(pF%k*pF%DeltaT)*pF%D,"erg s^-1" ! *k*DT/(D^2*rho) [W/kg] * D^3*rho (mass mantle) [kg] -> W = 10^7 erg/s
if (pH%use_magn_H.ge.4.0_dp) Q = Q/2.0_dp ! actually is angle-dependent, but divide here by 2 to get correct mantle-averaged Q
pH%Q_eddy = Q
!write(*,*) Q*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(rc(nr)**3-rc(1)**3))*1.0e+12
endif
! otherwise nothing to do here
......
......@@ -611,16 +611,16 @@ program CHIC
call exchange_bnd(field,cell,mc,pL)
if ((pD%nr_ph.gt.0.or.pW%vel_r.gt.0.0_dp).and.pD%set_tracers.eq.1) then
! trace surface and bottom material exchange -> 0 on top, 1 at bottom, 0.5 inbetween
! cell%o(:,:,:) = 0.5_dp
! if ((pD%nr_ph.gt.0.or.pW%vel_r.gt.0.0_dp).and.pD%set_tracers.eq.1) then
! ! trace surface and bottom material exchange -> 0 on top, 1 at bottom, 0.5 inbetween
!! cell%o(:,:,:) = 0.5_dp
!! cell%o(:,:,mesh%nr-nint(0.1*mesh%nr)+1:mesh%nr+1) = 0.0_dp
!! cell%o(:,:,0:nint(0.1*mesh%nr)) = 1.0_dp
! ! 0 at top, 0.5 at the bottom, 1 inbetween
! cell%o(:,:,mesh%nr-nint(0.1*mesh%nr)+1:mesh%nr+1) = 0.0_dp
! cell%o(:,:,0:nint(0.1*mesh%nr)) = 1.0_dp
! 0 at top, 0.5 at the bottom, 1 inbetween
cell%o(:,:,mesh%nr-nint(0.1*mesh%nr)+1:mesh%nr+1) = 0.0_dp
cell%o(:,:,nint(0.1*mesh%nr)+1:mesh%nr-nint(0.1*mesh%nr)) = 1.0_dp
cell%o(:,:,0:nint(0.1*mesh%nr)) = 0.5_dp
endif
! cell%o(:,:,nint(0.1*mesh%nr)+1:mesh%nr-nint(0.1*mesh%nr)) = 1.0_dp
! cell%o(:,:,0:nint(0.1*mesh%nr)) = 0.5_dp
! endif
! set initial compositional values to particles
if (pP%use_part.eq.1) then
......@@ -1021,7 +1021,8 @@ program CHIC
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,0.0_dp,&!Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,heat_av,0.0_dp
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*pI%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,&
&heat_m,heat_av,0.0_dp
endif
endif
endif
......@@ -1344,7 +1345,8 @@ program CHIC
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,0.0_dp,&!Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,heat_av,0.0_dp
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*pI%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,&
&heat_m,heat_av,0.0_dp
endif
endif
......@@ -2222,7 +2224,7 @@ program CHIC
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,0.0_dp,&!Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,heat_av,&
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*pI%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,heat_av,&
&timeS2-timeS1
endif
! write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
......@@ -2847,8 +2849,8 @@ program CHIC
write(60,'(F13.10,I3,I3,F13.10,F13.10,F13.8,ES13.6,F13.5,ES13.6)') t/pF%time_yr*1.0e-9_dp,magma_local,magma_global,depl_av,&
& maxval(cell%d),((pG%rmax*pF%D)**3-(pG%rmin*pF%D)**3)/(3.0_dp*(pG%rmax*pF%D)**2) * &
& pX%CO2_ini*pM%Chi_extr*depl_av*pI%rho*pF%rho*pI%gravity*pF%g*1.0e-5_dp,&
& pH%Q_eddy/(0.0001_dp)*(pF%k*pF%DeltaT)*pF%D*1000.0_dp,& ! Q eddy current in erg/s
& T_mean*pF%DeltaT+pF%Ts, radiogenic_heating(pI%H0,pI%lambda,t,pI)*1.0e-7_dp*(pF%k*pF%DeltaT)*pF%D &
& pH%Q_eddy/(0.0001_dp)*(pF%k*pF%DeltaT)*pF%D*1000.0_dp,& ! Q eddy current in erg/s ! ToDo: Check units, should it not be times D^2???
& T_mean*pF%DeltaT+pF%Ts, radiogenic_heating(pI%H0,pI%lambda,t,pI)*1.0e-7_dp*(pF%k*pF%DeltaT)*pF%D & !Todo: or better TW for Q_eddy and radiogenic?
& * 4.0_dp/3.0_dp*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3)
! & pM%Chi_CO2*pM%Chi_extr*depl_av*pI%rho*pF%rho*pI%gravity*pF%g*1.0e-5_dp,&
! & pH%Q_eddy/(0.0001_dp)*(pF%k*pF%DeltaT)*pF%D*1000.0_dp,& ! Q eddy current in erg/s
......
......@@ -1206,7 +1206,7 @@ contains
if (cell%fld_d) mc%solidus(:,:,:)=inival
if (cell%fld_d) mc%DeltaF(:,:,:)=inival
if (cell%fld_d) mc%serp(:,:,:)=inival
if (cell%fld_o) cell%o(:,:,:)=inival
!if (cell%fld_o) cell%o(:,:,:)=inival
if (cell%fld_e) cell%entropy(:,:,:)=inival
if (cell%fld_r) cell%FP(:,:,:)=inival
if (cell%fld_r) cell%FG(:,:,:)=inival
......@@ -1342,7 +1342,7 @@ contains
if (cell%fld_d) mc%solidus(cl,cy,cr) = part%s(m)
if (cell%fld_d) mc%DeltaF(cl,cy,cr) = part%df(m)
if (cell%fld_d) mc%serp(cl,cy,cr) = part%sp(m)
if (cell%fld_o) cell%o(cl,cy,cr) = part%o(m)
!if (cell%fld_o) cell%o(cl,cy,cr) = part%o(m)
if (cell%fld_e) cell%entropy(cl,cy,cr) = part%e(m)
if (cell%fld_r) cell%FP(cl,cy,cr) = part%fp(m)
if (cell%fld_r) cell%FG(cl,cy,cr) = part%fg(m)
......@@ -1416,7 +1416,7 @@ contains
if (cell%fld_d) mc%solidus(cl,cy,cr) = mc%solidus(cl,cy,cr)+weight/part%s(m)
if (cell%fld_d) mc%DeltaF(cl,cy,cr) = mc%DeltaF(cl,cy,cr)+weight/part%df(m)
if (cell%fld_d) mc%serp(cl,cy,cr) = mc%serp(cl,cy,cr)+weight/part%sp(m)
if (cell%fld_o) cell%o(cl,cy,cr) = cell%o(cl,cy,cr)+weight/part%o(m)
!if (cell%fld_o) cell%o(cl,cy,cr) = cell%o(cl,cy,cr)+weight/part%o(m)
if (cell%fld_e) cell%entropy(cl,cy,cr) = cell%entropy(cl,cy,cr)+weight/part%e(m)
if (cell%fld_r) cell%FP(cl,cy,cr) = cell%FP(cl,cy,cr)+weight/part%fp(m)
if (cell%fld_r) cell%FG(cl,cy,cr) = cell%FG(cl,cy,cr)+weight/part%fg(m)
......@@ -1449,7 +1449,7 @@ contains
if (cell%fld_d) mc%solidus(cl,cy,cr) = mc%solidus(cl,cy,cr)+part%s(m)*weight
if (cell%fld_d) mc%DeltaF(cl,cy,cr) = mc%DeltaF(cl,cy,cr)+part%df(m)*weight
if (cell%fld_d) mc%serp(cl,cy,cr) = mc%serp(cl,cy,cr)+part%sp(m)*weight
if (cell%fld_o) cell%o(cl,cy,cr) = cell%o(cl,cy,cr)+part%o(m)*weight
!if (cell%fld_o) cell%o(cl,cy,cr) = cell%o(cl,cy,cr)+part%o(m)*weight
if (cell%fld_e) cell%entropy(cl,cy,cr) = cell%entropy(cl,cy,cr)+part%e(m)*weight
if (cell%fld_r) cell%FP(cl,cy,cr) = cell%FP(cl,cy,cr)+part%fp(m)*weight
if (cell%fld_r) cell%FG(cl,cy,cr) = cell%FG(cl,cy,cr)+part%fg(m)*weight
......@@ -1482,7 +1482,7 @@ contains
if (cell%fld_d) mc%solidus(cl,cy,cr) = mc%solidus(cl,cy,cr)*part%s(m)**weight
if (cell%fld_d) mc%DeltaF(cl,cy,cr) = mc%DeltaF(cl,cy,cr)*part%df(m)**weight
if (cell%fld_d) mc%serp(cl,cy,cr) = mc%serp(cl,cy,cr)*part%sp(m)**weight
if (cell%fld_o) cell%o(cl,cy,cr) = cell%o(cl,cy,cr)*part%o(m)**weight
!if (cell%fld_o) cell%o(cl,cy,cr) = cell%o(cl,cy,cr)*part%o(m)**weight
if (cell%fld_e) cell%entropy(cl,cy,cr)= cell%entropy(cl,cy,cr)*part%e(m)**weight
if (cell%fld_r) cell%FP(cl,cy,cr) = cell%FP(cl,cy,cr)*part%fp(m)**weight
if (cell%fld_r) cell%FG(cl,cy,cr) = cell%FG(cl,cy,cr)*part%fg(m)**weight
......@@ -1537,7 +1537,7 @@ contains
if (cell%fld_d) mc%solidus(i,j,k) = dummy_s(i,j,k)
if (cell%fld_d) mc%DeltaF(i,j,k) = dummy_df(i,j,k)
if (cell%fld_d) mc%serp(i,j,k) = dummy_sp(i,j,k)
if (cell%fld_o) cell%o(i,j,k) = dummy_o(i,j,k)
!if (cell%fld_o) cell%o(i,j,k) = dummy_o(i,j,k)
if (cell%fld_e) cell%entropy(i,j,k) = dummy_e(i,j,k)
if (cell%fld_r) cell%FP(i,j,k) = dummy_fp(i,j,k)
if (cell%fld_r) cell%FG(i,j,k) = dummy_fg(i,j,k)
......@@ -1595,7 +1595,7 @@ contains
if (cell%fld_d) mc%solidus(i,j,k) = vec_w(i,j,k)/mc%solidus(i,j,k)
if (cell%fld_d) mc%DeltaF(i,j,k) = vec_w(i,j,k)/mc%DeltaF(i,j,k)
if (cell%fld_d) mc%serp(i,j,k) = vec_w(i,j,k)/mc%serp(i,j,k)
if (cell%fld_o) cell%o(i,j,k) = vec_w(i,j,k)/cell%o(i,j,k)
!if (cell%fld_o) cell%o(i,j,k) = vec_w(i,j,k)/cell%o(i,j,k)
if (cell%fld_e) cell%entropy(i,j,k) = vec_w(i,j,k)/cell%entropy(i,j,k)
if (cell%fld_r) cell%FP(i,j,k) = vec_w(i,j,k)/cell%FP(i,j,k)
if (cell%fld_r) cell%FG(i,j,k) = vec_w(i,j,k)/cell%FG(i,j,k)
......@@ -1628,7 +1628,7 @@ contains
if (cell%fld_d) mc%solidus(i,j,k) = mc%solidus(i,j,k)/vec_w(i,j,k)
if (cell%fld_d) mc%DeltaF(i,j,k) = mc%DeltaF(i,j,k)/vec_w(i,j,k)
if (cell%fld_d) mc%serp(i,j,k) = mc%serp(i,j,k)/vec_w(i,j,k)
if (cell%fld_o) cell%o(i,j,k) = cell%o(i,j,k)/vec_w(i,j,k)
!if (cell%fld_o) cell%o(i,j,k) = cell%o(i,j,k)/vec_w(i,j,k)
if (cell%fld_e) cell%entropy(i,j,k) = cell%entropy(i,j,k)/vec_w(i,j,k)
if (cell%fld_r) cell%FP(i,j,k) = cell%FP(i,j,k)/vec_w(i,j,k)
if (cell%fld_r) cell%FG(i,j,k) = cell%FG(i,j,k)/vec_w(i,j,k)
......@@ -1661,7 +1661,7 @@ contains
if (cell%fld_d) mc%solidus(i,j,k) = mc%solidus(i,j,k)**vec_w(i,j,k)
if (cell%fld_d) mc%DeltaF(i,j,k) = mc%DeltaF(i,j,k)**vec_w(i,j,k)
if (cell%fld_d) mc%serp(i,j,k) = mc%serp(i,j,k)**vec_w(i,j,k)
if (cell%fld_o) cell%o(i,j,k) = cell%o(i,j,k)**vec_w(i,j,k)
!if (cell%fld_o) cell%o(i,j,k) = cell%o(i,j,k)**vec_w(i,j,k)
if (cell%fld_e) cell%entropy(i,j,k) = cell%entropy(i,j,k)**vec_w(i,j,k)
if (cell%fld_r) cell%FP(i,j,k) = cell%FP(i,j,k)**vec_w(i,j,k)
if (cell%fld_r) cell%FG(i,j,k) = cell%FG(i,j,k)**vec_w(i,j,k)
......@@ -1726,8 +1726,8 @@ contains
if (cell%fld_d) mc%DeltaF(0:,:,nr+1) = mc%DeltaF(0:,:,nr)
if (cell%fld_d) mc%serp(0:,:,0) = mc%serp(0:,:,1)
if (cell%fld_d) mc%serp(0:,:,nr+1) = mc%serp(0:,:,nr)
if (cell%fld_o) cell%o(0:,:,0) = cell%o(0:,:,1)
if (cell%fld_o) cell%o(0:,:,nr+1) = cell%o(0:,:,nr)
!if (cell%fld_o) cell%o(0:,:,0) = cell%o(0:,:,1)
!if (cell%fld_o) cell%o(0:,:,nr+1) = cell%o(0:,:,nr)
if (cell%fld_e) cell%entropy(0:,:,0) = cell%entropy(0:,:,1)
if (cell%fld_e) cell%entropy(0:,:,nr+1) = cell%entropy(0:,:,nr)
if (cell%fld_r) cell%FP(0:,:,0) = cell%FP(0:,:,1)
......@@ -1792,8 +1792,8 @@ contains
if (cell%fld_d) mc%DeltaF(nl+1,:,:) = mc%DeltaF(1,:,:)
if (cell%fld_d) mc%serp(0,:,:) = mc%serp(nl,:,:)
if (cell%fld_d) mc%serp(nl+1,:,:) = mc%serp(1,:,:)
if (cell%fld_o) cell%o(0,:,:) = cell%o(nl,:,:)
if (cell%fld_o) cell%o(nl+1,:,:) = cell%o(1,:,:)
!if (cell%fld_o) cell%o(0,:,:) = cell%o(nl,:,:)
!if (cell%fld_o) cell%o(nl+1,:,:) = cell%o(1,:,:)
if (cell%fld_e) cell%entropy(0,:,:) = cell%entropy(nl,:,:)
if (cell%fld_e) cell%entropy(nl+1,:,:) = cell%entropy(1,:,:)
if (cell%fld_r) cell%FP(0,:,:) = cell%FP(nl,:,:)
......@@ -1859,8 +1859,8 @@ contains
if (cell%fld_d) mc%DeltaF(:,nl+1,:) = mc%DeltaF(:,1,:)
if (cell%fld_d) mc%serp(:,0,:) = mc%serp(:,nl,:)
if (cell%fld_d) mc%serp(:,nl+1,:) = mc%serp(:,1,:)
if (cell%fld_o) cell%o(:,0,:) = cell%o(:,nl,:)
if (cell%fld_o) cell%o(:,nl+1,:) = cell%o(:,1,:)
!if (cell%fld_o) cell%o(:,0,:) = cell%o(:,nl,:)
!if (cell%fld_o) cell%o(:,nl+1,:) = cell%o(:,1,:)
if (cell%fld_e) cell%entropy(:,0,:) = cell%entropy(:,nl,:)
if (cell%fld_e) cell%entropy(:,nl+1,:) = cell%entropy(:,1,:)
if (cell%fld_r) cell%FP(:,0,:) = cell%FP(:,nl,:)
......@@ -1927,8 +1927,8 @@ contains
if (cell%fld_d) mc%DeltaF(nl+1,:,0:) = mc%DeltaF(nl,:,0:)
if (cell%fld_d) mc%serp(0,:,0:) = mc%serp(1,:,0:)
if (cell%fld_d) mc%serp(nl+1,:,0:) = mc%serp(nl,:,0:)
if (cell%fld_o) cell%o(0,:,0:) = cell%o(1,:,0:)
if (cell%fld_o) cell%o(nl+1,:,0:) = cell%o(nl,:,0:)
!if (cell%fld_o) cell%o(0,:,0:) = cell%o(1,:,0:)
!if (cell%fld_o) cell%o(nl+1,:,0:) = cell%o(nl,:,0:)
if (cell%fld_e) cell%entropy(0,:,0:) = cell%entropy(1,:,0:)
if (cell%fld_e) cell%entropy(nl+1,:,0:) = cell%entropy(nl,:,0:)
if (cell%fld_r) cell%FP(0,:,0:) = cell%FP(1,:,0:)
......@@ -1994,8 +1994,8 @@ contains
if (cell%fld_d) mc%DeltaF(:,nl+1,:) = mc%DeltaF(:,nl,:)
if (cell%fld_d) mc%serp(:,0,:) = mc%serp(:,1,:)
if (cell%fld_d) mc%serp(:,nl+1,:) = mc%serp(:,nl,:)
if (cell%fld_o) cell%o(:,0,:) = cell%o(:,1,:)
if (cell%fld_o) cell%o(:,nl+1,:) = cell%o(:,nl,:)
!if (cell%fld_o) cell%o(:,0,:) = cell%o(:,1,:)
!if (cell%fld_o) cell%o(:,nl+1,:) = cell%o(:,nl,:)
if (cell%fld_e) cell%entropy(:,0,:) = cell%entropy(:,1,:)
if (cell%fld_e) cell%entropy(:,nl+1,:) = cell%entropy(:,nl,:)
if (cell%fld_r) cell%FP(:,0,:) = cell%FP(:,1,:)
......
......@@ -570,7 +570,7 @@ contains
type(visc),intent(inout)::eta
type(variables_unknowns),intent(in)::field
type(mesh_cp),intent(in)::mesh
type(cell_properties),intent(in)::cell
type(cell_properties),intent(inout)::cell ! out only needed for cell%o to trace buoyancy
type(melt_properties),intent(in)::mc
integer,intent(in)::mom_DuDt
real(dp),intent(in)::dt
......@@ -1871,6 +1871,7 @@ contains
endif
else
! add lithostatic pressure variations, body force term without alpha/temperature since already local density (no ref profile used)
! ToDo: Check Vorzeichen, sollte es nicht "+" sein???
BoussT = -RaB * (refR+refRu) * (g_t+g_u) * (FPt+FPu) / 8.0_dp + RaB * (cell%ref_p(i,j,k)-cell%ref_p(i,j,k-1))*dr
! new formulation, where actual fields are used for rho and p, and no reference profiles are used (hence artifical reference profile created here with rho=rho_ref*(1+alpha T):
! BoussT = Ra*( refR*g_t*refA*FPt*(field%T(i,j,k)+T0)*mesh%dVi(i,j,k)/(1.0_dp+refA*(field%T(i,j,k)+T0)) &
......@@ -1916,9 +1917,13 @@ contains
endif
if (Ra.ne.0.0_dp) BoussC = BoussC * Ra ! Ra can be zero for benchmarks, i.e. van Keken et al., 1997
!cell%o(i,j,k) = av_eta * (BoussC - BoussT)
! WEITER
!write(*,*) i, j ,k, BoussC, -BoussT, pI%drho, D_t, C, C_ref, pI%drho*D_t, C-C_ref-pI%drho*D_t
!if (D_t>0.0_dp) then
! write(*,*) i, j ,k, BoussC, BoussT, pI%drho, D_t, C, C_ref, pI%drho*D_t, C-C_ref-pI%drho*D_t
!endif
! 29 1 56 -24720775.145425387 -21375563.027024981 0.15151515151515152 0.10919216793520699 1.0000000000000000 1.0000000000000000 1.6544267868970758E-002 -1.6544267868970758E-002
! 16 1 56 -18147597.726186022 -21134632.330163386 0.15151515151515152 0.0000000000000000 1.0000000000000000 1.0000000000000000 0.0000000000000000 0.0000000000000000
! 16 1 55 -25416875.878874149 -27778944.839950610 0.15151515151515152 0.29999999999999999 1.0000000000000000 1.0000000000000000 4.5454545454545456E-002 -4.5454545454545456E-002
......
......@@ -703,13 +703,21 @@ contains
! Eddy current heating induced from the star
! if (pH%use_magn_H.ne.0.0_dp) source = source + pH%vec_magn_H(k)
if (pH%use_magn_H.ne.0.0_dp) then
if (pH%use_magn_H.le.2.0_dp) then
if (pH%use_magn_H.le.3.0_dp) then
! average heating independent of lattitude
source = source + pH%vec_magn_H(k)
cell%o(i,j,k) = pH%vec_magn_H(k)
else
! cos curve multiplied to have max value at equator, so that heating becomes 0 at poles (+90°/-90°)
! only works for spherical application, since in 2D box no angle given
source = source + pH%vec_magn_H(k)*cos(mesh%lc(i))
cell%o(i,j,k) = pH%vec_magn_H(k)*cos(mesh%lc(i))
!write(*,*) "i=",i,", lc=",mesh%lc(i),", cos(lc)=",cos(mesh%lc(i)),", H=",pH%vec_magn_H(k)*cos(mesh%lc(i))
! Check correct: for half sphere, lc goes from -i/2 to +pi/2; heating goes from 0 to max value to 0 (hence from pole
! via equ. to pole)
endif
! dimensionalize cell%o here already (anyway just output field:
cell%o(i,j,k) = cell%o(i,j,k)*pF%k*pF%DeltaT/(pF%rho*pF%D**2)
endif
! two-phase flow (set LatHeat!=0!!!)
......@@ -1462,13 +1470,17 @@ contains
if (pH%TidHeat.ne.0) source = source + tidal_heating(pH,pI,mesh%rc(k),eta(i,j,k),cell%shear(i,j,k),mesh%lc(i),pG)
! Eddy current heating induced from the star
if (pH%use_magn_H.ne.0.0_dp) then
if (pH%use_magn_H.le.2.0_dp) then
if (pH%use_magn_H.le.3.0_dp) then
! average heating independent of lattitude
source = source + pH%vec_magn_H(k)
source = source + pH%vec_magn_H(k)
cell%o(i,j,k) = pH%vec_magn_H(k)
else
! sin curve multiplied to max value at equator, so that heating becomes 0 at poles (+90°/-90°)
source = source + pH%vec_magn_H(k)*cos(mesh%lc(i))
cell%o(i,j,k) = pH%vec_magn_H(k)*cos(mesh%lc(i))
endif
! dimensionalize cell%o here already (anyway just output field:
cell%o(i,j,k) = cell%o(i,j,k)*pF%k*pF%DeltaT/(pF%rho*pF%D**2)
endif
!------------ Melt inclusion --------------
......@@ -2541,13 +2553,17 @@ contains
! Eddy current heating induced from the star
! if (pH%use_magn_H.ne.0.0_dp) source = source + pH%vec_magn_H(k)
if (pH%use_magn_H.ne.0.0_dp) then
if (pH%use_magn_H.le.2.0_dp) then
if (pH%use_magn_H.le.3.0_dp) then
! average heating independent of lattitude
source = source + pH%vec_magn_H(k)
cell%o(i,j,k) = pH%vec_magn_H(k)
else
! sin curve multiplied to max value at equator, so that heating becomes 0 at poles (+90°/-90°)
source = source + pH%vec_magn_H(k)*cos(mesh%lc(i))
cell%o(i,j,k) = pH%vec_magn_H(k)*cos(mesh%lc(i))
endif
! dimensionalize cell%o here already (anyway just output field:
cell%o(i,j,k) = cell%o(i,j,k)*pF%k*pF%DeltaT/(pF%rho*pF%D**2)
endif
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment