-
Notifications
You must be signed in to change notification settings - Fork 0
/
read_params.f90
368 lines (336 loc) · 10 KB
/
read_params.f90
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
38
39
40
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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
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
265
266
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
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
! Reads problem parameters from input file
subroutine read_params(inputfile)
include 'precision.inc'
include 'params.inc'
include 'arrays.inc'
character*200 inputfile
iu =4
open( iu, file=inputfile )
! MESH
call AdvanceToNextInputLine( 4 )
read(4,*) nx, nz
open(11,file='nxnz.0')
write(11,*) nx, nz
close(11)
if((nx.gt.mnx) .or. (nz.gt.mnz)) then
write(*,*) '# of elements exceed maximum. Increase mnx and mnz in "arrays.inc".'
stop 1
endif
nq = nx*nz
nx = nx+1
nz = nz+1
nmarkers = 9*(nz-1)*(nx-1)
call AdvanceToNextInputLine( 4 )
read(4,*) x0, z0
call AdvanceToNextInputLine( 4 )
read(4,*) rxbo, rzbo
call AdvanceToNextInputLine( 4 )
read(4,*) ircoord, coordfile
if (ircoord .gt. 0) then
! ignore next two lines
call AdvanceToNextInputLine( 4 )
read(4,*) nzonx
if (nzonx .ne. 0) stop 'ircoord is set, but nzonx is not 0!'
call AdvanceToNextInputLine( 4 )
read(4,*) nzony
if (nzony .ne. 0) stop 'ircoord is set, but nzony is not 0!'
go to 177
endif
call AdvanceToNextInputLine( 4 )
read(4,*) nzonx
if (nzonx .eq.0) then
nzonx = 1
nelz_x(1) = nx - 1
sizez_x(1) = 1.
go to 166
endif
do i = 1, nzonx
call AdvanceToNextInputLine( 4 )
read(4,*) nelz_x(i), sizez_x(i)
end do
166 continue
call AdvanceToNextInputLine( 4 )
read(4,*) nzony
if (nzony .eq.0) then
nzony = 1
nelz_y(1) = nz - 1
sizez_y(1) = 1.
go to 177
endif
do i = 1,nzony
sizez_y(i) = 1.
call AdvanceToNextInputLine( 4 )
read(4,*) nelz_y(i), sizez_y(i)
end do
177 continue
call AdvanceToNextInputLine( 4 )
read(4,*) iint_marker, iint_tracer
call AdvanceToNextInputLine( 4 )
read(4,*) nzone_marker
call AdvanceToNextInputLine( 4 )
do 253 i = 1, nzone_marker
read(iu,*) imx1(i),imy1(i),imx2(i),imy2(i)
253 continue
call AdvanceToNextInputLine( 4 )
read(4,*) nzone_tracer, dt_outtracer
dt_outtracer = dt_outtracer * 1000 * sec_year
call AdvanceToNextInputLine( 4 )
do 233 i = 1, nzone_tracer
read(iu,*) itx1(i),ity1(i),itx2(i),ity2(i)
233 continue
! MECHANICAL CONDITIONS
call AdvanceToNextInputLine( 4 )
read(4,*) ynstressbc,ydrsides
call AdvanceToNextInputLine( 4 )
read(4,*) nofbc
call AdvanceToNextInputLine( 4 )
do 21 i = 1,nofbc
read(iu,*) nofside(i),nbc1(i),nbc2(i),nbc(i), &
bca(i),bcb(i),bcc(i), &
bcd(i),bce(i),bcf(i),bcg(i),bch(i),bci(i)
21 continue
call AdvanceToNextInputLine( 4 )
! hydrostatic pressure applied at the bottom
call AdvanceToNextInputLine( 4 )
read(4,*) nyhydro,pisos,iphsub,drosub,damp_vis
! gravity
call AdvanceToNextInputLine( 4 )
read(4,*) g
! THERMAL CONDITIONS
call AdvanceToNextInputLine( 4 )
read(4,*) i_prestress
call AdvanceToNextInputLine( 4 )
read(4,*) itherm
call AdvanceToNextInputLine( 4 )
read(4,*) istress_therm
call AdvanceToNextInputLine( 4 ) ! thermal stresses
read (4,*) ishearh ! shear heating
call AdvanceToNextInputLine( 4 )
read (4,*) t_top
call AdvanceToNextInputLine( 4 )
read (4,*) t_bot
call AdvanceToNextInputLine( 4 )
read (4,*) hs, hr
! boundary conditions at the bottom (1-T,2-Flux)
call AdvanceToNextInputLine( 4 )
read (4,*) itemp_bc, bot_bc
if( itemp_bc.eq.2 ) bot_bc = bot_bc/1000 ! convert in W/m3
! temperature pertrubation (rectangular)
call AdvanceToNextInputLine( 4 )
read (4,*) temp_per, ix1t, ix2t, iy1t, iy2t
! Predefined distributions
call AdvanceToNextInputLine( 4 )
read(4,*) irtemp
if ( irtemp .gt. 0 ) then
call AdvanceToNextInputLine( 4 )
read(4,*) tempfile
else
call AdvanceToNextInputLine( 4 )
read(4,*)
endif
! time scale
call AdvanceToNextInputLine( 4 )
read (4,*) time_scale
! temp structure
call AdvanceToNextInputLine( 4 )
read (4,*) iynts,tbos
call AdvanceToNextInputLine( 4 )
read (4,*) iax1,iay1,ibx1,iby1,icx1,icy1,idx1,idy1
call AdvanceToNextInputLine( 4 )
read (4,*) nzone_age
call AdvanceToNextInputLine( 4 )
do i = 1, nzone_age
read (4,*) age_1(i),hc1(i),hc2(i),hc3(i),hc4(i),iph_col1(i),iph_col2(i), &
iph_col3(i),iph_col4(i),iph_col5(i),ixtb1(i),ixtb2(i)
enddo
! check smooth nzone_age
print *, 'B ixtb1 = ', ixtb1, 'B ixtb2 = ', ixtb2
do i = 1, nzone_age-1
iph_col_trans(i) = 0
if (ixtb2(i) == -1) then
if (ixtb1(i+1) /= -1) then
print *, 'Error: ixtb1 is not -1 at', i+1, 'column!'
endif
if (iph_col1(i) /= iph_col1(i+1)) then
print *, 'Error: iph_col1 at', i, i+1, 'columns are not equal!'
endif
if (iph_col2(i) /= iph_col2(i+1)) then
print *, 'Error: iph_col2 at', i, i+1, 'columns are not equal!'
endif
if (iph_col3(i) /= iph_col3(i+1)) then
print *, 'Error: iph_col3 at', i, i+1, 'columns are not equal!'
endif
if (iph_col4(i) /= iph_col4(i+1)) then
print *, 'Error: iph_col4 at', i, i+1, 'columns are not equal!'
endif
if (iph_col5(i) /= iph_col5(i+1)) then
print *, 'Error: iph_col5 at', i, i+1, 'columns are not equal!'
endif
! flag indicates smotth transition
iph_col_trans(i) = 1
! remove '-1' for
ixtb2(i) = ixtb2(i+1)
ixtb1(i+1) = ixtb1(i)
endif
enddo
print *, 'A ixtb1 = ', ixtb1, 'A ixtb2 = ', ixtb2
! RHEOLOGY
call AdvanceToNextInputLine( 4 )
read(4,*) nphase
do i = 1, nphase
call AdvanceToNextInputLine( 4 )
read(4,*) irheol(i),visc(i),den(i),alfa(i),beta(i),pln(i),acoef(i),eactiv(i),rl(i),rm(i), &
plstrain1(i),plstrain2(i),fric1(i),fric2(i),cohesion1(i),cohesion2(i), &
dilat1(i),dilat2(i), &
conduct(i),cp(i),ts(i),tl(i),tk(i),fk(i)
print *, i
end do
! Flag to take initial phase distribution from a file
call AdvanceToNextInputLine( 4 )
read(4,*) irphase
!write(*,*) irphase
if ( irphase .gt. 0 ) then
call AdvanceToNextInputLine( 4 )
read(4,*) phasefile
else
call AdvanceToNextInputLine( 4 )
read(4,*)
endif
! main phase
call AdvanceToNextInputLine( 4 )
read(4,*) mphase
! number of horizontal layers
call AdvanceToNextInputLine( 4 )
read(4,*) nphasl
! layers
do i = 1, nphasl
call AdvanceToNextInputLine( 4 )
read(4,*) ltop(i), lbottom(i), lphase(i)
end do
! inclusions
call AdvanceToNextInputLine( 4 )
read(4,*) inhom
if( inhom .gt. 50 ) then
call SysMsg('Read_params: Increase arrays for inhomogenities')
stop 26
endif
do i = 1, inhom
call AdvanceToNextInputLine( 4 )
read(4,*) ix1(i), ix2(i), iy1(i), iy2(i), inphase(i), igeom(i), xinitaps(i)
end do
! Tension cut-off
call AdvanceToNextInputLine( 4 )
read(4,*) ten_off
!linear healing parameter
call AdvanceToNextInputLine( 4 )
read(4,*) tau_heal
! viscosity limits
call AdvanceToNextInputLine( 4 )
read(4,*) v_min, v_max, ivis_shape,efoldc
call AdvanceToNextInputLine( 4 )
read(4,*)igeotherm,g_x0,g_y0c, g_amplitude,g_width
call AdvanceToNextInputLine( 4 )
!read(4,*) ny_inject, nelem_inject, nelem_inject1, nelem_inject2, rate_inject1, rate_inject2
!read(4,*) ny_inject, nelem_inject, rate_inject_brittle, rate_inject_ductile
read(4,*) ny_inject, nelem_inject, rate_inject_brittle, rate_inject_ductile
call AdvanceToNextInputLine( 4 )
read(4,*)rate_inject_ductile_e, rate_inject_ductile_s
call AdvanceToNextInputLine( 4 )
read(4,*) iinj1, iinj2, jinj1, jinj2, xlatheat, ratfac, fnu
call AdvanceToNextInputLine( 4 ) ! add to the last line of the input file and param.inc
!it=1 time-independent it=2 time-dependent; y = fa*sine(2pi/fb*time)+fc
!ip=1 turn on the hydromineral lame change; ip=2 turn off
!fsr full spreading rate
read(4,*) it, fa, fb, fc, fsr, ip
! REMESHING
call AdvanceToNextInputLine( 4 )
read(4,*) ny_rem, mode_rem, ntest_rem, angle_rem
if ( mode_rem.ne.1 .and. mode_rem.ne.11 .and. mode_rem.ne.3 ) then
call SysMsg('Illegal remeshing mode! Allowable - 1, 3 or 11')
stop
endif
! dx_rem - remeshing criteria for mode_rem=11
call AdvanceToNextInputLine( 4 )
read(4,*) dx_rem
! diffusion of topography
call AdvanceToNextInputLine( 4 )
read(4,*) topo_kappa, fac_kappa
! PROCESS CONTROL
! inertial mass scaling
call AdvanceToNextInputLine( 4 )
read(4,*) idt_scale
call AdvanceToNextInputLine( 4 )
read(4,*) dt_scale, strain_inert
call AdvanceToNextInputLine( 4 )
read(4,*) i_rey,xReyn
call AdvanceToNextInputLine( 4 )
read(4,*) ifreq_rmasses
call AdvanceToNextInputLine( 4 )
read(4,*) ifreq_imasses
call AdvanceToNextInputLine( 4 )
read(4,*) ifreq_visc
call AdvanceToNextInputLine( 4 )
read(4,*) ifreq_avgsr
! acceleration parameters
call AdvanceToNextInputLine( 4 )
read(4,*) amul, ratl, ratu
call AdvanceToNextInputLine( 4 )
read(4,*) frac, fracm
call AdvanceToNextInputLine( 4 )
read(4,*) n_boff_cutoff
call AdvanceToNextInputLine( 4 )
read(4,*) movegrid,ndim
call AdvanceToNextInputLine( 4 )
read(4,*) demf, mix_strain, mix_stress
! OUTPUT PARAMETERS
call AdvanceToNextInputLine( 4 )
read(4,*) time_max
time_max = time_max * 1000 * sec_year
call AdvanceToNextInputLine( 4 )
read (4,*) dtout_screen
dtout_screen = dtout_screen * 1000 * sec_year
call AdvanceToNextInputLine( 4 )
read(4,*) dtout_file
dtout_file = dtout_file * 1000 * sec_year
call AdvanceToNextInputLine( 4 )
!read(4,*) nprofil
!do 135 i=1,nprofil
!read(4,*) hv_out(i:i), iprof_out(i)
!135 continue
!call AdvanceToNextInputLine( 4 )
read(4,*) io_vel,io_srII,io_eII,io_aps,io_sII,io_sxx,io_szz,io_sxz,io_pres, &
io_temp,io_melt,io_visc,io_phas,io_mark,io_src,io_diss,io_forc,io_hfl,io_topo
call AdvanceToNextInputLine( 4 )
read(4,*) lastout
call AdvanceToNextInputLine( 4 )
read(4,*) dtsave_file
dtsave_file = dtsave_file * 1000 * sec_year
call AdvanceToNextInputLine( 4 )
read(4,*) lastsave
close (iu)
! ADDITIONAL PARTICULAR INPUT
!call ReadMoreParams()
return
end
!========================================================
subroutine AdvanceToNextInputLine( iu )
character*1 buf
10 read(iu, '(A1)') buf
if( buf(1:1).eq.';' ) then
goto 10
else
backspace( iu )
return
endif
print *, 'AdvanceToNextInputLine: EOF reached!'
stop
print *, 'AdvanceToNextInputLine: Error reading file!'
stop
return
end
subroutine ReadMoreParams()
call ReadIntrusions() ! - see user_ab.f90
call ReadHydro() ! - see user_luc.f90
!call ReadHeatinject() ! - see user_Lu.f90 (6/28/18)!
return
end