parameters.f90 22.6 KB
Newer Older
Lena Noack's avatar
Lena Noack committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
!>@brief Contains the type "parameters" , which contains an array of all parameters (read-in by input.f90)
!>@author Lena Noack
!>@date April 23, 2013
module parameters
  INCLUDE 'use_mpi.f90'
  use precision
  implicit none

  ! parameter structures pX, call of variable Z in structure pX via: pX%Z

  ! nondimensionalization factors
  type param_F
    real(dp)::D=1.0_dp
    real(dp)::time_yr=1.0_dp
    real(dp)::g=1.0_dp
    real(dp)::rho=1.0_dp
    real(dp)::kappa=1.0_dp
    real(dp)::Cp=1.0_dp
    real(dp)::alpha=1.0_dp
    real(dp)::k=1.0_dp
    real(dp)::grain=1.0_dp
    real(dp)::DeltaT=1.0_dp
    real(dp)::eta_ref=0.0_dp
    real(dp)::Ts=0.0_dp
    real(dp)::R=1.0_dp
  end type param_F

  type param_T
    real(dp)::ampl=0.1_dp !<amplitude of the spherical harmonics disturbation of the initial temperature field, if ampl=0 no disturbation
    integer::sph_l=1,sph_y=1 !<mode of spherical harmonics, if sph_l=0 random disturbation, sph_y only used for cartesian in 3D
    integer::linear=0 !< linear=0: constant initial temperature in mantle, linear=1: linear increase of temperature from surface to CMB
    real(dp)::isotherm=0.8_dp
    real(dp)::Tbottom=1.0_dp !<Temperature at bottom
    real(dp)::Ttop=0.0_dp  !< Temperature on top
    real(dp)::Tini=0.5_dp,s_bot=0.1_dp,s_top=0.1_dp !< values needed for initialization of temperature profile 3 (convective profile with TBLs)
    real(dp)::T0=1.0_dp !< non-dimensional surface temperature
    real(dp)::adiabat=0.0_dp
38
39
    real(dp)::Tmin=0.0_dp,Tmax=0.0_dp !< cut temperature if below/above these values
    real(dp)::DeltaTc=0.0_dp
Lena HPC's avatar
Lena HPC committed
40
    logical::CutSolidus=.false.
Lena Noack's avatar
Lena Noack committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
  end type param_T

  type param_M
    integer::use_melt=0,use_melt_old=0,use_lat_heat=0,update_solidus=0,add_crust=0,add_crust_incr=0 !< melt parameters, note: use_melt_old is not read-in, but defined by snapshot!
    integer::melt_extract=1
    real(dp)::Psurf=0.0_dp
    real(dp)::max_depl=1.0_dp
    real(dp)::entropy=0.0_dp
    real(dp)::lat_heat=0.0_dp
    real(dp)::melt_depth=1.0_dp
    real(dp)::melt_z_change=1.0_dp
    real(dp)::dehydr_melt=0.0_dp
    real(dp)::serp_instable_m=0.0_dp,serp_instable_T=0.0_dp
    real(dp)::serp_instable_P1=0.0_dp,serp_instable_P2=0.0_dp
    real(dp)::serp_instable_w=0.0_dp,serp_instable_w_T=0.0_dp,serp_instable_w_max=0.0_dp
    real(dp)::satur_w1=0.0_dp,satur_w2=0.0_dp,satur_exp=0.0_dp
    real(dp)::sol_hyd=0.0_dp,sol_hyd_exp=1.0_dp
    real(dp)::melt1s_0=0.0_dp,melt1s_1=0.0_dp,melt1s_2=0.0_dp,melt1s_3=0.0_dp
    real(dp)::melt1l_0=0.0_dp,melt1l_1=0.0_dp,melt1l_2=0.0_dp,melt1l_3=0.0_dp
    real(dp)::melt2s_0=0.0_dp,melt2s_1=0.0_dp,melt2s_2=0.0_dp,melt2s_3=0.0_dp,melt2s_4=0.0_dp
    real(dp)::melt2l_0=0.0_dp,melt2l_1=0.0_dp,melt2l_2=0.0_dp,melt2l_3=0.0_dp,melt2l_4=0.0_dp
    real(dp)::LatH=0.0_dp
    real(dp)::Dxi=0.0_dp
    real(dp)::Lambda=1.0_dp
    real(dp)::FCLambda=1.0_dp,CCLambda=1.0_dp
    real(dp)::DT_plume=0.0_dp,Per_plume=0.0_dp
