melt.f90 63.5 KB
Newer Older
Lena Noack's avatar
Lena Noack committed
1
2
3
4
5
6
7
8
!>@author Lena Noack
!>@date July 2013
module melt
use precision
use mesher, only:mesh_cp
use variables
use parameters
use particles, only:particle
9
use heating, only:radiogenic_nuclei
Lena Noack's avatar
Lena Noack committed
10
11
12
13
14
15
16
17
18
19
20
21
22
implicit none
contains

  ! solidus reduction by water
  function solidus_hydrated(sol_hyd,sol_hyd_exp,w)
!    type(param),intent(in)::p
    real(dp),intent(in)::sol_hyd,sol_hyd_exp,w
    real(dp)::solidus_hydrated

    solidus_hydrated = sol_hyd*w**sol_hyd_exp 

  end function solidus_hydrated

23
  subroutine init_melt_curves(mesh,mc,pF,pM,di,cell,X_FeM) !ToDo maybe trace Fe via particles and deplete mantle in iron, for now: use just average mantle value
Lena Noack's avatar
Lena Noack committed
24
25
26
27
28
29
      type(melt_properties),intent(inout)::mc
      type(mesh_cp),intent(in)::mesh
      type(param_F),intent(in)::pF
      type(param_M),intent(in)::pM
      integer,intent(in)::di
      type(cell_properties),intent(in)::cell
30
      real(dp),intent(in)::X_FeM
Lena Noack's avatar
Lena Noack committed
31
32

      integer::i,j,k,nl,ny,nr,jmin,jmax
33
      real(dp)::z,z_Fe_cut,DeltaTs,Pnond
Lena Noack's avatar
Lena Noack committed
34
35
36
37
38
39
40
41
42
43
44
45
46

      nl = mesh%nl
      ny = mesh%ny
      nr = mesh%nr

      if (di.eq.2) then
        jmin = 1
        jmax = 1
      else
        jmin = 0
        jmax = ny+1
      endif

47
48
49
50
51
52
53
54
      ! add iron-related solidus change below, two different formulations for pressures above/below 12 GPa
      Pnond = 1.0e+9/(pF%rho*pF%g*pF%D)
      z_Fe_cut = 12.0_dp*Pnond

      do k=0,nr+1
        ! get non-dimensional pressure, i.e. depth
        if (LBOUND(cell%ref_p,3).eq.0) then
          if (cell%ref_p(1,1,1).ne.cell%ref_p(1,1,nr)) then ! actual pressure value, not zero or dummy value
Lena Noack's avatar
Lena Noack committed
55
56
            z = cell%ref_p(1,1,k)
          else
57
            z = mesh%rmax - mesh%rc(k) + pM%Psurf ! Psurf is nondimensionalized to yield non-dimensional depth
Lena Noack's avatar
Lena Noack committed
58
          endif
59
60
61
62
63
64
65
66
67
68
        else
          z = mesh%rmax - mesh%rc(k) + pM%Psurf
        endif

        ! calculate DeltaT
        if (z.lt.z_Fe_cut) then ! iron content between 0 and 1
          DeltaTs = (0.1_dp-X_FeM)*(102.3_dp + 64.1_dp*z/Pnond - 3.62_dp*(z/Pnond)**2)/pF%DeltaT
        else
          DeltaTs = (0.1_dp-X_FeM)*360.0_dp/pF%DeltaT
        endif
69
        !write(*,*) DeltaTs
Lena Noack's avatar
Lena Noack committed
70

71
72
73
        ! calculate solidus
        if (pM%melt_z_change.eq.0) then ! one curve for the whole mantle
          mc%solidus_ref(:,:,k) = pM%melt1s_0 + pM%melt1s_1*z + pM%melt1s_2*z**2 + pM%melt1s_3*z**3 + DeltaTs
Lena Noack's avatar
Lena Noack committed
74
75
76
77
78
          mc%liquidus(:,:,k) = pM%melt1l_0 + pM%melt1l_1*z + pM%melt1l_2*z**2 + pM%melt1l_3*z**3

          ! derivative in radial direction ( ds_dr = ds_dz * dz_dr = - ds_dz ):
          mc%ds_dz(:,:,k) = -(pM%melt1s_1 + 2.0_dp*pM%melt1s_2*z + 3.0_dp*pM%melt1s_3*z**2)
          mc%dl_dz(:,:,k) = -(pM%melt1l_1 + 2.0_dp*pM%melt1l_2*z + 3.0_dp*pM%melt1l_3*z**2)
79
80
        else        
          ! one curve for upper mantle, one for lower mantle
Lena Noack's avatar
Lena Noack committed
81
          if (z.lt.pM%melt_z_change) then ! upper mantle
82
            mc%solidus_ref(:,:,k) = pM%melt1s_0 + pM%melt1s_1*z + pM%melt1s_2*z**2 + pM%melt1s_3*z**3 + DeltaTs
Lena Noack's avatar
Lena Noack committed
83
84
85
86
87
            mc%liquidus(:,:,k) = pM%melt1l_0 + pM%melt1l_1*z + pM%melt1l_2*z**2 + pM%melt1l_3*z**3

            mc%ds_dz(:,:,k) = -(pM%melt1s_1 + 2.0_dp*pM%melt1s_2*z + 3.0_dp*pM%melt1s_3*z**2)
            mc%dl_dz(:,:,k) = -(pM%melt1l_1 + 2.0_dp*pM%melt1l_2*z + 3.0_dp*pM%melt1l_3*z**2)
          else ! lower mantle
88
            mc%solidus_ref(:,:,k) = pM%melt2s_0 + pM%melt2s_1*z + pM%melt2s_2*z**2 + pM%melt2s_3*z**3 + pM%melt2s_4*z**4 + DeltaTs
Lena Noack's avatar
Lena Noack committed
89
90
91
92
93
94
            mc%liquidus(:,:,k) = pM%melt2l_0 + pM%melt2l_1*z + pM%melt2l_2*z**2 + pM%melt2l_3*z**3 + pM%melt2l_4*z**4

            mc%ds_dz(:,:,k) = -(pM%melt2s_1 + 2.0_dp*pM%melt2s_2*z + 3.0_dp*pM%melt2s_3*z**2 + 4.0_dp*pM%melt2s_4*z**3)
            mc%dl_dz(:,:,k) = -(pM%melt2l_1 + 2.0_dp*pM%melt2l_2*z + 3.0_dp*pM%melt2l_3*z**2 + 4.0_dp*pM%melt2l_4*z**3)
          endif

95
          !write(*,*) k,z*pF%D,z*pF%D*pF%rho*pF%g*1.0e-9_dp,mc%solidus_ref(1,1,k)*pF%DeltaT+pF%Ts
96
97
        endif
      enddo
Lena Noack's avatar
Lena Noack committed
98

99
100
101
      mc%solidus = 0.0_dp ! increase of solidus due to depletion; together with fixed solidus_ref and solidus_hydrated(w) this yields the final solidus curve
      mc%DeltaF = 0.0_dp
      mc%serp = 1.0_dp
Lena Noack's avatar
Lena Noack committed
102
103
104

  end subroutine init_melt_curves

Lena Noack's avatar
Lena Noack committed
105

106
  ! initial temperature profile cut down to solidus if necessary
Lena Noack's avatar
Lena Noack committed
107
108
109
110
  subroutine cut_T_solidus(field,mc)
    type(variables_unknowns),intent(inout)::field
    type(melt_properties),intent(in)::mc

