Commit 907346b2 authored by Lena Noack's avatar Lena Noack
Browse files

Separate recycling and mineral phase property (new property mp added)

parent 02b0c809
......@@ -97,7 +97,8 @@ module initialize
cell%fa = pI%mantle_fr_angl ! friction angle
cell%fc = pI%mantle_fr_coh ! cohesion
cell%g = pI%grain_size
cell%r = 0.0_dp ! recycling value for tracing
cell%r = 0.0_dp ! recycling value for tracing crust evolution
cell%MP = 0.0_dp ! mineral phase group
cell%cp= pI%Cp
cell%k = pI%k
cell%h = pI%H0
......
......@@ -684,7 +684,7 @@ program CHIC
if (pD%nr_ph.gt.0.and.pD%set_tracers.ge.2) then
if (pD%set_tracers.eq.2) then
! trace material exchange -> set to initial material value, will not be updated later on
part%o = part%r
part%o = part%mp
else
! trace material only from later time step on (with set_tracers being the i-th time step)
part%o = 0.0_dp
......@@ -888,7 +888,7 @@ program CHIC
if (pD%nr_ph.gt.0) then ! ignore everything above and set viscosity on each phase layer seperately
if ((pD%A_dif(1).ne.0.0_dp).or.(pD%A_dis(1).ne.0.0_dp)) then
call SetPhaseVisc(jmin_2D,eta%CV,field%strain_rate,cell%r,cell%ref_p,mesh%rc,field%T,pT%T0,pD,pF,pV,Lbound(cell%w,1),&
call SetPhaseVisc(jmin_2D,eta%CV,field%strain_rate,cell%MP,cell%ref_p,mesh%rc,field%T,pT%T0,pD,pF,pV,Lbound(cell%w,1),&
&Lbound(cell%g,1),cell%w,cell%wr,cell%g)
endif
else
......@@ -1380,7 +1380,7 @@ program CHIC
if (pD%nr_ph.gt.0.and.pD%set_tracers.gt.2.and.pD%set_tracers.eq.i) then
! trace material exchange -> set to initial material value, will not be updated later on
part%o = part%r
part%o = part%mp
endif
......@@ -1423,7 +1423,8 @@ program CHIC
if (cell%fld_c) write(*,*) "C",minval(cell%C_vOld),maxval(cell%C_vOld)
if (cell%fld_w) write(*,*) "w",minval(cell%w),maxval(cell%w)
if (cell%fld_g) write(*,*) "g",minval(cell%g),maxval(cell%g)
if (cell%fld_m) write(*,*) "r",minval(cell%r),maxval(cell%r)
if (cell%fld_r) write(*,*) "r",minval(cell%r),maxval(cell%r)
if (cell%fld_m) write(*,*) "r",minval(cell%MP),maxval(cell%MP)
if (cell%fld_k) write(*,*) "k",minval(cell%k),maxval(cell%k)
if (cell%fld_h) write(*,*) "h",minval(cell%h),maxval(cell%h)
if (cell%fld_v) write(*,*) "v",minval(cell%v),maxval(cell%v)
......@@ -1463,7 +1464,7 @@ program CHIC
!!!!!!!!!!!!!!!!!!!!!!!
if (pD%nr_ph.ne.0) then ! more than one phase material
if ((pD%A_dif(1).ne.0.0_dp).or.(pD%A_dis(1).ne.0.0_dp))&
& call SetPhaseVisc(jmin_2D,eta%CV,field%strain_rate,cell%r,cell%ref_p,mesh%rc,field%T,pT%T0,&
& call SetPhaseVisc(jmin_2D,eta%CV,field%strain_rate,cell%MP,cell%ref_p,mesh%rc,field%T,pT%T0,&
&pD,pF,pV,Lbound(cell%w,1),Lbound(cell%g,1),cell%w,cell%wr,cell%g)
call Visc_ph(jmin_2D,eta%CV,mesh%rc,field%T,pD,pF) ! SetPhaseVisc sets viscocity with rheology law, here multiply in addition viscosity jump if needed
else
......@@ -2050,7 +2051,7 @@ program CHIC
if(rank.eq.0) call print_time(pN%pr_time,timeD1,timeD2,"Part/Post-C4P")
!if ((abs(pI%compress).ge.3).or.(pD%nr_ph.ne.0)) &
if ((pW%weak_angle.eq.0.0_dp).and.(pW%BenchSubdFreeS.eq.0).and.(pD%nr_ph.ne.0)) &
& call get_phases(part,cell,jmin_2D,field%T,field%w,pD,pF,pI,mesh%rmax,pE%Psurf,pI%compress,pN%debug,t) !SetMatProp(part,jmin_2D,cell%r,mesh%rc,field%T,pE%Psurf,pD,pF,pI%compress)
& call get_phases(part,cell,jmin_2D,field%T,field%w,pD,pF,pI,mesh%rmax,pE%Psurf,pI%compress,pN%debug,t) !SetMatProp(part,jmin_2D,cell%MP,mesh%rc,field%T,pE%Psurf,pD,pF,pI%compress)
call part_part2cell(part,mc,mesh,cell,pG%geom,pN%debug,pG%periodic,pN%average,pN%p_average,pN%p_average2,&
&pN%pr_time,pP%moveRef) ! transfer particle properties to cells
endif
......@@ -2069,7 +2070,7 @@ program CHIC
! Calculate melt fraction for T^(n-1) (to be used in energy solver for T^(n))
if (pP%use_part.eq.1) then
call update_solidus_p(mesh,mc,part,pM) ! needs to be after particle_move to get correct delta_solidus
call melt_comp_p(field,cell,mesh,mc,part,pF,pM) ! melt fraction and depletion
call melt_comp_p(field,cell,mesh,mc,part,pF,pM,pV,pX) ! melt fraction and depletion
!call dehydration(part,mesh,cell,mc%mass_H2O,mc%total_mass_H2O,pG%dim,pM%dehydr_melt)
!call atmos_outgassing(part,mesh,cell,mc,mass_H2O,total_mass_H2O,di,mass_CO2,total_mass_CO2,X_CO2,fO2,CO2,fwm)
call update_solidus_p(mesh,mc,part,pM) ! just for plots
......
This diff is collapsed.
......@@ -49,6 +49,7 @@ module particles
real(dp),allocatable::e(:) !< entropy of the material
real(dp),allocatable::fp(:) !< factor phases for thermal and stokes solver when using phase transitions
real(dp),allocatable::fg(:) !< factor phases times clapeyron for thermal solver when using phase transitions
real(dp),allocatable::mp(:) !< mineral phases
real(dp),allocatable::cl(:) !< clapeyron slope for phase transitions, used in thermal solver
real(dp),allocatable::CO2(:) !< CO2 content of cell
real(dp),allocatable::X_CO2(:) !< CO2 concentration in melt (needed to have cell value, not stored in snapshots)
......@@ -114,6 +115,7 @@ contains
if (cell%fld_e) then ; allocate( part%e(1:part%n)) ; else ; allocate( part%e(1) ) ; endif
if (cell%fld_r) then ; allocate( part%fp(1:part%n)) ; else ; allocate( part%fp(1) ) ; endif
if (cell%fld_r) then ; allocate( part%fg(1:part%n)) ; else ; allocate( part%fg(1) ) ; endif
if (cell%fld_r) then ; allocate( part%mp(1:part%n)) ; else ; allocate( part%mp(1) ) ; endif
if (cell%fld_r) then ; allocate( part%cl(1:part%n)) ; else ; allocate( part%cl(1) ) ; endif
if (cell%fld_x) then ; allocate( part%CO2(1:part%n)) ; else ; allocate( part%CO2(1) ) ; endif
if (cell%fld_x) then ; allocate( part%X_CO2(1:part%n)) ; else ; allocate( part%X_CO2(1) ) ; endif
......@@ -131,8 +133,8 @@ contains
deallocate( part%pos_l, part%pos_y, part%pos_r, part%cell_l, part%cell_y, part%cell_r, part%inactive )
deallocate( part%c, part%w, part%g, part%r, part%o, part%cp, part%a, part%rho, part%t, part%k, part%h, &
& part%v, part%wr, part%fa,part%fc,part%d,part%s,part%df, part%sp, part%e, part%fp, part%fg, part%cl, &
& part%CO2, part%X_CO2, part%X_H2O, part%X_U238, part%X_U235, part%X_Th232, part%X_K40 )
& part%v, part%wr, part%fa,part%fc,part%d,part%s,part%df, part%sp, part%e, part%fp, part%fg, part%mp,&
& part%cl, part%CO2, part%X_CO2, part%X_H2O, part%X_U238, part%X_U235, part%X_Th232, part%X_K40 )
end subroutine part_deallocate
......@@ -218,6 +220,7 @@ contains
part%e = 1.0_dp
part%fp = 1.0_dp
part%fg = 0.0_dp
part%mp = 0.0_dp
part%cl = 0.0_dp
part%CO2 = 0.0_dp
part%X_CO2 = 0.0_dp
......@@ -256,6 +259,7 @@ contains
if (cell%fld_e) part%e(m) = cell%entropy(i,j,k)
if (cell%fld_r) part%fp(m) = cell%FP(i,j,k)
if (cell%fld_r) part%fg(m) = cell%FG(i,j,k)
if (cell%fld_r) part%mp(m) = cell%MP(i,j,k)
if (cell%fld_r) part%cl(m) = cell%CL(i,j,k)
if (cell%fld_x) part%CO2(m) = cell%CO2(i,j,k)
if (cell%fld_x) part%X_CO2(m) = cell%X_CO2(i,j,k)
......@@ -293,6 +297,7 @@ contains
write(*,*) "e :",part%e(1)
write(*,*) "fp:",part%fp(1)
write(*,*) "fg:",part%fg(1)
write(*,*) "mp:",part%mp(1)
write(*,*) "cl:",part%cl(1)
write(*,*) "CO:",part%CO2(1)
write(*,*) "XC:",part%X_CO2(1)
......@@ -340,6 +345,7 @@ contains
if (LBOUND(part%e,1).ne.UBOUND(part%e,1)) size_arr = size_arr+1
if (LBOUND(part%fp,1).ne.UBOUND(part%fp,1)) size_arr = size_arr+1
if (LBOUND(part%fg,1).ne.UBOUND(part%fg,1)) size_arr = size_arr+1
if (LBOUND(part%mp,1).ne.UBOUND(part%mp,1)) size_arr = size_arr+1
if (LBOUND(part%cl,1).ne.UBOUND(part%cl,1)) size_arr = size_arr+1
if (LBOUND(part%CO2,1).ne.UBOUND(part%CO2,1)) size_arr = size_arr+1
! part%X_CO2/part%X_H2O not stored in particle snapshot
......@@ -398,6 +404,7 @@ contains
if (LBOUND(part%e,1).ne.UBOUND(part%e,1)) then ; values(ind) = part%e(m) ; ind=ind+1 ; endif
if (LBOUND(part%fp,1).ne.UBOUND(part%fp,1)) then ; values(ind) = part%fp(m) ; ind=ind+1 ; endif
if (LBOUND(part%fg,1).ne.UBOUND(part%fg,1)) then ; values(ind) = part%fg(m) ; ind=ind+1 ; endif
if (LBOUND(part%mp,1).ne.UBOUND(part%mp,1)) then ; values(ind) = part%mp(m) ; ind=ind+1 ; endif
if (LBOUND(part%cl,1).ne.UBOUND(part%cl,1)) then ; values(ind) = part%cl(m) ; ind=ind+1 ; endif
if (LBOUND(part%CO2,1).ne.UBOUND(part%CO2,1)) then ; values(ind) = part%CO2(m) ; ind=ind+1 ; endif
if (LBOUND(part%X_U238,1).ne.UBOUND(part%X_U238,1)) then ; values(ind) = part%X_U238(m) ; ind=ind+1 ; endif
......@@ -459,6 +466,7 @@ contains
if (LBOUND(part%e,1).ne.UBOUND(part%e,1)) then ; values(ind) = part%e(m) ; ind=ind+1 ; endif
if (LBOUND(part%fp,1).ne.UBOUND(part%fp,1)) then ; values(ind) = part%fp(m) ; ind=ind+1 ; endif
if (LBOUND(part%fg,1).ne.UBOUND(part%fg,1)) then ; values(ind) = part%fg(m) ; ind=ind+1 ; endif
if (LBOUND(part%mp,1).ne.UBOUND(part%mp,1)) then ; values(ind) = part%mp(m) ; ind=ind+1 ; endif
if (LBOUND(part%cl,1).ne.UBOUND(part%cl,1)) then ; values(ind) = part%cl(m) ; ind=ind+1 ; endif
if (LBOUND(part%CO2,1).ne.UBOUND(part%CO2,1)) then ; values(ind) = part%CO2(m) ; ind=ind+1 ; endif
if (LBOUND(part%X_U238,1).ne.UBOUND(part%X_U238,1)) then ; values(ind) = part%X_U238(m) ; ind=ind+1 ; endif
......@@ -515,6 +523,7 @@ contains
if (LBOUND(part%e,1).ne.UBOUND(part%e,1)) size_arr = size_arr+1
if (LBOUND(part%fp,1).ne.UBOUND(part%fp,1)) size_arr = size_arr+1
if (LBOUND(part%fg,1).ne.UBOUND(part%fg,1)) size_arr = size_arr+1
if (LBOUND(part%mp,1).ne.UBOUND(part%mp,1)) size_arr = size_arr+1
if (LBOUND(part%cl,1).ne.UBOUND(part%cl,1)) size_arr = size_arr+1
if (LBOUND(part%CO2,1).ne.UBOUND(part%CO2,1)) size_arr = size_arr+1
if (LBOUND(part%X_U238,1).ne.UBOUND(part%X_U238,1)) size_arr = size_arr+1
......@@ -571,6 +580,7 @@ contains
if (LBOUND(part%e,1).ne.UBOUND(part%e,1)) then ; part%e(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%fp,1).ne.UBOUND(part%fp,1)) then ; part%fp(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%fg,1).ne.UBOUND(part%fg,1)) then ; part%fg(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%mp,1).ne.UBOUND(part%mp,1)) then ; part%mp(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%cl,1).ne.UBOUND(part%cl,1)) then ; part%cl(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%CO2,1).ne.UBOUND(part%CO2,1)) then ; part%CO2(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%X_U238,1).ne.UBOUND(part%X_U238,1)) then ; part%X_U238(m) = values(ind) ; ind=ind+1 ; endif
......@@ -622,6 +632,7 @@ contains
if (LBOUND(part%e,1).ne.UBOUND(part%e,1)) then ; part%e(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%fp,1).ne.UBOUND(part%fp,1)) then ; part%fp(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%fg,1).ne.UBOUND(part%fg,1)) then ; part%fg(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%mp,1).ne.UBOUND(part%mp,1)) then ; part%mp(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%cl,1).ne.UBOUND(part%cl,1)) then ; part%cl(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%CO2,1).ne.UBOUND(part%CO2,1)) then ; part%CO2(m) = values(ind) ; ind=ind+1 ; endif
if (LBOUND(part%X_U238,1).ne.UBOUND(part%X_U238,1)) then ; part%X_U238(m) = values(ind) ; ind=ind+1 ; endif
......@@ -1016,6 +1027,7 @@ contains
if (cell%fld_e) part%e(m) = cell%entropy(i,j,k)
if (cell%fld_r) part%fp(m) = cell%FP(i,j,k)
if (cell%fld_r) part%fg(m) = cell%FG(i,j,k)
if (cell%fld_r) part%mp(m) = cell%MP(i,j,k)
if (cell%fld_r) part%cl(m) = cell%CL(i,j,k)
if (cell%fld_x) part%CO2(m) = cell%CO2(i,j,k)
if (cell%fld_x) part%X_CO2(m) = cell%X_CO2(i,j,k)
......@@ -1060,7 +1072,7 @@ contains
real(dp),allocatable::vec_w(:,:,:),vec_w2(:,:,:),vec_w3(:,:,:),dummy_c(:,:,:),dummy_w(:,:,:),dummy_wr(:,:,:),dummy_fa(:,:,:), &
& dummy_fc(:,:,:),dummy_g(:,:,:),dummy_r(:,:,:),dummy_cp(:,:,:),dummy_a(:,:,:),dummy_rho(:,:,:),dummy_t(:,:,:),&
& dummy_k(:,:,:),dummy_h(:,:,:),dummy_v(:,:,:),dummy_d(:,:,:),dummy_s(:,:,:),dummy_df(:,:,:),dummy_o(:,:,:), &
& dummy_e(:,:,:),dummy_fp(:,:,:),dummy_fg(:,:,:),dummy_cl(:,:,:),dummy_CO2(:,:,:),dummy_X_CO2(:,:,:), &
& dummy_e(:,:,:),dummy_fp(:,:,:),dummy_fg(:,:,:),dummy_mp(:,:,:),dummy_cl(:,:,:),dummy_CO2(:,:,:),dummy_X_CO2(:,:,:), &
& dummy_X_H2O(:,:,:), dummy_X_U238(:,:,:),dummy_X_U235(:,:,:),dummy_X_Th232(:,:,:),dummy_X_K40(:,:,:)
real(dp),allocatable::dummy_sp(:,:,:)
real(dp),allocatable::weight_vec(:)
......@@ -1109,6 +1121,7 @@ contains
if (cell%fld_e) then ; allocate( dummy_e(0:nl+1,jmin:jmax,0:nr+1) ) ; else ; allocate( dummy_e(1,1,1) ) ; endif
if (cell%fld_r) then ; allocate( dummy_fp(0:nl+1,jmin:jmax,0:nr+1) ) ; else ; allocate( dummy_fp(1,1,1) ) ; endif
if (cell%fld_r) then ; allocate( dummy_fg(0:nl+1,jmin:jmax,0:nr+1) ) ; else ; allocate( dummy_fg(1,1,1) ) ; endif
if (cell%fld_r) then ; allocate( dummy_mp(0:nl+1,jmin:jmax,0:nr+1) ) ; else ; allocate( dummy_mp(1,1,1) ) ; endif
if (cell%fld_r) then ; allocate( dummy_cl(0:nl+1,jmin:jmax,0:nr+1) ) ; else ; allocate( dummy_cl(1,1,1) ) ; endif
if (cell%fld_x) then ; allocate( dummy_CO2(0:nl+1,jmin:jmax,0:nr+1) ) ; else ; allocate( dummy_CO2(1,1,1) ) ; endif
if (cell%fld_x) then ; allocate( dummy_X_CO2(0:nl+1,jmin:jmax,0:nr+1) ) ; else ; allocate( dummy_X_CO2(1,1,1) ) ; endif
......@@ -1148,6 +1161,7 @@ contains
if (cell%fld_e) dummy_e(:,:,:)=cell%entropy(:,:,:)
if (cell%fld_r) dummy_fp(:,:,:)=cell%FP(:,:,:)
if (cell%fld_r) dummy_fg(:,:,:)=cell%FG(:,:,:)
if (cell%fld_r) dummy_mp(:,:,:)=cell%MP(:,:,:)
if (cell%fld_r) dummy_cl(:,:,:)=cell%CL(:,:,:)
if (cell%fld_x) dummy_CO2(:,:,:)=cell%CO2(:,:,:)
if (cell%fld_x) dummy_X_CO2(:,:,:)=0.0_dp !cell%X_CO2(:,:,:)
......@@ -1196,6 +1210,7 @@ contains
if (cell%fld_e) cell%entropy(:,:,:)=inival
if (cell%fld_r) cell%FP(:,:,:)=inival
if (cell%fld_r) cell%FG(:,:,:)=inival
if (cell%fld_r) cell%MP(:,:,:)=inival
if (cell%fld_r) cell%CL(:,:,:)=inival
if (cell%fld_x) cell%CO2(:,:,:)=inival
if (cell%fld_x) cell%X_CO2(:,:,:)=inival
......@@ -1331,6 +1346,7 @@ contains
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)
if (cell%fld_r) cell%MP(cl,cy,cr) = part%mp(m)
if (cell%fld_r) cell%CL(cl,cy,cr) = part%cl(m)
if (cell%fld_x) cell%CO2(cl,cy,cr) = part%CO2(m)
if (cell%fld_x) cell%X_CO2(cl,cy,cr) = part%X_CO2(m)
......@@ -1404,6 +1420,7 @@ contains
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)
if (cell%fld_r) cell%MP(cl,cy,cr) = cell%MP(cl,cy,cr)+weight/part%mp(m)
if (cell%fld_r) cell%CL(cl,cy,cr) = cell%CL(cl,cy,cr)+weight/part%cl(m)
if (cell%fld_x) cell%CO2(cl,cy,cr) = cell%CO2(cl,cy,cr)+weight/part%CO2(m)
if (cell%fld_x) cell%X_CO2(cl,cy,cr) = cell%X_CO2(cl,cy,cr)+weight/part%X_CO2(m)
......@@ -1436,6 +1453,7 @@ contains
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
if (cell%fld_r) cell%MP(cl,cy,cr) = cell%MP(cl,cy,cr)+part%mp(m)*weight
if (cell%fld_r) cell%CL(cl,cy,cr) = cell%CL(cl,cy,cr)+part%cl(m)*weight
if (cell%fld_x) cell%CO2(cl,cy,cr) = cell%CO2(cl,cy,cr)+part%CO2(m)*weight
if (cell%fld_x) cell%X_CO2(cl,cy,cr) = cell%X_CO2(cl,cy,cr)+part%X_CO2(m)*weight
......@@ -1467,7 +1485,8 @@ contains
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%FP(cl,cy,cr) = cell%FG(cl,cy,cr)*part%fg(m)**weight
if (cell%fld_r) cell%FG(cl,cy,cr) = cell%FG(cl,cy,cr)*part%fg(m)**weight
if (cell%fld_r) cell%MP(cl,cy,cr) = cell%MP(cl,cy,cr)*part%mp(m)**weight
if (cell%fld_r) cell%CL(cl,cy,cr) = cell%CL(cl,cy,cr)*part%cl(m)**weight
if (cell%fld_x) cell%CO2(cl,cy,cr) = cell%CO2(cl,cy,cr)*part%CO2(m)**weight
if (cell%fld_x) cell%X_CO2(cl,cy,cr) = cell%X_CO2(cl,cy,cr)*part%X_CO2(m)**weight
......@@ -1522,6 +1541,7 @@ contains
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)
if (cell%fld_r) cell%MP(i,j,k) = dummy_mp(i,j,k)
if (cell%fld_r) cell%CL(i,j,k) = dummy_cl(i,j,k)
if (cell%fld_x) cell%CO2(i,j,k) = dummy_CO2(i,j,k)
if (cell%fld_x) cell%X_CO2(i,j,k) = dummy_X_CO2(i,j,k)
......@@ -1579,6 +1599,7 @@ contains
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)
if (cell%fld_r) cell%MP(i,j,k) = vec_w(i,j,k)/cell%MP(i,j,k)
if (cell%fld_r) cell%CL(i,j,k) = vec_w(i,j,k)/cell%CL(i,j,k)
if (cell%fld_x) cell%CO2(i,j,k) = vec_w(i,j,k)/cell%CO2(i,j,k)
if (cell%fld_x) cell%X_CO2(i,j,k) = vec_w(i,j,k)/cell%X_CO2(i,j,k)
......@@ -1611,6 +1632,7 @@ contains
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)
if (cell%fld_r) cell%MP(i,j,k) = cell%MP(i,j,k)/vec_w(i,j,k)
if (cell%fld_r) cell%CL(i,j,k) = cell%CL(i,j,k)/vec_w(i,j,k)
if (cell%fld_x) cell%CO2(i,j,k) = cell%CO2(i,j,k)/vec_w(i,j,k)
if (cell%fld_x) cell%X_CO2(i,j,k) = cell%X_CO2(i,j,k)/vec_w(i,j,k)
......@@ -1643,6 +1665,7 @@ contains
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)
if (cell%fld_r) cell%MP(i,j,k) = cell%MP(i,j,k)**vec_w(i,j,k)
if (cell%fld_r) cell%CL(i,j,k) = cell%CL(i,j,k)**vec_w(i,j,k)
if (cell%fld_x) cell%CO2(i,j,k) = cell%CO2(i,j,k)**vec_w(i,j,k)
if (cell%fld_x) cell%X_CO2(i,j,k) = cell%X_CO2(i,j,k)**vec_w(i,j,k)
......@@ -1711,6 +1734,8 @@ contains
if (cell%fld_r) cell%FP(0:,:,nr+1) = cell%FP(0:,:,nr)
if (cell%fld_r) cell%FG(0:,:,0) = cell%FG(0:,:,1)
if (cell%fld_r) cell%FG(0:,:,nr+1) = cell%FG(0:,:,nr)
if (cell%fld_r) cell%MP(0:,:,0) = cell%MP(0:,:,1)
if (cell%fld_r) cell%MP(0:,:,nr+1) = cell%MP(0:,:,nr)
if (cell%fld_r) cell%CL(0:,:,0) = cell%CL(0:,:,1)
if (cell%fld_r) cell%CL(0:,:,nr+1) = cell%CL(0:,:,nr)
if (cell%fld_x) cell%CO2(0:,:,0) = cell%CO2(0:,:,1)
......@@ -1775,6 +1800,8 @@ contains
if (cell%fld_r) cell%FP(nl+1,:,:) = cell%FP(1,:,:)
if (cell%fld_r) cell%FG(0,:,:) = cell%FG(nl,:,:)
if (cell%fld_r) cell%FG(nl+1,:,:) = cell%FG(1,:,:)
if (cell%fld_r) cell%MP(0,:,:) = cell%MP(nl,:,:)
if (cell%fld_r) cell%MP(nl+1,:,:) = cell%MP(1,:,:)
if (cell%fld_r) cell%CL(0,:,:) = cell%CL(nl,:,:)
if (cell%fld_r) cell%CL(nl+1,:,:) = cell%CL(1,:,:)
if (cell%fld_x) cell%CO2(0,:,:) = cell%CO2(nl,:,:)
......@@ -1840,6 +1867,8 @@ contains
if (cell%fld_r) cell%FP(:,nl+1,:) = cell%FP(:,1,:)
if (cell%fld_r) cell%FG(:,0,:) = cell%FG(:,nl,:)
if (cell%fld_r) cell%FG(:,nl+1,:) = cell%FG(:,1,:)
if (cell%fld_r) cell%MP(:,0,:) = cell%MP(:,nl,:)
if (cell%fld_r) cell%MP(:,nl+1,:) = cell%MP(:,1,:)
if (cell%fld_r) cell%CL(:,0,:) = cell%CL(:,nl,:)
if (cell%fld_r) cell%CL(:,nl+1,:) = cell%CL(:,1,:)
if (cell%fld_x) cell%CO2(:,0,:) = cell%CO2(:,nl,:)
......@@ -1906,6 +1935,8 @@ contains
if (cell%fld_r) cell%FP(nl+1,:,0:) = cell%FP(nl,:,0:)
if (cell%fld_r) cell%FG(0,:,0:) = cell%FG(1,:,0:)
if (cell%fld_r) cell%FG(nl+1,:,0:) = cell%FG(nl,:,0:)
if (cell%fld_r) cell%MP(0,:,0:) = cell%MP(1,:,0:)
if (cell%fld_r) cell%MP(nl+1,:,0:) = cell%MP(nl,:,0:)
if (cell%fld_r) cell%CL(0,:,0:) = cell%CL(1,:,0:)
if (cell%fld_r) cell%CL(nl+1,:,0:) = cell%CL(nl,:,0:)
if (cell%fld_x) cell%CO2(0,:,0:) = cell%CO2(1,:,0:)
......@@ -1971,6 +2002,8 @@ contains
if (cell%fld_r) cell%FP(:,nl+1,:) = cell%FP(:,nl,:)
if (cell%fld_r) cell%FG(:,0,:) = cell%FG(:,1,:)
if (cell%fld_r) cell%FG(:,nl+1,:) = cell%FG(:,nl,:)
if (cell%fld_r) cell%MP(:,0,:) = cell%MP(:,1,:)
if (cell%fld_r) cell%MP(:,nl+1,:) = cell%MP(:,nl,:)
if (cell%fld_r) cell%CL(:,0,:) = cell%CL(:,1,:)
if (cell%fld_r) cell%CL(:,nl+1,:) = cell%CL(:,nl,:)
if (cell%fld_x) cell%CO2(:,0,:) = cell%CO2(:,1,:)
......@@ -2000,7 +2033,7 @@ contains
deallocate(vec_w,vec_w2,vec_w3,dummy_c,dummy_w,dummy_wr,dummy_fa,dummy_fc,dummy_g,dummy_r,dummy_o,dummy_cp,&
&dummy_a,dummy_rho,dummy_t,dummy_k,dummy_h,dummy_v,dummy_d,dummy_s,dummy_df,dummy_sp,dummy_e,dummy_fp,&
&dummy_fg,dummy_cl,dummy_CO2,dummy_X_CO2,dummy_X_H2O,dummy_X_U238,dummy_X_U235,dummy_X_Th232,dummy_X_K40)
&dummy_fg,dummy_mp,dummy_cl,dummy_CO2,dummy_X_CO2,dummy_X_H2O,dummy_X_U238,dummy_X_U235,dummy_X_Th232,dummy_X_K40)
deallocate(weight_vec)
......@@ -2639,6 +2672,7 @@ contains
if (cell%fld_e) part%e(m) = cell%entropy(part%cell_l(m),part%cell_y(m),part%cell_r(m))
if (cell%fld_r) part%fp(m) = cell%FP(part%cell_l(m),part%cell_y(m),part%cell_r(m))
if (cell%fld_r) part%fg(m) = cell%FG(part%cell_l(m),part%cell_y(m),part%cell_r(m))
if (cell%fld_r) part%mp(m) = cell%MP(part%cell_l(m),part%cell_y(m),part%cell_r(m))
if (cell%fld_r) part%cl(m) = cell%CL(part%cell_l(m),part%cell_y(m),part%cell_r(m))
if (cell%fld_x) part%CO2(m) = cell%CO2(part%cell_l(m),part%cell_y(m),part%cell_r(m))
if (cell%fld_x) part%X_CO2(m) = cell%X_CO2(part%cell_l(m),part%cell_y(m),part%cell_r(m))
......
......@@ -13,6 +13,7 @@ contains
! combine all functions and routines in one main subroutine
! ToDo: exchange depth rmax-r+Psurf by cell%ref_p
! ToDo: adapt initialize_rheol_vec to phase layers from here!!!
subroutine get_phases(part,cell,jmin2D,T,wr,pD,pF,pI,rmax,Psurf,compr,debug,time)
type(particle),intent(inout)::part
type(cell_properties),intent(inout)::cell
......@@ -30,7 +31,7 @@ contains
real(dp),allocatable::ClSl(:),counter(:),RhoSl(:),PSl(:),TSl(:),mass_tr(:)
!re-set field to zero, to be set below
part%r = 0.0_dp ! r = material value, is set below by phase transitions
part%mp = 0.0_dp ! mp = mineral phase value, is set below by phase transitions
ph = 0 ! ph is phase (i.e. first phase with material 3)
allocate(ClSl(pD%nr_ph),counter(pD%nr_ph),RhoSl(pD%nr_ph),PSl(pD%nr_ph),TSl(pD%nr_ph),mass_tr(pD%nr_ph))
......@@ -120,7 +121,7 @@ contains
! ((pD%m_ph(p_i).eq.6).and.((rmax-r+Psurf)*pF%rho*pF%g*pF%D.lt.1.0e+11_dp))))) then ! if pressure below 100 GPa, phase 7 (ppv) not stable
! found the right phase
part%r(m) = pD%m_ph(p_i)
part%mp(m) = pD%m_ph(p_i)
ph = p_i
! write(*,*) "Particle ",m," has phase ",ph," (",pD%m_ph(ph),"), Cl=",Cl_ph*pF%rho*pF%g*pF%D/pF%DeltaT*1.0e-6_dp,"MPa/K"&
......@@ -250,7 +251,7 @@ contains
enddo
if (ph.eq.0) then !last phase
part%r(m) = pD%m_ph(pD%nr_ph+1)
part%mp(m) = pD%m_ph(pD%nr_ph+1)
ph = pD%nr_ph+1
! write(*,*) "Particle ",m," has phase ",ph," (",pD%m_ph(ph),")"
......@@ -331,13 +332,13 @@ contains
part%a(m) = pD%alpha_ph(ph)
part%cp(m) = 1.0_dp
else
values1 = material_prop_ND(rmax-r+Psurf,T(i,j,k),nint(part%r(m)),error,pF,pI) ! needed, otherwise wrong phase values are set, not clear why though...
values1 = material_prop_ND(rmax-r+Psurf,T(i,j,k),nint(part%mp(m)),error,pF,pI) ! needed, otherwise wrong phase values are set, not clear why though...
part%a(m) = values1(3)
part%rho(m) = values1(2)
part%cp(m) = values1(4)
part%t(m) = T(i,j,k) ! local evaluation of rho, hence refT = T
!write(*,*) "Ol: rho=",part%rho(m)*pF%rho,", cp=",part%cp(m)*pF%Cp,&
! &", alpha=",part%a(m)*pF%alpha," for particle ",m," and phase ",part%r(m)
! &", alpha=",part%a(m)*pF%alpha," for particle ",m," and phase ",part%mp(m)
endif
! write(*,'("m ",i6," T=",F7.2," p=",ES10.3," rho=",F7.2," Tr=",F7.2," a=",ES10.3," Cp=",F7.2," k=",F7.2," h=",F7.2)') &
......@@ -363,7 +364,7 @@ contains
enddo
!write(*,*) "Phases range from ",minval(part%r)," to ",maxval(part%r)
!write(*,*) "Phases range from ",minval(part%mp)," to ",maxval(part%mp)
deallocate(ClSl,counter,RhoSl,PSl,TSl,mass_tr)
......@@ -616,7 +617,7 @@ contains
subroutine SetMatProp(part,jmin2D,mat,rc,T,Psurf,pD,pF,pI,comp)
type(particle),intent(inout)::part
integer,intent(in)::jmin2D
real(dp),intent(in)::mat(0:,jmin2D:,0:) ! mat = cell%r, this field normally convects with particle flow, but here always set back to transitions
real(dp),intent(in)::mat(0:,jmin2D:,0:) ! mat = cell%MP, this field normally convects with particle flow, but here always set back to transitions
real(dp),intent(in)::rc(0:),T(0:,jmin2D:,0:)
real(dp),intent(in)::Psurf
real(dp),intent(in)::comp
......@@ -629,7 +630,7 @@ contains
real(dp)::values(9)
real(dp)::r_ph,r_max
part%r = pD%m_ph(pD%nr_ph+1) ! set all material to lower-most shell, shells above are set to correct value below
part%mp= pD%m_ph(pD%nr_ph+1) ! set all material to lower-most shell, shells above are set to correct value below
part%k = pD%k_ph(pD%nr_ph+1)
part%h = pD%h_ph(pD%nr_ph+1) ! ToDo: problem with heat source depletion by melting: always overwritten again by initial phase value...
part%a = pD%alpha_ph(pD%nr_ph+1)
......@@ -645,7 +646,7 @@ contains
j = part%cell_y(m) ! equals 1 in 2D
k = part%cell_r(m)
values = material_prop_ND(r_max-rc(k)+Psurf,T(i,j,k),int(part%r(m)),error,pF,pI)
values = material_prop_ND(r_max-rc(k)+Psurf,T(i,j,k),int(part%mp(m)),error,pF,pI)
part%e(m) = values(9)
! add phase transition jumps from surface downwards
......@@ -653,7 +654,7 @@ contains
p_i = pD%nr_ph+1-p_inv ! from bottom phase to upper phase, to not over-write values
r_ph = pD%r_ph(p_i)-pD%Cl_ph(p_i)*(T(i,j,k)-pD%T_ph(p_i))
if (rc(k).gt.r_ph) then
part%r(m) = pD%m_ph(p_i)
part%mp(m)= pD%m_ph(p_i)
part%k(m) = pD%k_ph(p_i)
part%h(m) = pD%h_ph(p_i)
......@@ -663,18 +664,18 @@ contains
part%cp(m) = 1.0_dp
!part%t(m) = T(i,j,k) ! Di-dependent reference profile, no update
!part%e(m) = 1.0_dp
values = material_prop_ND(r_max-rc(k)+Psurf,T(i,j,k),int(part%r(m)),error,pF,pI)
values = material_prop_ND(r_max-rc(k)+Psurf,T(i,j,k),int(part%mp(m)),error,pF,pI)
part%e(m) = values(9)
endif
endif
enddo
if (int(part%r(m)).eq.0) write(*,*) i,j,k,int(part%r(m)),mat(i,j,k)
if (int(part%mp(m)).eq.0) write(*,*) i,j,k,int(part%mp(m)),mat(i,j,k)
if (comp.ge.3.0_dp) then
values = material_prop_ND(r_max-rc(k)+Psurf,T(i,j,k),int(part%r(m)),error,pF,pI)
values = material_prop_ND(r_max-rc(k)+Psurf,T(i,j,k),int(part%mp(m)),error,pF,pI)
part%a(m) = values(3)
part%rho(m) = values(2) !material_density(rc(k),0.0_dp,int(part%r(m)),pF) ! values(2); no temp influence here, this is what alpha-term is for in buocancy term in stoks/energy
part%rho(m) = values(2) !material_density(rc(k),0.0_dp,int(part%mp(m)),pF) ! values(2); no temp influence here, this is what alpha-term is for in buocancy term in stoks/energy
part%cp(m) = values(4)
part%t(m) = T(i,j,k) ! local evaluation of rho, hence refT = T
part%e(m) = values(9) ! entropy
......@@ -686,7 +687,7 @@ contains
subroutine SetPhaseVisc(j2D,eta,str_rate,material,ref_p,r,T,T0,pD,pF,pV,useW,useG,w_H2O,w_r,grain)
integer,intent(in)::j2D,useW,useG
real(dp),intent(inout)::eta(0:,j2D:,0:) ! eta%CV
real(dp),intent(in)::material(0:,j2D:,0:),str_rate(0:,j2D:,0:),ref_p(0:,j2D:,0:),w_H2O(useW:,min(useW+j2D,1):,useW:),& ! material = cell%r
real(dp),intent(in)::material(0:,j2D:,0:),str_rate(0:,j2D:,0:),ref_p(0:,j2D:,0:),w_H2O(useW:,min(useW+j2D,1):,useW:),& ! material = cell%MP
& w_r(useW:,min(useW+j2D,1):,useW:),grain(useG:,min(useG+j2D,1):,useG:)
real(dp),intent(in)::r(0:),T(0:,j2D:,0:),T0
type(param_D),intent(in)::pD
......@@ -712,7 +713,7 @@ contains
do p_i=1,pD%nr_ph+1
!write(*,*) i,j,k,int(pD%m_ph(p_i)),pD%m_ph(p_i),int(material(i,j,k)),material(i,j,k)
if (nint(material(i,j,k)).eq.int(pD%m_ph(p_i))) then ! use int values since cell%r values are average of particles -> int leads to round values
if (nint(material(i,j,k)).eq.int(pD%m_ph(p_i))) then ! use int values since cell%MP values are average of particles -> int leads to round values
! if prefactors are not used, set to 1:
if (pD%A_dif(p_i).eq.0.0_dp) then ; A_dif=1.0_dp ; else ; A_dif=pD%A_dif(p_i) ; endif
if (pD%A_dis(p_i).eq.0.0_dp) then ; A_dis=1.0_dp ; else ; A_dis=pD%A_dis(p_i) ; endif
......
......@@ -125,6 +125,7 @@ module variables
real(dp),allocatable::entropy(:,:,:)!e !< entropy of material (for phase transitions)
real(dp),allocatable::FP(:,:,:) !r !< factor phases
real(dp),allocatable::FG(:,:,:) !r !< phases derivative times Clapeyron slope
real(dp),allocatable::MP(:,:,:) !r !< mineral phases
real(dp),allocatable::CL(:,:,:) !r !< Clapeyron slope
real(dp),allocatable::CO2(:,:,:) !x !< CO2 concentration in solid (reduced with time due to melt depletion)
real(dp),allocatable::X_CO2(:,:,:) !x !< CO2 concentration in melt (depends on local conditions and oxygen fugacity)
......@@ -365,6 +366,7 @@ contains
if (cell%fld_e) then ; allocate( cell%entropy(0:n1+1,jmin:jmax,0:n3+1) ) ; else ; allocate( cell%entropy(1,1,1) ) ; endif
if (cell%fld_r) then ; allocate( cell%FP(0:n1+1,jmin:jmax,0:n3+1) ) ; else ; allocate( cell%FP(1,1,1) ) ; endif
if (cell%fld_r) then ; allocate( cell%FG(0:n1+1,jmin:jmax,0:n3+1) ) ; else ; allocate( cell%FG(1,1,1) ) ; endif
if (cell%fld_r) then ; allocate( cell%MP(0:n1+1,jmin:jmax,0:n3+1) ) ; else ; allocate( cell%MP(1,1,1) ) ; endif
if (cell%fld_r) then ; allocate( cell%CL(0:n1+1,jmin:jmax,0:n3+1) ) ; else ; allocate( cell%CL(1,1,1) ) ; endif
if (cell%fld_x) then ; allocate( cell%CO2(0:n1+1,jmin:jmax,0:n3+1) ) ; else ; allocate( cell%CO2(1,1,1) ) ; endif
if (cell%fld_x) then ; allocate( cell%X_CO2(0:n1+1,jmin:jmax,0:n3+1) ) ; else ; allocate( cell%X_CO2(1,1,1) ) ; endif
......@@ -380,7 +382,7 @@ contains
cell%wr = 0.0_dp ; cell%fa = 0.0_dp ; cell%fc = 0.0_dp ; cell%d = 0.0_dp
cell%ref_T = 0.0_dp ; cell%ref_r = 0.0_dp ; cell%ref_a = 0.0_dp ; cell%o = 0.0_dp ; cell%shear = 0.0_dp
cell%ref_g = 1.0_dp ; cell%ref_p = 0.0_dp ; cell%bulk = 0.0_dp ; cell%entropy = 0.0_dp ; cell%FP = 1.0_dp ; cell%FG = 0.0_dp
cell%CL = 0.0_dp ; cell%CO2 = 0.0_dp ; cell%X_CO2 = 0.0_dp ; cell%X_H2O = 0.0_dp
cell%MP = 0.0_dp ; cell%CL = 0.0_dp ; cell%CO2 = 0.0_dp ; cell%X_CO2 = 0.0_dp ; cell%X_H2O = 0.0_dp
cell%X_U238 = 0.0_dp ; cell%X_U235 = 0.0_dp ; cell%X_Th232 = 0.0_dp ; cell%X_K40 = 0.0_dp
if (mc%fld_d) then ; allocate( mc%solidus(0:n1+1,jmin:jmax,0:n3+1) ) ; else ; allocate( mc%solidus(1,1,1) ) ; endif
......@@ -462,7 +464,8 @@ contains
deallocate( field%v_old )
deallocate( field%w_old )
deallocate( cell%c, cell%w, cell%wr, cell%fa, cell%fc, cell%g, cell%r, cell%cp, cell%k, cell%h, cell%v, cell%d, cell%ref_T, &
& cell%ref_r, cell%ref_a, cell%ref_g, cell%ref_p, cell%o, cell%shear, cell%bulk, cell%entropy, cell%FP, cell%FG, cell%CL )
& cell%ref_r, cell%ref_a, cell%ref_g, cell%ref_p, cell%o, cell%shear, cell%bulk, cell%entropy, cell%FP, cell%FG, &
& cell%MP, cell%CL )