source: trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 @ 3532

Last change on this file since 3532 was 3498, checked in by jbclement, 7 weeks ago

PEM:

  • Correction of the variable name for the ice table depth in "pemetat0.F90". So it is now got as intended from the "startpem.nc" file;
  • Renaming of the tendencies in the PEM with the prefix 'd_' instead of 'tend_';
  • Modification of the PEM time step type from integer to real. As a consequence, all time variables are now of real type. This change adds the possibility to consider fractions of year as time step.

JBC

File size: 30.7 KB
Line 
1MODULE layering_mod
2
3!=======================================================================
4! Purpose: Manage layered deposits in the PEM with a linked list data structure
5!
6! Author: JBC
7! Date: 08/03/2024
8!=======================================================================
9
10
11implicit none
12
13!---- Declarations
14! Numerical parameter
15real, parameter :: eps = epsilon(1.), tol = 1.e2*eps
16integer         :: nb_str_max
17logical         :: layering_algo
18
19! Physical parameters
20real, parameter :: d_dust = 5.78e-2             ! Tendency of dust [kg.m-2.y-1]
21real, parameter :: dry_lag_porosity = 0.5       ! Porosity of dust lag
22real, parameter :: dry_regolith_porosity = 0.45 ! Porosity of regolith
23real, parameter :: co2ice_porosity = 0.3        ! Porosity of CO2 ice
24real, parameter :: h2oice_porosity = 0.3        ! Porosity of H2O ice
25real, parameter :: rho_co2ice = 1650.           ! Density of CO2 ice [kg.m-3]
26real, parameter :: rho_h2oice = 920.            ! Density of H2O ice [kg.m-3]
27real, parameter :: rho_dust = 2500.             ! Density of dust [kg.m-3]
28real, parameter :: rho_rock = 3200.             ! Density of rock [kg.m-3]
29
30! Lag layer parameters -> see Levrard et al. 2007
31real, parameter :: hmin_lag = 0.5  ! Minimal height of the lag deposit to reduce the sblimation rate
32real, parameter :: fred_subl = 0.1 ! Reduction factor of sublimation rate
33
34! Stratum = node of the linked list
35type :: stratum
36    real                   :: thickness      ! Layer thickness [m]
37    real                   :: top_elevation  ! Layer top_elevation (top height from the surface) [m]
38    real                   :: co2ice_volfrac ! CO2 ice volumetric fraction
39    real                   :: h2oice_volfrac ! H2O ice volumetric fraction
40    real                   :: dust_volfrac   ! Dust volumetric fraction
41    real                   :: air_volfrac    ! Air volumetric fraction inside pores
42    type(stratum), pointer :: up => null()   ! Upper stratum (next node)
43    type(stratum), pointer :: down => null() ! Lower stratum (previous node)
44end type stratum
45
46! Layering = linked list
47type layering
48    integer                :: nb_str = 0       ! Number of strata
49    type(stratum), pointer :: bottom => null() ! First stratum (head of the list)
50    type(stratum), pointer :: top => null()    ! Last stratum (tail of the list)
51end type layering
52
53! Array of pointers to "stratum"
54type ptrarray
55    type(stratum), pointer :: p
56end type ptrarray
57
58!=======================================================================
59contains
60! Procedures to manipulate the data structure:
61!     > print_stratum
62!     > add_stratum
63!     > insert_stratum
64!     > rm_stratum
65!     > get_stratum
66!     > ini_layering
67!     > del_layering
68!     > print_layering
69!     > get_nb_str_max
70!     > stratif2array
71!     > array2stratif
72!     > print_stratif
73! Procedures for the algorithm to evolve the layering:
74!     > make_layering
75!     > subl_co2ice_layering
76!     > subl_h2oice_layering
77!     > lift_dust_lags
78!=======================================================================
79! To display a stratum
80SUBROUTINE print_stratum(this)
81
82implicit none
83
84!---- Arguments
85type(stratum), intent(in) :: this
86
87!---- Code
88write(*,'(a20,es13.6)') 'Thickness      : ', this%thickness
89write(*,'(a20,es13.6)') 'Top elevation  : ', this%top_elevation
90write(*,'(a20,es13.6)') 'CO2 ice (m3/m3): ', this%co2ice_volfrac
91write(*,'(a20,es13.6)') 'H2O ice (m3/m3): ', this%h2oice_volfrac
92write(*,'(a20,es13.6)') 'Dust (m3/m3)   : ', this%dust_volfrac
93write(*,'(a20,es13.6)') 'Air (m3/m3)    : ', this%air_volfrac
94
95END SUBROUTINE print_stratum
96
97!=======================================================================
98! To add a stratum to the top of a layering
99SUBROUTINE add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,air)
100
101implicit none
102
103!---- Arguments
104type(layering), intent(inout)  :: this
105real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, air
106
107!---- Local variables
108type(stratum), pointer :: str
109
110!---- Code
111! Creation of the stratum
112allocate(str)
113nullify(str%up,str%down)
114str%thickness = 1.
115str%top_elevation = 1.
116str%co2ice_volfrac = 0.
117str%h2oice_volfrac = 0.
118str%dust_volfrac = 0.
119str%air_volfrac = 1.
120if (present(thickness)) str%thickness = thickness
121if (present(top_elevation)) str%top_elevation = top_elevation
122if (present(co2ice)) str%co2ice_volfrac = co2ice
123if (present(h2oice)) str%h2oice_volfrac = h2oice
124if (present(dust)) str%dust_volfrac = dust
125if (present(air)) str%air_volfrac = air
126
127! Verification of volume fraction
128if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) &
129     error stop 'add_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!'
130
131! Increment the number of strata
132this%nb_str = this%nb_str + 1
133
134! Add the stratum to the layering
135if (.not. associated(this%bottom)) then ! If it is the very first one
136    this%bottom => str
137    this%top => str
138else
139    str%down => this%top
140    this%top%up => str
141    this%top => str
142endif
143
144END SUBROUTINE add_stratum
145
146!=======================================================================
147! To insert a stratum after another one in a layering
148SUBROUTINE insert_stratum(this,str_prev,thickness,top_elevation,co2ice,h2oice,dust,air)
149
150implicit none
151
152!---- Arguments
153type(layering),         intent(inout) :: this
154type(stratum), pointer, intent(inout) :: str_prev
155real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, air
156
157!---- Local variables
158type(stratum), pointer :: str, current
159
160!---- Code
161if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering
162    call add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,air)
163else
164    ! Creation of the stratum
165    allocate(str)
166    nullify(str%up,str%down)
167    str%thickness = 1.
168    str%top_elevation = 1.
169    str%co2ice_volfrac = 0.
170    str%h2oice_volfrac = 0.
171    str%dust_volfrac = 0.
172    str%air_volfrac = 1.
173    if (present(thickness)) str%thickness = thickness
174    if (present(top_elevation)) str%top_elevation = top_elevation
175    if (present(co2ice)) str%co2ice_volfrac = co2ice
176    if (present(h2oice)) str%h2oice_volfrac = h2oice
177    if (present(dust)) str%dust_volfrac = dust
178    if (present(air)) str%air_volfrac = air
179
180    ! Verification of volume fraction
181    if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) &
182         error stop 'insert_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!'
183
184    ! Increment the number of strata
185    this%nb_str = this%nb_str + 1
186
187    ! Insert the stratum into the layering at the right place
188    str%down => str_prev
189    str%up => str_prev%up
190    str_prev%up => str
191    str%up%down => str
192
193    ! Correct the value of top_elevation for the upper strata
194    current => str%up
195    do while (associated(current))
196        current%top_elevation = current%down%top_elevation + current%thickness
197        current => current%up
198    enddo
199endif
200
201END SUBROUTINE insert_stratum
202
203!=======================================================================
204! To remove a stratum in a layering
205SUBROUTINE rm_stratum(this,str)
206
207implicit none
208
209!---- Arguments
210type(layering),         intent(inout) :: this
211type(stratum), pointer, intent(inout) :: str
212
213!---- Local variables
214type(stratum), pointer :: new_str
215
216!---- Code
217! Protect the 3 sub-surface strata from removing
218if (str%top_elevation <= 0.) return
219
220! Decrement the number of strata
221this%nb_str = this%nb_str - 1
222
223! Remove the stratum from the layering
224str%down%up => str%up
225if (associated(str%up)) then ! If it is not the last one at the top of the layering
226    str%up%down => str%down
227else
228    this%top => str%down
229endif
230new_str => str%down
231nullify(str%up,str%down)
232deallocate(str)
233nullify(str)
234str => new_str
235
236END SUBROUTINE rm_stratum
237
238!=======================================================================
239! To get a specific stratum in a layering
240FUNCTION get_stratum(lay,i) RESULT(str)
241
242implicit none
243
244!---- Arguments
245type(layering), intent(in) :: lay
246integer,        intent(in) :: i
247type(stratum), pointer     :: str
248
249!---- Local variables
250integer :: k
251
252!---- Code
253if (i < 1 .or. i > lay%nb_str) error stop 'get_stratum: requested number is beyond the number of strata of the layering!'
254k = 1
255str => lay%bottom
256do while (k /= i)
257    str => str%up
258    k = k + 1
259enddo
260
261END FUNCTION get_stratum
262
263!=======================================================================
264! To initialize the layering
265SUBROUTINE ini_layering(this)
266
267implicit none
268
269!---- Arguments
270type(layering), intent(inout) :: this
271
272!---- Code
273! Creation of three strata at the bottom of the layering
274! 1) Dry porous deep regolith
275!call add_stratum(this,18.,-12.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
276! 2) Ice-cemented regolith (ice table)
277!call add_stratum(this,10.,-2.,0.,1. - dry_regolith_porosity,dry_regolith_porosity,0.)
278! 3) Dry porous regolith
279call add_stratum(this,2.,0.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
280
281END SUBROUTINE ini_layering
282
283!=======================================================================
284! To delete the entire layering
285SUBROUTINE del_layering(this)
286
287implicit none
288
289!---- Arguments
290type(layering), intent(inout) :: this
291
292!---- Local variables
293type(stratum), pointer :: current, next
294
295!---- Code
296current => this%bottom
297do while (associated(current))
298    next => current%up
299    deallocate(current)
300    nullify(current)
301    current => next
302enddo
303this%nb_str = 0
304
305END SUBROUTINE del_layering
306
307!=======================================================================
308! To display a layering
309SUBROUTINE print_layering(this)
310
311implicit none
312
313!---- Arguments
314type(layering), intent(in) :: this
315
316!---- Local variables
317type(stratum), pointer :: current
318integer                :: i
319
320!---- Code
321current => this%top
322
323i = this%nb_str
324do while (associated(current))
325    write(*,'(a8,i4)') 'Stratum ', i
326    call print_stratum(current)
327    write(*,*)
328    current => current%down
329    i = i - 1
330enddo
331
332END SUBROUTINE print_layering
333
334!=======================================================================
335! To get the maximum number of "stratum" in the stratification (layerings)
336FUNCTION get_nb_str_max(stratif,ngrid,nslope) RESULT(nb_str_max)
337
338implicit none
339
340!---- Arguments
341type(layering), dimension(ngrid,nslope), intent(in) :: stratif
342integer,                                 intent(in) :: ngrid, nslope
343integer :: nb_str_max
344
345!---- Local variables
346integer :: ig, islope
347
348!---- Code
349nb_str_max = 0
350do islope = 1,nslope
351    do ig = 1,ngrid
352        nb_str_max = max(stratif(ig,islope)%nb_str,nb_str_max)
353    enddo
354enddo
355
356END FUNCTION get_nb_str_max
357
358!=======================================================================
359! To convert the stratification data structure (layerings) into an array able to be outputted
360SUBROUTINE stratif2array(stratif,ngrid,nslope,stratif_array)
361
362implicit none
363
364!---- Arguments
365integer,                                 intent(in) :: ngrid, nslope
366type(layering), dimension(ngrid,nslope), intent(in) :: stratif
367real, dimension(:,:,:,:), allocatable, intent(inout) :: stratif_array
368
369!---- Local variables
370type(stratum), pointer :: current
371integer                :: ig, islope, k
372
373!---- Code
374stratif_array = 0.
375do islope = 1,nslope
376    do ig = 1,ngrid
377        current => stratif(ig,islope)%bottom
378        k = 1
379        do while (associated(current))
380            stratif_array(ig,islope,k,1) = current%thickness
381            stratif_array(ig,islope,k,2) = current%top_elevation
382            stratif_array(ig,islope,k,3) = current%co2ice_volfrac
383            stratif_array(ig,islope,k,4) = current%h2oice_volfrac
384            stratif_array(ig,islope,k,5) = current%dust_volfrac
385            stratif_array(ig,islope,k,6) = current%air_volfrac
386            current => current%up
387            k = k + 1
388        enddo
389    enddo
390enddo
391
392END SUBROUTINE stratif2array
393
394!=======================================================================
395! To convert the stratification array into the stratification data structure (layerings)
396SUBROUTINE array2stratif(stratif_array,ngrid,nslope,stratif)
397
398implicit none
399
400!---- Arguments
401integer,                               intent(in) :: ngrid, nslope
402real, dimension(:,:,:,:), allocatable, intent(in) :: stratif_array
403type(layering), dimension(ngrid,nslope), intent(inout) :: stratif
404
405!---- Local variables
406integer :: ig, islope, k
407
408!---- Code
409do islope = 1,nslope
410    do ig = 1,ngrid
411        stratif(ig,islope)%bottom%thickness      = stratif_array(ig,islope,1,1)
412        stratif(ig,islope)%bottom%top_elevation  = stratif_array(ig,islope,1,2)
413        stratif(ig,islope)%bottom%co2ice_volfrac = stratif_array(ig,islope,1,3)
414        stratif(ig,islope)%bottom%h2oice_volfrac = stratif_array(ig,islope,1,4)
415        stratif(ig,islope)%bottom%dust_volfrac   = stratif_array(ig,islope,1,5)
416        stratif(ig,islope)%bottom%air_volfrac    = stratif_array(ig,islope,1,6)
417        do k = 2,size(stratif_array,3)
418            if (all(abs(stratif_array(ig,islope,k,:)) < eps)) exit
419            call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6))
420        enddo
421    enddo
422enddo
423
424END SUBROUTINE array2stratif
425
426!=======================================================================
427! To display a stratification (layerings)
428SUBROUTINE print_stratif(this,ngrid,nslope)
429
430implicit none
431
432!---- Arguments
433integer,                                 intent(in) :: ngrid, nslope
434type(layering), dimension(ngrid,nslope), intent(in) :: this
435
436!---- Local variables
437integer :: ig, islope
438
439!---- Code
440do ig = 1,ngrid
441    write(*,'(a10,i4)') 'Grid cell ', ig
442    do islope = 1,nslope
443        write(*,'(a13,i1)') 'Slope number ', islope
444        call print_layering(this(ig,islope))
445        write(*,*) '~~~~~~~~~~'
446    enddo
447    write(*,*) '--------------------'
448enddo
449
450END SUBROUTINE print_stratif
451
452!=======================================================================
453! Layering algorithm
454SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)
455
456implicit none
457
458!---- Arguments
459type(layering),         intent(inout) :: this
460type(stratum), pointer, intent(inout) :: current1, current2
461logical,                intent(inout) :: new_str, new_lag1, new_lag2
462integer,                intent(inout) :: stopPEM
463real,     intent(in) :: d_co2ice, d_h2oice, d_dust
464
465!---- Local variables
466real                   :: dh_co2ice, dh_h2oice, dh_dust
467real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
468real                   :: thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
469real                   :: h_toplag
470type(stratum), pointer :: currentloc
471
472!---- Code
473dh_co2ice = d_co2ice/rho_co2ice
474dh_h2oice = d_h2oice/rho_h2oice
475dh_dust = d_dust/rho_dust
476
477if (dh_dust < 0.) then ! Dust lifting only
478    if (abs(dh_co2ice) > eps .or. abs(dh_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
479    write(*,'(a)') ' Stratification -> Dust lifting'
480    h2lift_tot = abs(dh_dust)
481    do while (h2lift_tot > 0. .and. associated(current1))
482        ! Is the considered stratum a dust lag?
483        if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we can lift dust
484            h_dust_old = current1%dust_volfrac*current1%thickness
485            ! How much dust can be lifted in the considered dust lag?
486            if (h2lift_tot - h_dust_old <= eps) then ! There is enough to lift
487                h2lift = h2lift_tot
488                h2lift_tot = 0.
489                call lift_dust_lags(this,current1,h2lift)
490            else ! Only a fraction can be lifted and so we move to the underlying stratum
491                h2lift = h_dust_old
492                h2lift_tot = h2lift_tot - h2lift
493                call lift_dust_lags(this,current1,h2lift)
494                current1 => current1%down
495            endif
496        else ! No, we stop
497            write(*,'(a)') ' Stratification -> Dust lifting: no dust lag!'
498            stopPEM = 8
499            exit
500        endif
501    enddo
502    if (h2lift_tot > 0.) then
503        write(*,'(a)') ' Stratification -> Not enough dust in the lag layers to complete the lifting!'
504        stopPEM = 8
505    endif
506
507else
508
509!------ Dust sedimentation only
510    if (abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then
511        write(*,'(a)') ' Stratification -> Dust sedimentation'
512        ! New stratum at the layering top by sedimentation of dust
513        thickness = dh_dust/(1. - dry_regolith_porosity)
514        if (thickness > 0.) then ! Only if the stratum is building up
515             if (new_str) then
516                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity)
517                 new_str = .false.
518             else
519                 this%top%thickness = this%top%thickness + thickness
520                 this%top%top_elevation = this%top%top_elevation + thickness
521             endif
522        endif
523
524!------ Condensation of CO2 ice + H2O ice
525    else if ((dh_co2ice >= 0. .and. dh_h2oice > 0.) .or. (dh_co2ice > 0. .and. dh_h2oice >= 0.)) then
526        write(*,'(a)') ' Stratification -> CO2 and H2O ice condensation'
527        ! New stratum at the layering top by condensation of CO2 and H2O ice
528        thickness = dh_co2ice/(1. - co2ice_porosity) + dh_h2oice/(1. - h2oice_porosity) + dh_dust
529        if (thickness > 0.) then ! Only if the stratum is building up
530             if (new_str) then
531                 call add_stratum(this,thickness,this%top%top_elevation + thickness,dh_co2ice/thickness,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_co2ice/thickness + dh_h2oice/thickness + dh_dust/thickness))
532                 new_str = .false.
533             else
534                 this%top%thickness = this%top%thickness + thickness
535                 this%top%top_elevation = this%top%top_elevation + thickness
536             endif
537        endif
538
539!------ Sublimation of CO2 ice + Condensation of H2O ice
540    else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then
541        write(*,'(a)') ' Stratification -> Sublimation of CO2 ice + Condensation of H2O ice'
542        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
543        h2subl_tot = abs(dh_co2ice)
544        do while (h2subl_tot > 0. .and. associated(current1))
545            h_co2ice_old = current1%co2ice_volfrac*current1%thickness
546            h_h2oice_old = current1%h2oice_volfrac*current1%thickness
547            h_dust_old = current1%dust_volfrac*current1%thickness
548
549            ! How much is CO2 ice sublimation inhibited by the top dust lag?
550            h_toplag = 0.
551            if (associated(current1%up)) then ! If there is an upper stratum
552                currentloc => current1%up
553                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
554                    h_toplag = h_toplag + currentloc%thickness
555                    if (.not. associated(currentloc%up)) exit
556                    currentloc => currentloc%up
557                enddo
558            endif
559            h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
560
561            ! How much CO2 ice can sublimate in the considered stratum?
562            if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
563                h2subl = h2subl_tot
564                h2subl_tot = 0.
565                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
566            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
567                ! Is the considered stratum a dust lag?
568                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
569                    current1 => current1%down
570                    new_lag1 = .true.
571                else ! No, so we stop here
572                    exit
573                endif
574            else ! Only a fraction can sublimate and so we move to the underlying stratum
575                h2subl = h_co2ice_old
576                h2subl_tot = h2subl_tot - h2subl
577                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
578                current1 => current1%down
579                new_lag1 = .true.
580            endif
581        enddo
582        if (h2subl_tot > 0.) then
583            write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
584            stopPEM = 8
585        endif
586        ! New stratum at the layering top by condensation of H2O ice
587        thickness = dh_h2oice/(1. - h2oice_porosity) + dh_dust
588        if (thickness > 0.) then ! Only if the stratum is building up
589             if (new_str) then
590                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness))
591                 new_str = .false.
592             else
593                 this%top%thickness = this%top%thickness + thickness
594                 this%top%top_elevation = this%top%top_elevation + thickness
595             endif
596        endif
597
598!------ Sublimation of CO2 ice + H2O ice
599    else if ((dh_co2ice <= 0. .and. dh_h2oice < 0.) .or. (dh_co2ice < 0. .and. dh_h2oice <= 0.)) then
600        write(*,'(a)') ' Stratification -> Sublimation of CO2 and H2O ice'
601        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
602        h2subl_tot = abs(dh_co2ice)
603        do while (h2subl_tot > 0. .and. associated(current1))
604            h_co2ice_old = current1%co2ice_volfrac*current1%thickness
605            h_h2oice_old = current1%h2oice_volfrac*current1%thickness
606            h_dust_old = current1%dust_volfrac*current1%thickness
607
608            ! How much is CO2 ice sublimation inhibited by the top dust lag?
609            h_toplag = 0.
610            if (associated(current1%up)) then ! If there is an upper stratum
611                currentloc => current1%up
612                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
613                    h_toplag = h_toplag + currentloc%thickness
614                    if (.not. associated(currentloc%up)) exit
615                    currentloc => currentloc%up
616                enddo
617            endif
618            h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
619
620            ! How much CO2 ice can sublimate in the considered stratum?
621            if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
622                h2subl = h2subl_tot
623                h2subl_tot = 0.
624                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
625            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
626                ! Is the considered stratum a dust lag?
627                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
628                    current1 => current1%down
629                    new_lag1 = .true.
630                else ! No, so we stop here
631                    exit
632                endif
633            else ! Only a fraction can sublimate and so we move to the underlying stratum
634                h2subl = h_co2ice_old
635                h2subl_tot = h2subl_tot - h2subl
636                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
637                current1 => current1%down
638                new_lag1 = .true.
639            endif
640        enddo
641        if (h2subl_tot > 0.) then
642            write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
643            stopPEM = 8
644        endif
645        ! H2O ice sublimation in the considered stratum + New stratum for dust lag
646        h2subl_tot = abs(dh_h2oice)
647        do while (h2subl_tot > 0. .and. associated(current2))
648            h_co2ice_old = current2%co2ice_volfrac*current2%thickness
649            h_h2oice_old = current2%h2oice_volfrac*current2%thickness
650            h_dust_old = current2%dust_volfrac*current2%thickness
651
652            ! How much is H2O ice sublimation inhibited by the top dust lag?
653            h_toplag = 0.
654            if (associated(current2%up)) then ! If there is an upper stratum
655                currentloc => current2%up
656                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
657                    h_toplag = h_toplag + currentloc%thickness
658                    if (.not. associated(currentloc%up)) exit
659                    currentloc => currentloc%up
660                enddo
661            endif
662            h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
663
664            ! How much H2O ice can sublimate in the considered stratum?
665            if (h2subl_tot - h_h2oice_old <= eps) then ! There is enough to sublimate
666                h2subl = h2subl_tot
667                h2subl_tot = 0.
668                call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag2)
669            else if (h_h2oice_old < eps) then ! There is nothing to sublimate
670                ! Is the considered stratum a dust lag?
671                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
672                    current1 => current1%down
673                    new_lag1 = .true.
674                else ! No, so we stop here
675                    exit
676                endif
677            else ! Only a fraction can sublimate and so we move to the underlying stratum
678                h2subl = h_h2oice_old
679                h2subl_tot = h2subl_tot - h2subl
680                call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag2)
681                current2 => current2%down
682                new_lag2 = .true.
683            endif
684        enddo
685        if (h2subl_tot > 0.) then
686            write(*,'(a)') ' Stratification -> Not enough H2O ice in the layering to complete the sublimation!'
687            stopPEM = 8
688        endif
689
690!------ Condensation of CO2 ice + Sublimation of H2O ice
691    else if (dh_co2ice > 0. .and. dh_h2oice < 0.) then
692        error stop 'Impossible situation: CO2 ice condensation + H2O ice sublimation!'
693    endif
694endif
695
696END SUBROUTINE make_layering
697
698!=======================================================================
699! To sublimate CO2 ice in the layering
700SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag)
701
702implicit none
703
704!---- Arguments
705type(layering),         intent(inout) :: this
706type(stratum), pointer, intent(inout) :: str
707logical,                intent(inout) :: new_lag
708real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
709
710!---- Local variables
711real                   :: r_subl, hsubl_dust, h_lag
712type(stratum), pointer :: current
713
714!---- Code
715! How much dust does CO2 ice sublimation involve to create the lag?
716r_subl = h2subl/(1. - co2ice_porosity) &
717        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
718hsubl_dust = r_subl*str%dust_volfrac*str%thickness
719
720str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust
721if (str%thickness < tol) then
722    ! Remove the sublimating stratum if there is no ice anymore
723    call rm_stratum(this,str)
724else
725    ! Update of properties in the sublimating stratum
726    str%top_elevation = str%top_elevation - h2subl/(1. - co2ice_porosity) - hsubl_dust
727    str%co2ice_volfrac = (h_co2ice_old - h2subl)/str%thickness
728    str%h2oice_volfrac = h_h2oice_old/str%thickness
729    str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
730    str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
731    ! Correct the value of top_elevation for the upper strata
732    current => str%up
733    do while (associated(current))
734        current%top_elevation = current%down%top_elevation + current%thickness
735        current => current%up
736    enddo
737endif
738
739! New stratum = dust lag
740h_lag = hsubl_dust/(1. - dry_lag_porosity)
741if (h_lag > 0.) then ! Only if the dust lag is building up
742    if (new_lag) then
743        call insert_stratum(this,str,h_lag,str%top_elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
744        new_lag = .false.
745    else
746        str%up%thickness = str%up%thickness + h_lag
747        str%up%top_elevation = str%up%top_elevation + h_lag
748    endif
749endif
750
751END SUBROUTINE subl_co2ice_layering
752
753!=======================================================================
754! To sublimate H2O ice in the layering
755SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag)
756
757implicit none
758
759!---- Arguments
760type(layering),         intent(inout) :: this
761type(stratum), pointer, intent(inout) :: str
762logical,                intent(inout) :: new_lag
763real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
764
765!---- Local variables
766real                   :: r_subl, hsubl_dust, h_lag
767type(stratum), pointer :: current
768
769!---- Code
770! How much dust does CO2 ice sublimation involve to create the lag?
771r_subl = h2subl/(1. - h2oice_porosity) &
772        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
773hsubl_dust = r_subl*str%dust_volfrac*str%thickness
774
775str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust
776if (str%thickness < tol) then
777    ! Remove the sublimating stratum if there is no ice anymore
778    call rm_stratum(this,str)
779else
780    ! Update of properties in the sublimating stratum
781    str%top_elevation = str%top_elevation - h2subl/(1. - h2oice_porosity) - hsubl_dust
782    str%co2ice_volfrac = h_co2ice_old/str%thickness
783    str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
784    str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
785    str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
786    ! Correct the value of top_elevation for the upper strata
787    current => str%up
788    do while (associated(current))
789        current%top_elevation = current%down%top_elevation + current%thickness
790        current => current%up
791    enddo
792endif
793
794! New stratum = dust lag
795h_lag = hsubl_dust/(1. - dry_lag_porosity)
796if (h_lag > 0.) then ! Only if the dust lag is building up
797    if (new_lag) then
798        call insert_stratum(this,str,h_lag,str%top_elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
799        new_lag = .false.
800    else
801        str%up%thickness = str%up%thickness + h_lag
802        str%up%top_elevation = str%up%top_elevation + h_lag
803    endif
804endif
805
806END SUBROUTINE subl_h2oice_layering
807
808!=======================================================================
809! To lift dust in the layering
810SUBROUTINE lift_dust_lags(this,str,h2lift)
811
812implicit none
813
814!---- Arguments
815type(layering),         intent(inout) :: this
816type(stratum), pointer, intent(inout) :: str
817real, intent(in) :: h2lift
818
819!---- Code
820! Update of properties in the eroding dust lag
821str%thickness = str%thickness - h2lift/(1. - str%air_volfrac)
822str%top_elevation = str%top_elevation - h2lift/(1. - str%air_volfrac)
823! Remove the eroding dust lag if there is no dust anymore
824if (str%thickness < tol) call rm_stratum(this,str)
825
826END SUBROUTINE lift_dust_lags
827
828END MODULE layering_mod
Note: See TracBrowser for help on using the repository browser.