111
    where(field%T.gt.mc%solidus_ref) field%T = mc%solidus_ref ! ToDo: add hydration effect
Lena Noack's avatar
Lena Noack committed
112
113
114
  end subroutine


115
  subroutine melt_comp_p(field,cell,mesh,mc,part,pF,pM)  
Lena Noack's avatar
Lena Noack committed
116
      type(variables_unknowns),intent(in)::field
Lena Noack's avatar
Lena Noack committed
117
      type(cell_properties),intent(in)::cell
Lena Noack's avatar
Lena Noack committed
118
119
120
      type(mesh_cp),intent(in)::mesh
      type(melt_properties),intent(in)::mc
      type(particle),intent(inout)::part
121
      type(param_F),intent(in)::pF
Lena Noack's avatar
Lena Noack committed
122
123
124
      type(param_M),intent(in)::pM

      integer::i,j,k,m
Lena Noack's avatar
Lena Noack committed
125
      real(dp)::z,solidus,liquidus
Lena Noack's avatar
Lena Noack committed
126
127
128
129
130
131
132
133
134
135
136
137
138

      ! change depletion, TODO later on: add change in heat sources etc. (maybe in extra function like dehydration?)
      do m=1,part%n
        i = part%cell_l(m)
        j = part%cell_y(m) ! equals 1 in 2D
        k = part%cell_r(m)
        z = mesh%rmax - mesh%rc(k) + pM%Psurf
        solidus = mc%solidus_ref(i,j,k)+part%s(m)-solidus_hydrated(pM%sol_hyd,pM%sol_hyd_exp,part%w(m))   ! solidus is the initial reference solidus + solidus increase (part%s) - hydration
        if (solidus.lt.0.0_dp) then
          print *," +++++ SERIOUS PROBLEM WITH T/SOLIDUS +++++ "
          print *,solidus,mc%solidus_ref(i,j,k),part%s(m),-solidus_hydrated(pM%sol_hyd,pM%sol_hyd_exp,part%w(m))
        endif

139
140
141
142
143
144
145
146
!        write(*,*) m,k,z*pF%D*pF%rho*pF%g*1.0e-9_dp,mc%solidus_ref(i,j,k)*pF%DeltaT+pF%Ts,field%T(i,j,k)*pF%DeltaT+pF%Ts,&
!                   &mc%solidus_ref(i,j,k+1)*pF%DeltaT+pF%Ts,&
!                    &(pM%melt1s_0 + pM%melt1s_1*z + pM%melt1s_2*z**2 + pM%melt1s_3*z**3)*pF%DeltaT+pF%Ts
!        write(*,*) m,z*pF%D*pF%rho*pF%g*1.0e-9_dp,mc%solidus_ref(i,j,k)*pF%DeltaT+pF%Ts,field%T(i,j,k)*pF%DeltaT+pF%Ts,&
!                    &(solidus-field%T(i,j,k))*pF%DeltaT,mc%solidus_ref(i,j,k+1)*pF%DeltaT+pF%Ts,&
!                    &(pM%melt1s_0 + pM%melt1s_1*z + pM%melt1s_2*z**2 + pM%melt1s_3*z**3)*pF%DeltaT+pF%Ts,&
!                    &(pM%melt2s_0 + pM%melt2s_1*z + pM%melt2s_2*z**2 + pM%melt2s_3*z**3 + pM%melt2s_4*z**4)*pF%DeltaT+pF%Ts,&
!                    &mesh%rmax,mesh%rb(mesh%nr+1),mesh%dr_vec(mesh%nr),z
Lena Noack's avatar
Lena Noack committed
147
148
149
        if (pM%Rt.eq.0.0_dp) then !simple calculation of depletion and melt fraction
!          if ((z.lt.pM%melt_depth).and.(field%T(i,j,k).gt.solidus)) then
          if ((z.lt.pM%melt_depth).and.(field%T(i,j,k).gt.mc%solidus_ref(i,j,k))) then !negativ melt fractions allowed if melt is kept in system and freezes out
150
151
            if ((pM%melt_extract.eq.1).and.(field%T(i,j,k).gt.solidus)) then
              part%df(m) = (field%T(i,j,k)-solidus)/(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k))
152
153
154
              !write(*,*) "Melting occurred",i,k,field%T(i,j,k),solidus,mc%liquidus(i,j,k),part%df(m),part%d(m)
              !write(*,*) "Solidus=",solidus,"=",mc%solidus_ref(i,j,k),"+",part%s(m),"+",&
              !         &-solidus_hydrated(pM%sol_hyd,pM%sol_hyd_exp,part%w(m))
155
156
157
158
159
            else if (pM%melt_extract.ne.1) then
              part%df(m) = (field%T(i,j,k)-solidus)/(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k))
            else
              part%df(m) = 0.0_dp
            endif
Lena Noack's avatar
Lena Noack committed
160
161
162
163
164
          else
!            write(*,*) i,j,k,solidus,field%T(i,j,k),mc%liquidus(i,j,k),mc%solidus_ref(i,j,k),z,pM%melt_depth
            part%df(m) = 0.0_dp
          endif

Lena Noack's avatar
Lena Noack committed
165
166
167
          if ((i.eq.1).and.(j.eq.1).and.(k.eq.mesh%nr+1)) &
            & write(*,*) "Melt p 2:",part%df(m),part%w(m)

Lena Noack's avatar
Lena Noack committed
168
169
170
171
172
173
174
!          if ((part%df(m).gt.(pM%max_depl-part%d(m))).and.(pM%max_depl.ne.0.0_dp)) then ! more melt than defined in pM%max_depl is not allowed
!            part%df(m) = pM%max_depl-part%d(m)
!          endif
          if (pM%max_depl.ne.0.0_dp) part%df(m) = min( part%df(m),pM%max_depl-part%d(m) )
 
          if (pM%melt_extract.eq.1) part%df(m) = max(part%df(m),0.0_dp) ! no melt kept
          
Lena Noack's avatar
Lena Noack committed
175
          ! ToDo below has to be checked, does not really make sense, since depletion increases with melt???
Lena Noack's avatar
Lena Noack committed
176
177
178
179
180
181
182
183
184
185
          ! change in the whole code (especially energy solver): d increases from 0 until max_depl is reached...
          ! depletion value counts from hydrated rock on! Only if all water is gone, then normal melting from solidus (including increase due to depletion) starts 
          ! - therefore depletion has to include hydration effects
          !part%d(m) = part%d(m) - part%df(m)
          !if (part%d(m).lt.0.0_dp) then
          !  part%d(m) = 0.0_dp
          !endif

          part%d(m) = max(part%d(m) + part%df(m),0.0_dp) ! shoudln't go under 0.0 anyway, but just in case
        else
Lena Noack's avatar
Lena Noack committed
186
187
188
189
190
191
192
193
194
195
196
          ! two-phase flow benchmark
          ! assume here that liquidus temperature increases as well as solidus temperature!!!
          solidus = mc%solidus_ref(i,j,k) + (cell%d(i,j,k)-mc%porosity(i,j,k))*(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k)) ! no hydration considered here
          liquidus = mc%liquidus(i,j,k) + (cell%d(i,j,k)-mc%porosity(i,j,k))*(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k))
          if ((field%T(i,j,k).gt.solidus).and.(field%T(i,j,k).lt.liquidus)) then
            part%df(m) = (field%T(i,j,k)-solidus)/(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k)) ! part%d is set later, after particles moved
          else if (field%T(i,j,k).ge.liquidus) then
            part%df(m) = 1.0_dp
          else
            part%df(m) = 0.0_dp
          endif
