Ignore:
Timestamp:
May 26, 2025, 10:34:57 AM (4 weeks ago)
Author:
jbclement
Message:

PEM:

  • New subroutine to detect whether there is subsurface ice or not
  • Rework of the initialization/update/finalization of the situation regarding the layering data structure
  • Introduction of a threshold 'h_patchy_dust' under which the top dust layer is not considered as a stratum
  • 'deposits' is renamed as 'layerings_map'
  • Few cleanings

JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3770 r3778  
    22
    33!=======================================================================
    4 ! Purpose: Manage layered deposits in the PEM with a linked list data structure
     4! Purpose: Manage layered layerings_map in the PEM with a linked list data structure
    55!
    66! Author: JBC
    … …  
    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]
     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 :: dry_regolith_porosity = 0.4     ! Porosity of regolith
     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]
     29real, parameter :: h_patchy_dust         = 0.05    ! Height under which a dust statum is considered patchy [m]
    2930
    3031! Lag layer parameters -> see Levrard et al. 2007
    31 real, parameter :: hmin_lag = 0.5 ! Minimal height of the lag deposit to reduce the sublimation rate
     32real, parameter :: hmin_lag  = 0.5 ! Minimal height of the lag deposit to reduce the sublimation rate
    3233real, parameter :: fred_subl = 0.1 ! Reduction factor of sublimation rate
    3334
    … …  
    4950    type(stratum), pointer :: bottom => null() ! First stratum (head of the list)
    5051    type(stratum), pointer :: top => null()    ! Last stratum (tail of the list)
     52    !real                   :: h_toplag         ! Height of dust layer at the top [m]
     53    !type(stratum), pointer :: subice => null() ! Shallower stratum containing ice
     54    !type(stratum), pointer :: surf => null()   ! Surface stratum
    5155end type layering
    5256
    … …  
    7074!     > stratif2array
    7175!     > array2stratif
    72 !     > print_deposits
     76!     > print_layerings_map
    7377! Procedures to get information about a stratum:
    7478!     > thickness
     79!     > thickness_toplag
    7580!     > is_co2ice_str
    7681!     > is_h2oice_str
    7782!     > is_dust_str
    7883!     > is_dust_lag
    79 !     > compute_h_toplag
     84!     > subsurface_ice_layering
    8085! Procedures for the algorithm to evolve the layering:
    8186!     > make_layering
    … …  
    118123allocate(str)
    119124nullify(str%up,str%down)
    120 str%top_elevation = 1.
    121 str%h_co2ice = 0.
    122 str%h_h2oice = 0.
    123 str%h_dust = 1.
    124 str%h_pore = 0.
     125str%top_elevation   = 1.
     126str%h_co2ice        = 0.
     127str%h_h2oice        = 0.
     128str%h_dust          = 1.
     129str%h_pore          = 0.
    125130str%poreice_volfrac = 0.
    126 if (present(top_elevation)) str%top_elevation = top_elevation
    127 if (present(co2ice)) str%h_co2ice = co2ice
    128 if (present(h2oice)) str%h_h2oice = h2oice
    129 if (present(dust)) str%h_dust = dust
    130 if (present(pore)) str%h_pore = pore
    131 if (present(poreice)) str%poreice_volfrac = poreice
     131if (present(top_elevation)) str%top_elevation   = top_elevation
     132if (present(co2ice))        str%h_co2ice        = co2ice
     133if (present(h2oice))        str%h_h2oice        = h2oice
     134if (present(dust))          str%h_dust          = dust
     135if (present(pore))          str%h_pore          = pore
     136if (present(poreice))       str%poreice_volfrac = poreice
    132137
    133138! Increment the number of strata
    … …  
    167172    allocate(str)
    168173    nullify(str%up,str%down)
    169     str%top_elevation = 1.
    170     str%h_co2ice = 0.
    171     str%h_h2oice = 0.
    172     str%h_dust = 1.
    173     str%h_pore = 0.
     174    str%top_elevation   = 1.
     175    str%h_co2ice        = 0.
     176    str%h_h2oice        = 0.
     177    str%h_dust          = 1.
     178    str%h_pore          = 0.
    174179    str%poreice_volfrac = 0.
    175     if (present(top_elevation)) str%top_elevation = top_elevation
    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
    180     if (present(poreice)) str%poreice_volfrac = poreice
     180    if (present(top_elevation)) str%top_elevation   = top_elevation
     181    if (present(co2ice))        str%h_co2ice        = co2ice
     182    if (present(h2oice))        str%h_h2oice        = h2oice
     183    if (present(dust))          str%h_dust          = dust
     184    if (present(pore))          str%h_pore          = pore
     185    if (present(poreice))       str%poreice_volfrac = poreice
    181186
    182187    ! Increment the number of strata
    … …  
    353358!=======================================================================
    354359! To get the maximum number of "stratum" in the stratification (layerings)
    355 FUNCTION get_nb_str_max(stratif,ngrid,nslope) RESULT(nb_str_max)
    356 
    357 implicit none
    358 
    359 !---- Arguments
    360 type(layering), dimension(ngrid,nslope), intent(in) :: stratif
     360FUNCTION get_nb_str_max(layerings_map,ngrid,nslope) RESULT(nb_str_max)
     361
     362implicit none
     363
     364!---- Arguments
     365type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map
    361366integer,                                 intent(in) :: ngrid, nslope
    362367integer :: nb_str_max
    … …  
    369374do islope = 1,nslope
    370375    do ig = 1,ngrid
    371         nb_str_max = max(stratif(ig,islope)%nb_str,nb_str_max)
     376        nb_str_max = max(layerings_map(ig,islope)%nb_str,nb_str_max)
    372377    enddo
    373378enddo
    … …  
    376381
    377382!=======================================================================
    378 ! To convert the stratification data structure (layerings) into an array able to be outputted
    379 SUBROUTINE stratif2array(stratif,ngrid,nslope,stratif_array)
     383! To convert the layerings map into an array able to be outputted
     384SUBROUTINE stratif2array(layerings_map,ngrid,nslope,stratif_array)
    380385
    381386implicit none
    … …  
    383388!---- Arguments
    384389integer,                                 intent(in) :: ngrid, nslope
    385 type(layering), dimension(ngrid,nslope), intent(in) :: stratif
     390type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map
    386391real, dimension(:,:,:,:), allocatable, intent(inout) :: stratif_array
    387392
    … …  
    394399do islope = 1,nslope
    395400    do ig = 1,ngrid
    396         current => stratif(ig,islope)%bottom
     401        current => layerings_map(ig,islope)%bottom
    397402        k = 1
    398403        do while (associated(current))
    … …  
    412417
    413418!=======================================================================
    414 ! To convert the stratification array into the stratification data structure (layerings)
    415 SUBROUTINE array2stratif(stratif_array,ngrid,nslope,stratif)
     419! To convert the stratification array into the layerings map
     420SUBROUTINE array2stratif(stratif_array,ngrid,nslope,layerings_map)
    416421
    417422implicit none
    … …  
    420425integer,                               intent(in) :: ngrid, nslope
    421426real, dimension(:,:,:,:), allocatable, intent(in) :: stratif_array
    422 type(layering), dimension(ngrid,nslope), intent(inout) :: stratif
     427type(layering), dimension(ngrid,nslope), intent(inout) :: layerings_map
    423428
    424429!---- Local variables
    … …  
    429434    do ig = 1,ngrid
    430435        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))
     436            call add_stratum(layerings_map(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6))
    432437        enddo
    433438    enddo
    … …  
    437442
    438443!=======================================================================
    439 ! To display deposits
    440 SUBROUTINE print_deposits(this,ngrid,nslope)
     444! To display the map of layerings
     445SUBROUTINE print_layerings_map(this,ngrid,nslope)
    441446
    442447implicit none
    … …  
    460465enddo
    461466
    462 END SUBROUTINE print_deposits
     467END SUBROUTINE print_layerings_map
    463468
    464469!=======================================================================
    … …  
    477482
    478483END FUNCTION thickness
    479 
    480 !=======================================================================
    481 ! To know whether the stratum can be considered as a CO2 ice layer
    482 FUNCTION is_co2ice_str(str) RESULT(co2ice_str)
    483 
    484 implicit none
    485 
    486 !---- Arguments
    487 type(stratum), intent(in) :: str
    488 logical :: co2ice_str
    489 
    490 !---- Code
    491 co2ice_str = .false.
    492 if (str%h_co2ice > str%h_h2oice .and. str%h_co2ice > str%h_dust) co2ice_str = .true.
    493 
    494 END FUNCTION is_co2ice_str
    495 
    496 !=======================================================================
    497 ! To know whether the stratum can be considered as a H2O ice layer
    498 FUNCTION is_h2oice_str(str) RESULT(h2oice_str)
    499 
    500 implicit none
    501 
    502 !---- Arguments
    503 type(stratum), intent(in) :: str
    504 logical :: h2oice_str
    505 
    506 !---- Code
    507 h2oice_str = .false.
    508 if (str%h_h2oice > str%h_co2ice .and. str%h_h2oice > str%h_dust) h2oice_str = .true.
    509 
    510 END FUNCTION is_h2oice_str
    511 
    512 !=======================================================================
    513 ! To know whether the stratum can be considered as a dust layer
    514 FUNCTION is_dust_str(str) RESULT(dust_str)
    515 
    516 implicit none
    517 
    518 !---- Arguments
    519 type(stratum), intent(in) :: str
    520 logical :: dust_str
    521 
    522 !---- Code
    523 dust_str = .false.
    524 if (str%h_dust > str%h_co2ice .and. str%h_dust > str%h_h2oice) dust_str = .true.
    525 
    526 END FUNCTION is_dust_str
    527 
    528 !=======================================================================
    529 ! To know whether the stratum can be considered as a dust lag
    530 FUNCTION is_dust_lag(str) RESULT(dust_lag)
    531 
    532 implicit none
    533 
    534 !---- Arguments
    535 type(stratum), intent(in) :: str
    536 logical :: dust_lag
    537 
    538 !---- Code
    539 dust_lag = .false.
    540 if (str%h_co2ice < tol .and. str%h_h2oice < tol .and. str%poreice_volfrac < tol) dust_lag = .true.
    541 
    542 END FUNCTION is_dust_lag
    543484
    544485!=======================================================================
    … …  
    559500current => this%top
    560501! Is the considered stratum a dust lag?
    561 do while (current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < tol .and. associated(current))
     502do while (is_dust_lag(current) .and. associated(current))
    562503    h_toplag = h_toplag + thickness(current)
    563504    current => current%down
    … …  
    565506
    566507END FUNCTION thickness_toplag
     508
     509!=======================================================================
     510! To know whether the stratum can be considered as a CO2 ice layer
     511FUNCTION is_co2ice_str(str) RESULT(co2ice_str)
     512
     513implicit none
     514
     515!---- Arguments
     516type(stratum), intent(in) :: str
     517logical :: co2ice_str
     518
     519!---- Code
     520co2ice_str = .false.
     521if (str%h_co2ice > str%h_h2oice .and. str%h_co2ice > str%h_dust) co2ice_str = .true.
     522
     523END FUNCTION is_co2ice_str
     524
     525!=======================================================================
     526! To know whether the stratum can be considered as a H2O ice layer
     527FUNCTION is_h2oice_str(str) RESULT(h2oice_str)
     528
     529implicit none
     530
     531!---- Arguments
     532type(stratum), intent(in) :: str
     533logical :: h2oice_str
     534
     535!---- Code
     536h2oice_str = .false.
     537if (str%h_h2oice > str%h_co2ice .and. str%h_h2oice > str%h_dust) h2oice_str = .true.
     538
     539END FUNCTION is_h2oice_str
     540
     541!=======================================================================
     542! To know whether the stratum can be considered as a dust layer
     543FUNCTION is_dust_str(str) RESULT(dust_str)
     544
     545implicit none
     546
     547!---- Arguments
     548type(stratum), intent(in) :: str
     549logical :: dust_str
     550
     551!---- Code
     552dust_str = .false.
     553if (str%h_dust > str%h_co2ice .and. str%h_dust > str%h_h2oice) dust_str = .true.
     554
     555END FUNCTION is_dust_str
     556
     557!=======================================================================
     558! To know whether the stratum can be considered as a dust lag
     559FUNCTION is_dust_lag(str) RESULT(dust_lag)
     560
     561implicit none
     562
     563!---- Arguments
     564type(stratum), intent(in) :: str
     565logical :: dust_lag
     566
     567!---- Code
     568dust_lag = .false.
     569if (str%h_co2ice < tol .and. str%h_h2oice < tol .and. str%poreice_volfrac < tol) dust_lag = .true.
     570
     571END FUNCTION is_dust_lag
     572
     573!=======================================================================
     574! To get data about possible subsurface ice in a layering
     575SUBROUTINE subsurface_ice_layering(this,is_subsurface_ice,h_toplag,h2o_ice)
     576
     577implicit none
     578
     579!---- Arguments
     580type(layering), intent(in) :: this
     581logical, intent(out) :: is_subsurface_ice
     582real,    intent(out) :: h2o_ice, h_toplag
     583
     584!---- Local variables
     585type(stratum), pointer :: current
     586
     587!---- Code
     588h_toplag = 0.
     589h2o_ice = 0.
     590is_subsurface_ice = .false.
     591current => this%top
     592! Is the considered stratum a dust lag?
     593if (is_dust_lag(current)) return
     594do while (is_dust_lag(current) .and. associated(current))
     595    h_toplag = h_toplag + thickness(current)
     596    current => current%down
     597enddo
     598! Is the stratum below the dust lag made of ice?
     599if (is_h2oice_str(current)) then
     600    if (h_toplag >= h_patchy_dust) then
     601        is_subsurface_ice = .true.
     602    else
     603        h_toplag = 0.
     604        h2o_ice = current%h_h2oice
     605    endif
     606endif
     607
     608END SUBROUTINE subsurface_ice_layering
    567609
    568610!=======================================================================
    … …  
    597639zlag = 0.
    598640
     641! NOTHING HAPPENS
     642if (abs(dh_co2ice) < tol .and. abs(dh_h2oice) < tol .and. abs(dh_dust) < tol) then
     643    ! So we do nothing
     644!-----------------------------------------------------------------------
    599645! DUST SEDIMENTATION with possibly ice condensation
    600 if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice .and. dh_co2ice >= 0. .and. dh_h2oice >= 0.) then
     646else if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice .and. dh_co2ice >= 0. .and. dh_h2oice >= 0.) then
    601647    write(*,*) '> Layering: dust sedimentation'
    602648    h_pore = dh_dust*dry_regolith_porosity/(1. - dry_regolith_porosity)
    … …  
    608654    else
    609655        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
     656        this%top%h_co2ice      = this%top%h_co2ice      + dh_co2ice
     657        this%top%h_h2oice      = this%top%h_h2oice      + dh_h2oice
     658        this%top%h_dust        = this%top%h_dust        + dh_dust
     659        this%top%h_pore        = this%top%h_pore        + h_pore
    614660    endif
    615661!-----------------------------------------------------------------------
    616 ! DUST LIFTING
    617 !(with possibly ice sublimation???)
     662! DUST LIFTING with possibly ice sublimation
    618663!else if (dh_dust < dh_co2ice .and. dh_dust < dh_h2oice .and. dh_co2ice <= 0. .and. dh_h2oice <= 0.) then
    619664else if (dh_dust < 0. .and. abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then
    … …  
    651696    else
    652697        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
     698        this%top%h_co2ice      = this%top%h_co2ice      + dh_co2ice
     699        this%top%h_h2oice      = this%top%h_h2oice      + dh_h2oice
     700        this%top%h_dust        = this%top%h_dust        + dh_dust
     701        this%top%h_pore        = this%top%h_pore        + h_pore
    657702    endif
    658703!-----------------------------------------------------------------------
    … …  
    710755
    711756!-----------------------------------------------------------------------
    712 ! NOTHING HAPPENS
    713 else if (abs(dh_co2ice) < tol .and. abs(dh_h2oice) < tol .and. abs(dh_dust) < tol) then
    714     ! We do nothing
    715 !-----------------------------------------------------------------------
    716757! UNUSUAL (WEIRD) SITUATION
    717758else
    … …  
    741782! Update of properties in the eroding dust lag
    742783str%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust)
    743 str%h_pore = str%h_pore - h2lift*str%h_pore/str%h_dust
    744 str%h_dust = str%h_dust - h2lift
     784str%h_pore        = str%h_pore        - h2lift*str%h_pore/str%h_dust
     785str%h_dust        = str%h_dust        - h2lift
    745786
    746787! Remove the eroding dust lag if there is no dust anymore
    … …  
    775816! Update of properties in the sublimating stratum
    776817str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hsubl_dust
    777 str%h_co2ice = str%h_co2ice - h2subl_co2ice
    778 str%h_h2oice = str%h_h2oice - h2subl_h2oice
    779 str%h_dust = str%h_dust - hsubl_dust
    780 str%h_pore = str%h_pore - h2subl*h_pore/h_ice
     818str%h_co2ice      = str%h_co2ice      - h2subl_co2ice
     819str%h_h2oice      = str%h_h2oice      - h2subl_h2oice
     820str%h_dust        = str%h_dust        - hsubl_dust
     821str%h_pore        = str%h_pore        - h2subl*h_pore/h_ice
    781822
    782823! Correct the value of top_elevation for the upper strata
Note: See TracChangeset for help on using the changeset viewer.