Lena Noack's avatar
Lena Noack committed
67
68
    real(dp)::Ra_melt=0.0_dp,Ra_depl=0.0_dp,rho_melt=1.0_dp,rho_depl=1.0_dp,Rt=0.0_dp
    real(dp)::rho_f=1.0_dp,rate_freezing=10.0_dp,n_perm=3 !for Schmeling-Melt benchmark
Lena Noack's avatar
Lena Noack committed
69
70
    integer::No_del=0
    real(dp)::add_crust_depth=0.0_dp
71
    real(dp)::Chi_CO2=0.0_dp,Chi_extr=0.0_dp ! outgassing parameters
Lena Noack's avatar
Lena Noack committed
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
  end type param_M

  type param_V
    real(dp)::T_ref=0.0_dp  !< reference temperature for viscosity
    real(dp)::z_ref=0.0_dp  !< reference depth for viscosity
    real(dp)::gamma=0.0_dp  !< T-coefficient in FK approximation of the Arrhenius law
    real(dp)::gammaP=0.0_dp  !< p-coefficient in FK approximation
    real(dp)::R=8.314472_dp
    real(dp)::E=0.0_dp,V=0.0_dp !< for linear Arrhenius viscosity
    real(dp)::E_dis=0.0_dp,V_dis=0.0_dp,n_fac=0.0_dp,n_fac_wet=0.0_dp,A_dif=0.0_dp,A_dis=0.0_dp,p_fac=0.0_dp,p_fac_wet=0.0_dp !< for Non-newtonian / mixed versions
    real(dp)::E_wet=0.0_dp,V_wet=0.0_dp,E_dis_wet=0.0_dp,V_dis_wet=0.0_dp,A_dif_wet=0.0_dp,A_dis_wet=0.0_dp !< for wet mantle/lithosphere
    integer::StartNN=0 !< for Non-newtonian / mixed versions
    real(dp)::water_r_dif=0.0_dp,water_r_dis=0.0_dp
    integer::AdVisc=0 !< adapt viscosity to "real" boundaries: use avergae viscosity at upper and lower boundary
    real(dp)::eta_min=0.0_dp,eta_max=0.0_dp
    integer::rheol_vec=0
88
    real(dp)::DeltaVisc=1.0_dp,DeltaVisc_LM=1.0_dp
Lena Noack's avatar
Lena Noack committed
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
  end type param_V

  ! crust (only pre-defined, development of new crust in melt parameter structure?)
  type param_C
    integer::cont_nr=0
    real(dp)::cont_w=0.0_dp,cont_d=0.0_dp
    real(dp)::cont_CC_rho=0.0_dp,cont_CC_rhoRef=0.0_dp,cont_CC_TRef=0.0_dp,cont_CC_eta=0.0_dp
    real(dp)::cont_CC_cp=0.0_dp,cont_CC_k=0.0_dp,cont_CC_h=0.0_dp,cont_CC_alpha=1.0_dp
    real(dp)::cont_CC_w=0.0_dp,cont_CC_wr=0.0_dp,cont_CC_fr_angl=0.0_dp,cont_CC_fr_angl_m=0.0_dp,cont_CC_fr_coh=0.0_dp
    real(dp)::cont_FC_rho=0.0_dp,cont_FC_rhoRef=0.0_dp,cont_FC_TRef=0.0_dp,cont_FC_eta=0.0_dp
    real(dp)::cont_FC_cp=0.0_dp,cont_FC_k=0.0_dp,cont_FC_h=0.0_dp,cont_FC_alpha=1.0_dp
    real(dp)::cont_FC_w=0.0_dp,cont_FC_wr=0.0_dp,cont_FC_fr_angl=0.0_dp,cont_FC_fr_angl_m=0.0_dp,cont_FC_fr_coh=0.0_dp
    real(dp)::cont_OC_rho=0.0_dp,cont_OC_rhoRef=0.0_dp,cont_OC_TRef=0.0_dp,cont_OC_eta=0.0_dp
    real(dp)::cont_OC_cp=0.0_dp,cont_OC_k=0.0_dp,cont_OC_h=0.0_dp,cont_OC_alpha=1.0_dp
    real(dp)::cont_OC_w=0.0_dp,cont_OC_wr=0.0_dp,cont_OC_fr_angl=0.0_dp,cont_OC_fr_angl_m=0.0_dp,cont_OC_fr_coh=0.0_dp
    real(dp)::LHB_time=0.0_dp,crust_ini=0.0_dp
    real(dp)::weak_angle_l=0.0_dp,weak_angle_r=0.0_dp,weak_w=0.0_dp,weak_d=0.0_dp,weak_l=0.0_dp,weak_r=0.0_dp,ridge_spread=0.0_dp
  end type param_C

  type param_I
    real(dp)::alpha=1.0_dp,alpha_c=1.0_dp  !< termal expansion coefficient
    real(dp)::Cp=1.0_dp,Cp_c=1.0_dp,Cp_cr=1.0_dp,Cp_r=1.0_dp !< heat capacity, is connected to density via: Cp = k/(rho*kappa)
    real(dp)::gravity=1.0_dp,gravity_c=1.0_dp  !< acceleration of gravity
    real(dp)::kappa=1.0_dp !< thermal DIFFUSION constant
    real(dp)::Ra=0.0_dp  !< Rayleigh number