Lena Noack's avatar
Lena Noack committed
197
        endif
198

199
200
201
202
203
204
!        if (part%df(m).gt.0.0_dp) then
!        if (mc%solidus_ref(i,j,k) + part%d(m)*(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k)).gt.field%T(i,j,k)) then
!          write(*,*) m,part%d(m),part%df(m),solidus,field%T(i,j,k),&
!             & mc%solidus_ref(i,j,k) + part%d(m)*(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k))
!        endif
!        endif
Lena Noack's avatar
Lena Noack committed
205
206
      enddo

207
!      if (sum(part%df).gt.0.0_dp) write(*,*) "Melting occured, average DeltaF=",sum(part%df)/real(part%n)
208

Lena Noack's avatar
Lena Noack committed
209
210
  end subroutine melt_comp_p

Boissinot Pierre's avatar
Boissinot Pierre committed
211

Lena Noack's avatar
Lena Noack committed
212
 ! not needed anymore, since now part of atmosphere.f90
Boissinot Pierre's avatar
Boissinot Pierre committed
213
 subroutine dehydration(part,mesh,cell,vol_H2O,total_H2O,di,dehydr_melt)
Lena Noack's avatar
Lena Noack committed
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
      type(particle),intent(inout)::part
      type(mesh_cp),intent(in)::mesh
      type(cell_properties),intent(inout)::cell
      real(dp),intent(inout)::vol_H2O,total_H2O,dehydr_melt
      integer,intent(in)::di

      integer::i,j,k,m,jmin,jmax

      if (di.eq.2) then
        jmin = 1
        jmax = 1
      else
        jmin = 0
        jmax = mesh%ny+1
      endif

      do m=1,part%n
Lena Noack's avatar
Lena Noack committed
231

Lena Noack's avatar
Lena Noack committed
232
233
234
235
!        if (part%d(m).lt.dehydr_melt) then ! immediate dehydration when melt occurs if dehydr_melt = depl_max; no dehydration at all if dehydr_melt = 0
        if ((dehydr_melt.gt.0.0_dp).and.(part%d(m).gt.dehydr_melt)) then ! dehydr_melt is melt in percent where all water is dehydrated, e.g. 1-5%
          part%w(m) = 0.0_dp ! Complete dehydration
        endif
Lena Noack's avatar
Lena Noack committed
236

Lena Noack's avatar
Lena Noack committed
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
      enddo

      vol_H2O = 0.0_dp
      do i=0,mesh%nl+1
        do j=jmin,jmax
          do k=0,mesh%nr+1
            if ((dehydr_melt.gt.0.0_dp).and.(cell%d(i,j,k).gt.dehydr_melt)) then
!            if (cell%d(i,j,k).lt.dehydr_melt) then
              vol_H2O = vol_H2O + cell%w(i,j,k)*mesh%dVi(i,j,k)
              cell%w(i,j,k) = 0.0_dp
            endif
          enddo
        enddo
      enddo
      total_H2O = total_H2O + vol_H2O

  end subroutine dehydration


Boissinot Pierre's avatar
Boissinot Pierre committed
256
257
 

Lena Noack's avatar
Lena Noack committed
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
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
  ! at a specific depth/temperature, water is released since subducted serpentine becomes instable
  ! shallow dewatering is considered by limiting the amount of water that can be subducted (everything else is assumed to be released at very shallow depths and not traced here)
  ! this max water limit depends on thickness of sediment layer, porosity, etc. and is applied by defining accurate water inclusion values of the lithosphere/crust in the input file 
  subroutine release_water(field,mesh,mc,cell,part,pM,rhogD,di)
      type(variables_unknowns),intent(in)::field
      type(mesh_cp),intent(in)::mesh
      type(melt_properties),intent(inout)::mc
      type(cell_properties),intent(inout)::cell
      type(particle),intent(inout)::part
      type(param_M),intent(in)::pM
      real(dp),intent(in)::rhogD
      integer,intent(in)::di

      integer::i,j,k,m,jmin,jmax,cl,cy,cr,k_up,instable
      real(dp)::free_w,satur_w,old_w,keep_H2O

      if (di.eq.2) then
        jmin = 1
        jmax = 1
      else
        jmin = 0
        jmax = mesh%ny+1
      endif

      do i=0,mesh%nl+1
       do j=jmin,jmax
        do k=0,mesh%nr+1
          instable = 0
          if (cell%w(i,j,k).gt.0.0_dp) then ! water that may be released
           if (mc%serp(i,j,k).ge.0.5_dp) then ! no serpentine instability has occured yet in main rock amount (>50%) of cell
            
            if (pM%serp_instable_m.ne.0.0_dp) then
              ! fixed depth at which serpentine gets instable
              if ((mesh%rmax-mesh%rc(k)+pM%Psurf).gt.pM%serp_instable_m) instable = 1
            elseif (pM%serp_instable_T.ne.0.0_dp) then
              ! fixed temperature
              if (field%T(i,j,k).gt.pM%serp_instable_T) instable = 1
            elseif ((mesh%rmax-mesh%rc(k)+pM%Psurf).gt.(pM%serp_instable_P1-pM%serp_instable_P2*field%T(i,j,k))) then 
              ! instability function (follows Irifune et al., 1998)
                instable = 1
            endif

            if (instable.eq.1) then
              ! how much water will be kept and how much released?
              keep_H2O = pM%serp_instable_w - pM%serp_instable_w_T*field%T(i,j,k)  ! in %
              keep_H2O = min(pM%serp_instable_w_max,max(0.0_dp,keep_H2O)) ! truncate between 0 and max%
              keep_H2O = keep_H2O/100.0_dp ! 0=all water is lost, 1=all water is kept

              ! determine amount of free water
              free_w = cell%w(i,j,k)*mesh%dVi(i,j,k)*(1.0_dp-keep_H2O)

              ! set water to zero in field and particles
              cell%w(i,j,k) = keep_H2O !0.0_dp
              mc%serp(i,j,k) = 0.0_dp
              do m=1,part%n
                cl = part%cell_l(m)
                cy = part%cell_y(m)
                cr = part%cell_r(m)
                if ((cl.eq.i).and.(cy.eq.j).and.(cr.eq.k)) then
                  part%w(m) = keep_H2O !0.0_dp
                  part%sp(m) = 0.0_dp
                endif
              enddo

              print *,"Free-water for (",i,",",j,",",k,"): ",free_w,", keep H2O=",&
                 &keep_H2O*100.0_dp,"%",", T=",field%T(i,j,k),", z=",mesh%rmax-mesh%rc(k)+pM%Psurf

              ! distribute water to regions above the destruction zone
              ! water is reduced until saturation point is reached, then goes to next level
              do k_up=k+1,mesh%nr+1
                satur_w = pM%satur_w1 * ((mesh%rmax-mesh%rc(k_up)+pM%Psurf)*rhogD)**pM%satur_exp &
                 &+ pM%satur_w2 * ((mesh%rmax-mesh%rc(k_up)+pM%Psurf)*rhogD) ! ToDo: "/rhogD" instead of "*rhogD" ???
                if (k_up.eq.mesh%nr+1) then ! satur_w is NaN with formular above
                  satur_w = 1e+10_dp ! actually infinity, water goes back into ocean...
                endif
                if ((free_w/mesh%dVi(i,j,k_up) + cell%w(i,j,k_up)).gt.satur_w) then
                  old_w = cell%w(i,j,k_up)
                  free_w = free_w - (satur_w-cell%w(i,j,k_up))*mesh%dVi(i,j,k_up)
                  cell%w(i,j,k_up) = satur_w
                  do m=1,part%n
                    cl = part%cell_l(m)
                    cy = part%cell_y(m)
                    cr = part%cell_r(m)
                    if ((cl.eq.i).and.(cy.eq.j).and.(cr.eq.k_up)) then
                      part%w(m) = satur_w
                    endif
                  enddo
                  print *," goes to (:",i,",",j,",",k_up,"), residual: ",free_w,"; old w: ",old_w,", new w:",&
                           &cell%w(i,j,k_up), ", satur: ",satur_w
                else if (free_w.gt.0.0_dp) then !else if ((free_w/mesh%dVi(i,j,k_up) + cell%w(i,j,k_up)).gt.0.0_dp) then ! rest of the free water
                  !free_w = 0.0_dp
                  old_w = cell%w(i,j,k_up)
                  cell%w(i,j,k_up) = cell%w(i,j,k_up) + free_w/mesh%dVi(i,j,k_up)
                  free_w = max( 0.0_dp , free_w - cell%w(i,j,k_up)*mesh%dVi(i,j,k_up) )
                  do m=1,part%n
                    cl = part%cell_l(m)
                    cy = part%cell_y(m)
                    cr = part%cell_r(m)
                    if ((cl.eq.i).and.(cy.eq.j).and.(cr.eq.k_up)) then
                      part%w(m) = cell%w(i,j,k_up)
                    endif
                  enddo
  write(6,'(" residual goes to (",i3,",",i3,",",i3,"),residual: ",ES10.3,"; old w: ",ES10.3,", new w:",ES10.3,", satur: ",ES10.3)')&
                       & i,j,k_up,free_w,old_w,cell%w(i,j,k_up),satur_w
                endif
              enddo

            endif
           endif
          endif
        enddo
       enddo
      enddo

  end subroutine release_water

  ! if mantel material has more water than can be stored (above saturation value), the water is transported to upper regions and to surface
  subroutine saturate_water(mesh,cell,part,pM,rhogD,di)
      type(mesh_cp),intent(in)::mesh
      type(cell_properties),intent(inout)::cell
      type(particle),intent(inout)::part
      type(param_M),intent(in)::pM
      real(dp),intent(in)::rhogD
      integer,intent(in)::di

      integer::i,j,k,m,cl,cy,cr,k_up,jmin,jmax
      real(dp)::free_w,satur_w

