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

Last change on this file since 3831 was 3821, checked in by jbclement, 7 days ago

PEM:
Addition of the flag 'impose_dust_ratio' (default: .false.) to impose a dust-to-ice ratio instead of a dust tendency for the layering algorithm. The ratio value can be set in the "run_PEM.def" by 'dust_to_ice_ratio' (default: 0.01).
JBC

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