114
    real(dp)::RaB=1.0_dp !< ALA/compress Rayleigh number Ra*B
Lena Noack's avatar
Lena Noack committed
115
116
117
118
119
120
121
122
    real(dp)::eta_ref=1.0_dp !< pre-factor for viscosity, different from pF%eta_ref, which stores input value
    real(dp)::rho=1.0_dp,rho_c=1.0_dp,rho_cr=1.0_dp,rho_r=1.0_dp,rho_mc=1.0_dp !< density
    real(dp)::rhogD=1.0_dp !< rho * g * D
    real(dp)::p_surf=0.0_dp !< surface pressure (default 0), TODO: implement p_surf in code, but where to add? seems like cancels out everywhere...
    real(dp)::D=1.0_dp
    real(dp)::Di=0.0_dp,DiT=0.0_dp !< actually Di (for rho) and DiT (for T) should both be the standard dissipation number, but DiT can be set to zero if only density variations with depth are considered...
    real(dp)::C_ref=1.0_dp !< reference composition (should correspond to reference density)
    real(dp)::B=0.0_dp !< Buoyancy ratio for chemical diffusion
123
    real(dp)::drho=220.0_dp !< density difference between depleted and undepleted material, where depleted is scaled to 100% depletion
Lena Noack's avatar
Lena Noack committed
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
    real(dp)::Le=0.0_dp,LeV=0.0_dp !< Lewis number for chemical diffusion and viscosity diffsuion
    real(dp)::comp_init=0.0_dp !< initialization of composition: +-1 = linear, +-2 = rho value; positive: add cylinder, negative: add random perturbation
    real(dp)::comp_ampl=0.0_dp !< amplitude of initial perturbation
    real(dp)::lambda=0.0_dp !< ln(2)/half life
    real(dp)::H0=0.0_dp !< heat generation rate
    real(dp)::k=1.0_dp,k_max=0.0_dp,k_l=1.0_dp,k_cr=1.0_dp,k_r=1.0_dp !< thermal conductivity, if k_max ne 0 then use linear conductivity profile (only for benchmark purposes...)
    real(dp)::mantle_w_H2O=0.0_dp,mantle_w_r=1.0_dp,mantle_fr_angl=0.0_dp,mantle_fr_angl_m=0.0_dp,mantle_fr_coh=0.0_dp
    real(dp)::Tm_ref=0.0_dp,rho_ref=1.0_dp
    real(dp)::grain_size=1.0_dp
    real(dp)::H0U238=0.0_dp,H0U235=0.0_dp,H0Th232=0.0_dp,H0K40=0.0_dp
    real(dp)::t12U238=0.0_dp,t12U235=0.0_dp,t12Th232=0.0_dp,t12K40=0.0_dp
    real(dp)::X0U238=0.0_dp,X0U235=0.0_dp,X0Th232=0.0_dp,X0K40=0.0_dp
    real(dp)::C_U=0.0_dp,C_Th=0.0_dp,C_K=0.0_dp
    real(dp)::Rc=1.0_dp,Rp=2.0_dp
    integer::CCool=0
    real(dp)::compress=0.0_dp ! set to 1.0 (TALA) or 2.0 (ALA) if compressible
    real(dp)::a1=1.0_dp,a2=0.0_dp,a3=0.0_dp,a4=0.0_dp ! compressibility parameters fro depth-dependent alpha-function
    real(dp)::Gr=1.0_dp ! Grueneisen Parameter
    integer::TinclTref=1 ! TinclTref = 1: simulation T is T' + Tref, TinclTref = 0: simulation T is only T', -1: like 0 but Work for T instead of T', 2: like 1 but with Work with T' instead of T
143
    real(dp)::Psurf=0.0_dp,mu=0.6_dp ! only for stress induced by overlying pressure