!      print *,"Pre-Cut: ",sum(cell%w)/(mesh%nl+2)/(mesh%nr+2)

      if (di.eq.2) then
        jmin = 1
        jmax = 1
      else
        jmin = 0
        jmax = mesh%ny+1
      endif

      if (rhogD.eq.0.0_dp) then
        print *,"++++ Set rhogD = rho*g*D [in GPa] in input.txt! ++++"
      endif

      do i=0,mesh%nl+1
       do j=jmin,jmax
        do k=0,mesh%nr+1

          !write(*,*) i,j,k
          if (k.le.mesh%nr) then
            satur_w = pM%satur_w1 * ((mesh%rmax-mesh%rc(k)+pM%Psurf)*rhogD)**pM%satur_exp &
                  & + pM%satur_w2 * ((mesh%rmax-mesh%rc(k)+pM%Psurf)*rhogD)
          else
            satur_w = 0.0_dp
          endif
!          if (k.eq.mesh%nr+1) then
!            satur_w = 1e+10_dp !0.0_dp
!            !cell%w(i,j,k) = 0.0_dp ! boundary condition, no water can be stored at crust surface, is automatically released in form of oceans
!          endif
!          if (k.eq.10) then
!            write(*,*) k,satur_w
!          endif

          if (cell%w(i,j,k).gt.satur_w) then ! water that has to be released
              free_w = cell%w(i,j,k)*mesh%dVi(i,j,k) - satur_w*mesh%dVi(i,j,k)

              ! set water to saturation value in field and particles
              cell%w(i,j,k) = satur_w
              do m=1,part%n
                cl = part%cell_l(m)
                cy = part%cell_y(m)
                cr = part%cell_r(m)
                if ((cl.eq.i).and.(cy.eq.j).and.(cr.eq.k)) then
                  part%w(m) = satur_w
                endif
              enddo

              ! distribute water to regions above the destruction zone
              ! water is reduced until saturation point is reached, then goes to next level
              ! if water is left after k_up=nr+1, then water goes directly back into ocean (shallow dewatering, hydrothermal vents)
              if (k.le.mesh%nr) then
               do k_up=k+1,mesh%nr+1
!                write(*,*) "k_up=",k_up," in interval [",k+1,",",mesh%nr+1,"]"
                if (k_up.lt.mesh%nr+1) then
                  satur_w = pM%satur_w1 * ((mesh%rmax-mesh%rc(k_up)+pM%Psurf)*rhogD)**pM%satur_exp &
                        & + pM%satur_w2 * ((mesh%rmax-mesh%rc(k_up)+pM%Psurf)*rhogD)
                else
                  ! at surface, rc(k_up) larger than rmax and error would occur here
                  satur_w = 0.0_dp
                endif
                !write(*,*) k,k_up,mesh%dVi(i,j,k_up),satur_w
                if ((free_w/mesh%dVi(i,j,k_up) + cell%w(i,j,k_up)).gt.satur_w) then
!                  write(*,*) "Case 1"
                  free_w = free_w - (satur_w-cell%w(i,j,k_up))*mesh%dVi(i,j,k_up)
                  cell%w(i,j,k_up) = satur_w
                  do m=1,part%n
                    cl = part%cell_l(m)
                    cy = part%cell_y(m)
                    cr = part%cell_r(m)
                    if ((cl.eq.i).and.(cy.eq.j).and.(cr.eq.k_up)) then
                      part%w(m) = satur_w
                    endif
                  enddo
                else if ((free_w/mesh%dVi(i,j,k_up) + cell%w(i,j,k_up)).gt.0.0_dp) then ! rest of the free water
!                  write(*,*) "Case 2"
                  free_w = 0.0_dp
                  cell%w(i,j,k_up) = cell%w(i,j,k_up) + free_w/mesh%dVi(i,j,k_up)
                  do m=1,part%n
                    cl = part%cell_l(m)
                    cy = part%cell_y(m)
                    cr = part%cell_r(m)
                    if ((cl.eq.i).and.(cy.eq.j).and.(cr.eq.k_up)) then
                      part%w(m) = cell%w(i,j,k_up)
                    endif
                  enddo
                endif
!                write(*,*) "Finished for k_up=",k_up
               enddo
              endif

          endif
        enddo
       enddo
      enddo

