Changeset 3159


Ignore:
Timestamp:
Dec 14, 2023, 6:56:40 PM (11 months ago)
Author:
jbclement
Message:

PEM:

  • A number of corrections related to r3149 for CO2 ice management:

    'co2_ice' was not transferred to 'perennial_co2ice' at the end of the PEM run;
    In the Mars PCM, 'perennial_co2ice' is now correctly handled regarding the frost in case of CO2 sublimation/condensation.

  • Addition of (CO2 and H2O) ice metamorphism: if the frost is above a given value, then the excess frost is transferred to the perennial ice.
  • Thresholds value for ice management can now be set in the "run_PEM.def".

JBC

Location:
trunk
Files:
11 edited

Legend:

Unmodified
Added
Removed

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

    r3152 r3159  
    172172- The PEM deftank folder is now in LMDZ.COMMON/libf/evolution/ (not anymore in the Mars PCM deftank).
    173173- Small corrections related to r3149.
     174
     175== 14/12/2023 == JBC
     176- A number of corrections related to r3149 for CO2 ice management:
     177      > 'co2_ice' was not transferred to 'perennial_co2ice' at the end of the PEM run;
     178      > In the Mars PCM, 'perennial_co2ice' is now correctly handled regarding the frost in case of CO2 sublimation/condensation.
     179- Addition of (CO2 and H2O) ice metamorphism: if the frost is above a given value, then the excess frost is transferred to the perennial ice.
     180- Thresholds value for ice management can now be set in the "run_PEM.def".

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

    r3149 r3159  
    2323#endif
    2424
    25 use time_evol_mod,    only: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
     25use time_evol_mod,         only: year_bp_ini, dt_pem, h2o_ice_crit, co2_ice_crit, ps_criterion, Max_iter_pem, &
    2626                          evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
    27 use comsoil_h_pem,    only: soil_pem, fluxgeo, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, reg_thprop_dependp
    28 use adsorption_mod,   only: adsorption_pem
    29 use glaciers_mod,     only: co2_ice_flow, h2o_ice_flow
    30 use ice_table_mod,    only: icetable_equilibrium, icetable_dynamic
    31 use time_phylmdz_mod, only: ecritphy
     27use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, reg_thprop_dependp
     28use adsorption_mod,        only: adsorption_pem
     29use glaciers_mod,          only: co2_ice_flow, h2o_ice_flow
     30use ice_table_mod,         only: icetable_equilibrium, icetable_dynamic
     31use time_phylmdz_mod,      only: ecritphy
     32use constants_marspem_mod, only: inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold
    3233
    3334implicit none
    … …  
    5253close(100)
    5354
    54 ! #---------- Output parameters
     55!#---------- Output parameters ----------#
    5556! Frequency of outputs for the PEM
    5657ecritphy = 1 ! Default value: every year
    5758call getin('ecritpem',ecritphy)
    5859
    59 ! #---------- Orbital parameters
     60!#---------- Orbital parameters ----------#
    6061evol_orbit_pem = .false.
    6162call getin('evol_orbit_pem',evol_orbit_pem)
    … …  
    7778write(*,*) 'Does Ls peri vary ?',var_lsp
    7879
    79 ! #---------- Stopping criteria parameters
     80!#---------- Stopping criteria parameters ----------#
    8081Max_iter_pem = 100000000
    8182call getin('Max_iter_pem',Max_iter_pem)
    8283
    83 water_ice_criterion = 0.2
    84 call getin('water_ice_criterion',water_ice_criterion)
     84h2o_ice_crit = 0.2
     85call getin('h2o_ice_crit',h2o_ice_crit)
    8586
    86 co2_ice_criterion = 0.2
    87 call getin('co2_ice_criterion',co2_ice_criterion)
     87co2_ice_crit = 0.2
     88call getin('co2_ice_crit',co2_ice_crit)
    8889
    8990ps_criterion = 0.15
    … …  
    9394call getin('dt_pem', dt_pem)
    9495
    95 ! #---------- Subsurface parameters
     96!#---------- Subsurface parameters ----------#
    9697soil_pem = .true.
    9798call getin('soil_pem',soil_pem)
    … …  
    110111write(*,*)  'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
    111112
    112 ! #---------- Layering parameters
     113!#---------- Layering parameters ----------#
    113114fluxgeo = 0.
    114115call getin('fluxgeo',fluxgeo)
    … …  
    156157endif
    157158
     159!#---------- Ice management parameters ----------#
    158160ini_h2o_bigreservoir = 1.e4
    159161call getin('ini_h2o_bigreservoir',ini_h2o_bigreservoir)
     162
     163inf_h2oice_threshold = 2.e3
     164call getin('inf_h2oice_threshold',inf_h2oice_threshold)
     165
     166metam_h2oice_threshold = 5.e-2
     167call getin('metam_h2oice_threshold',metam_h2oice_threshold)
     168
     169metam_co2ice_threshold = 16.e3
     170call getin('metam_co2ice_threshold',metam_co2ice_threshold)
    160171
    161172END SUBROUTINE conf_pem

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

    r3149 r3159  
    3737real, parameter :: Lco2 =  5.71e5 ! Pilorget and Forget 2016
    3838
    39 ! Threshold to consider the amount of H2O ice as an infinite reservoir
    40 real, parameter :: inf_h2oice_threshold = 2000.  !~ 1 m
     39! Thresholds for ice management
     40real, save :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir
     41real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice
     42real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice
    4143
    4244END MODULE constants_marspem_mod

  • trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

    r3149 r3159  
    2828#  Max_iter_pem=100000000
    2929
    30 # Acceptance rate of sublimating water ice surface change? Default = 0.2
    31 # water_ice_criterion=0.2
     30# Acceptance rate of sublimating H2o ice surface change? Default = 0.2
     31# h2o_ice_crit=0.2
    3232
    3333# Acceptance rate of sublimating CO2 ice surface change? Default = 0.2
    34 # co2_ice_criterion=0.2
     34# co2_ice_crit=0.2
    3535
    3636# Acceptance rate of pressure surface change? Default = 0.15
    … …  
    6969# icetable_dynamic=.false.
    7070
    71 #---------- H2O reservoir parameters ----------#
    72 # Amount of h2o ice to initialize the big reservoir if the variable is not present in startfi_PEM.nc? Default = 1.e4
    73 # ini_h2o_bigreservoir=1e4
     71#---------- Ice management parameters ----------#
     72# Amount of h2o ice to initialize the big reservoir if the variable is not present in "startfi_PEM.nc"? Default = 1.e4
     73# ini_h2o_bigreservoir=1.e4
     74
     75# Threshold to consider the amount of H2O ice as an infinite reservoir? Default = 2.e3
     76# inf_h2oice_threshold=2.e3
     77
     78# Threshold to consider frost is becoming perennial H2O ice? Default = 5.e-2
     79# metam_h2oice_threshold=5.e-2
     80
     81# Threshold to consider frost is becoming perennial CO2 ice? Default = 16.e3
     82# metam_co2ice_threshold=16.e3
    7483
    7584# Some definitions for the physics, in file 'callphys.def'

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

    r3149 r3159  
    4242use stopping_crit_mod,           only: stopping_crit_h2o_ice, stopping_crit_co2
    4343use constants_marspem_mod,       only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
    44                                        m_noco2, inf_h2oice_threshold
     44                                       m_noco2, inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold
    4545use evol_ice_mod,                only: evol_co2_ice, evol_h2o_ice
    4646use comsoil_h_PEM,               only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
    … …  
    189189real, dimension(:,:,:), allocatable :: co2_ice_PCM        ! Physics x NSLOPE x Times field: co2 ice given by the PCM [kg/m^2]
    190190real, dimension(:,:),   allocatable :: co2_ice_ini_sublim ! physical point field: Logical array indicating sublimating point of co2 ice
    191 real, dimension(:,:),   allocatable :: initial_h2o_ice    ! physical point field: Logical array indicating if there is water ice at initial state
     191real, dimension(:,:),   allocatable :: initial_h2o_ice    ! physical point field: Logical array indicating if there is h2o ice at initial state
    192192real, dimension(:,:),   allocatable :: tend_co2_ice       ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
    193193real, dimension(:,:),   allocatable :: tend_co2_ice_ini   ! physical point x slope field x nslope: Tendency of evolution of perennial co2 ice over a year in the PCM
    … …  
    550550!    I_g Save initial PCM situation
    551551!------------------------
    552 ! We save the places where water ice is sublimating
    553 ! We compute the surface of water ice sublimating
     552! We save the places where h2o ice is sublimating
     553! We compute the surface of h2o ice sublimating
    554554co2_surf_ini = 0.
    555555h2o_surf_ini = 0.
    … …  
    848848! II_d.5 Update the mass of the regolith adsorbed
    849849        if (adsorption_pem) then
    850             call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice,       &
     850            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice,                   &
    851851                                     h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
    852852                                     h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
    … …  
    858858                    do l = 1,nsoilmx_PEM - 1
    859859                        totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
    860                         subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
     860                                               subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    861861                        totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
    862                         subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
     862                                               subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    863863                    enddo
    864864                enddo
    … …  
    907907        select case (stopPEM)
    908908            case(1)
    909                 write(*,*) "STOPPING because surface of water ice sublimating is too low:", stopPEM, "See message above."
     909                write(*,*) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above."
    910910            case(2)
    911                 write(*,*) "STOPPING because tendencies on water ice = 0:", stopPEM, "See message above."
     911                write(*,*) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above."
    912912            case(3)
    913913                write(*,*) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
    … …  
    940940! III_a.1 Ice update (for startfi)
    941941
    942 ! H2O ice
    943942watercap = 0.
     943perennial_co2ice = co2_ice
    944944do ig = 1,ngrid
     945    ! H2O ice metamorphism
     946    if (sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
     947        h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
     948        qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
     949    endif
     950
     951    ! Is H2O ice still considered as an infinite reservoir for the PCM?
    945952    if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
    946     ! If there is enough ice to be considered as an infinite reservoir
     953        ! There is enough ice to be considered as an infinite reservoir
    947954        watercaptag(ig) = .true.
    948955    else
    949     ! If there too little ice to be considered as an infinite reservoir
    950     ! Ice is transferred to the frost
     956        ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
    951957        watercaptag(ig) = .false.
    952958        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
    953959        h2o_ice(ig,:) = 0.
     960    endif
     961
     962    ! CO2 ice metamorphism
     963    if (sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
     964        perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
     965        qsurf(ig,igcm_co2,:) = metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
    954966    endif
    955967enddo

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

    r3149 r3159  
    8282#ifdef CPP_STD
    8383    logical, dimension(ngrid) :: watercaptag
    84     watercaptag(:) = .false.
     84    watercaptag = .false.
    8585#endif
    8686
    … …  
    152152                            TI_PEM(ig,iloop,islope) = TI_breccia
    153153                        enddo
    154                     else ! we keep the high ti values
     154                    else ! we keep the high TI values
    155155                        do iloop = index_breccia + 1,index_bedrock
    156156                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)

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

    r3149 r3159  
    1616SUBROUTINE stopping_crit_h2o_ice(cell_area,surf_ini,h2o_ice,stopPEM,ngrid)
    1717
    18 use time_evol_mod, only: water_ice_criterion
     18use time_evol_mod, only: h2o_ice_crit
    1919use comslope_mod,  only: subslope_dist, nslope
    2020
    … …  
    3232integer,                       intent(in) :: ngrid     ! # of physical grid points
    3333real, dimension(ngrid),        intent(in) :: cell_area ! Area of the cells
    34 real, dimension(ngrid,nslope), intent(in) :: h2o_ice   ! Actual density of water ice
     34real, dimension(ngrid,nslope), intent(in) :: h2o_ice   ! Actual density of h2o ice
    3535real,                          intent(in) :: surf_ini  ! Initial surface of h2o ice that was sublimating
    3636!   OUTPUT
    … …  
    5151
    5252! Check of the criterion
    53 if (surf_now < surf_ini*(1. - water_ice_criterion)) then
     53if (surf_now < surf_ini*(1. - h2o_ice_crit)) then
    5454    stopPEM = 1
    55     write(*,*) "Reason of stopping: the surface of water ice sublimating reaches the threshold"
    56     write(*,*) "surf_now < surf_ini*(1. - water_ice_criterion)", surf_now < surf_ini*(1. - water_ice_criterion)
    57     write(*,*) "Current surface of water ice sublimating =", surf_now
    58     write(*,*) "Initial surface of water ice sublimating =", surf_ini
    59     write(*,*) "Percentage of change accepted =", water_ice_criterion*100
    60 else if (surf_now > surf_ini*(1. + water_ice_criterion)) then
     55    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
     56    write(*,*) "surf_now < surf_ini*(1. - h2o_ice_crit)", surf_now < surf_ini*(1. - h2o_ice_crit)
     57    write(*,*) "Current surface of h2o ice sublimating =", surf_now
     58    write(*,*) "Initial surface of h2o ice sublimating =", surf_ini
     59    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
     60else if (surf_now > surf_ini*(1. + h2o_ice_crit)) then
    6161    stopPEM = 1
    62     write(*,*) "Reason of stopping: the surface of water ice sublimating reaches the threshold"
    63     write(*,*) "surf_now > surf_ini*(1. + water_ice_criterion)", surf_now > surf_ini*(1. + water_ice_criterion)
    64     write(*,*) "Current surface of water ice sublimating =", surf_now
    65     write(*,*) "Initial surface of water ice sublimating =", surf_ini
    66     write(*,*) "Percentage of change accepted =", water_ice_criterion*100
     62    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
     63    write(*,*) "surf_now > surf_ini*(1. + h2o_ice_crit)", surf_now > surf_ini*(1. + h2o_ice_crit)
     64    write(*,*) "Current surface of h2o ice sublimating =", surf_now
     65    write(*,*) "Initial surface of h2o ice sublimating =", surf_ini
     66    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
    6767endif
    6868
    … …  
    7575SUBROUTINE stopping_crit_co2(cell_area,surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
    7676
    77 use time_evol_mod, only: co2_ice_criterion, ps_criterion
     77use time_evol_mod, only: co2_ice_crit, ps_criterion
    7878use comslope_mod,  only: subslope_dist
    7979
    … …  
    9191integer,                       intent(in) :: ngrid, nslope        ! # of grid physical grid points
    9292real, dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
    93 real, dimension(ngrid,nslope), intent(in) :: co2_ice              ! Actual density of water ice
     93real, dimension(ngrid,nslope), intent(in) :: co2_ice              ! Actual density of h2o ice
    9494real,                          intent(in) :: surf_ini             ! Initial surface of co2 ice that was sublimating
    9595real,                          intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files
    … …  
    113113
    114114! Check of the criterion
    115 if (surf_now < surf_ini*(1. - co2_ice_criterion)) then
     115if (surf_now < surf_ini*(1. - co2_ice_crit)) then
    116116    stopPEM = 3
    117117    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
    118     write(*,*) "surf_now < surf_ini*(1. - co2_ice_criterion)", surf_now < surf_ini*(1. - co2_ice_criterion)
     118    write(*,*) "surf_now < surf_ini*(1. - co2_ice_crit)", surf_now < surf_ini*(1. - co2_ice_crit)
    119119    write(*,*) "Current surface of co2 ice sublimating =", surf_now
    120120    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
    121     write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
    122 else if (surf_now > surf_ini*(1. + co2_ice_criterion)) then
     121    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
     122else if (surf_now > surf_ini*(1. + co2_ice_crit)) then
    123123    stopPEM = 3
    124124    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
    125     write(*,*) "surf_now > surf_ini*(1. + co2_ice_criterion)", surf_now > surf_ini*(1. + co2_ice_criterion)
     125    write(*,*) "surf_now > surf_ini*(1. + co2_ice_crit)", surf_now > surf_ini*(1. + co2_ice_crit)
    126126    write(*,*) "Current surface of co2 ice sublimating =", surf_now
    127127    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
    128     write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
     128    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
    129129endif
    130130

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

    r3149 r3159  
    33implicit none
    44
    5 integer :: year_bp_ini         ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
    6 integer :: dt_pem              ! Time step used by the PEM in Planetary years
    7 real    :: convert_years       ! Conversion ratio from Planetary years to Earth years
    8 real    :: water_ice_criterion ! Percentage of change of the surface of water ice sublimating before stopping the PEM
    9 real    :: co2_ice_criterion   ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
    10 real    :: ps_criterion        ! Percentage of change of averaged surface pressure before stopping the PEM
    11 integer :: Max_iter_pem        ! Maximal number of iteration when converging to a steady state, read in evol.def
    12 logical :: evol_orbit_pem      ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
    13 logical :: var_obl             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
    14 logical :: var_ecc             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for eccenticity
    15 logical :: var_lsp             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
     5integer :: year_bp_ini    ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
     6integer :: dt_pem         ! Time step used by the PEM in Planetary years
     7real    :: convert_years  ! Conversion ratio from Planetary years to Earth years
     8real    :: h2o_ice_crit   ! Percentage of change of the surface of h2o ice sublimating before stopping the PEM
     9real    :: co2_ice_crit   ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
     10real    :: ps_criterion   ! Percentage of change of averaged surface pressure before stopping the PEM
     11integer :: Max_iter_pem   ! Maximal number of iteration when converging to a steady state, read in evol.def
     12logical :: evol_orbit_pem ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
     13logical :: var_obl        ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
     14logical :: var_ecc        ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for eccenticity
     15logical :: var_lsp        ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
    1616
    1717END MODULE time_evol_mod

  • trunk/LMDZ.MARS/libf/phymars/albedocaps.F90

    r3130 r3159  
    105105      psolaralb(ig,2)=psolaralb(ig,1)
    106106    else
    107       psolaralb(ig,1)=albedice(icap)
    108       psolaralb(ig,2)=albedice(icap)
    109       if(paleoclimate) then
    110          if((piceco2_peren(ig).gt.0.).and.(piceco2(ig).lt.piceco2_peren(ig))) then
    111              psolaralb(ig,1) = albedo_perennialco2
    112              psolaralb(ig,2) = albedo_perennialco2
    113              piceco2_peren(ig) = piceco2(ig)
    114          endif
    115       endif
     107      psolaralb(ig,:)=albedice(icap)
     108      if (paleoclimate .and. piceco2_peren(ig) > 0. .and. abs(piceco2(ig)) < 1.e-10) psolaralb(ig,:) = albedo_perennialco2
    116109    endif
    117110  else if (watercaptag(ig) .and. water) then

  • trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F

    r3142 r3159  
    55      logical, save :: scavco2cond = .false. ! flag for using scavenging_by_co2
    66!$OMP THREADPRIVATE(scavco2cond)
    7      
     7
    88      CONTAINS
    99
    … …  
    186186      REAL   :: zcondices_tmp(ngrid)          ! temporary condensation rate [kg/m^2/s]
    187187      REAL   :: piceco2_tmp(ngrid)            ! temporary amount of CO2 frost on the surface [kg/m^2]
    188       REAL   :: perennial_co2ice_tmp(ngrid)    ! temporary amount of perennial CO2 frost on the surface [kg/m^2]
     188      REAL   :: perennial_co2ice_tmp(ngrid)   ! temporary amount of perennial CO2 frost on the surface [kg/m^2]
    189189      LOGICAL :: condsub_tmp(ngrid)           ! Boolean to check if CO2 ice is condensing / sublimating on the sub grid surface [1]
    190190      REAL :: zfallice_tmp(ngrid,nlayer+1)    ! temporary amount of ice falling from layer l for a specific sub-grid surface [kg/m^2/s]
    191191      REAL :: condens_layer_tmp(ngrid,nlayer) ! temporary condensation rate in layer l (co2 cloud) for a specific sub-grid surface [kg/m^2/s]
    192192      INTEGER :: islope                       ! index for loop variables
    193       REAL :: pdqsc_tmp(ngrid,nq)         ! tendency on surface tracers (grid-mesh average) after scavenging
     193      REAL :: pdqsc_tmp(ngrid,nq)             ! tendency on surface tracers (grid-mesh average) after scavenging
    194194c----------------------------------------------------------------------
    195195
    … …  
    227227      ENDIF ! of IF (firstcall)
    228228      zcpi=1./cpp
     229      if (paleoclimate) piceco2 = piceco2 + perennial_co2ice
    229230
    230231c
    … …  
    564565c  ********************************************************************
    565566
     567!     Redistribute piceco2 into piceco2 (frost) and perennial_co2ice
     568!     --------------------------------------------------------------
     569      if (paleoclimate) then
     570          where (piceco2 > perennial_co2ice) ! Perennial co2 ice has not been affected
     571              ! It means:
     572              !     - In case of sublimation, only frost is lost
     573              !     - In case of condensation, only frost accumulates new ice
     574              piceco2 = piceco2 - perennial_co2ice
     575          else where ! Perennial co2 ice has been affected
     576              ! It means that frost disappeared with sublimation and perennial ice is being lost
     577              perennial_co2ice = piceco2
     578              piceco2 = 0.
     579          end where
     580      endif
     581
    566582!     Check that amont of CO2 ice is not problematic
    567       DO ig=1,ngrid
     583!     ----------------------------------------------
     584      DO ig = 1,ngrid
    568585         DO islope = 1,nslope
    569            if(.not.piceco2(ig,islope).ge.0.) THEN
    570               if(piceco2(ig,islope).le.-5.e-8) print*,
    571      $         'WARNING co2condens piceco2(',ig,')=', piceco2(ig,islope)
    572                piceco2(ig,islope)=0.
     586           if (piceco2(ig,islope) < 0.) then
     587              if (piceco2(ig,islope) <= -5.e-8) print*,
     588     $         'WARNING co2condens piceco2(',ig,',',islope,
     589     $         ') =',piceco2(ig,islope)
     590               piceco2(ig,islope) = 0.
    573591           endif
    574592        ENDDO
    575593      ENDDO
     594      if (paleoclimate) then
     595      DO ig = 1,ngrid
     596         DO islope = 1,nslope
     597           if (perennial_co2ice(ig,islope) < 0.) then
     598              if (perennial_co2ice(ig,islope) <= -5.e-8) print*,
     599     $         'WARNING co2condens perennial_co2ice(',ig,',',islope,
     600     $         ') =',perennial_co2ice(ig,islope)
     601               perennial_co2ice(ig,islope) = 0.
     602           endif
     603        ENDDO
     604      ENDDO
     605      endif
    576606     
    577607!     Set albedo and emissivity of the surface
    578608!     ----------------------------------------
    579609      DO islope = 1,nslope
    580         piceco2_tmp(:) = piceco2(:,islope)
    581         alb_tmp(:,:) = psolaralb(:,:,islope)
    582         emisref_tmp(:) = 0.
    583         perennial_co2ice_tmp(:) =  perennial_co2ice(:,islope)
     610        piceco2_tmp = piceco2(:,islope)
     611        alb_tmp = psolaralb(:,:,islope)
     612        emisref_tmp = 0.
     613        perennial_co2ice_tmp =  perennial_co2ice(:,islope)
    584614        CALL albedocaps(zls,ngrid,piceco2_tmp,perennial_co2ice_tmp,
    585615     &                  alb_tmp,emisref_tmp)
    586         perennial_co2ice(:,islope) = perennial_co2ice_tmp(:)
     616        perennial_co2ice(:,islope) = perennial_co2ice_tmp
    587617        psolaralb(:,1,islope) =  alb_tmp(:,1)
    588618        psolaralb(:,2,islope) =  alb_tmp(:,2)
    589         emisref(:,islope) = emisref_tmp(:)
     619        emisref(:,islope) = emisref_tmp
    590620      ENDDO
    591621
    592 ! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
    593       DO ig=1,ngrid
    594         DO islope = 1,nslope
    595           if (piceco2(ig,islope).eq.0) then
    596             pemisurf(ig,islope)=emissiv
    597           endif
    598         ENDDO
    599       ENDDO
     622!     Set pemisurf() to emissiv when there is bare surface (needed for co2snow)
     623!     -------------------------------------------------------------------------
     624      where (piceco2 == 0.) pemisurf = emissiv
     625
    600626
    601627!         firstcall2=.false.

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r3157 r3159  
    24732473      IF (refill_watercap) THEN
    24742474
    2475         DO ig=1,ngrid
     2475        DO ig = 1,ngrid
    24762476         DO islope = 1,nslope
    2477            if (watercaptag(ig).and. (qsurf(ig,igcm_h2o_ice,islope)
    2478      &        .gt.frost_metam_threshold)) then
    2479 
    2480                 watercap(ig,islope)=watercap(ig,islope)
    2481      &          +qsurf(ig,igcm_h2o_ice,islope)
    2482      &         - frost_metam_threshold
     2477           if (watercaptag(ig) .and. (qsurf(ig,igcm_h2o_ice,islope)
     2478     &        > frost_metam_threshold)) then
     2479
     2480                watercap(ig,islope) = watercap(ig,islope)
     2481     &          + qsurf(ig,igcm_h2o_ice,islope)
     2482     &          - frost_metam_threshold
    24832483                qsurf(ig,igcm_h2o_ice,islope) = frost_metam_threshold
    2484            endif ! (watercaptag(ig).and.
     2484           endif ! watercaptag
    24852485          ENDDO
    24862486        ENDDO
    24872487
    2488       ENDIF ! (refill_watercap) THEN
     2488      ENDIF ! refill_watercap
    24892489
    24902490
Note: See TracChangeset for help on using the changeset viewer.