Lena Noack's avatar
Lena Noack committed
144
    integer::read_profs=0 ! read-in of profiles for density, alpha etc.
145
    real(dp)::X_FeM=0.1_dp ! 10 wt-% iron / Mg# 0.9
146
    real(dp)::X_M0=0.0_dp ! depression of liquidus melt temp by other minerals to solidus after Stixrude et al., 2014 (e.g. 0.21 for MgO and SiO2)
Lena HPC's avatar
Lena HPC committed
147
    real(dp)::X_H2OM=0.0_dp ! 0 wt-% water in mantle
148
    integer::prof_nm=0 ! size of mantle part in profs
Lena Noack's avatar
Lena Noack committed
149
    integer::ec_model=4 ! model to be used for upper mantle for electrical conductivity in model1D.f90
Lena Noack's avatar
Lena Noack committed
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
  end type param_I

  type param_H
    real(dp)::eccentricity=0.0_dp
    real(dp)::shear_modulus=0.0_dp
    real(dp)::shear_modulus_c=0.0_dp  
    real(dp)::rho_av=0.0_dp
    real(dp)::diss_fact=0.0_dp
    real(dp)::orb_period=0.0_dp    
    integer::TidHeat=0
    integer::TidHeatAdapt=0
    real(dp)::k2=0.0_dp
    real(dp)::h2=0.0_dp
    real(dp)::l2=0.0_dp
    real(dp)::eta_ref_c=1.0_dp
    real(dp)::G_fac=0.0_dp
    real(dp)::mu_fac=0.0_dp
    real(dp)::TidH_fac=0.0_dp
    real(dp)::tidalH_global=0.0_dp
    complex(dp)::FA_=0.0_dp    
    complex(dp)::FB_=0.0_dp  
    complex(dp)::FC_=0.0_dp
    complex(dp)::Hmu_=0.0_dp    
173
    real(dp),allocatable::vec_magn_H(:),vals_H(:,:)
174
175
    real(dp)::use_magn_H=0
    integer::stop_magm_oc=0
Lena HPC's avatar
Lena HPC committed
176
    real(dp)::Q_eddy