!      print *,"Post-Cut: ",sum(cell%w)/(mesh%nl+2)/(mesh%nr+2)

  end subroutine saturate_water

  ! if pM%use_lat_heat.eq.-1 then subtr. lat heat from thermal field after energy solver used 
  ! if value is -2 then T is just set to T_solidus
  ! ToDo: add hydrated solidus?
  subroutine T_lat_heat(j2D,T,depl,ind_r,ref_r,entropy,lat_heat,T0,use_lat_heat,mc,maxD)
      integer,intent(in)::j2D,ind_r ! index j2 (jmin=0/1), index for reference density
      real(dp),intent(inout)::T(0:,j2D:,0:) !< T corrected: subtract latent heat from T field
      real(dp),intent(in)::depl(0:,j2D:,0:)
      real(dp),intent(in)::ref_r(ind_r:,j2D:,ind_r:)
      real(dp),intent(in)::entropy,lat_heat,T0
      integer,intent(in)::use_lat_heat
      type(melt_properties),intent(in)::mc
      real(dp),intent(in)::maxD
 
      integer::i,j,k,nl,ny,nr,jmin,jmax
      real(dp)::solidus,p,q,T_new,refR

      nl = UBound(T,1)-1
      ny = UBound(T,2)-1
      nr = UBound(T,3)-1
      if (ny.eq.0) then
        jmin = 1
        jmax = 1
      else
        jmin = 0
        jmax = ny+1
      endif

      do i=0,nl+1
       do j=jmin,jmax
        do k=0,nr+1
          ! use either entropy (=L/cp) or entropyT (=Delta S), because in the second version (T+T0) is multiplied to entropyT (now only if -2)

          solidus = mc%solidus_ref(i,j,k)+mc%solidus(i,j,k) ! ToDO: - hydrated solidus??? or since anyway only used for simple benchmarks, hydration effect not needed here?

          T_new = T(i,j,k) ! if for some reason none of the function below is used -> keep old T value

          if (ind_r.eq.1) then ! only one reference density value
            refR = ref_r(i,j,1)
          else
            refR = ref_r(i,j,k)
          endif

          if (use_lat_heat.eq.-2) then
            T_new = min(T(i,j,k),solidus) ! cut down T to solidus, also if DeltaF=0, which can be if depletion already reached max depletion...
          else if (use_lat_heat.eq.-3) then
            if ((T(i,j,k).gt.solidus).and.((depl(i,j,k).lt.maxD).or.(maxD.eq.0.0_dp))) then ! DeltaF not yet known, depends on corrected temperature
              if (entropy.eq.0.0_dp) then ! use latent heat (term: + rhoR*L*DF/Dt)
                T_new = ( T(i,j,k) + refR*lat_heat*solidus/(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k)) ) / &
                       & ( 1.0_dp + refR*lat_heat/(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k)) )
              else ! use entropy formulation (term: + (T+T0)*rhoR*entropy*DF/Dt)
                ! here term T^n appears also squared -> use p-q formula
                p = (mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k))/(refR*entropy) - solidus + T0
                q = -T0*( (mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k))/(refR*entropy) + solidus ) ! with this q, root is definitely positive, therefore use: p/2 + sqrt( (p/2)^2 - q )
                T_new = 0.5_dp*p + sqrt( (0.5_dp*p)**2 - q )

                ! if for some reason p-q formula gives a strange result, use instead demi-explicit version (T^nc + T0) with non-corrected temperature, could also be old temperature from last TS...
                if (.not.((T_new.ge.0.0_dp).and.(T_new.le.T(i,j,k)))) then
                  T_new = ( T(i,j,k) + refR*(T(i,j,k)+T0)*entropy*solidus/(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k)) ) / &
                       & ( 1.0_dp + refR*(T(i,j,k)+T0)*entropy/(mc%liquidus(i,j,k)-mc%solidus_ref(i,j,k)) )
                endif
              endif
546
547
            else if (T(i,j,k).gt.solidus) then ! max allowed depletion already reached, temp is not allowed to become larger than solidus
              T_new = min(T(i,j,k),solidus)
Lena Noack's avatar
Lena Noack committed
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
            endif
          else ! i.e. use_lat_heat = -1
            if (entropy.eq.0.0_dp) then ! use latent heat (term: + rhoR*L*DF/Dt)
              T_new = T(i,j,k) - refR*lat_heat*mc%DeltaF(i,j,k) !ToDo: where is "/Dt" ???
            else ! use entropy formulation (term: + (T+T0)*rhoR*entropy*DF/Dt)
              T_new = (T(i,j,k) - T0*refR*entropy*mc%DeltaF(i,j,k)) / (1.0_dp + refR*entropy*mc%DeltaF(i,j,k)) ! is T(i,j,k) if DeltaF = 0
            endif
          endif

          if (T_new.gt.T(i,j,k)) then
            write(*,*) "Strange, corrected temperature larger than original temperature from energy solver",i,j,k
            write(*,*) T_new,T(i,j,k),refR,lat_heat,mc%DeltaF(i,j,k)
            write(*,*) lat_heat,entropy
          endif

          T(i,j,k) = T_new
        enddo
       enddo
      enddo

  end subroutine T_lat_heat

