Ignore:
Timestamp:
May 30, 2025, 5:57:10 PM (3 weeks ago)
Author:
jbclement
Message:

PEM:

  • Rework of ice sublimation in the layering algorithm to sublimate ice and to lift dust at the same time
  • Addition of the possibility to make H2O ice or mixed dust/H2O ice lag when CO2 ice is sublimating
  • Addition of a lag resistance for CO2 ice sublimation
  • Consideration of a 'wet_lag_porosity' for H2O ice lag layer
  • Consideration of a threshold under which H2O ice lag layer is considered patchy
  • Deletion of unused function 'thickness_toplag'

JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3784 r3785  
    681681== 30/05/2025 == JBC
    682682Correction of 'h2o_ice_depth' update in the main loop.
     683
     684== 30/05/2025 == JBC
     685- Rework of ice sublimation in the layering algorithm to sublimate ice and to lift dust at the same time
     686- Addition of the possibility to make H2O ice or mixed dust/H2O ice lag when CO2 ice is sublimating
     687- Addition of a lag resistance for CO2 ice sublimation
     688- Consideration of a 'wet_lag_porosity' for H2O ice lag layer
     689- Consideration of a threshold under which H2O ice lag layer is considered patchy
     690- Deletion of unused function 'thickness_toplag'

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

    r3782 r3785  
    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
    23 real, parameter :: dry_regolith_porosity = 0.4     ! Porosity of regolith
    24 real, parameter :: breccia_porosity      = 0.1     ! Porosity of breccia
    25 real, parameter :: co2ice_porosity       = 0.      ! Porosity of CO2 ice
    26 real, parameter :: h2oice_porosity       = 0.      ! Porosity of H2O ice
    27 real, parameter :: rho_dust              = 2500.   ! Density of dust [kg.m-3]
    28 real, parameter :: rho_rock              = 3200.   ! Density of rock [kg.m-3]
    29 real, parameter :: h_patchy_dust         = 0.05    ! Height under which a dust statum is considered patchy [m]
     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
     23real, parameter :: wet_lag_porosity  = 0.1     ! Porosity of water ice lag
     24real, parameter :: regolith_porosity = 0.4     ! Porosity of regolith
     25real, parameter :: breccia_porosity  = 0.1     ! Porosity of breccia
     26real, parameter :: co2ice_porosity   = 0.      ! Porosity of CO2 ice
     27real, parameter :: h2oice_porosity   = 0.      ! Porosity of H2O ice
     28real, parameter :: rho_dust          = 2500.   ! Density of dust [kg.m-3]
     29real, parameter :: rho_rock          = 3200.   ! Density of rock [kg.m-3]
     30real, parameter :: h_patchy_dust     = 0.05    ! Height under which a dust lag is considered patchy [m]
     31real, parameter :: h_patchy_ice      = 0.05    ! Height under which a H2O ice lag is considered patchy [m]
    3032
    3133! Parameters for CO2 flux correction (inspired by Levrard et al. 2007)
    … …  
    7779! Procedures to get information about a stratum:
    7880!     > thickness
    79 !     > thickness_toplag
    8081!     > is_co2ice_str
    8182!     > is_h2oice_str
    … …  
    8586! Procedures for the algorithm to evolve the layering:
    8687!     > make_layering
    87 !     > lift_dust_lag
     88!     > lift_dust
    8889!     > subl_ice_str
    8990!=======================================================================
    … …  
    294295    do i = index_breccia - 1,2,-1
    295296        h_tot = layer_PEM(i) - layer_PEM(i - 1)
    296         h_pore = h_tot*dry_regolith_porosity
     297        h_pore = h_tot*regolith_porosity
    297298        h_soil = h_tot - h_pore
    298299        call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.)
    299300    enddo
    300     h_pore = layer_PEM(1)*dry_regolith_porosity
     301    h_pore = layer_PEM(1)*regolith_porosity
    301302    h_soil = layer_PEM(1) - h_pore
    302303    call add_stratum(this,0.,0.,0.,h_soil,h_pore,0.)
    … …  
    482483
    483484END FUNCTION thickness
    484 
    485 !=======================================================================
    486 ! To get the top lag height
    487 FUNCTION thickness_toplag(this) RESULT(h_toplag)
    488 
    489 implicit none
    490 
    491 !---- Arguments
    492 type(layering), intent(in) :: this
    493 real :: h_toplag
    494 
    495 !---- Local variables
    496 type(stratum), pointer :: current
    497 
    498 !---- Code
    499 h_toplag = 0.
    500 current => this%top
    501 ! Is the considered stratum a dust lag?
    502 do while (is_dust_lag(current) .and. associated(current))
    503     h_toplag = h_toplag + thickness(current)
    504     current => current%down
    505 enddo
    506 
    507 END FUNCTION thickness_toplag
    508485
    509486!=======================================================================
    … …  
    610587!=======================================================================
    611588!=======================================================================
     589! To lift dust in a stratum
     590SUBROUTINE lift_dust(this,str,h2lift)
     591
     592implicit none
     593
     594!---- Arguments
     595type(layering),         intent(inout) :: this
     596type(stratum), pointer, intent(inout) :: str
     597real, intent(in) :: h2lift
     598
     599!---- Code
     600! Update of properties in the eroding dust lag
     601str%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust)
     602str%h_pore        = str%h_pore        - h2lift*str%h_pore/str%h_dust
     603str%h_dust        = str%h_dust        - h2lift
     604
     605! Remove the eroding dust lag if there is no dust anymore
     606if (str%h_dust < tol) call del_stratum(this,str)
     607
     608END SUBROUTINE lift_dust
     609
     610!=======================================================================
     611! To sublimate ice in a stratum
     612SUBROUTINE sublimate_ice(this,str,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)
     613
     614implicit none
     615
     616!---- Arguments
     617type(layering),         intent(inout) :: this
     618type(stratum), pointer, intent(inout) :: str
     619logical,                intent(inout) :: new_lag
     620real,                   intent(inout) :: zlag
     621real, intent(in) :: h2subl_co2ice, h2subl_h2oice
     622
     623!---- Local variables
     624real                   :: h_ice, h2subl, h_pore, h_tot
     625real                   :: hlag_dust, hlag_h2oice
     626type(stratum), pointer :: current
     627
     628!---- Code
     629h_ice = str%h_co2ice + str%h_h2oice
     630h2subl = h2subl_co2ice + h2subl_h2oice
     631
     632! How much dust and H2O ice does ice sublimation involve to create the lag?
     633hlag_dust = str%h_dust*h2subl/h_ice
     634hlag_h2oice = 0.
     635if (h2subl_co2ice > 0. .and. h2subl_h2oice < tol) hlag_h2oice = str%h_h2oice*h2subl/h_ice
     636
     637! Update of properties in the sublimating stratum
     638str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hlag_dust - hlag_h2oice
     639str%h_co2ice      = str%h_co2ice      - h2subl_co2ice
     640str%h_h2oice      = str%h_h2oice      - h2subl_h2oice - hlag_h2oice
     641str%h_dust        = str%h_dust        - hlag_dust
     642str%h_pore        = str%h_pore        - h2subl*h_pore/h_ice
     643
     644! Correct the value of top_elevation for the upper strata
     645current => str%up
     646do while (associated(current))
     647    current%top_elevation = current%down%top_elevation + thickness(current)
     648    current => current%up
     649enddo
     650
     651! Remove the sublimating stratum if there is no ice anymore
     652if (str%h_co2ice < tol .and. str%h_h2oice < tol) call del_stratum(this,str)
     653
     654! Which porosity is considered?
     655if (hlag_dust >= hlag_h2oice) then
     656    h_pore = hlag_dust*dry_lag_porosity/(1. - dry_lag_porosity)
     657else
     658    h_pore = hlag_h2oice*wet_lag_porosity/(1. - wet_lag_porosity)
     659endif
     660h_tot = hlag_dust + hlag_h2oice + h_pore
     661
     662! New stratum above the current stratum as a lag layer
     663if (new_lag) then
     664    call insert_stratum(this,str,str%top_elevation + h_tot,0.,hlag_h2oice,hlag_dust,h_pore,0.)
     665    new_lag = .false.
     666else
     667    str%up%top_elevation = str%up%top_elevation + h_tot
     668    str%up%h_h2oice      = str%up%h_h2oice      + hlag_h2oice
     669    str%up%h_dust        = str%up%h_dust        + hlag_dust
     670    str%up%h_pore        = str%up%h_pore        + h_pore
     671endif
     672
     673zlag = zlag + h_tot
     674
     675END SUBROUTINE sublimate_ice
     676
     677!=======================================================================
    612678! Layering algorithm
    613679SUBROUTINE make_layering(this,d_co2ice,d_h2oice,new_str,zshift_surf,new_lag,zlag,current)
    … …  
    626692
    627693!---- Local variables
    628 real    :: dh_co2ice, dh_h2oice, dh_dust
    629 real    :: h_co2ice_old, h_h2oice_old, h_dust_old
    630 real    :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
    631 real    :: h2lift, h2lift_tot, h_pore, h_tot
    632 logical :: unexpected
     694real                   :: dh_co2ice, dh_h2oice, dh_dust
     695real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
     696real                   :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
     697real                   :: h2lift, h2lift_tot, h_pore, h_tot, h_toplag, R_subl
     698logical                :: unexpected
     699type(stratum), pointer :: currentloc
    633700
    634701!---- Code
    … …  
    653720        h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity)
    654721    else if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice) then ! Dust is dominant
    655         h_pore = dh_dust*dry_regolith_porosity/(1. - dry_regolith_porosity)
     722        h_pore = dh_dust*regolith_porosity/(1. - regolith_porosity)
    656723    else
    657724        unexpected = .true.
    … …  
    670737    endif
    671738!-----------------------------------------------------------------------
    672 else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. dh_dust <= 0.) then ! A stratum is retreating
    673     ! DUST LIFTING
    674     if (abs(dh_co2ice) < tol .and. abs(dh_co2ice) < tol) then
    675         write(*,*) '> Layering: dust lifting'
    676         h2lift_tot = abs(dh_dust)
    677         do while (h2lift_tot > 0. .and. associated(current))
    678             ! Is the considered stratum a dust lag?
    679             if (current%h_co2ice > 0. .or. current%h_h2oice > 0. .or. current%poreice_volfrac > 0.) then ! No, there is ice cementing the dust
    680                 dh_dust = 0. ! The tendency is set to 0
    681                 exit
    682             else ! Yes, we can lift dust
    683                 h2lift = min(h2lift_tot,current%h_dust) ! How much dust can be lifted?
    684                 h2lift_tot = h2lift_tot - h2lift
    685                 call lift_dust_lag(this,current,h2lift)
    686                 if (h2lift_tot > 0.) current => current%down ! Still dust to be lifted so we move to the underlying stratum
    687             endif
    688         enddo
    689         if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
    690 !-----------------------------------------------------------------------
    691     ! ICE SUBLIMATION
    692     else if (abs(dh_dust) < tol) then
    693         write(*,*) '> Layering: ice sublimation'
     739 ! RETREATING A STRATUM (ice sublimation and/or dust lifting)
     740else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. dh_dust <= 0.) then
    694741        h2subl_co2ice_tot = abs(dh_co2ice)
    695742        h2subl_h2oice_tot = abs(dh_h2oice)
    696         do while (h2subl_co2ice_tot > 0. .and. h2subl_h2oice_tot > 0. .and. associated(current))
     743        h2lift_tot = abs(dh_dust)
     744        do while (h2subl_co2ice_tot > 0. .and. h2subl_h2oice_tot > 0. .and. h2lift_tot > 0. .and. associated(current))
    697745            h_co2ice_old = current%h_co2ice
    698746            h_h2oice_old = current%h_h2oice
    699    
     747
     748            ! How much is CO2 ice sublimation inhibited by the top dust lag?
     749            R_subl = 1.
     750            if (associated(current%up)) then ! If there is an upper stratum
     751                h_toplag = 0.
     752                currentloc => current%up
     753                do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
     754                    h_toplag = h_toplag + thickness(currentloc)
     755                    currentloc => currentloc%up
     756                enddo
     757                if (currentloc%h_h2oice > h_patchy_ice) then
     758                    R_subl = 0.
     759                else if (0. < currentloc%h_h2oice .and. currentloc%h_h2oice <= h_patchy_ice) then
     760                    h_toplag = h_toplag + thickness(currentloc)
     761                    R_subl = fred_subl**(h_toplag/hmin_lag)
     762                else
     763                    R_subl = fred_subl**(h_toplag/hmin_lag)
     764                endif
     765            endif
     766           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
     767
    700768            ! Is there CO2 ice to sublimate?
    701769            h2subl_co2ice = 0.
    … …  
    710778                h2subl_h2oice_tot = h2subl_h2oice_tot - h2subl_h2oice
    711779            endif
    712             ! Sublimation
    713             if (h2subl_co2ice > 0. .or. h2subl_h2oice > 0.) call subl_ice_str(this,current,h2subl_co2ice,h2subl_h2oice,zlag,new_lag)
    714             ! Still ice to be sublimated and no more ice in the current stratum so we move to the underlying stratum
    715             if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0.) .and. current%h_co2ice < tol .and. current%h_h2oice < tol) then
     780            ! Ice sublimation
     781            if (h2subl_co2ice > 0. .or. h2subl_h2oice > 0.) call sublimate_ice(this,current,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)
     782
     783            ! Is there dust to lift?
     784            h2lift = 0.
     785            if (is_dust_lag(current) .and. h2lift_tot > 0.) then
     786                h2lift = min(h2lift_tot,current%h_dust)
     787                h2lift_tot = h2lift_tot - h2lift
     788                ! Dust lifting
     789                if (h2lift > 0.) call lift_dust(this,current,h2lift)
     790            else if (associated(current%up)) then
     791                if (is_dust_lag(current%up) .and. h2lift_tot > 0.) then
     792                    h2lift = min(h2lift_tot,current%up%h_dust)
     793                    h2lift_tot = h2lift_tot - h2lift
     794                    ! Dust lifting
     795                    if (h2lift > 0.) call lift_dust(this,current%up,h2lift)
     796                endif
     797            endif
     798
     799            ! Still ice to be sublimated or dust to be lifted and no more ice in the current stratum so we move to the underlying stratum
     800            if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0. .or. h2lift_tot > 0.) .and. current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < tol) then
    716801                current => current%down
    717802                new_lag = .true.
    … …  
    722807        if (h2subl_co2ice_tot > 0.) dh_co2ice = 0. ! No CO2 ice available anymore so the tendency is set to 0
    723808        if (h2subl_h2oice_tot > 0.) dh_h2oice = 0. ! No H2O ice available anymore so the tendency is set to 0
    724     else
    725         unexpected = .true.
    726     endif
     809        if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
    727810!-----------------------------------------------------------------------
    728811! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION
    729812else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then
    730 !~     if (abs(dh_co2ice) > dh_h2oice .and. dh_co2ice < dh_dust .and. dh_dust <= 0.) then
    731 !~         h2subl_co2ice_tot = abs(dh_co2ice)
    732 !~         do while (h2subl_co2ice_tot > 0. associated(current))
    733 !~             h_co2ice_old = current%h_co2ice
    734 !~             h_h2oice_old = current%h_h2oice
    735    
    736 !~             ! Is there CO2 ice to sublimate?
    737 !~             h2subl_co2ice = 0.
    738 !~             if (h_co2ice_old > 0. .and. h2subl_co2ice_tot > 0.) then
    739 !~                 h2subl_co2ice = min(h2subl_co2ice_tot,h_co2ice_old)
    740 !~                 h2subl_co2ice_tot = h2subl_co2ice_tot - h2subl_co2ice
    741 !~             endif
    742 !~             ! Sublimation
    743 !~             if (h2subl_co2ice > 0.) call subl_ice_str(this,current,h2subl_co2ice,0.,zlag,new_lag)
    744 !~             ! Still ice to be sublimated and no more ice in the current stratum so we move to the underlying stratum
    745 !~             if (h2subl_co2ice_tot > 0. .and. current%h_co2ice < tol .and. current%h_h2oice < tol) then
    746 !~                 current => current%down
    747 !~                 new_lag = .true.
    748 !~             else
    749 !~                 exit
    750 !~             endif
    751 !~         enddo
    752 !~         if (h2subl_co2ice_tot > 0.) dh_co2ice = 0. ! No CO2 ice available anymore so the tendency is set to 0
    753 !~     else if (dh_h2oice > abs(dh_co2ice) .and. 0 <= dh_dust .and. dh_dust < dh_h2oice) then
     813!~     if (abs(dh_co2ice) > dh_h2oice .and. dh_co2ice < dh_dust .and. dh_dust <= 0.) then ! CO2 sublimation is dominant
     814
     815!~     else if (dh_h2oice > abs(dh_co2ice) .and. 0 <= dh_dust .and. dh_dust < dh_h2oice) then ! H2O condensation is dominant
    754816
    755817!~     else
    … …  
    774836END SUBROUTINE make_layering
    775837
    776 !=======================================================================
    777 ! To lift dust in a stratum
    778 SUBROUTINE lift_dust_lag(this,str,h2lift)
    779 
    780 implicit none
    781 
    782 !---- Arguments
    783 type(layering),         intent(inout) :: this
    784 type(stratum), pointer, intent(inout) :: str
    785 real, intent(in) :: h2lift
    786 
    787 !---- Code
    788 ! Update of properties in the eroding dust lag
    789 str%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust)
    790 str%h_pore        = str%h_pore        - h2lift*str%h_pore/str%h_dust
    791 str%h_dust        = str%h_dust        - h2lift
    792 
    793 ! Remove the eroding dust lag if there is no dust anymore
    794 if (str%h_dust < tol) call del_stratum(this,str)
    795 
    796 END SUBROUTINE lift_dust_lag
    797 
    798 !=======================================================================
    799 ! To sublimate ice in a stratum
    800 SUBROUTINE subl_ice_str(this,str,h2subl_co2ice,h2subl_h2oice,zlag,new_lag)
    801 
    802 implicit none
    803 
    804 !---- Arguments
    805 type(layering),         intent(inout) :: this
    806 type(stratum), pointer, intent(inout) :: str
    807 logical,                intent(inout) :: new_lag
    808 real,                   intent(inout) :: zlag
    809 real, intent(in) :: h2subl_co2ice, h2subl_h2oice
    810 
    811 !---- Local variables
    812 real                   :: hsubl_dust, h_pore, h_ice, h_tot, h2subl
    813 type(stratum), pointer :: current
    814 
    815 !---- Code
    816 h_ice = str%h_co2ice + str%h_h2oice
    817 h2subl = h2subl_co2ice + h2subl_h2oice
    818 
    819 ! How much dust does ice sublimation involve to create the lag?
    820 hsubl_dust = str%h_dust*h2subl/h_ice
    821 
    822 ! Update of properties in the sublimating stratum
    823 str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hsubl_dust
    824 str%h_co2ice      = str%h_co2ice      - h2subl_co2ice
    825 str%h_h2oice      = str%h_h2oice      - h2subl_h2oice
    826 str%h_dust        = str%h_dust        - hsubl_dust
    827 str%h_pore        = str%h_pore        - h2subl*h_pore/h_ice
    828 
    829 ! Correct the value of top_elevation for the upper strata
    830 current => str%up
    831 do while (associated(current))
    832     current%top_elevation = current%down%top_elevation + thickness(current)
    833     current => current%up
    834 enddo
    835 
    836 ! Remove the sublimating stratum if there is no ice anymore
    837 if (str%h_co2ice < tol .and. str%h_h2oice < tol) call del_stratum(this,str)
    838 
    839 ! New stratum above the current stratum as a dust lag
    840 if (hsubl_dust > 0.) then ! Only if the dust lag is building up
    841     h_pore = hsubl_dust*dry_lag_porosity/(1. - dry_lag_porosity)
    842     h_tot = hsubl_dust + h_pore
    843     if (new_lag) then
    844         call insert_stratum(this,str,str%top_elevation + h_tot,0.,0.,hsubl_dust,h_pore)
    845         new_lag = .false.
    846     else
    847         str%up%top_elevation = str%up%top_elevation + h_tot
    848     endif
    849 endif
    850 zlag = zlag + h_tot
    851 
    852 END SUBROUTINE subl_ice_str
    853 
    854838END MODULE layering_mod

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

    r3784 r3785  
    6969use writediagpem_mod,           only: writediagpem, writediagsoilpem
    7070use co2condens_mod,             only: CO2cond_ps
    71 use layering_mod,               only: ptrarray, stratum, layering, del_layering, make_layering, layering_algo, &
    72                                       get_nb_str_max, nb_str_max, thickness_toplag, subsurface_ice_layering,   &
    73                                       is_dust_lag, is_co2ice_str, is_h2oice_str
     71use layering_mod,               only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, &
     72                                      ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str
    7473use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
    7574use version_info_mod,           only: print_version_info
Note: See TracChangeset for help on using the changeset viewer.