Ignore:
Timestamp:
May 26, 2025, 10:34:57 AM (3 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

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

Legend:

Unmodified
Added
Removed

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

    r3777 r3778  
    663663== 23/05/2025 == JBC
    664664Big improvements of Python scripts to output the stratifications (user-friendly prompt, error checks, code modularity, nice visualization, etc).
     665
     666== 26/05/2025 == JBC
     667- New subroutine to detect whether there is subsurface ice or not
     668- Rework of the initialization/update/finalization of the situation regarding the layering data structure
     669- Introduction of a threshold 'h_patchy_dust' under which the top dust layer is not considered as a stratum
     670- 'deposits' is renamed as 'layerings_map'
     671- Few cleanings

  • 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

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

    r3770 r3778  
    6969use writediagpem_mod,           only: writediagpem, writediagsoilpem
    7070use co2condens_mod,             only: CO2cond_ps
    71 use 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
     71use 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
    7374use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
    7475use version_info_mod,           only: print_version_info
    … …  
    163164real                            :: total_surface ! Total surface of the planet
    164165
    165 ! Variables for h2o_ice evolution
     166! Variables for h2o ice evolution
    166167real, dimension(:,:),    allocatable  :: h2o_ice              ! h2o ice in the PEM
    167168real, dimension(:,:),    allocatable  :: d_h2oice             ! physical point x slope field: Tendency of evolution of perennial h2o ice
    … …  
    182183real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
    183184
    184 ! Variables for co2_ice evolution
     185! Variables for co2 ice evolution
    185186real,    dimension(:,:), allocatable :: co2_ice                ! co2 ice in the PEM
    186187real,    dimension(:,:), allocatable :: d_co2ice               ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
    … …  
    194195real,    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]
    195196
    196 ! Variables for the evolution of layered deposits
    197 type(layering), dimension(:,:), allocatable :: deposits          ! Layering for each grid point and slope
     197! Variables for the evolution of layered layerings_map
     198type(layering), dimension(:,:), allocatable :: layerings_map     ! Layering for each grid point and slope
    198199type(ptrarray), dimension(:,:), allocatable :: current           ! Current active stratum in the layering
    199200logical,        dimension(:,:), allocatable :: new_str, new_lag  ! Flags for the layering algorithm
    200201real                                        :: h2o_ice_depth_old ! Old depth of subsurface ice layer
     202logical                                     :: is_subsurface_ice ! Boolean to know if there is subsurface ice
    201203
    202204! Variables for slopes
    … …  
    645647!------------------------
    646648write(*,*) '> Reading "startpem.nc"'
    647 allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),deposits(ngrid,nslope))
     649allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),layerings_map(ngrid,nslope))
    648650allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid))
    649651delta_h2o_icetablesublim = 0.
    650652call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
    651653              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,            &
    652               watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,deposits)
     654              watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,layerings_map)
    653655deallocate(tsurf_avg_yr1)
    654656
    … …  
    658660co2ice_sublim_surf_ini = 0.
    659661h2oice_ini_surf = 0.
    660 h2o_ice_depth = 0.
    661662is_co2ice_sublim_ini = .false.
    662663is_h2oice_sublim_ini = .false.
    … …  
    666667    do ig = 1,ngrid
    667668        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
     669            co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
     670            h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
    670671        enddo
    671672    enddo
    … …  
    681682            is_h2oice_sublim_ini(i,islope) = .true.
    682683            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))
    684684        endif
    685685    enddo
    … …  
    687687write(*,*) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini
    688688write(*,*) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf
     689
     690! We put the sublimating tendency coming from subsurface ice
     691where (zdqsdif_ssi_tot > 0.)
     692    d_h2oice = -zdqsdif_ssi_tot
     693end where
    689694
    690695if (adsorption_pem) then
    … …  
    742747    do islope = 1,nslope
    743748        do ig = 1,ngrid
    744             current(ig,islope)%p => deposits(ig,islope)%top
     749            current(ig,islope)%p => layerings_map(ig,islope)%top
    745750        enddo
    746751    enddo
    … …  
    822827        do islope = 1,nslope
    823828            do ig = 1,ngrid
    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
     829                call make_layering(layerings_map(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)
     830                co2_ice = 0.
     831                h2o_ice = 0.
     832                h2o_ice_depth = 0.
     833                if (is_co2ice_str(layerings_map(ig,islope)%top)) then
     834                    co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
     835                else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
     836                    h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
     837                else
     838                    call subsurface_ice_layering(layerings_map(ig,islope),is_subsurface_ice,h2o_ice_depth(ig,islope),h2o_ice(ig,islope))
     839                endif
    827840            enddo
    828841        enddo
    … …  
    844857        do islope = 1,nslope
    845858            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)
     859                layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
     860                layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
    848861            enddo
    849862        enddo
    … …  
    10281041                if (is_h2oice_sublim_ini(ig,islope)) then
    10291042                    h2o_ice_depth_old = h2o_ice_depth(ig,islope)
    1030                     h2o_ice_depth(ig,islope) = thickness_toplag(deposits(ig,islope))
     1043                    h2o_ice_depth(ig,islope) = thickness_toplag(layerings_map(ig,islope))
    10311044                    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))
    10321045                endif
    … …  
    10971110!    III_a Update surface values for the PCM start files
    10981111!------------------------
    1099 ! III_a.1 Ice update for start file
    11001112write(*,*)
    11011113write(*,*) '********* PEM finalization *********'
     1114! III_a.1 Ice update for start file
    11021115write(*,*) '> Reconstructing perennial ice and frost for the PCM'
    11031116watercap = 0.
    … …  
    11271140    !endif
    11281141enddo
    1129 
    1130 ! III_a.2 Subsurface-surface interaction update for start file
    1131 zdqsdif_ssi_tot = 0.
    1132 h2o_ice_depth = 0.
    1133 if (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
    1143 endif
    11441142
    11451143! III.a.3. Tsurf update for start file
    … …  
    13131311!------------------------
    13141312write(*,*) '> Writing "restartpem.nc"'
    1315 if (layering_algo) nb_str_max = get_nb_str_max(deposits,ngrid,nslope) ! Get the maximum number of "stratum" in the depositsication (layerings)
     1313if (layering_algo) nb_str_max = get_nb_str_max(layerings_map,ngrid,nslope) ! Get the maximum number of "stratum" in the layerings_mapication (layerings)
    13161314call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
    13171315call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
    1318              co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,deposits)
     1316             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,layerings_map)
    13191317
    13201318call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
    … …  
    13311329    do islope = 1,nslope
    13321330        do i = 1,ngrid
    1333             call del_layering(deposits(i,islope))
     1331            call del_layering(layerings_map(i,islope))
    13341332        enddo
    13351333    enddo
    13361334endif
    13371335deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
    1338 deallocate(co2_ice,h2o_ice,deposits)
     1336deallocate(co2_ice,h2o_ice,layerings_map)
    13391337!----------------------------- END OUTPUT ------------------------------
    13401338

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

    r3770 r3778  
    99SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
    1010                    tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global,d_h2oice,d_co2ice,co2_ice,h2o_ice,                         &
    11                     watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
     11                    watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,m_h2o_regolith_phys,deltam_h2o_regolith_phys,layerings_map)
    1212
    1313use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
    … …  
    5454real, dimension(ngrid,nslope),           intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
    5555real, dimension(ngrid,nslope),           intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
    56 type(layering), dimension(ngrid,nslope), intent(inout) :: stratif             ! stratification (layerings)
     56type(layering), dimension(ngrid,nslope), intent(inout) :: layerings_map             ! Layerings
    5757real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
    5858real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
    … …  
    7171character(2)                            :: num               ! intermediate string to read PEM start sloped variables
    7272logical                                 :: startpem_file     ! boolean to check if we read the startfile or not
    73 real, dimension(:,:,:,:), allocatable   :: stratif_array     ! Array for stratification (layerings)
     73real, dimension(:,:,:,:), allocatable   :: stratif_array     ! Array for layerings
    7474
    7575#ifdef CPP_STD
    … …  
    132132
    133133!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    134     ! Stratification (layerings)
     134    ! Layerings
    135135    nb_str_max = 68
    136136    if (layering_algo) then
    … …  
    178178        enddo ! islope
    179179        if (.not. found) then
    180             write(*,*) 'So the deposits are initialized with sub-surface strata.'
     180            write(*,*) 'So the layerings are initialized with sub-surface strata.'
    181181            write(*,*) 'Ice is added with ''ini_huge_h2oice'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.'
    182182            do ig = 1,ngrid
    183183                if (watercaptag(ig)) then
    184184                    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.)
     185                        call ini_layering(layerings_map(ig,islope))
     186                        call add_stratum(layerings_map(ig,islope),ini_huge_h2oice,0.,ini_huge_h2oice,0.,0.,0.)
    187187                    enddo
    188188                else
    189189                    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.)
     190                        call ini_layering(layerings_map(ig,islope))
     191                        if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),perennial_co2ice(ig,islope),perennial_co2ice(ig,islope),0.,0.,0.,0.)
    192192                    enddo
    193193                endif
    194194            enddo
    195195        else
    196             call array2stratif(stratif_array,ngrid,nslope,stratif)
     196            call array2stratif(stratif_array,ngrid,nslope,layerings_map)
    197197        endif
    198198        deallocate(stratif_array)
    … …  
    383383    co2_ice = perennial_co2ice
    384384
    385     ! Stratification (layerings)
     385    ! Layerings
    386386    nb_str_max = 68
    387387    if (layering_algo) then
    388         write(*,*)'So the deposits are initialized with sub-surface strata.'
     388        write(*,*)'So the layerings are initialized with sub-surface strata.'
    389389        write(*,*)'Ice is added with ''ini_huge_h2oice'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.'
    390390        do ig = 1,ngrid
    391391            if (watercaptag(ig)) then
    392392                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.)
     393                    call ini_layering(layerings_map(ig,islope))
     394                    call add_stratum(layerings_map(ig,islope),ini_huge_h2oice,0.,ini_huge_h2oice,0.,0.,0.)
    395395                enddo
    396396            else
    397397                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.)
     398                    call ini_layering(layerings_map(ig,islope))
     399                    if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),perennial_co2ice(ig,islope),perennial_co2ice(ig,islope),0.,0.,0.,0.)
    400400                enddo
    401401            endif

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

    r3770 r3778  
    5858
    5959SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, &
    60                    icetable_depth,icetable_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif)
     60                   icetable_depth,icetable_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,layerings_map)
    6161
    6262! write time-dependent variable to restart file
    … …  
    7878real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith
    7979real, dimension(ngrid,nslope),           intent(in) :: h2o_ice
    80 type(layering), dimension(ngrid,nslope), intent(in) :: stratif ! Stratification (layerings)
     80type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map ! Layerings
    8181
    8282integer                               :: islope
    … …  
    9696if (layering_algo) then
    9797    allocate(stratif_array(ngrid,nslope,nb_str_max,6))
    98     call stratif2array(stratif,ngrid,nslope,stratif_array)
     98    call stratif2array(layerings_map,ngrid,nslope,stratif_array)
    9999    do islope = 1,nslope
    100100        write(num,fmt='(i2.2)') islope
Note: See TracChangeset for help on using the changeset viewer.