570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
  ! temp not allowed to be above solidus if max depletion reached!
  subroutine cut_temp_solidus(j2D,T,depl,mc,maxD)
      integer,intent(in)::j2D               !< index j2 (jmin=0/1)
      real(dp),intent(inout)::T(0:,j2D:,0:) !< T cut to solidus if max_depl reached
      real(dp),intent(inout)::depl(0:,j2D:,0:)
      type(melt_properties),intent(in)::mc
      real(dp),intent(in)::maxD

      integer::i,j,k,nl,ny,nr,jmin,jmax
      real(dp)::solidus

      nl = UBound(T,1)-1
      ny = UBound(T,2)-1
      nr = UBound(T,3)-1
      if (ny.eq.0) then
        jmin = 1
        jmax = 1
      else
        jmin = 0
        jmax = ny+1
      endif

      do i=0,nl+1
       do j=jmin,jmax
        do k=0,nr+1
          solidus = mc%solidus_ref(i,j,k)+mc%solidus(i,j,k) ! ToDO: - hydrated solidus??? or since anyway only used for simple benchmarks, hydration effect not needed here?
                                    ! -solidus_hydrated(pM%sol_hyd,pM%sol_hyd_exp,cell%w(i,j,k)

          if ((T(i,j,k).gt.solidus).and.(depl(i,j,k).ge.maxD)) then
            T(i,j,k) = min(T(i,j,k),solidus) ! cut down T to solidus, if depletion already reached max depletion...
            depl(i,j,k) = maxD
          endif
        enddo
       enddo
      enddo

  end subroutine


Lena Noack's avatar
Lena Noack committed
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
  subroutine update_solidus_p(mesh,mc,part,pM)
      type(melt_properties),intent(in)::mc
      type(mesh_cp),intent(in)::mesh
      type(particle),intent(inout)::part
      type(param_M),intent(in)::pM

      integer::m,np,cl,cy,cr
      real(dp)::z

      np = part%n

      ! update solidus difference (saved in part%s / mc%solidus)
      ! hydration effect here not considered: water is lost when melting occurs. If T>solidus then T>hydrated solidus, T<solidus, then only w reduced, not the solidus functions
      if (pM%update_solidus.eq.1) then
        do m=1,np
          cl = part%cell_l(m)
          cy = part%cell_y(m)
          cr = part%cell_r(m)
          z = mesh%rmax - mesh%rc(cr) + pM%Psurf
          !if ((z.lt.pM%melt_depth).and.(part%d(m).gt.0.0_dp)) then ! not at larger depths, Delta_F=0.0 anyway...
              part%s(m) = part%d(m)*(mc%liquidus(cl,cy,cr)-mc%solidus_ref(cl,cy,cr))
              !if (part%d(m).gt.0.15_dp) write(*,*) cl,cr,part%d(m),part%df(m),part%s(m),mc%liquidus(cl,cy,cr),mc%solidus_ref(cl,cy,cr)
          !endif
!          if (part%s(m).gt.(mc%liquidus(cl,cy,cr)-mc%solidus_ref(cl,cy,cr))) then
!              part%s(m) = mc%liquidus(cl,cy,cr)-mc%solidus_ref(cl,cy,cr)
!          endif
        enddo
      endif
  end subroutine update_solidus_p


  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!                                                                                             !!
  !!                                     Treatment of crust                                      !!
  !!                                                                                             !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


647
  subroutine set_init_crust(mc,cell,mesh,pC,pI,pM,pX)
Lena Noack's avatar
Lena Noack committed
648
649
650
651
652
      type(melt_properties),intent(inout)::mc
      type(cell_properties),intent(inout)::cell
      type(mesh_cp),intent(in)::mesh
      type(param_C),intent(in)::pC
      type(param_I),intent(in)::pI
653
654
      type(param_M),intent(in)::pM
      type(param_X),intent(in)::pX
Lena Noack's avatar
Lena Noack committed
655
656
657
658
659

      integer::k

      mc%oc_vol_tot = 0.0_dp

660
661
!      do k=1,mesh%nr+1
      do k=0,mesh%nr
Lena Noack's avatar
Lena Noack committed
662
663
664
665
666
667
668
        if (mesh%rc(k).ge.(mesh%rmax-pC%crust_ini)) then !crust
          cell%r(:,:,k) = 7.0_dp
          cell%c(:,:,k) = pC%cont_OC_rho
          if (cell%fld_t) cell%ref_T(:,:,k) = pC%cont_OC_TRef
          if (cell%fld_v) cell%v(:,:,k) = pC%cont_OC_eta
          if (cell%fld_p) cell%cp(:,:,k)= pC%cont_OC_cp
          if (cell%fld_k) cell%k(:,:,k) = pC%cont_OC_k
669
670
671
672
673
674
675
676
677
678
679
          if (pX%dU238.eq.0.0_dp) then
            ! no redistribution of radioactive nuclei over time, instead constant pre-factor
            if (cell%fld_h) cell%h(:,:,k) = pC%cont_OC_h 
          else
            ! mantle initial value in initialize.f90, here set initial crustal values using crustal enrichment factor Crust_Lambda
            cell%h = 0.0_dp ! not used
            cell%X_U238(:,:,k) = cell%X_U238(1,1,1)*pM%Lambda ! at core-mantle boundary definitly mantle, not crust -> multiply with initial crustal enrichmant factor lambda
            cell%X_U235(:,:,k) = cell%X_U235(1,1,1)*pM%Lambda
            cell%X_Th232(:,:,k) = cell%X_Th232(1,1,1)*pM%Lambda
            cell%X_K40(:,:,k) = cell%X_K40(1,1,1)*pM%Lambda
          endif
Lena Noack's avatar
Lena Noack committed
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
          if (cell%fld_w) cell%w(:,:,k) = pC%cont_OC_w
          if (cell%fld_w) cell%wr(:,:,k) = pC%cont_OC_wr
          if (cell%fld_f) cell%fa(:,:,k) = pC%cont_OC_fr_angl
          if (cell%fld_f) cell%fc(:,:,k) = pC%cont_OC_fr_coh
          if (cell%fld_a) cell%ref_a(:,:,k) = pC%cont_OC_alpha
          if (mesh%rb(k).gt.(mesh%rmax-pC%crust_ini)) then ! whole cell in crust
            mc%oc_vol_tot(:,:) = mc%oc_vol_tot(:,:) + mesh%dVi(:,:,k)
          else
            mc%oc_vol_tot(:,:) = mc%oc_vol_tot(:,:) + mesh%dVi(:,:,k)*(mesh%rb(k+1)-(mesh%rmax-pC%crust_ini))/mesh%dr_vec(k)
          endif
        else if (mesh%rb(k+1).ge.(mesh%rmax-pC%crust_ini)) then !upper part of cell still in crust
          ! main part of cell in mantle, do not change properties of cell
          mc%oc_vol_tot(:,:) = mc%oc_vol_tot(:,:) + mesh%dVi(:,:,k)*(mesh%rb(k+1)-(mesh%rmax-pC%crust_ini))/mesh%dr_vec(k)
        endif
      enddo

696
697
      write(*,*) "+++Crust+++",pM%Lambda,cell%X_U238(1,1,1),cell%X_U238(1,1,mesh%nr)

Lena Noack's avatar
Lena Noack committed
698
699
700
701
  end subroutine


  ! add crust depending on new melt fraction
702
703
  ! includes redistribution of incompatible elements (only H2O and CO2 are treated seperately in atmosphere.f90)
  subroutine add_crust(mesh,cell,mc,part,t,dt,pC,pX)
Lena Noack's avatar
Lena Noack committed
704
705
706
707
708
709
      type(melt_properties),intent(inout)::mc
      type(mesh_cp),intent(in)::mesh
      type(cell_properties),intent(inout)::cell
      type(particle),intent(inout)::part
      real(dp),intent(in)::t,dt
      type(param_C),intent(in)::pC
710
      type(param_X),intent(in)::pX
Lena Noack's avatar
Lena Noack committed
711
712

      integer::i,j,k,m,nl,ny,nr,min_crust,shift_crust
713
714
      real(dp)::new_vol,dummy_oc,dummy_fc,dummy_cc,crust_old,crust_new
      real(dp)::X_U238_old,X_U235_old,X_Th232_old,X_K40_old
Lena Noack's avatar
Lena Noack committed
715
716
717
718
719
720
721
722
723
724

      nl = mesh%nl
      ny = mesh%ny
      nr = mesh%nr

      ! set these values to 0 to have only the new added crust, whereas *c_vol_tot carries old crustal volumes
      mc%oc_vol = 0.0_dp
      mc%fc_vol = 0.0_dp
      mc%cc_vol = 0.0_dp

725
726
727
728
729
      mc%X_U238_add = 0.0_dp
      mc%X_U235_add = 0.0_dp
      mc%X_Th232_add = 0.0_dp
      mc%X_K40_add = 0.0_dp

Lena Noack's avatar
Lena Noack committed
730
731
732
733
734
735
      do i=1,nl
       do j=1,ny
        do k=1,nr
          if (mc%DeltaF(i,j,k).ne.0.0_dp) then ! bring new crust to surface
            mc%cr_age(i,j) = t ! always assume that newest crust goes (partly) to the surface

736
737
738
739
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            !! Update crust thickness and crust composition !!
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Lena Noack's avatar
Lena Noack committed
740
741
742
743
744
            ! Primordial material -> basaltic crust
            if (cell%r(i,j,k).lt.1.5_dp) then ! is (in average) primordial mantle material or wedge material
              mc%oc_vol(i,j) = mc%oc_vol(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k) ! is added to basaltic crust
              mc%oc_vol_tot(i,j) = mc%oc_vol_tot(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)
            elseif ((cell%r(i,j,k).ge.3.5_dp).and.(cell%r(i,j,k).lt.4.5_dp)) then ! is (in average) primordial mantle material (TBL of OP)
745
              ! only for predefined subduction zones
Lena Noack's avatar
Lena Noack committed
746
747
              mc%oc_vol(i,j) = mc%oc_vol(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k) ! is added to basaltic crust
              mc%oc_vol_tot(i,j) = mc%oc_vol_tot(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)
748
              ! only for predefined subduction zones
Lena Noack's avatar
Lena Noack committed
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
            elseif ((cell%r(i,j,k).ge.6.5_dp).and.(cell%r(i,j,k).lt.7.5_dp)) then ! is (in average) primordial mantle material (TBL of SP)
              mc%oc_vol(i,j) = mc%oc_vol(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k) ! is added to basaltic crust
              mc%oc_vol_tot(i,j) = mc%oc_vol_tot(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)
  
            ! Felsic or continental crust -> continental crust
            elseif (cell%r(i,j,k).ge.8.5_dp) then ! is (in average) continental crust, surf r will not change; or (in average) felsic material, that now becomes continental crust
              if (pC%cont_CC_rho.ne.0.0_dp) then
                mc%cc_vol(i,j) = mc%cc_vol(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)
                mc%cc_vol_tot(i,j) = mc%cc_vol_tot(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)
              endif

            ! Basaltic crust -> Felsic crust
            else ! is (in average) basaltic material, that has been recycled (delamination or surface mobilization) or later on formed basaltic crust, and now becomes felsic crust
              if (pC%cont_FC_rho.ne.0.0_dp) then
                mc%fc_vol(i,j) = mc%fc_vol(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k) ! is added to felsic crust
                mc%fc_vol_tot(i,j) = mc%fc_vol_tot(i,j) + mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)
              endif
            endif