Lena Noack's avatar
Lena Noack committed
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
  end type param_H

  type param_G
    integer::geom=0
    integer::dim=2
    integer::nl=0   !< discretization size for lateral dimension (would be nx for cartesian or nphi for polar)
    integer::ny=0   !< discretization size for perpendicular dimension (exists only in cartesian domain)
    integer::nr=0   !< discretization size for radial dimension (would be ny for cartesian or nr for polar)
    real(dp)::lmin=0.0_dp !< minimum abscissae of the domain, not used in polar version 
    real(dp)::lmax=0.0_dp !< maximum abscissae of the domain, not used in polar version
    real(dp)::ymin=0.0_dp !< minimum depth of the 3D box domain (not used in spherical geom. or 2D box version)
    real(dp)::ymax=0.0_dp !< maximum depth of the 3D box domain (not used in spherical geom. or 2D box version)
    real(dp)::rmin=0.0_dp !< minimum ordinate of the domain
    real(dp)::rmax=0.0_dp !< maximum ordinate of the domain
    integer::noslip_t=0,noslip_b=0,noslip_w=0
    integer::inner_bound=0,outer_bound=0 !< isolation of the boundary when set to 1
    integer::open_b_bound=0 !< open boundary at the bottom: mass flus is possible if value is set to 1
    integer::periodic=0
    real(dp)::phi_factor=0.0_dp !< 1=full annulus, 0.5=half annulus, ...
    real(dp)::theta_factor=0.0_dp
    real(dp)::save_lmin=0.0_dp,save_lmax=0.0_dp,save_ymin=0.0_dp,save_ymax=0.0_dp,save_rmin=0.0_dp,save_rmax=0.0_dp
    integer::nl_s=0,ny_s=0,nr_s=0 !ToDo: To be deleted, also in stokes and other solvers...
  end type param_G

  type param_W
    real(dp)::weak_x0=0.0_dp,weak_x1=0.0_dp,weak_angle=0.0_dp,weak_d=0.0_dp,weak_w=0.0_dp,weak_visc=1.0_dp,weak_rho=1.0_dp
    real(dp)::weak_cp=1.0_dp,weak_k=1.0_dp,weak_h=0.0_dp,weak_T=0.0_dp,weak_fr_coh=0.0_dp,weak_fr_angl=0.0_dp
    real(dp)::weak_fr_angl_m=0.0_dp,weak_w_H2O=0.0_dp,weak_w_r=1.0_dp !< predefined weak zone in lithosphere

    real(dp)::OP_BOC=0.0_dp,OP_BOC_rho=1.0_dp,OP_BOC_eta=1.0_dp,OP_BOC_cp=1.0_dp,OP_BOC_k=1.0_dp
    real(dp)::OP_SHB=0.0_dp,OP_SHB_rho=1.0_dp,OP_SHB_eta=1.0_dp,OP_SHB_cp=1.0_dp,OP_SHB_k=1.0_dp
    real(dp)::OP_ThL=0.0_dp,OP_ThL_rho=1.0_dp,OP_ThL_eta=1.0_dp,OP_ThL_age=0.0_dp,OP_ThL_cp=1.0_dp,OP_ThL_k=1.0_dp !< overriding plate

    real(dp)::SP_BOC=0.0_dp,SP_BOC_rho=1.0_dp,SP_BOC_eta=1.0_dp,SP_BOC_cp=1.0_dp,SP_BOC_k=1.0_dp
    real(dp)::SP_SHB=0.0_dp,SP_SHB_rho=1.0_dp,SP_SHB_eta=1.0_dp,SP_SHB_cp=1.0_dp,SP_SHB_k=1.0_dp
    real(dp)::SP_ThL=0.0_dp,SP_ThL_rho=1.0_dp,SP_ThL_eta=1.0_dp,SP_ThL_age=0.0_dp,SP_ThL_cp=1.0_dp,SP_ThL_k=1.0_dp !< subducting plate

    real(dp)::OP_BOC_T=0.0_dp,OP_SHB_T=0.0_dp,OP_ThL_T=0.0_dp,SP_BOC_T=0.0_dp,SP_SHB_T=0.0_dp,SP_ThL_T=0.0_dp
    real(dp)::OP_BOC_h=0.0_dp,OP_SHB_h=0.0_dp,OP_ThL_h=0.0_dp,SP_BOC_h=0.0_dp,SP_SHB_h=0.0_dp,SP_ThL_h=0.0_dp

    real(dp)::OP_BOC_fr_coh=0.0_dp,OP_BOC_fr_angl=0.0_dp,OP_BOC_fr_angl_m=0.0_dp,OP_BOC_w_H2O=0.0_dp,OP_BOC_w_r=1.0_dp
    real(dp)::OP_SHB_fr_coh=0.0_dp,OP_SHB_fr_angl=0.0_dp,OP_SHB_fr_angl_m=0.0_dp,OP_SHB_w_H2O=0.0_dp,OP_SHB_w_r=1.0_dp
    real(dp)::OP_ThL_fr_coh=0.0_dp,OP_ThL_fr_angl=0.0_dp,OP_ThL_fr_angl_m=0.0_dp,OP_ThL_w_H2O=0.0_dp,OP_ThL_w_r=1.0_dp
    real(dp)::SP_BOC_fr_coh=0.0_dp,SP_BOC_fr_angl=0.0_dp,SP_BOC_fr_angl_m=0.0_dp,SP_BOC_w_H2O=0.0_dp,SP_BOC_w_r=1.0_dp
    real(dp)::SP_SHB_fr_coh=0.0_dp,SP_SHB_fr_angl=0.0_dp,SP_SHB_fr_angl_m=0.0_dp,SP_SHB_w_H2O=0.0_dp,SP_SHB_w_r=1.0_dp
    real(dp)::SP_ThL_fr_coh=0.0_dp,SP_ThL_fr_angl=0.0_dp,SP_ThL_fr_angl_m=0.0_dp,SP_ThL_w_H2O=0.0_dp,SP_ThL_w_r=1.0_dp

    real(dp)::weak_alpha=1.0_dp,OP_BOC_alpha=1.0_dp,OP_SHB_alpha=1.0_dp,OP_ThL_alpha=1.0_dp,SP_BOC_alpha=1.0_dp
    real(dp)::SP_SHB_alpha=1.0_dp,SP_ThL_alpha=1.0_dp

    real(dp)::OS_SuL=0.0_dp,OS_VelRed=1.0_dp,OS_SuL_rho=1.0_dp,OS_SuL_T=0.0_dp,OS_SuL_eta=1.0_dp,OS_SuL_cp=1.0_dp,OS_SuL_k=1.0_dp
    real(dp)::OS_SuL_alpha=1.0_dp,OS_SuL_h=0.0_dp,OS_SuL_fr_coh=0.0_dp,OS_SuL_fr_angl=0.0_dp,OS_SuL_fr_angl_m=0.0_dp
    real(dp)::OS_SuL_w_H2O=0.0_dp,OS_SuL_w_r=1.0_dp,compr_w=0.0_dp,rhogDw=1.0_dp

    real(dp)::vel_z_in=0.0_dp,vel_z_out=0.0_dp,vel_r=0.0_dp !< values for imposed plate velocity

    integer::T_init_plate=0,UseRefDensity=0
    integer::BenchSubdFreeS=0
    real(dp)::eta_ref_m=1.0_dp,eta_ref_s=1.0_dp,n_m=1.0_dp,n_s=1.0_dp
  end type param_W

  type param_N
    integer::test=0 !< Used only for test issues, default is zero for normal simulations
    real(dp)::p_penalty=1.0e-7_dp !< Penalty term for pressure calculation in the mass equation (adding the penalty term simulates an almost-incompressible convection), 0 neglects penalty term
    real(dp)::beta=1.0_dp !< beta-factor for flux-limiter
    integer::average=3,p_average=2,p_average2=2,pm_average=2 !< averaging scheme to be used for viscosity/particles: 1=harmonic, 2=arithmetic, 3=geometric
    integer::debug=0 !< set debug between 1 and 5 to get detailled information during the simulation
    integer::noBouss=0 !< if set to one then no Boussinesq is used in momentum solver
    integer::e_solver=3 !< choose between the different energy solvers: 0 = first-oder explicit, 1 = second-order explicit, 2 = second-order implicit
    integer::c_solver=3 !< choose between the double-diffusive solvers: 0 = first-oder explicit, 1 = second-order explicit, 2 = second-order implicit
    integer::m_solver=1 !< 0 (default) in one matrix, 1 decoupld in u-v-w and p matricex, 2 decoupled in u,v,w and p matrices
    integer::t_method=1,s_method=1 !< 0 = FDM, 1 = FVM1 (with FDM for second derivatives), 2 = FVM2 (using twice FVM for second-order derivatives), t=thermal, s=stokes
    integer::itmax_bcg=1000,rout_solver=0,itol_bcg=1 ! max iterations to be used in linear bcg solver, solver to be used (0:Pardiso, 1:linBCG, -1:linBG in thermal only)
    real(dp)::tol_bcg=1.0e-9_dp
    integer::mom_DuDt=0 !< if 1 then include Du/Dt term in momentum equation, default is zero
    integer::no_m_solver=0 !< if one then velocity is not solved, otherwise standard stokes approach
    integer::CorrRot=0
    integer::iter_pard=2 !< number of iterations for pardiso solver
    integer::pr_time=0
    integer::vel_extr=0 !< velocity boundary extrapolation method if MPI is used
    integer::PreStokes=0 !< ToDo: Not used anymore???
    integer::fieldA=0,fieldB=0,fieldC=0,fieldD=0,fieldF=0,fieldG=0,fieldH=0,fieldK=0,fieldL=0,fieldM=0,fieldO=0,fieldP=0
    integer::fieldR=0,fieldS=0,fieldT=0,fieldV=0,fieldW=0,fields_All=0 
    logical::NoMMSolveA=.false.,NoMMSolveBC=.false.,NoStDRho=.false.,NoVD=.false.,NoEnEBrho=.false.
    real(dp)::relax=0.8_dp,relax_v=1.0_dp,relax_tol=1.2_dp ! relaxation factor for pressure correction, should be inbetween (but excluding) 0 and 1
    real(dp)::CN_theta=-1.0_dp ! if set to 0.5 then Crank-Nicolson method in energy solver, if between 0 and 1 then generalized Euler method, if -1 (default) then 3rd-order time with upwind (e_solver parameter)
    integer::UseUpwExpl=0,UseUpwImpl=0
    integer::aver_v_str=0 ! use laterally-averaged profiles for temperature and strain-rate in viscosity calculation
