Ignore:
Timestamp:
May 21, 2025, 3:57:33 PM (4 weeks ago)
Author:
jbclement
Message:

PEM:
Revision of the layering structure and algorithm:

  • 'stratum' components are now expressed in height
  • deletion of the (redundant) thickness feature of 'stratum'
  • 'stratif' in the PEM main program is renamed 'deposits'
  • addition of several procedures to get useful information about a stratum (major component or thickness)
  • all subsurface layers are now represented in the layering structure by strata with negative top elevation
  • simplification of the different situations arising from the input tendencies
  • porosity is determined for the entire stratum (and not anymore for each component) based on dominant component
  • sublimation of CO2 and H2O ice is now handled simultaneously (more realistic) in a stratum
  • linking the layering algorithm with the PEM initilization/finalization regarding PCM data and with the PEM stopping criteria
  • making separate cases for glaciers vs layering management
  • H2O sublimation flux correction is now handled with the layering when a dust lag layer layer is growing
  • update of 'h2o_ice_depth' and 'zdqsdif' accordingly at the PEM end for the PCM

JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90

    r3704 r3770  
    1919
    2020! Physical parameters
    21 real, parameter :: d_dust = 5.78e-2             ! Tendency of dust [kg.m-2.y-1]
    22 real, parameter :: dry_lag_porosity = 0.4       ! Porosity of dust lag
     21real, parameter :: d_dust = 5.78e-2            ! Tendency of dust [kg.m-2.y-1]
     22real, parameter :: dry_lag_porosity = 0.4      ! Porosity of dust lag
    2323real, parameter :: dry_regolith_porosity = 0.4 ! Porosity of regolith
    24 real, parameter :: co2ice_porosity = 0        ! Porosity of CO2 ice
    25 real, parameter :: h2oice_porosity = 0        ! Porosity of H2O ice
    26 real, parameter :: rho_dust = 2500.             ! Density of dust [kg.m-3]
    27 real, parameter :: rho_rock = 3200.             ! Density of rock [kg.m-3]
     24real, parameter :: breccia_porosity = 0.1      ! Porosity of breccia
     25real, parameter :: co2ice_porosity = 0.        ! Porosity of CO2 ice
     26real, parameter :: h2oice_porosity = 0.        ! Porosity of H2O ice
     27real, parameter :: rho_dust = 2500.            ! Density of dust [kg.m-3]
     28real, parameter :: rho_rock = 3200.            ! Density of rock [kg.m-3]
    2829
    2930! Lag layer parameters -> see Levrard et al. 2007
    30 real, parameter :: hmin_lag = 0  ! Minimal height of the lag deposit to reduce the sublimation rate
    31 real, parameter :: fred_subl = 1 ! Reduction factor of sublimation rate
     31real, parameter :: hmin_lag = 0.5  ! Minimal height of the lag deposit to reduce the sublimation rate
     32real, parameter :: fred_subl = 0.1 ! Reduction factor of sublimation rate
    3233
    3334! Stratum = node of the linked list
    3435type :: stratum
    35     real                   :: thickness       ! Layer thickness [m]
    3636    real                   :: top_elevation   ! Layer top_elevation (top height from the surface) [m]
    37     real                   :: co2ice_volfrac  ! CO2 ice volumetric fraction
    38     real                   :: h2oice_volfrac  ! H2O ice volumetric fraction
    39     real                   :: dust_volfrac    ! Dust volumetric fraction
    40     real                   :: pore_volfrac    ! Pore volumetric fraction = pores
     37    real                   :: h_co2ice        ! Height of CO2 ice [m]
     38    real                   :: h_h2oice        ! Height of H2O ice [m]
     39    real                   :: h_dust          ! Height of dust [m]
     40    real                   :: h_pore          ! Height of pores [m]
    4141    real                   :: poreice_volfrac ! Pore-filled ice volumetric fraction
    4242    type(stratum), pointer :: up => null()    ! Upper stratum (next node)
    … …  
    7070!     > stratif2array
    7171!     > array2stratif
    72 !     > print_stratif
     72!     > print_deposits
     73! Procedures to get information about a stratum:
     74!     > thickness
     75!     > is_co2ice_str
     76!     > is_h2oice_str
     77!     > is_dust_str
     78!     > is_dust_lag
     79!     > compute_h_toplag
    7380! Procedures for the algorithm to evolve the layering:
    7481!     > make_layering
    75 !     > subl_co2ice_layering
    76 !     > subl_h2oice_layering
    77 !     > lift_dust_lags
     82!     > lift_dust_lag
     83!     > subl_ice_str
    7884!=======================================================================
    7985! To display a stratum
    … …  
    8692
    8793!---- Code
    88 write(*,'(a,es13.6)') 'Thickness              : ', this%thickness
    89 write(*,'(a,es13.6)') 'Top elevation          : ', this%top_elevation
    90 write(*,'(a,es13.6)') 'CO2 ice (m3/m3)        : ', this%co2ice_volfrac
    91 write(*,'(a,es13.6)') 'H2O ice (m3/m3)        : ', this%h2oice_volfrac
    92 write(*,'(a,es13.6)') 'Dust (m3/m3)           : ', this%dust_volfrac
    93 write(*,'(a,es13.6)') 'Porosity (m3/m3)       : ', this%pore_volfrac
    94 write(*,'(a,es13.6)') 'Pore-filled ice (m3/m3): ', this%poreice_volfrac
     94write(*,'(a,es13.6)') 'Top elevation       [m]: ', this%top_elevation
     95write(*,'(a,es13.6)') 'Height of CO2 ice   [m]: ', this%h_co2ice
     96write(*,'(a,es13.6)') 'Height of H2O ice   [m]: ', this%h_h2oice
     97write(*,'(a,es13.6)') 'Height of dust      [m]: ', this%h_dust
     98write(*,'(a,es13.6)') 'Height of pores     [m]: ', this%h_pore
     99write(*,'(a,es13.6)') 'Pore-filled ice [m3/m3]: ', this%poreice_volfrac
    95100
    96101END SUBROUTINE print_stratum
    … …  
    98103!=======================================================================
    99104! To add a stratum to the top of a layering
    100 SUBROUTINE add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice)
     105SUBROUTINE add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice)
    101106
    102107implicit none
    … …  
    104109!---- Arguments
    105110type(layering), intent(inout)  :: this
    106 real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, pore, poreice
     111real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice
    107112
    108113!---- Local variables
    … …  
    113118allocate(str)
    114119nullify(str%up,str%down)
    115 str%thickness = 1.
    116120str%top_elevation = 1.
    117 str%co2ice_volfrac = 0.
    118 str%h2oice_volfrac = 0.
    119 str%dust_volfrac = 0.
    120 str%pore_volfrac = 1.
     121str%h_co2ice = 0.
     122str%h_h2oice = 0.
     123str%h_dust = 1.
     124str%h_pore = 0.
    121125str%poreice_volfrac = 0.
    122 if (present(thickness)) str%thickness = thickness
    123126if (present(top_elevation)) str%top_elevation = top_elevation
    124 if (present(co2ice)) str%co2ice_volfrac = co2ice
    125 if (present(h2oice)) str%h2oice_volfrac = h2oice
    126 if (present(dust)) str%dust_volfrac = dust
    127 if (present(pore)) str%pore_volfrac = pore
     127if (present(co2ice)) str%h_co2ice = co2ice
     128if (present(h2oice)) str%h_h2oice = h2oice
     129if (present(dust)) str%h_dust = dust
     130if (present(pore)) str%h_pore = pore
    128131if (present(poreice)) str%poreice_volfrac = poreice
    129 
    130 ! Verification of volume fraction
    131 if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%pore_volfrac)) > tol) then
    132     call print_stratum(str)
    133     error stop 'add_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!'
    134 endif
    135132
    136133! Increment the number of strata
    … …  
    151148!=======================================================================
    152149! To insert a stratum after another one in a layering
    153 SUBROUTINE insert_stratum(this,str_prev,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice)
     150SUBROUTINE insert_stratum(this,str_prev,top_elevation,co2ice,h2oice,dust,pore,poreice)
    154151
    155152implicit none
    … …  
    158155type(layering),         intent(inout) :: this
    159156type(stratum), pointer, intent(inout) :: str_prev
    160 real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, pore, poreice
     157real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice
    161158
    162159!---- Local variables
    … …  
    165162!---- Code
    166163if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering
    167     call add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice)
     164    call add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice)
    168165else
    169166    ! Creation of the stratum
    170167    allocate(str)
    171168    nullify(str%up,str%down)
    172     str%thickness = 1.
    173169    str%top_elevation = 1.
    174     str%co2ice_volfrac = 0.
    175     str%h2oice_volfrac = 0.
    176     str%dust_volfrac = 0.
    177     str%pore_volfrac = 1.
     170    str%h_co2ice = 0.
     171    str%h_h2oice = 0.
     172    str%h_dust = 1.
     173    str%h_pore = 0.
    178174    str%poreice_volfrac = 0.
    179     if (present(thickness)) str%thickness = thickness
    180175    if (present(top_elevation)) str%top_elevation = top_elevation
    181     if (present(co2ice)) str%co2ice_volfrac = co2ice
    182     if (present(h2oice)) str%h2oice_volfrac = h2oice
    183     if (present(dust)) str%dust_volfrac = dust
    184     if (present(pore)) str%pore_volfrac = pore
     176    if (present(co2ice)) str%h_co2ice = co2ice
     177    if (present(h2oice)) str%h_h2oice = h2oice
     178    if (present(dust)) str%h_dust = dust
     179    if (present(pore)) str%h_pore = pore
    185180    if (present(poreice)) str%poreice_volfrac = poreice
    186 
    187     ! Verification of volume fraction
    188     if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%pore_volfrac)) > tol) then
    189         call print_stratum(str)
    190         error stop 'add_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!'
    191     endif
    192181
    193182    ! Increment the number of strata
    … …  
    203192    current => str%up
    204193    do while (associated(current))
    205         current%top_elevation = current%down%top_elevation + current%thickness
     194        current%top_elevation = current%down%top_elevation + thickness(current)
    206195        current => current%up
    207196    enddo
    … …  
    224213
    225214!---- Code
    226 ! Protect the 3 sub-surface strata from removing
     215! Protect the sub-surface strata from removing
    227216if (str%top_elevation <= 0.) return
    228217
    … …  
    274263SUBROUTINE ini_layering(this)
    275264
     265use comsoil_h_PEM, only: soil_pem, nsoilmx_PEM, layer_PEM, index_breccia, index_bedrock
     266
    276267implicit none
    277268
    … …  
    279270type(layering), intent(inout) :: this
    280271
    281 !---- Code
    282 ! Creation of three strata at the bottom of the layering
    283 ! 1) Dry porous deep regolith
    284 !call add_stratum(this,18.,-12.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
    285 ! 2) Ice-cemented regolith (ice table)
    286 !call add_stratum(this,10.,-2.,0.,1. - dry_regolith_porosity,dry_regolith_porosity,0.)
    287 ! 3) Dry porous regolith
    288 call add_stratum(this,2.,0.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
     272!---- Local variables
     273integer :: i
     274real    :: h_soil, h_pore, h_tot
     275
     276!---- Code
     277! Creation of strata at the bottom of the layering to describe the sub-surface
     278if (soil_pem) then
     279    do i = nsoilmx_PEM,index_bedrock,-1
     280        h_soil = layer_PEM(i) - layer_PEM(i - 1) ! No porosity
     281        call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,0.,0.)
     282    enddo
     283    do i = index_bedrock - 1,index_breccia,-1
     284        h_tot = layer_PEM(i) - layer_PEM(i - 1)
     285        h_pore = h_tot*breccia_porosity
     286        h_soil = h_tot - h_pore
     287        call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.)
     288    enddo
     289    do i = index_breccia - 1,2,-1
     290        h_tot = layer_PEM(i) - layer_PEM(i - 1)
     291        h_pore = h_tot*dry_regolith_porosity
     292        h_soil = h_tot - h_pore
     293        call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.)
     294    enddo
     295    h_pore = layer_PEM(1)*dry_regolith_porosity
     296    h_soil = layer_PEM(1) - h_pore
     297    call add_stratum(this,0.,0.,0.,h_soil,h_pore,0.)
     298endif
    289299
    290300END SUBROUTINE ini_layering
    … …  
    387397        k = 1
    388398        do while (associated(current))
    389             stratif_array(ig,islope,k,1) = current%thickness
    390             stratif_array(ig,islope,k,2) = current%top_elevation
    391             stratif_array(ig,islope,k,3) = current%co2ice_volfrac
    392             stratif_array(ig,islope,k,4) = current%h2oice_volfrac
    393             stratif_array(ig,islope,k,5) = current%dust_volfrac
    394             stratif_array(ig,islope,k,6) = current%pore_volfrac
    395             stratif_array(ig,islope,k,7) = current%poreice_volfrac
     399            stratif_array(ig,islope,k,1) = current%top_elevation
     400            stratif_array(ig,islope,k,2) = current%h_co2ice
     401            stratif_array(ig,islope,k,3) = current%h_h2oice
     402            stratif_array(ig,islope,k,4) = current%h_dust
     403            stratif_array(ig,islope,k,5) = current%h_pore
     404            stratif_array(ig,islope,k,6) = current%poreice_volfrac
    396405            current => current%up
    397406            k = k + 1
    … …  
    419428do islope = 1,nslope
    420429    do ig = 1,ngrid
    421         stratif(ig,islope)%bottom%thickness       = stratif_array(ig,islope,1,1)
    422         stratif(ig,islope)%bottom%top_elevation   = stratif_array(ig,islope,1,2)
    423         stratif(ig,islope)%bottom%co2ice_volfrac  = stratif_array(ig,islope,1,3)
    424         stratif(ig,islope)%bottom%h2oice_volfrac  = stratif_array(ig,islope,1,4)
    425         stratif(ig,islope)%bottom%dust_volfrac    = stratif_array(ig,islope,1,5)
    426         stratif(ig,islope)%bottom%pore_volfrac    = stratif_array(ig,islope,1,6)
    427         stratif(ig,islope)%bottom%poreice_volfrac = stratif_array(ig,islope,1,7)
    428         do k = 2,size(stratif_array,3)
    429             if (all(abs(stratif_array(ig,islope,k,:)) < eps)) exit
    430             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),stratif_array(ig,islope,k,7))
     430        do k = 1,size(stratif_array,3)
     431            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))
    431432        enddo
    432433    enddo
    … …  
    436437
    437438!=======================================================================
    438 ! To display a stratification (layerings)
    439 SUBROUTINE print_stratif(this,ngrid,nslope)
     439! To display deposits
     440SUBROUTINE print_deposits(this,ngrid,nslope)
    440441
    441442implicit none
    … …  
    459460enddo
    460461
    461 END SUBROUTINE print_stratif
    462 
     462END SUBROUTINE print_deposits
     463
     464!=======================================================================
     465!=======================================================================
     466! To compute the thickness of a stratum
     467FUNCTION thickness(this) RESULT(h_str)
     468
     469implicit none
     470
     471!---- Arguments
     472type(stratum), intent(in) :: this
     473real :: h_str
     474
     475!---- Code
     476h_str = this%h_co2ice + this%h_h2oice + this%h_dust + this%h_pore
     477
     478END FUNCTION thickness
     479
     480!=======================================================================
     481! To know whether the stratum can be considered as a CO2 ice layer
     482FUNCTION is_co2ice_str(str) RESULT(co2ice_str)
     483
     484implicit none
     485
     486!---- Arguments
     487type(stratum), intent(in) :: str
     488logical :: co2ice_str
     489
     490!---- Code
     491co2ice_str = .false.
     492if (str%h_co2ice > str%h_h2oice .and. str%h_co2ice > str%h_dust) co2ice_str = .true.
     493
     494END FUNCTION is_co2ice_str
     495
     496!=======================================================================
     497! To know whether the stratum can be considered as a H2O ice layer
     498FUNCTION is_h2oice_str(str) RESULT(h2oice_str)
     499
     500implicit none
     501
     502!---- Arguments
     503type(stratum), intent(in) :: str
     504logical :: h2oice_str
     505
     506!---- Code
     507h2oice_str = .false.
     508if (str%h_h2oice > str%h_co2ice .and. str%h_h2oice > str%h_dust) h2oice_str = .true.
     509
     510END FUNCTION is_h2oice_str
     511
     512!=======================================================================
     513! To know whether the stratum can be considered as a dust layer
     514FUNCTION is_dust_str(str) RESULT(dust_str)
     515
     516implicit none
     517
     518!---- Arguments
     519type(stratum), intent(in) :: str
     520logical :: dust_str
     521
     522!---- Code
     523dust_str = .false.
     524if (str%h_dust > str%h_co2ice .and. str%h_dust > str%h_h2oice) dust_str = .true.
     525
     526END FUNCTION is_dust_str
     527
     528!=======================================================================
     529! To know whether the stratum can be considered as a dust lag
     530FUNCTION is_dust_lag(str) RESULT(dust_lag)
     531
     532implicit none
     533
     534!---- Arguments
     535type(stratum), intent(in) :: str
     536logical :: dust_lag
     537
     538!---- Code
     539dust_lag = .false.
     540if (str%h_co2ice < tol .and. str%h_h2oice < tol .and. str%poreice_volfrac < tol) dust_lag = .true.
     541
     542END FUNCTION is_dust_lag
     543
     544!=======================================================================
     545! To get the top lag height
     546FUNCTION thickness_toplag(this) RESULT(h_toplag)
     547
     548implicit none
     549
     550!---- Arguments
     551type(layering), intent(in) :: this
     552real :: h_toplag
     553
     554!---- Local variables
     555type(stratum), pointer :: current
     556
     557!---- Code
     558h_toplag = 0.
     559current => this%top
     560! Is the considered stratum a dust lag?
     561do while (current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < tol .and. associated(current))
     562    h_toplag = h_toplag + thickness(current)
     563    current => current%down
     564enddo
     565
     566END FUNCTION thickness_toplag
     567
     568!=======================================================================
    463569!=======================================================================
    464570! Layering algorithm
    465 SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,zshift_surf,new_lag1,new_lag2,zlag,stopPEM,current1,current2)
     571SUBROUTINE make_layering(this,d_co2ice,d_h2oice,new_str,zshift_surf,new_lag,zlag,current)
    466572
    467573implicit none
    … …  
    469575!---- Arguments
    470576! Inputs
    471 real, intent(in) :: d_co2ice, d_h2oice, d_dust
    472577
    473578! Outputs
    474579type(layering),         intent(inout) :: this
    475 type(stratum), pointer, intent(inout) :: current1, current2
    476 logical,                intent(inout) :: new_str, new_lag1, new_lag2
    477 integer,                intent(inout) :: stopPEM
     580type(stratum), pointer, intent(inout) :: current
     581real,                   intent(inout) :: d_co2ice, d_h2oice
     582logical,                intent(inout) :: new_str, new_lag
    478583real, intent(out) :: zshift_surf, zlag
    479584
    … …  
    481586real                   :: dh_co2ice, dh_h2oice, dh_dust
    482587real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
    483 real                   :: thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
    484 real                   :: h_toplag
     588real                   :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
     589real                   :: h2lift, h2lift_tot, h_pore, h_tot, h_toplag
    485590type(stratum), pointer :: currentloc
    486591
    … …  
    489594dh_h2oice = d_h2oice/rho_h2oice
    490595dh_dust = d_dust/rho_dust
    491 zshift_surf = 0.
     596zshift_surf = this%top%top_elevation
    492597zlag = 0.
    493598
    494 if (dh_dust < 0.) then ! Dust lifting only
    495     if (abs(dh_co2ice) > eps .or. abs(dh_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
    496     write(*,'(a)') ' Stratification -> Dust lifting'
     599! DUST SEDIMENTATION with possibly ice condensation
     600if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice .and. dh_co2ice >= 0. .and. dh_h2oice >= 0.) then
     601    write(*,*) '> Layering: dust sedimentation'
     602    h_pore = dh_dust*dry_regolith_porosity/(1. - dry_regolith_porosity)
     603    h_tot = dh_co2ice + dh_h2oice + dh_dust + h_pore
     604    ! New stratum at the top of the layering
     605    if (new_str) then
     606        call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,0.)
     607        new_str = .false.
     608    else
     609        this%top%top_elevation = this%top%top_elevation + h_tot
     610        this%top%h_co2ice = this%top%h_co2ice + dh_co2ice
     611        this%top%h_h2oice = this%top%h_h2oice + dh_h2oice
     612        this%top%h_dust = this%top%h_dust + dh_dust
     613        this%top%h_pore = this%top%h_pore + h_pore
     614    endif
     615!-----------------------------------------------------------------------
     616! DUST LIFTING
     617!(with possibly ice sublimation???)
     618!else if (dh_dust < dh_co2ice .and. dh_dust < dh_h2oice .and. dh_co2ice <= 0. .and. dh_h2oice <= 0.) then
     619else if (dh_dust < 0. .and. abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then
     620    write(*,*) '> Layering: dust lifting'
    497621    h2lift_tot = abs(dh_dust)
    498     do while (h2lift_tot > 0. .and. associated(current1))
     622    do while (h2lift_tot > 0. .and. associated(current))
    499623        ! Is the considered stratum a dust lag?
    500         if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we can lift dust
    501             h_dust_old = current1%dust_volfrac*current1%thickness
    502             ! How much dust can be lifted in the considered dust lag?
    503             if (h2lift_tot - h_dust_old <= eps) then ! There is enough to lift
    504                 h2lift = h2lift_tot
    505                 h2lift_tot = 0.
    506                 call lift_dust_lags(this,current1,h2lift)
    507             else ! Only a fraction can be lifted and so we move to the underlying stratum
    508                 h2lift = h_dust_old
    509                 h2lift_tot = h2lift_tot - h2lift
    510                 call lift_dust_lags(this,current1,h2lift)
    511                 current1 => current1%down
    512             endif
    513         else ! No, we stop
    514             write(*,'(a)') ' Stratification -> Dust lifting: no dust lag!'
    515             stopPEM = 8
     624        !if ((current%h_co2ice > eps .and. abs(dh_co2ice) < tol) .or. (current%h_h2oice > eps .and. abs(dh_h2oice) < tol)) then ! No, there is ice cementing the dust
     625        if (current%h_co2ice > 0. .or. current%h_h2oice > 0. .or. current%poreice_volfrac > 0.) then ! No, there is ice cementing the dust
     626            dh_dust = 0. ! The tendency is set to 0
     627            exit
     628        else ! Yes, we can lift dust
     629            h2lift = min(h2lift_tot,current%h_dust) ! How much dust can be lifted?
     630            h2lift_tot = h2lift_tot - h2lift
     631            call lift_dust_lag(this,current,h2lift)
     632            if (h2lift_tot > 0.) current => current%down ! Still dust to be lifted so we move to the underlying stratum
     633        endif
     634    enddo
     635    if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
     636!-----------------------------------------------------------------------
     637! ICE CONDENSATION
     638else if ((dh_co2ice > dh_dust .or. dh_h2oice > dh_dust) .and. dh_dust >= 0.) then
     639    write(*,*) '> Layering: ice condensation'
     640    ! Which porosity is considered?
     641    if (dh_co2ice > dh_h2oice) then ! CO2 ice is dominant
     642        h_pore = dh_co2ice*co2ice_porosity/(1. - co2ice_porosity)
     643    else ! H2O ice is dominant
     644        h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity)
     645    endif
     646    h_tot = dh_co2ice + dh_h2oice + dh_dust + h_pore
     647    ! New stratum at the top of the layering
     648    if (new_str) then
     649        call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,0.)
     650        new_str = .false.
     651    else
     652        this%top%top_elevation = this%top%top_elevation + h_tot
     653        this%top%h_co2ice = this%top%h_co2ice + dh_co2ice
     654        this%top%h_h2oice = this%top%h_h2oice + dh_h2oice
     655        this%top%h_dust = this%top%h_dust + dh_dust
     656        this%top%h_pore = this%top%h_pore + h_pore
     657    endif
     658!-----------------------------------------------------------------------
     659! ICE SUBLIMATION
     660else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. abs(dh_dust) < tol) then
     661    write(*,*) '> Layering: ice sublimation'
     662    h2subl_co2ice_tot = abs(dh_co2ice)
     663    h2subl_h2oice_tot = abs(dh_h2oice)
     664    do while (h2subl_co2ice_tot > 0. .and. h2subl_h2oice_tot > 0. .and. associated(current))
     665        h_co2ice_old = current%h_co2ice
     666        h_h2oice_old = current%h_h2oice
     667        h_dust_old = current%h_dust
     668
     669!~         ! How much is CO2 ice sublimation inhibited by the top dust lag?
     670!~         h_toplag = 0.
     671!~         if (associated(current%up)) then ! If there is an upper stratum
     672!~             currentloc => current%up
     673!~             ! Is the considered stratum a dust lag?
     674!~             do while (.not. (current%h_co2ice > 0. .or. current%h_h2oice > 0. .or. current%poreice_volfrac > 0.) .and. associated(currentloc%up))
     675!~                 h_toplag = h_toplag + thickness(currentloc)
     676!~                 currentloc => currentloc%up
     677!~             enddo
     678!~         endif
     679!~         h2subl_co2ice_tot = h2subl_co2ice_tot*fred_subl**(h_toplag/hmin_lag)
     680
     681        ! Is there CO2 ice to sublimate?
     682        h2subl_co2ice = 0.
     683        if (h_co2ice_old > 0. .and. h2subl_co2ice_tot > 0.) then
     684            h2subl_co2ice = min(h2subl_co2ice_tot,h_co2ice_old)
     685            h2subl_co2ice_tot = h2subl_co2ice_tot - h2subl_co2ice
     686        endif
     687        ! Is there H2O ice to sublimate?
     688        h2subl_h2oice = 0.
     689        if (h_h2oice_old > 0. .and. h2subl_h2oice_tot > 0.) then
     690            h2subl_h2oice = min(h2subl_h2oice_tot,h_h2oice_old)
     691            h2subl_h2oice_tot = h2subl_h2oice_tot - h2subl_h2oice
     692        endif
     693        ! Sublimation
     694        if (h2subl_co2ice > 0. .or. h2subl_h2oice >0.) call subl_ice_str(this,current,h2subl_co2ice,h2subl_h2oice,zlag,new_lag)
     695        ! Still ice to be sublimated and no more ice in the current stratum so we move to the underlying stratum
     696        if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0.) .and. current%h_co2ice < tol .and. current%h_h2oice < tol) then
     697            current => current%down
     698            new_lag = .true.
     699        else
    516700            exit
    517701        endif
    518702    enddo
    519     if (h2lift_tot > 0.) then
    520         write(*,'(a)') ' Stratification -> Not enough dust in the lag layers to complete the lifting!'
    521         stopPEM = 8
    522     endif
    523     zshift_surf = dh_dust + h2lift_tot
    524 
     703    if (h2subl_co2ice_tot > 0.) dh_co2ice = 0. ! No CO2 ice available anymore so the tendency is set to 0
     704    if (h2subl_h2oice_tot > 0.) dh_h2oice = 0. ! No H2O ice available anymore so the tendency is set to 0
     705!-----------------------------------------------------------------------
     706!~ ! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION
     707!~ else if (dh_co2ice < 0. .and. dh_h2oice > 0. .and. (abs(dh_dust) < abs(dh_co2ice) .or. abs(dh_dust) < abs(dh_h2oice))) then
     708
     709! TO DO????
     710
     711!-----------------------------------------------------------------------
     712! NOTHING HAPPENS
     713else if (abs(dh_co2ice) < tol .and. abs(dh_h2oice) < tol .and. abs(dh_dust) < tol) then
     714    ! We do nothing
     715!-----------------------------------------------------------------------
     716! UNUSUAL (WEIRD) SITUATION
    525717else
    526 
    527 !------ Dust sedimentation only
    528     if (abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then
    529         write(*,'(a)') ' Stratification -> Dust sedimentation'
    530         ! New stratum at the layering top by sedimentation of dust
    531         thickness = dh_dust/(1. - dry_regolith_porosity)
    532         if (thickness > 0.) then ! Only if the stratum is building up
    533              if (new_str) then
    534                  call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity,0.)
    535                  new_str = .false.
    536              else
    537                  this%top%thickness = this%top%thickness + thickness
    538                  this%top%top_elevation = this%top%top_elevation + thickness
    539              endif
    540         endif
    541         zshift_surf = thickness
    542 
    543 !------ Condensation of CO2 ice + H2O ice
    544     else if ((dh_co2ice >= 0. .and. dh_h2oice > 0.) .or. (dh_co2ice > 0. .and. dh_h2oice >= 0.)) then
    545         write(*,'(a)') ' Stratification -> CO2 and H2O ice condensation'
    546         ! New stratum at the layering top by condensation of CO2 and H2O ice
    547         thickness = dh_co2ice/(1. - co2ice_porosity) + dh_h2oice/(1. - h2oice_porosity) + dh_dust
    548         if (thickness > 0.) then ! Only if the stratum is building up
    549              if (new_str) then
    550                  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),0.)
    551                  new_str = .false.
    552              else
    553                  this%top%thickness = this%top%thickness + thickness
    554                  this%top%top_elevation = this%top%top_elevation + thickness
    555              endif
    556         endif
    557         zshift_surf = thickness
    558 
    559 !------ Sublimation of CO2 ice + Condensation of H2O ice
    560     else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then
    561         write(*,'(a)') ' Stratification -> Sublimation of CO2 ice + Condensation of H2O ice'
    562         ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
    563         h2subl_tot = abs(dh_co2ice)
    564         do while (h2subl_tot > 0. .and. associated(current1))
    565             h_co2ice_old = current1%co2ice_volfrac*current1%thickness
    566             h_h2oice_old = current1%h2oice_volfrac*current1%thickness
    567             h_dust_old = current1%dust_volfrac*current1%thickness
    568 
    569             ! How much is CO2 ice sublimation inhibited by the top dust lag?
    570             h_toplag = 0.
    571             if (associated(current1%up)) then ! If there is an upper stratum
    572                 currentloc => current1%up
    573                 do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol)
    574                     h_toplag = h_toplag + currentloc%thickness
    575                     if (.not. associated(currentloc%up)) exit
    576                     currentloc => currentloc%up
    577                 enddo
    578             endif
    579             h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
    580 
    581             ! How much CO2 ice can sublimate in the considered stratum?
    582             if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
    583                 h2subl = h2subl_tot
    584                 h2subl_tot = 0.
    585                 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag)
    586             else if (h_co2ice_old < eps) then ! There is nothing to sublimate
    587                 ! Is the considered stratum a dust lag?
    588                 if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum
    589                     current1 => current1%down
    590                     new_lag1 = .true.
    591                 else ! No, so we stop here
    592                     exit
    593                 endif
    594             else ! Only a fraction can sublimate and so we move to the underlying stratum
    595                 h2subl = h_co2ice_old
    596                 h2subl_tot = h2subl_tot - h2subl
    597                 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag)
    598                 current1 => current1%down
    599                 new_lag1 = .true.
    600             endif
    601         enddo
    602         if (h2subl_tot > 0.) then
    603             write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
    604             stopPEM = 8
    605         endif
    606         ! New stratum at the layering top by condensation of H2O ice
    607         thickness = dh_h2oice/(1. - h2oice_porosity) + dh_dust
    608         if (thickness > 0.) then ! Only if the stratum is building up
    609              if (new_str) then
    610                  call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness),0.)
    611                  new_str = .false.
    612              else
    613                  this%top%thickness = this%top%thickness + thickness
    614                  this%top%top_elevation = this%top%top_elevation + thickness
    615              endif
    616         endif
    617         zshift_surf = zshift_surf + thickness
    618 
    619 !------ Sublimation of CO2 ice + H2O ice
    620     else if ((dh_co2ice <= 0. .and. dh_h2oice < 0.) .or. (dh_co2ice < 0. .and. dh_h2oice <= 0.)) then
    621         write(*,'(a)') ' Stratification -> Sublimation of CO2 and H2O ice'
    622         ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
    623         h2subl_tot = abs(dh_co2ice)
    624         do while (h2subl_tot > 0. .and. associated(current1))
    625             h_co2ice_old = current1%co2ice_volfrac*current1%thickness
    626             h_h2oice_old = current1%h2oice_volfrac*current1%thickness
    627             h_dust_old = current1%dust_volfrac*current1%thickness
    628 
    629             ! How much is CO2 ice sublimation inhibited by the top dust lag?
    630             h_toplag = 0.
    631             if (associated(current1%up)) then ! If there is an upper stratum
    632                 currentloc => current1%up
    633                 do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol)
    634                     h_toplag = h_toplag + currentloc%thickness
    635                     if (.not. associated(currentloc%up)) exit
    636                     currentloc => currentloc%up
    637                 enddo
    638             endif
    639             h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
    640 
    641             ! How much CO2 ice can sublimate in the considered stratum?
    642             if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
    643                 h2subl = h2subl_tot
    644                 h2subl_tot = 0.
    645                 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag)
    646             else if (h_co2ice_old < eps) then ! There is nothing to sublimate
    647                 ! Is the considered stratum a dust lag?
    648                 if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum
    649                     current1 => current1%down
    650                     new_lag1 = .true.
    651                 else ! No, so we stop here
    652                     exit
    653                 endif
    654             else ! Only a fraction can sublimate and so we move to the underlying stratum
    655                 h2subl = h_co2ice_old
    656                 h2subl_tot = h2subl_tot - h2subl
    657                 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag)
    658                 current1 => current1%down
    659                 new_lag1 = .true.
    660             endif
    661         enddo
    662         if (h2subl_tot > 0.) then
    663             write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
    664             stopPEM = 8
    665         endif
    666         ! H2O ice sublimation in the considered stratum + New stratum for dust lag
    667         h2subl_tot = abs(dh_h2oice)
    668         do while (h2subl_tot > 0. .and. associated(current2))
    669             h_co2ice_old = current2%co2ice_volfrac*current2%thickness
    670             h_h2oice_old = current2%h2oice_volfrac*current2%thickness
    671             h_dust_old = current2%dust_volfrac*current2%thickness
    672 
    673             ! How much is H2O ice sublimation inhibited by the top dust lag?
    674             h_toplag = 0.
    675             if (associated(current2%up)) then ! If there is an upper stratum
    676                 currentloc => current2%up
    677                 do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol)
    678                     h_toplag = h_toplag + currentloc%thickness
    679                     if (.not. associated(currentloc%up)) exit
    680                     currentloc => currentloc%up
    681                 enddo
    682             endif
    683             h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
    684 
    685             ! How much H2O ice can sublimate in the considered stratum?
    686             if (h2subl_tot - h_h2oice_old <= eps) then ! There is enough to sublimate
    687                 h2subl = h2subl_tot
    688                 h2subl_tot = 0.
    689                 call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag2,zlag)
    690             else if (h_h2oice_old < eps) then ! There is nothing to sublimate
    691                 ! Is the considered stratum a dust lag?
    692                 if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum
    693                     current1 => current1%down
    694                     new_lag1 = .true.
    695                 else ! No, so we stop here
    696                     exit
    697                 endif
    698             else ! Only a fraction can sublimate and so we move to the underlying stratum
    699                 h2subl = h_h2oice_old
    700                 h2subl_tot = h2subl_tot - h2subl
    701                 call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag2,zlag)
    702                 current2 => current2%down
    703                 new_lag2 = .true.
    704             endif
    705         enddo
    706         if (h2subl_tot > 0.) then
    707             write(*,'(a)') ' Stratification -> Not enough H2O ice in the layering to complete the sublimation!'
    708             stopPEM = 8
    709         endif
    710 
    711 !------ Condensation of CO2 ice + Sublimation of H2O ice
    712     else if (dh_co2ice > 0. .and. dh_h2oice < 0.) then
    713         error stop 'Impossible situation: CO2 ice condensation + H2O ice sublimation!'
    714     endif
     718    write(*,'(a)') 'Error: situation for the layering construction not managed!'
     719    write(*,'(a,es12.3)') '    > dh_co2ice [m/y] = ', dh_co2ice
     720    write(*,'(a,es12.3)') '    > dh_h2oice [m/y] = ', dh_h2oice
     721    write(*,'(a,es12.3)') '    > dh_dust   [m/y] = ', dh_dust, tol
     722!~     error stop
    715723endif
    716724
     725zshift_surf = this%top%top_elevation - zshift_surf
     726
    717727END SUBROUTINE make_layering
    718728
    719729!=======================================================================
    720 ! To sublimate CO2 ice in the layering
    721 SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag)
     730! To lift dust in a stratum
     731SUBROUTINE lift_dust_lag(this,str,h2lift)
     732
     733implicit none
     734
     735!---- Arguments
     736type(layering),         intent(inout) :: this
     737type(stratum), pointer, intent(inout) :: str
     738real, intent(in) :: h2lift
     739
     740!---- Code
     741! Update of properties in the eroding dust lag
     742str%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust)
     743str%h_pore = str%h_pore - h2lift*str%h_pore/str%h_dust
     744str%h_dust = str%h_dust - h2lift
     745
     746! Remove the eroding dust lag if there is no dust anymore
     747if (str%h_dust < tol) call rm_stratum(this,str)
     748
     749END SUBROUTINE lift_dust_lag
     750
     751!=======================================================================
     752! To sublimate ice in a stratum
     753SUBROUTINE subl_ice_str(this,str,h2subl_co2ice,h2subl_h2oice,zlag,new_lag)
    722754
    723755implicit none
    … …  
    727759type(stratum), pointer, intent(inout) :: str
    728760logical,                intent(inout) :: new_lag
    729 real,                   intent(inout) :: zshift_surf, zlag
    730 real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
    731 
    732 !---- Local variables
    733 real                   :: r_subl, hsubl_dust, h_lag
     761real,                   intent(inout) :: zlag
     762real, intent(in) :: h2subl_co2ice, h2subl_h2oice
     763
     764!---- Local variables
     765real                   :: hsubl_dust, h_pore, h_ice, h_tot, h2subl
    734766type(stratum), pointer :: current
    735767
    736768!---- Code
    737 ! How much dust does CO2 ice sublimation involve to create the lag?
    738 r_subl = h2subl/(1. - co2ice_porosity) &
    739         /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
    740 hsubl_dust = r_subl*str%dust_volfrac*str%thickness
    741 
    742 str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust
    743 zshift_surf = zshift_surf - h2subl/(1. - co2ice_porosity) - hsubl_dust
    744 if (str%thickness < tol) then
    745     ! Remove the sublimating stratum if there is no ice anymore
    746     call rm_stratum(this,str)
    747 else
    748     ! Update of properties in the sublimating stratum
    749     str%top_elevation = str%top_elevation - h2subl/(1. - co2ice_porosity) - hsubl_dust
    750     str%co2ice_volfrac = (h_co2ice_old - h2subl)/str%thickness
    751     str%h2oice_volfrac = h_h2oice_old/str%thickness
    752     str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
    753     str%pore_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
    754     str%poreice_volfrac = 0.
    755     ! Correct the value of top_elevation for the upper strata
    756     current => str%up
    757     do while (associated(current))
    758         current%top_elevation = current%down%top_elevation + current%thickness
    759         current => current%up
    760     enddo
    761 endif
    762 
    763 ! New stratum = dust lag
    764 h_lag = hsubl_dust/(1. - dry_lag_porosity)
    765 if (h_lag > 0.) then ! Only if the dust lag is building up
     769h_ice = str%h_co2ice + str%h_h2oice
     770h2subl = h2subl_co2ice + h2subl_h2oice
     771
     772! How much dust does ice sublimation involve to create the lag?
     773hsubl_dust = str%h_dust*h2subl/h_ice
     774
     775! Update of properties in the sublimating stratum
     776str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hsubl_dust
     777str%h_co2ice = str%h_co2ice - h2subl_co2ice
     778str%h_h2oice = str%h_h2oice - h2subl_h2oice
     779str%h_dust = str%h_dust - hsubl_dust
     780str%h_pore = str%h_pore - h2subl*h_pore/h_ice
     781
     782! Correct the value of top_elevation for the upper strata
     783current => str%up
     784do while (associated(current))
     785    current%top_elevation = current%down%top_elevation + thickness(current)
     786    current => current%up
     787enddo
     788
     789! Remove the sublimating stratum if there is no ice anymore
     790if (str%h_co2ice < tol .and. str%h_h2oice < tol) call rm_stratum(this,str)
     791
     792! New stratum above the current stratum as a dust lag
     793if (hsubl_dust > 0.) then ! Only if the dust lag is building up
     794    h_pore = hsubl_dust*dry_lag_porosity/(1. - dry_lag_porosity)
     795    h_tot = hsubl_dust + h_pore
    766796    if (new_lag) then
    767         call insert_stratum(this,str,h_lag,str%top_elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
     797        call insert_stratum(this,str,str%top_elevation + h_tot,0.,0.,hsubl_dust,h_pore)
    768798        new_lag = .false.
    769799    else
    770         str%up%thickness = str%up%thickness + h_lag
    771         str%up%top_elevation = str%up%top_elevation + h_lag
     800        str%up%top_elevation = str%up%top_elevation + h_tot
    772801    endif
    773802endif
    774 zlag = zlag + h_lag
    775 
    776 END SUBROUTINE subl_co2ice_layering
    777 
    778 !=======================================================================
    779 ! To sublimate H2O ice in the layering
    780 SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag)
    781 
    782 implicit none
    783 
    784 !---- Arguments
    785 type(layering),         intent(inout) :: this
    786 type(stratum), pointer, intent(inout) :: str
    787 logical,                intent(inout) :: new_lag
    788 real,                   intent(inout) :: zshift_surf, zlag
    789 real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
    790 
    791 !---- Local variables
    792 real                   :: r_subl, hsubl_dust, h_lag
    793 type(stratum), pointer :: current
    794 
    795 !---- Code
    796 ! How much dust does CO2 ice sublimation involve to create the lag?
    797 r_subl = h2subl/(1. - h2oice_porosity) &
    798         /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
    799 hsubl_dust = r_subl*str%dust_volfrac*str%thickness
    800 
    801 str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust
    802 zshift_surf = zshift_surf - h2subl/(1. - h2oice_porosity) - hsubl_dust
    803 if (str%thickness < tol) then
    804     ! Remove the sublimating stratum if there is no ice anymore
    805     call rm_stratum(this,str)
    806 else
    807     ! Update of properties in the sublimating stratum
    808     str%top_elevation = str%top_elevation - h2subl/(1. - h2oice_porosity) - hsubl_dust
    809     str%co2ice_volfrac = h_co2ice_old/str%thickness
    810     str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
    811     str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
    812     str%pore_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
    813     str%poreice_volfrac = 0.
    814     ! Correct the value of top_elevation for the upper strata
    815     current => str%up
    816     do while (associated(current))
    817         current%top_elevation = current%down%top_elevation + current%thickness
    818         current => current%up
    819     enddo
    820 endif
    821 
    822 ! New stratum = dust lag
    823 h_lag = hsubl_dust/(1. - dry_lag_porosity)
    824 if (h_lag > 0.) then ! Only if the dust lag is building up
    825     if (new_lag) then
    826         call insert_stratum(this,str,h_lag,str%top_elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
    827         new_lag = .false.
    828     else
    829         str%up%thickness = str%up%thickness + h_lag
    830         str%up%top_elevation = str%up%top_elevation + h_lag
    831     endif
    832 endif
    833 zlag = zlag + h_lag
    834 
    835 END SUBROUTINE subl_h2oice_layering
    836 
    837 !=======================================================================
    838 ! To lift dust in the layering
    839 SUBROUTINE lift_dust_lags(this,str,h2lift)
    840 
    841 implicit none
    842 
    843 !---- Arguments
    844 type(layering),         intent(inout) :: this
    845 type(stratum), pointer, intent(inout) :: str
    846 real, intent(in) :: h2lift
    847 
    848 !---- Code
    849 ! Update of properties in the eroding dust lag
    850 str%thickness = str%thickness - h2lift/(1. - str%pore_volfrac)
    851 str%top_elevation = str%top_elevation - h2lift/(1. - str%pore_volfrac)
    852 ! Remove the eroding dust lag if there is no dust anymore
    853 if (str%thickness < tol) call rm_stratum(this,str)
    854 
    855 END SUBROUTINE lift_dust_lags
     803zlag = zlag + h_tot
     804
     805END SUBROUTINE subl_ice_str
    856806
    857807END MODULE layering_mod
Note: See TracChangeset for help on using the changeset viewer.