768
769
770
771
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            !! Deplete residual mantle !!
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!

772
773
774
775
776
777
778
779
780
            if ((pX%dU238.ne.0.0_dp).and.(cell%r(i,j,k).lt.0.5_dp)) then ! is (in average) primordial mantle material or wedge material (should be 0 if primordial mantle)
              !write(*,*) i,j,k,mc%DeltaF(i,j,k),(1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU238),cell%X_U238(i,j,k),&
              !              & cell%X_U238(i,j,k) * (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU238),cell%r(i,j,k)

              ! should not be applied to cell values but to particle values, otherwise cell values will be overwritten by particles with initial values from set_init_crust
              !cell%X_U238(i,j,k) = cell%X_U238(i,j,k) * (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU238)
              !cell%X_U235(i,j,k) = cell%X_U235(i,j,k) * (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU235)
              !cell%X_Th232(i,j,k) = cell%X_Th232(i,j,k) * (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dTh232)
              !cell%X_K40(i,j,k) = cell%X_K40(i,j,k) * (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dK40)
781
782

              ! sum up radioactive nuclei that go into crust, will be added to crust in set_new_crust
783
              mc%X_U238_add(i,j) = mc%X_U238_add(i,j) + mesh%dVi(i,j,k) * cell%X_U238(i,j,k) * &
784
                                            & (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU238))  ! ToDo: U238 in mol, weight or volume fraction? here used as volume fraction
785
              mc%X_U235_add(i,j) = mc%X_U235_add(i,j) + mesh%dVi(i,j,k) * cell%X_U235(i,j,k) * &
786
                                            & (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU235))
787
              mc%X_Th232_add(i,j) = mc%X_Th232_add(i,j) + mesh%dVi(i,j,k) * cell%X_Th232(i,j,k) * &
788
                                            & (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dTh232))
789
              mc%X_K40_add(i,j) = mc%X_K40_add(i,j) + mesh%dVi(i,j,k) * cell%X_K40(i,j,k) * &
790
                                            & (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dK40))
791
              !write(*,*) "    sum:",i,j,k,mc%X_U238_add(i,j),mc%X_U235_add(i,j),mc%X_Th232_add(i,j),mc%X_K40_add(i,j)
792
793
            endif