265
    logical::binary=.true.
266
    integer::ReduceMaxTS=1
Lena Noack's avatar
Lena Noack committed
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
  end type param_N

  type param_A
    real(dp)::ad_min_l=0.0_dp,ad_min_y=0.0_dp,ad_min_r=0.0_dp
    real(dp)::ad_max_l=0.0_dp,ad_max_y=0.0_dp,ad_max_r=0.0_dp
    real(dp)::res_ratio_l=0.0_dp,res_ratio_y=0.0_dp,res_ratio_r=0.0_dp
    integer::ad_factor=0
  end type param_A

  type param_O
    integer::max_outer=1,min_outer=1 !< number of outer iterations (velocity and temperature are solved alternately in each outer iteration until solution converges or max_outer is reached)
    integer::max_inner=1,max_inner_init=1,max_inner_TS=0
    integer::nmesures=0 !< number of iterations to be made while measuring physical data
    integer::nbiteration=0
    integer::iter_output=0,no_init_out=0,FirstOutput=1
    real(dp)::t_output=0.0_dp
    integer::use_snap=0,last_snap=0,start_snap=0
    real(dp)::convSS=0.0_dp !< convergence criteria asked to stop time iterations, ( i.e assumption of steaty state)
    real(dp)::convV=1.0e-4_dp !< convergence criteria for velocity (inner iteration convergence)
    real(dp)::convT=1.0e-6_dp !< thershold value for convergence of energy solver (outer iteration convergence)
    real(dp)::CourantCoeff=1.0_dp !< coefficient C in courant criterion to determine timestep dt.
    real(dp)::dt_min=1.0e-10_dp,dt_max=1.0e-2_dp,dt_ini=1.0e-6_dp !< maximal, minimal and initial dt, dt_ini is used as fixed dt if CourantCoeff=0 (but should be between dt_min and dt_max)
    real(dp)::tmax=0.0_dp !<maximal nondimensional simulation time
    integer::output_lid=0,OutputCompr=0
  end type param_O

  ! strain, strain rate, subduction (plate tectonics)
  type param_S
    real(dp)::iniStrRate=1.0_dp,minStrRate=0.0_dp,str_pre_fac=1.0_dp
    real(dp)::strain_min=0.0_dp,strain_max=0.0_dp !< nondimens. strain weakening values (for plate tectonics with plasticity=2)
    real(dp)::yield_stress=0.0_dp !< Parameters for plate tectonics simulations
    real(dp)::yield_stress_z=0.0_dp !< Parameters for plate tectonics simulations
