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

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
7 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

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

    r3742 r3770  
    2121
    2222! III Output
    23 !    III_a Update surface value for the PCM start files
     23!    III_a Update surface values for the PCM start files
    2424!    III_b Write the "restart.nc" and "restartfi.nc"
    2525!    III_c Write the "restartpem.nc"
    … …  
    6969use writediagpem_mod,           only: writediagpem, writediagsoilpem
    7070use co2condens_mod,             only: CO2cond_ps
    71 use layering_mod,               only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
     71use layering_mod,               only: ptrarray, stratum, layering, del_layering, make_layering, get_nb_str_max, &
     72                                      nb_str_max, layering_algo, thickness_toplag, is_dust_lag, is_h2oice_str
    7273use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
    7374use version_info_mod,           only: print_version_info
    … …  
    119120include "comgeom.h"
    120121include "iniprint.h"
    121 include "callkeys.h"
    122 
    123 integer ngridmx
    124 parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
    125122
    126123! Same variable names as in the PCM
     124integer, parameter :: ngridmx = 2 + (jjm - 1)*iim - 1/jjm
    127125integer, parameter :: nlayer = llm ! Number of vertical layer
    128126integer            :: ngrid        ! Number of physical grid points
    … …  
    196194real,    dimension(:,:), allocatable :: q_co2_PEM_phys         ! Grid points x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
    197195
    198 ! Variables for stratification (layering) evolution
    199 type(layering), dimension(:,:), allocatable :: stratif                     ! Layering (linked list of "stratum") for each grid point and slope
    200 type(ptrarray), dimension(:,:), allocatable :: current1, current2          ! Current active stratum in the layering
    201 logical,        dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm
     196! Variables for the evolution of layered deposits
     197type(layering), dimension(:,:), allocatable :: deposits          ! Layering for each grid point and slope
     198type(ptrarray), dimension(:,:), allocatable :: current           ! Current active stratum in the layering
     199logical,        dimension(:,:), allocatable :: new_str, new_lag  ! Flags for the layering algorithm
     200real                                        :: h2o_ice_depth_old ! Old depth of subsurface ice layer
    202201
    203202! Variables for slopes
    … …  
    646645!------------------------
    647646write(*,*) '> Reading "startpem.nc"'
    648 allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),stratif(ngrid,nslope))
     647allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),deposits(ngrid,nslope))
    649648allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid))
    650 if (layering_algo) then
    651     do islope = 1,nslope
    652         do i = 1,ngrid
    653             call ini_layering(stratif(i,islope))
    654         enddo
    655     enddo
    656 endif
    657 
    658649delta_h2o_icetablesublim = 0.
    659650call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
    660651              tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice,            &
    661               watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif)
     652              watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,deposits)
    662653deallocate(tsurf_avg_yr1)
    663654
    … …  
    667658co2ice_sublim_surf_ini = 0.
    668659h2oice_ini_surf = 0.
     660h2o_ice_depth = 0.
    669661is_co2ice_sublim_ini = .false.
    670662is_h2oice_sublim_ini = .false.
    671663is_co2ice_ini = .false.
    672664co2ice_disappeared = .false.
     665if (layering_algo) then
     666    do ig = 1,ngrid
     667        do islope = 1,nslope
     668            co2_ice(ig,islope) = deposits(ig,islope)%top%h_co2ice
     669            h2o_ice(ig,islope) = deposits(ig,islope)%top%h_h2oice
     670        enddo
     671    enddo
     672endif
    673673do i = 1,ngrid
    674674    do islope = 1,nslope
    … …  
    681681            is_h2oice_sublim_ini(i,islope) = .true.
    682682            h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
     683            if (layering_algo) h2o_ice_depth(ig,islope) = thickness_toplag(deposits(ig,islope))
    683684        endif
    684685    enddo
    … …  
    736737stopPEM = 0
    737738if (layering_algo) then
    738     allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope))
     739    allocate(new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
    739740    new_str = .true.
    740     new_lag1 = .true.
    741     new_lag2 = .true.
     741    new_lag = .true.
    742742    do islope = 1,nslope
    743743        do ig = 1,ngrid
    744             current1(ig,islope)%p => stratif(ig,islope)%top
    745             current2(ig,islope)%p => stratif(ig,islope)%top
     744            current(ig,islope)%p => deposits(ig,islope)%top
    746745        enddo
    747746    enddo
    … …  
    820819!------------------------
    821820    allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
    822     call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
    823     call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
    824821    if (layering_algo) then
    825822        do islope = 1,nslope
    826823            do ig = 1,ngrid
    827                 call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),zshift_surf(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),zlag(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
     824                call make_layering(deposits(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),new_str(ig,islope),zshift_surf(ig,islope),new_lag(ig,islope),zlag(ig,islope),current(ig,islope)%p)
     825                co2_ice(ig,islope) = deposits(ig,islope)%top%h_co2ice
     826                h2o_ice(ig,islope) = deposits(ig,islope)%top%h_h2oice
    828827            enddo
    829828        enddo
    830829    else
    831830        zlag = 0.
     831        call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
     832        call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
    832833    endif
    833834
    … …  
    840841                                                            ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
    841842    if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
     843    if (layering_algo) then
     844        do islope = 1,nslope
     845            do ig = 1,ngrid
     846                deposits(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
     847                deposits(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
     848            enddo
     849        enddo
     850    endif
    842851
    843852!------------------------
    … …  
    891900    if (soil_pem) then
    892901! II_d.2 Shifting soil temperature to surface
    893         !call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
     902        call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
    894903
    895904! II_d.3 Update soil temperature
    … …  
    10141023    call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_old,ps_avg_global_new)
    10151024    write(*,*) "> Updating the H2O sub-surface ice tendency due to lag layer"
    1016     do ig = 1,ngrid
    1017         do islope = 1,nslope
    1018             call recomp_tend_h2o(icetable_depth_old(ig,islope),icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope))
    1019         enddo
    1020     enddo
     1025    if (layering_algo) then
     1026        do ig = 1,ngrid
     1027            do islope = 1,nslope
     1028                if (is_h2oice_sublim_ini(ig,islope)) then
     1029                    h2o_ice_depth_old = h2o_ice_depth(ig,islope)
     1030                    h2o_ice_depth(ig,islope) = thickness_toplag(deposits(ig,islope))
     1031                    call recomp_tend_h2o(h2o_ice_depth_old,h2o_ice_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope))
     1032                endif
     1033            enddo
     1034        enddo
     1035!~     else
     1036!~         do ig = 1,ngrid
     1037!~             do islope = 1,nslope
     1038!~                 call recomp_tend_h2o(icetable_depth_old(ig,islope),icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope))
     1039!~             enddo
     1040!~         enddo
     1041    endif
    10211042    if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old)
    10221043    deallocate(vmr_co2_PEM_phys)
    10231044    write(*,*) "> Updating the H2O sub-surface ice depth"
    1024     do ig = 1,ngrid
    1025         do islope = 1,nslope
    1026            if (icetable_depth(ig,islope) > 0) then
    1027                icetable_depth(ig,islope)=icetable_depth(ig,islope) + d_h2oice(ig,islope)*0.01 !!! 0.01 is for 1 percent dust, needs to be updated !!!
    1028             endif
    1029         enddo
    1030     enddo
    1031 
    10321045
    10331046!------------------------
    … …  
    10361049!------------------------
    10371050    write(*,*) "> Checking the stopping criteria"
    1038     h2o_ice=sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.))
    1039     co2_ice=sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.))
    10401051    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
    10411052    call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope)
    … …  
    10621073            case(7)
    10631074                write(*,'(a,i0)') " **** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
    1064             case(8)
    1065                 write(*,'(a,i0)') " **** STOPPING because the layering algorithm met an hasty end: ", stopPEM
    10661075            case default
    10671076                write(*,'(a,i0)') " **** STOPPING with unknown stopping criterion code: ", stopPEM
    … …  
    10801089deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
    10811090deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
    1082 if (layering_algo) deallocate(new_str,new_lag1,new_lag2,current1,current2)
     1091if (layering_algo) deallocate(new_str,new_lag,current)
    10831092!------------------------------ END RUN --------------------------------
    10841093
    … …  
    10861095!------------------------
    10871096! III Output
    1088 !    III_a Update surface value for the PCM start files
     1097!    III_a Update surface values for the PCM start files
    10891098!------------------------
    10901099! III_a.1 Ice update for start file
    … …  
    10931102write(*,*) '> Reconstructing perennial ice and frost for the PCM'
    10941103watercap = 0.
    1095 perennial_co2ice =sum(stratif(ig,islope)%top%co2ice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.))
     1104perennial_co2ice = co2_ice
    10961105do ig = 1,ngrid
    10971106    ! H2O ice metamorphism
    … …  
    11021111
    11031112    ! Is H2O ice still considered as an infinite reservoir for the PCM?
    1104     !if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
    1105     if (sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then   
    1106        ! There is enough ice to be considered as an infinite reservoir
     1113    if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
     1114        ! There is enough ice to be considered as an infinite reservoir
    11071115        watercaptag(ig) = .true.
    11081116    else
    1109         ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
     1117        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
    11101118        watercaptag(ig) = .false.
    1111 
    1112        ! qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
    1113 
    1114        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.))
     1119        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
    11151120        h2o_ice(ig,:) = 0.
    11161121    endif
    … …  
    11231128enddo
    11241129
    1125 ! III.a.2. Tsurf update for start file
     1130! III_a.2 Subsurface-surface interaction update for start file
     1131zdqsdif_ssi_tot = 0.
     1132h2o_ice_depth = 0.
     1133if (layering_algo) then
     1134    write(*,*) '> Reconstructing ice sublimation flux from subsurface to surface for the PCM'
     1135    do ig = 1,ngrid
     1136        do islope = 1,nslope
     1137            if (is_dust_lag(deposits(ig,islope)%top) .and. is_h2oice_str(deposits(ig,islope)%top%down)) then
     1138                zdqsdif_ssi_tot(ig,islope) = d_h2oice(ig,islope)
     1139                h2o_ice_depth(ig,islope) = thickness_toplag(deposits(ig,islope))
     1140            endif
     1141        enddo
     1142    enddo
     1143endif
     1144
     1145! III.a.3. Tsurf update for start file
    11261146write(*,*) '> Reconstructing the surface temperature for the PCM'
    11271147tsurf = tsurf_avg + tsurf_dev
    11281148deallocate(tsurf_dev)
    11291149
    1130 ! III_a.3 Tsoil update for start file
     1150! III_a.4 Tsoil update for start file
    11311151if (soil_pem) then
    11321152    write(*,*) '> Reconstructing the soil temperature profile for the PCM'
    … …  
    11431163deallocate(tsurf_avg,tsoil_dev)
    11441164
    1145 ! III_a.4 Pressure update for start file
     1165! III_a.5 Pressure update for start file
    11461166write(*,*) '> Reconstructing the pressure for the PCM'
    11471167allocate(ps_start(ngrid))
    … …  
    11501170deallocate(ps_avg,ps_dev)
    11511171
    1152 ! III_a.5 Tracers update for start file
     1172! III_a.6 Tracers update for start file
    11531173write(*,*) '> Reconstructing the tracer VMR for the PCM'
    11541174allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1))
    … …  
    11871207                q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
    11881208                write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
    1189            endif
     1209            endif
    11901210            if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
    11911211        enddo
    … …  
    11941214deallocate(zplev_new)
    11951215
    1196 ! III_a.6 Albedo update for start file
     1216! III_a.7 Albedo update for start file
    11971217write(*,*) '> Reconstructing the albedo for the PCM'
    11981218do ig = 1,ngrid
    … …  
    12311251enddo
    12321252
    1233 ! III_a.7 Orbital parameters update for start file
     1253! III_a.8 Orbital parameters update for start file
    12341254write(*,*) '> Setting the new orbital parameters'
    12351255if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
    … …  
    12931313!------------------------
    12941314write(*,*) '> Writing "restartpem.nc"'
    1295 if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
     1315if (layering_algo) nb_str_max = get_nb_str_max(deposits,ngrid,nslope) ! Get the maximum number of "stratum" in the depositsication (layerings)
    12961316call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
    12971317call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
    1298              co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,stratif)
     1318             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,deposits)
    12991319
    13001320call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
    … …  
    13021322write(*,*)
    13031323write(*,*) '****** PEM final information *******'
    1304 write(*,'(a,f10.2,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
    1305 write(*,'(a,f10.2,a,f10.2,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myear*convert_years, " Earth years."
    1306 write(*,'(a,f10.2,a)') " + The reached date is now ", (year_bp_ini + i_myear)*convert_years, " Earth years."
     1324write(*,'(a,f16.4,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
     1325write(*,'(a,f16.4,a,f16.4,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myear*convert_years, " Earth years."
     1326write(*,'(a,f16.4,a)') " + The reached date is now ", (year_bp_ini + i_myear)*convert_years, " Earth years."
    13071327write(*,*) "+ PEM: so far, so good!"
    13081328write(*,*) '************************************'
    … …  
    13111331    do islope = 1,nslope
    13121332        do i = 1,ngrid
    1313             call del_layering(stratif(i,islope))
     1333            call del_layering(deposits(i,islope))
    13141334        enddo
    13151335    enddo
    13161336endif
    13171337deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
    1318 deallocate(co2_ice,h2o_ice,stratif)
     1338deallocate(co2_ice,h2o_ice,deposits)
    13191339!----------------------------- END OUTPUT ------------------------------
    13201340

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

    r3673 r3770  
    2121use abort_pem_mod,              only: abort_pem
    2222use compute_soiltemp_mod,       only: ini_tsoil_pem, compute_tsoil_pem
    23 use layering_mod,               only: layering, array2stratif, nb_str_max, layering_algo
     23use layering_mod,               only: layering, array2stratif, nb_str_max, layering_algo, add_stratum, ini_layering
     24use callkeys_mod,               only: startphy_file
    2425
    2526#ifndef CPP_STD
    … …  
    3132
    3233implicit none
    33 
    34 include "callkeys.h"
    3534
    3635character(*),                   intent(in) :: filename      ! name of the startfi_PEM.nc
    … …  
    134133!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    135134    ! Stratification (layerings)
    136     nb_str_max = 1
     135    nb_str_max = 68
    137136    if (layering_algo) then
    138137        found = inquire_dimension("nb_str_max")
    139138        if (.not. found) then
    140139            write(*,*) 'Pemetat0: failed loading <nb_str_max>!'
    141             write(*,*) '''nb_str_max'' is set to 1!'
     140            write(*,*) '''nb_str_max'' is set to 68 (sub-surface layers)!'
    142141        else
    143142            nb_str_max = int(inquire_dimension_length('nb_str_max'))
    144143        endif
    145         allocate(stratif_array(ngrid,nslope,nb_str_max,7))
     144        allocate(stratif_array(ngrid,nslope,nb_str_max,6))
    146145        stratif_array = 0.
    147146        do islope = 1,nslope
    148147            write(num,'(i2.2)') islope
    149             call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
    150             call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
    151             found = found .or. found2
    152             call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
    153             found = found .or. found2
    154             call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
    155             found = found .or. found2
    156             call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
    157             found = found .or. found2
    158             call get_field('stratif_slope'//num//'_pore_volfrac',stratif_array(:,islope,:,6),found2)
    159             found = found .or. found2
    160             call get_field('stratif_slope'//num//'_icepore_volfrac',stratif_array(:,islope,:,7),found2)
    161             found = found .or. found2
    162             if (.not. found) then
    163                 write(*,*) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
    164             endif ! not found
     148            call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,1),found)
     149            if (.not. found) then
     150                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_top_elevation>!'
     151                exit
     152            endif
     153            call get_field('stratif_slope'//num//'_h_co2ice',stratif_array(:,islope,:,2),found)
     154            if (.not. found) then
     155                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_h_co2ice>!'
     156                exit
     157            endif
     158            call get_field('stratif_slope'//num//'_h_h2oice',stratif_array(:,islope,:,3),found)
     159            if (.not. found) then
     160                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_h_h2oice>!'
     161                exit
     162            endif
     163            call get_field('stratif_slope'//num//'_h_dust',stratif_array(:,islope,:,4),found)
     164            if (.not. found) then
     165                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_h_dust>!'
     166                exit
     167            endif
     168            call get_field('stratif_slope'//num//'_h_pore',stratif_array(:,islope,:,5),found)
     169            if (.not. found) then
     170                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_h_pore>!'
     171                exit
     172            endif
     173            call get_field('stratif_slope'//num//'_icepore_volfrac',stratif_array(:,islope,:,6),found)
     174            if (.not. found) then
     175                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_icepore_volfrac>!'
     176                exit
     177            endif
    165178        enddo ! islope
    166         call array2stratif(stratif_array,ngrid,nslope,stratif)
     179        if (.not. found) then
     180            write(*,*) 'So the deposits are initialized with sub-surface strata.'
     181            write(*,*) 'Ice is added with ''ini_huge_h2oice'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.'
     182            do ig = 1,ngrid
     183                if (watercaptag(ig)) then
     184                    do islope = 1,nslope
     185                        call ini_layering(stratif(ig,islope))
     186                        call add_stratum(stratif(ig,islope),ini_huge_h2oice,0.,ini_huge_h2oice,0.,0.,0.)
     187                    enddo
     188                else
     189                    do islope = 1,nslope
     190                        call ini_layering(stratif(ig,islope))
     191                        if (perennial_co2ice(ig,islope) > 0.) call add_stratum(stratif(ig,islope),perennial_co2ice(ig,islope),perennial_co2ice(ig,islope),0.,0.,0.,0.)
     192                    enddo
     193                endif
     194            enddo
     195        else
     196            call array2stratif(stratif_array,ngrid,nslope,stratif)
     197        endif
    167198        deallocate(stratif_array)
    168199    endif
    … …  
    337368    call close_startphy
    338369
    339 else !No startpem, let's build all by hand
     370else ! No startpem, let's build all by hand
    340371
    341372    write(*,*)'There is no "startpem.nc"!'
    … …  
    349380
    350381    ! co2 ice
    351     write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
     382    write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice'' found in the PCM.'
    352383    co2_ice = perennial_co2ice
    353384
    354385    ! Stratification (layerings)
    355     nb_str_max = 1
    356     if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with 1 sub-surface stratum.'
     386    nb_str_max = 68
     387    if (layering_algo) then
     388        write(*,*)'So the deposits are initialized with sub-surface strata.'
     389        write(*,*)'Ice is added with ''ini_huge_h2oice'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.'
     390        do ig = 1,ngrid
     391            if (watercaptag(ig)) then
     392                do islope = 1,nslope
     393                    call ini_layering(stratif(ig,islope))
     394                    call add_stratum(stratif(ig,islope),ini_huge_h2oice,0.,ini_huge_h2oice,0.,0.,0.)
     395                enddo
     396            else
     397                do islope = 1,nslope
     398                    call ini_layering(stratif(ig,islope))
     399                    if (perennial_co2ice(ig,islope) > 0.) call add_stratum(stratif(ig,islope),perennial_co2ice(ig,islope),perennial_co2ice(ig,islope),0.,0.,0.,0.)
     400                enddo
     401            endif
     402        enddo
     403    endif
    357404
    358405    if (soil_pem) then

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

    r3673 r3770  
    6868implicit none
    6969
    70 #ifndef CPP_STD
    71     include "callkeys.h"
    72 #endif
    73 
    7470character(*),                            intent(in) :: filename
    7571integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope
    … …  
    9995
    10096if (layering_algo) then
    101     allocate(stratif_array(ngrid,nslope,nb_str_max,7))
     97    allocate(stratif_array(ngrid,nslope,nb_str_max,6))
    10298    call stratif2array(stratif,ngrid,nslope,stratif_array)
    10399    do islope = 1,nslope
    104100        write(num,fmt='(i2.2)') islope
    105         call put_field('stratif_slope'//num//'_thickness','Layering thickness',stratif_array(:,islope,:,1),Year)
    106         call put_field('stratif_slope'//num//'_top_elevation','Layering top elevation',stratif_array(:,islope,:,2),Year)
    107         call put_field('stratif_slope'//num//'_co2ice_volfrac','Layering CO2 ice volume fraction',stratif_array(:,islope,:,3),Year)
    108         call put_field('stratif_slope'//num//'_h2oice_volfrac','Layering H2O ice volume fraction',stratif_array(:,islope,:,4),Year)
    109         call put_field('stratif_slope'//num//'_dust_volfrac','Layering dust volume fraction',stratif_array(:,islope,:,5),Year)
    110         call put_field('stratif_slope'//num//'_pore_volfrac','Layering pore volume fraction',stratif_array(:,islope,:,6),Year)
    111         call put_field('stratif_slope'//num//'_icepore_volfrac','Layering ice pore volume fraction',stratif_array(:,islope,:,7),Year)
     101        call put_field('stratif_slope'//num//'_top_elevation','Layering top elevation',stratif_array(:,islope,:,1),Year)
     102        call put_field('stratif_slope'//num//'_h_co2ice','Layering CO2 ice height',stratif_array(:,islope,:,2),Year)
     103        call put_field('stratif_slope'//num//'_h_h2oice','Layering H2O ice height',stratif_array(:,islope,:,3),Year)
     104        call put_field('stratif_slope'//num//'_h_dust','Layering dust height',stratif_array(:,islope,:,4),Year)
     105        call put_field('stratif_slope'//num//'_h_pore','Layering pore height',stratif_array(:,islope,:,5),Year)
     106        call put_field('stratif_slope'//num//'_icepore_volfrac','Layering ice pore volume fraction',stratif_array(:,islope,:,6),Year)
    112107    enddo
    113108    deallocate(stratif_array)

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

    r3707 r3770  
    9292Rz_old = icetable_depth_old*zcdv/coef_diff ! Old resistance from PCM
    9393Rz_new = icetable_depth_new*zcdv/coef_diff ! New resistance based on new depth
    94 R_dec = Rz_new/Rz_old ! Decrease because of resistance
     94if (abs(Rz_old) < 1.e-10) then
     95    R_dec = 1.
     96else
     97    R_dec = Rz_new/Rz_old ! Decrease because of resistance
     98endif
    9599
    96100! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location
    … …  
    103107
    104108! Lower humidity due to growing lag layer (higher depth)
     109print*, psv_max_old,psv_max_new
    105110hum_dec = psv_max_old/psv_max_new ! Decrease because of lower water vapor pressure at the new depth
    106111
    107112! Flux correction (decrease)
     113print*, 'aaaaah', R_dec,hum_dec
    108114d_h2oice = d_h2oice/R_dec/hum_dec
    109115

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

    r3584 r3770  
    5757
    5858do iloop = 0,nsoil_PEM - 1
    59     if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then
     59    if (lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then
    6060        index_powerlaw = iloop
    6161        exit
    … …  
    6969! Build layer_PEM()
    7070do iloop = 1,nsoil_PEM - 1
    71 layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2.
     71    layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2.
    7272enddo
    7373layer_PEM(nsoil_PEM) = 2*mlayer_PEM(nsoil_PEM - 1) - mlayer_PEM(nsoil_PEM - 2)

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

    r3711 r3770  
    6969endif
    7070
    71 !if (abs(h2oice_ini_surf) < 1.e-5) stopPEM = 0
    72 
    7371END SUBROUTINE stopping_crit_h2o_ice
    7472
    … …  
    133131endif
    134132
    135 !if (abs(co2ice_sublim_surf_ini) < 1.e-5) stopPEM = 0
    136 
    137133if (ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)) then
    138134    stopPEM = 4
Note: See TracChangeset for help on using the changeset viewer.