Lena Noack's avatar
Lena Noack committed
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
!            ! Re-set values
!            cell%r(i,j,k) = -20.0_dp ! ! neg so it is known later that particle has to change, will then be 0 after set_new_crust
!            ! Take average value from surrounding cells, that fill up this cell after depletion
!            cell%c(i,j,k) = ( cell%c(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%c(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                            & cell%c(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%c(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                            & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) ) 
!            cell%v(i,j,k) = ( cell%v(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%v(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                            & cell%v(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%v(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                            & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) )
!            cell%cp(i,j,k) = ( cell%cp(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%cp(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                             & cell%cp(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%cp(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                             & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) )
!            cell%k(i,j,k) = ( cell%k(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%k(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                            & cell%k(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%k(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                            & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) )
!            cell%h(i,j,k) = ( cell%h(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%h(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                            & cell%h(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%h(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                            & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) )
!            cell%w(i,j,k) = ( cell%w(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%w(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                            & cell%w(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%w(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                            & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) )
!            cell%wr(i,j,k) = ( cell%wr(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%wr(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                             & cell%wr(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%wr(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                             & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) )
!            cell%fa(i,j,k) = ( cell%fa(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%fa(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                             & cell%fa(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%fa(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                             & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) )
!            cell%fc(i,j,k) = ( cell%fc(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%fc(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                             & cell%fc(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%fc(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                             & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) )
!            cell%ref_a(i,j,k) = ( cell%ref_a(i-1,j,k)*mesh%dVi(i-1,j,k) + cell%ref_a(i+1,j,k)*mesh%dVi(i+1,j,k) + &
!                             & cell%ref_a(i,j,k-1)*mesh%dVi(i,j,k-1) + cell%ref_a(i,j,k-1)*mesh%dVi(i,j,k-1) ) &
!                             & / (mesh%dVi(i-1,j,k)+mesh%dVi(i+1,j,k)+mesh%dVi(i,j,k-1)+mesh%dVi(i,j,k-1) )
          endif
        enddo
829
830
831
832
833
834

        ! concentrations mc%X_..._add were multiplied with new melt volume
        ! get new average concentration in crust (take here OC,FC,CC similarly for now), below in set_crust these values will be set for FC,CC,OC
        crust_new = mc%oc_vol_tot(i,j) + mc%fc_vol_tot(i,j) + mc%cc_vol_tot(i,j)
        crust_old = crust_new - mc%oc_vol(i,j) - mc%fc_vol(i,j) - mc%cc_vol(i,j)

835
836
        !if (crust_new.gt.0.0_dp) then
        if (crust_new.gt.crust_old) then
837
838
839
840
          !X_U238_old = cell%X_U238(i,j,nr) ! only works if start with initial crust, otherwise should be zero instead of mantle value, or not? ToDo CHECK
          !X_U235_old = cell%X_U235(i,j,nr)
          !X_Th232_old = cell%X_Th232(i,j,nr)
          !X_K40_old = cell%X_K40(i,j,nr)
841

842
843
          !write(*,*) "   New: ",crust_new,crust_old,i,j,k,X_U238_old,(mc%X_U238_add(i,j)+X_U238_old*crust_old)/crust_new

844
845
846
847
848
849
850
851
          !mc%X_U238_add(i,j) = ( mc%X_U238_add(i,j) + X_U238_old*crust_old ) / crust_new ! new average crust value for i,j
          !mc%X_U235_add(i,j) = ( mc%X_U235_add(i,j) + X_U235_old*crust_old ) / crust_new
          !mc%X_Th232_add(i,j) = ( mc%X_Th232_add(i,j) + X_Th232_old*crust_old ) / crust_new
          !mc%X_K40_add(i,j) = ( mc%X_K40_add(i,j) + X_K40_old*crust_old ) / crust_new
          mc%X_U238_add(i,j) = mc%X_U238_add(i,j) / crust_new ! what will be added to particles in new crust
          mc%X_U235_add(i,j) = mc%X_U235_add(i,j) / crust_new
          mc%X_Th232_add(i,j) = mc%X_Th232_add(i,j) / crust_new
          mc%X_K40_add(i,j) = mc%X_K40_add(i,j) / crust_new
852
853
        else
          ! no changes in local composition
854
855
856
857
858
859
860
861
          !mc%X_U238_add(i,j) = cell%X_U238(i,j,nr)
          !mc%X_U235_add(i,j) = cell%X_U235(i,j,nr)
          !mc%X_Th232_add(i,j) = cell%X_Th232(i,j,nr)
          !mc%X_K40_add(i,j) = cell%X_K40(i,j,nr)
          mc%X_U238_add(i,j) = 0.0_dp          
          mc%X_U235_add(i,j) = 0.0_dp
          mc%X_Th232_add(i,j) = 0.0_dp
          mc%X_K40_add(i,j) = 0.0_dp
862
863
        endif

864
865
866

        write(*,*) i,j,k,mc%X_U238_add(i,j),mc%X_U235_add(i,j),mc%X_Th232_add(i,j),mc%X_K40_add(i,j),crust_new,crust_old

Lena Noack's avatar
Lena Noack committed
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
       enddo
      enddo

      if ((pC%LHB_time.gt.0.0_dp).and.(pC%LHB_time.gt.t-dt).and.(pC%LHB_time.le.t)) then !re-set crust age and values at LHB
        write(*,*) "LHB, re-set crust to LHB age"
        mc%oc_vol_tot = mc%oc_vol_tot + mc%fc_vol_tot + mc%cc_vol_tot
        mc%fc_vol_tot = 0.0_dp
        mc%cc_vol_tot = 0.0_dp
        mc%cr_age = t
      endif

  end subroutine add_crust


  !ToDO: Adapt to 3D !!!
  !ToDo: until now assumes only uniform grid...
  subroutine move_crust(mc,mesh,field,dt)
      type(melt_properties),intent(inout)::mc
      type(mesh_cp),intent(in)::mesh
      type(variables_unknowns),intent(in)::field
      real(dp),intent(in)::dt

      integer::i,j,nl,shift_i
      real(dp)::shift
      real(dp),allocatable::dummy_oc(:,:),dummy_fc(:,:),dummy_cc(:,:)

      nl = mesh%nl
      j = 1

      allocate( dummy_oc(0:nl+1,1), dummy_fc(0:nl+1,1), dummy_cc(0:nl+1,1) )

      dummy_oc = 0.0_dp
      dummy_fc = 0.0_dp
      dummy_cc = 0.0_dp

      do i=1,nl
        ! determine shift of crustal material
        shift = field%u(i,j,mesh%nr)*dt ! can be positive or negative
        shift_i = 0
        if (abs(shift).gt.mesh%dl_vec(i)) then
         if (shift.gt.mesh%dl_vec(i)) then
          do while (shift.gt.mesh%dl_vec(i+shift_i))
            shift_i = shift_i + 1
            if ((i+shift_i).gt.nl+1) write(*,*) "Error Melt r: i+shift_i out of range with ",i+shift_i,", i=",i
            shift = shift - mesh%dl_vec(i+shift_i)
          enddo
         else if ((-shift).gt.mesh%dl_vec(i)) then
          do while ((-shift).gt.mesh%dl_vec(i+shift_i))
            shift_i = shift_i - 1
            if ((i+shift_i).lt.0) write(*,*) "Error Melt l: i+shift_i out of range with ",i+shift_i,", i=",i
            shift = shift + mesh%dl_vec(i+shift_i)
          enddo
         endif
        endif

        ! boundary treatment ;ToDo: add reflective treatment?
        ! can there be an index problem, i.e. i+shift_i not in [0,nl+1] ?
        ! set dummy[0] -> dummy[1] and dummy[nl+1] -> dummy[nl] ?
        ! or average dummy values at boundary?
        if (((i+shift_i).lt.0).or.((i+shift_i).gt.nl+1)) write(*,*) "Error Melt: i+shift_i out of range with ",i+shift_i,i,shift_i

        ! apply to crust via dummy vector
        dummy_oc(i+shift_i,j) = dummy_oc(i+shift_i,j) + mc%oc_vol_tot(i,j)
        dummy_fc(i+shift_i,j) = dummy_fc(i+shift_i,j) + mc%fc_vol_tot(i,j)
        dummy_cc(i+shift_i,j) = dummy_cc(i+shift_i,j) + mc%cc_vol_tot(i,j)
      enddo
      ! ToDo: bnd dummy values

      ! set new crust vectors
      mc%oc_vol_tot = dummy_oc
      mc%fc_vol_tot = dummy_fc
      mc%cc_vol_tot = dummy_cc

      deallocate( dummy_oc,dummy_fc,dummy_cc )

  end subroutine move_crust

  ! only called if pM%add_crust_depth.gt.0.0_dp
  subroutine distr_crust(mesh,mc,ac_depth,add_crust_incr)
      type(melt_properties),intent(inout)::mc
      type(mesh_cp),intent(in)::mesh
      real(dp),intent(inout)::ac_depth  ! add crust depth
      integer,intent(in)::add_crust_incr

      integer::i,j,d,i_2
      real(dp)::ac_vol,cr_vol,shift_vol,shift_vol_r,shift_loc

      if (mesh%geom.eq.1) then
        d=2
      else
        d=1
      endif

      ! check if total crust volume is more than allowed, then increase depth accordingly if pM%add_crust_incr=1
      ! ToDo: add for non-uniform lateral angles
      cr_vol = sum(mc%oc_vol_tot(1:mesh%nl,1:mesh%ny))/real(mesh%nl*mesh%ny) ! average crust volume
      if (mesh%geom.eq.0) then
        ! 2D or 3D box
        ac_vol = mesh%dl_vec(1)*mesh%dy_vec(1)*ac_depth
      else
        ! here only 2D version
        ac_vol = mesh%dl_vec(1)*((mesh%rmax-0.5_dp*ac_depth)**d)*ac_depth
      endif

      if ((add_crust_incr.eq.1).and.(cr_vol.gt.ac_vol)) then
        ! increase ac_vol / ac_depth
        if (mesh%geom.eq.0) then
          ! 2D or 3D box
          ac_depth = cr_vol / (mesh%dl_vec(1)*mesh%dy_vec(1))
        else
          ! here only 2D version, assume that new ac_depth is similar to old ac_depth (needed for radius)
          ac_depth = cr_vol / (mesh%dl_vec(1)*((mesh%rmax-0.5_dp*ac_depth)**d))
        endif
      endif

      do i=1,mesh%nl
        do j=1,mesh%ny
          if (mesh%geom.eq.0) then
            ! 2D or 3D box
            ac_vol = mesh%dl_vec(i)*mesh%dy_vec(j)*ac_depth
          else
            ! here only 2D version
            ac_vol = mesh%dl_vec(i)*((mesh%rmax-0.5_dp*ac_depth)**d)*ac_depth
          endif

          if (mc%oc_vol_tot(i,j).gt.ac_vol) then ! more crust than allowed
            shift_vol = mc%oc_vol_tot(i,j) - ac_vol
            mc%oc_vol_tot(i,j) = ac_vol
            mc%oc_vol(i,j) = mc%oc_vol(i,j) - shift_vol

            shift_vol = 0.5*shift_vol
            shift_vol_r = shift_vol

            ! assume here that ac_vol is the same for all columns
For faster browsing, not all history is shown. View entire blame