Lena Noack's avatar
Lena Noack committed
299
    real(dp)::ys_water_density=0.0_dp !< subtracted from density multiplied to yield_stress_z, set to 1000 to include water effect, 0 otherwise
Lena Noack's avatar
Lena Noack committed
300
301
302
303
304
305
306
307
308
    integer::plasticity=0 !< choose between different plasticity modules: 0-Bingham, 1-nonsmooth, 2-angle-dependent
    real(dp)::PY_eta_min=0.0_dp
    real(dp)::friction=0.0_dp !< nondimens. factor for friction for plasticity=2
  end type param_S

  ! density variations by phase transitions
  type param_D
    integer::nr_ph=0 ! for now max 4: either Ice VI-VII-X or olivine(alpha-beta-gamma)-pv-ppv
    real(dp),allocatable::T_ph(:),r_ph(:),Cl_ph(:),P_ph(:),eta_ph(:) ! phase temperature, radius (dim. pressure), Clapeyron slope, Phase parameter and pre-defined viscosity jump
Lena Noack's avatar
Lena Noack committed
309
    real(dp),allocatable::alpha_ph(:),k_ph(:),h_ph(:) ! layer properties
Lena Noack's avatar
Lena Noack committed
310
311
312
313
314
    integer,allocatable::m_ph(:) ! layer properties
    real(dp)::dr_ph=0.05_dp ! thickness of phase transition (needed for numerical stability: jump not possible, use smoother transition over dr)
    real(dp),allocatable::A_dif(:),A_dis(:),A_dif_w(:),A_dis_w(:),n_fac(:),p_fac(:),wr_dif(:),wr_dis(:)
    real(dp),allocatable::E_dif(:),V_dif(:),E_dis(:),V_dis(:),E_dif_w(:),V_dif_w(:),E_dis_w(:),V_dis_w(:)
    real(dp)::Layer=0.0_dp,LayerDensity=1.0_dp,LayerVisc=1.0_dp,LayerMaxV=1.0_dp,LayerMix=0.5_dp,LayerOutput=1.0_dp,LayerAmpl=0.0_dp
Lena Noack's avatar
Lena Noack committed
315
    integer::PhT_param=0
316
    integer::set_tracers=1
Lena Noack's avatar
Lena Noack committed
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
  end type param_D

  ! particles
  type param_P
    integer::use_part=0 !< use particles
    integer::part_n=0,part_n_cell=0,rk=4
    real(dp)::part_depth=1.0_dp !< depth of layer where initial particles should be set (1 in general)
    real(dp)::part_pref=1.0_dp
    integer::restart_p=1
    integer::bnd_l=1,bnd_y=1,bnd_r=1 !bnd_l=2 used for subduction zone (but then save_right_bnd needed, otherwise errors can occur)
    integer::max_part=0,min_part=0
    integer::moveRef=0 ! Normally, reference values like density, alpha, Cp, k and temp should not convect, but for subduction benchmark they do!
  end type param_P

  ! 1D evolution
  type param_E
    real(dp)::Tm=0.0_dp,Tcr=0.0_dp,Tc0=0.0_dp,DTc=0.0_dp
    real(dp)::Dl=0.0_dp,Dcr0=0.0_dp,Dreg=0.0_dp,Delta_u=0.0_dp,Delta_c=0.0_dp,Dlcr=0.0_dp,Delta_c_max,minDcr=0.0_dp
    real(dp)::eps_c=1.1_dp,eps_m=1.0_dp,eps_i=1.0_dp,eps_w=1.0_dp
    real(dp)::beta=0.33333333_dp
    real(dp)::fc=0.0_dp
    real(dp)::Ra_crit=450.0_dp
    real(dp)::Theta=2.9_dp
340
341
    real(dp)::gc=0.0_dp,Gr=0.0_dp !6.67384e-11_dp
    real(dp)::u0=2.0e-12_dp,DTsol=0.0_dp,Dref=0.0_dp
Lena Noack's avatar
Lena Noack committed
342
343
344
345
    real(dp)::Mpcr0U238=0.0_dp,Mpcr0U235=0.0_dp,Mpcr0Th232=0.0_dp,Mpcr0K40=0.0_dp
    integer::output_dim=1,H_from_surf=0,CondStop=0
    integer::core_freezing=-1,CheckCore=0
    real(dp)::XiS0=0.0_dp,Tm0=0.0_dp,Tm1=0.0_dp,Tm2=0.0_dp,Ta1=0.0_dp,Ta2=0.0_dp,alpha_S=0.0_dp
Lena Noack's avatar
Lena Noack committed
346
    real(dp)::P0=0.0_dp,Pc=0.0_dp,Psurf=0.0_dp !1.0e-4_dp
347
    real(dp)::L=0.0_dp,t_rc=0.0_dp,dV_V=0.03_dp,M_E=0.0_dp,X_Fe=0.0_dp,X_H2O=0.0_dp,X_Mg=0.0_dp
Lena Noack's avatar
Lena Noack committed
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
    real(dp)::Doc=0.0_dp,T_OMB=0.0_dp,DTmo=0.0_dp,DTmi=0.0_dp,DTbo=0.0_dp,Delta_o=0.0_dp,Delta_i=0.0_dp,Dli=0.0_dp ! ocean-layer
    integer::n_oc=100 ! number of grid points in ocean
    integer::TestIce=0,DimValues=0
    real(dp)::E_ice=0.0_dp,V_ice=0.0_dp,eta_ref_ice=0.0_dp,t_start_count=0.0_dp
    integer::IntStruct=0,PT=0
  end type param_E

  ! paraLlelism
  type param_L
    integer::nprocs_l=0,nprocs_y=0,nprocs_r=0 !< read-in in input.txt, number of processors per direction, e.g. (nprocs, 1, 1) for pure lateral split
    integer::shift_nl=0,shift_ny=0,shift_nr=0 !< shift of nl/ny/nr indices per domain needed for correct snapshots
    integer::nb_rank_r=-1,nb_rank_l=-1,nb_rank_f=-1,nb_rank_b=-1,nb_rank_t=-1,nb_rank_u=-1
    integer::save_nl=0,save_ny=0,save_nr=0
    integer,allocatable::nl_vec(:),ny_vec(:),nr_vec(:)
  end type param_L

364
365
366
  ! X mantle concentrations and outgassing
  type param_X
    integer::use_atmos=0
Boissinot Pierre's avatar
Boissinot Pierre committed
367
    real(dp)::CO2_ini=0.0_dp,fO2=0.0_dp,fwm=36.594_dp,dH2O=0.01_dp,dCO2=0.1_dp
368
    !real(dp),allocatable::CO2(:,:,:),X_CO2(:,:,:)
369
    real(dp)::dU238=0.0_dp,dU235=0.0_dp,dTh232=0.0_dp,dK40=0.0_dp
370
371
372
    ! below for gas speciation model from Gianluigi Ortenzi
    integer::gas_spec=0
    real(dp)::alpha_melt=3.0e-5_dp,Cp_melt=1793.0_dp,density_melt=3000.0_dp ! for melt adiabat
373
374
  end type param_X

Lena Noack's avatar
Lena Noack committed
375
376
377
378
379
380
381
382
383
! contains

  ! nothing





end module parameters