Ignore:
Timestamp:
May 15, 2024, 4:45:14 PM (7 months ago)
Author:
jbclement
Message:

PEM:

  • Update of some default parameters;
  • Correction of a bug when the ice table was too low for the PEM subsurface discretization;
  • Correction of the computation for the change of sublimating ice necessary to the stopping criteria (the mistake was introduced in r3149) + simplification of the algorithm + renaming of variables with more explicit names;
  • Cleaning of "soil_thermalproperties_mod.F90".

JBC

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

Legend:

Unmodified
Added
Removed

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

    r3324 r3327  
    300300Correction of a bug when setting the albedo / emissivity: it should be the PCM
    301301which sets the emissivity / albedo of the CO2 ice (because of CO2 snowfall) and
    302 not the PEM. I thus removed the lines in the PEM accordingly
     302not the PEM. I thus removed the lines in the PEM accordingly
     303
     304== 15/05/2024 == JBC
     305- Update of some default parameters;
     306- Correction of a bug when the ice table was too low for the PEM subsurface discretization;
     307- Correction of the computation for the change of sublimating ice necessary to the stopping criteria (the mistake was introduced in r3149) + simplification of the algorithm + renaming of variables with more explicit names;
     308- Cleaning of "soil_thermalproperties_mod.F90".

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

    r3320 r3327  
    1212#---------- Orbital parameters ----------#
    1313# Do you want to follow an orbital forcing read in "obl_ecc_lsp.asc"? Default = .false.
    14   evol_orbit_pem=.true.
     14# evol_orbit_pem=.false.
    1515
    1616# If evol_orbit_pem=.true., number of Earth years before present to start the PEM run? Default = 0
    17   year_earth_bp_ini=-1000000
     17# year_earth_bp_ini=0
    1818
    1919# Do you want to vary the obliquity when following "obl_ecc_lsp.asc"? Default = .true.
    20   var_obl=.true.
     20# var_obl=.false.
    2121
    2222# Do you want to vary the eccentricity when following "obl_ecc_lsp.asc"? Default = .true.
    23   var_ecc=.true.
     23# var_ecc=.false.
    2424
    2525# Do you want to vary the ls perihelie when following "obl_ecc_lsp.asc"? Default = .true.
    26   var_lsp=.true.
     26# var_lsp=.false.
    2727
    2828#---------- Stopping criteria parameters ----------#
    2929# If evol_orbit_pem=.false., maximal number of iterations if no stopping criterion is reached? Default=100000000
    30 #  Max_iter_pem=100000000
     30# Max_iter_pem=10
    3131
    3232# Acceptance rate of sublimating H2O ice surface change? Default = 0.2
    … …  
    4646# soil_pem=.true.
    4747
    48 # Do you want to run with adsoprtion in the PEM? Default = .true.
     48# Do you want to run with adsoprtion in the PEM? Default = .false.
    4949# adsorption_pem=.true.
    5050
    … …  
    9090
    9191# Do you want the CO2 ice to flow along subslope inside a cell? Default = .true.
    92 # co2ice_flow=.true.
     92# co2ice_flow=.false.
    9393
    9494#---------- Layering parameters ----------#

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

    r3321 r3327  
    1616
    1717!-----------------------------------------------------------------------
    18   contains
     18contains
    1919!-----------------------------------------------------------------------
    2020SUBROUTINE ini_ice_table_porefilling(ngrid,nslope)
    … …  
    5858implicit none
    5959
    60 !   Inputs
    61 ! --------
     60! Inputs
     61! ------
    6262integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM
    6363logical, dimension(ngrid),               intent(in) :: watercaptag              ! Boolean to check the presence of a perennial glacier
    … …  
    6565real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: rhowatersoil_ave         ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
    6666real, dimension(ngrid,nslope),           intent(in) :: regolith_inertia         ! Thermal inertia of the regolith layer [SI]
    67 !   Ouputs
    68 ! --------
     67! Ouputs
     68! ------
    6969real, dimension(ngrid,nslope), intent(out) :: ice_table_beg       ! ice table depth [m]
    7070real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m]
    71 !   Locals
    72 ! --------
     71! Locals
     72! ------
    7373integer                       :: ig, islope, isoil, isoilend ! loop variables
    7474real, dimension(nsoil_PEM)    :: diff_rho                    ! difference of water vapor density between the surface and at depth [kg/m^3]
    … …  
    7777real                          :: stretch                     ! stretch factor to improve the convergence of the ice table
    7878real                          :: wice_inertia                ! Water Ice thermal Inertia [USI]
    79 !   Code
    80 ! ------
     79! Code
     80! ----
    8181previous_icetable_depth = ice_table_beg
    8282do ig = 1,ngrid
    … …  
    147147use comslope_mod,          only: subslope_dist, def_slope_mean
    148148use constants_marspem_mod, only: porosity
    149 use vdifc_mod,             only: compute_Tice
    150149#ifndef CPP_STD
    151150    use comcstfi_h,   only: pi
    … …  
    156155implicit none
    157156
    158 !   Inputs
    159 ! --------
     157! Inputs
     158! ------
    160159integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
    161160real, dimension(ngrid,nslope),           intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m]
    … …  
    164163real, dimension(ngrid,nslope),           intent(in) :: tsurf                      ! Surface temperature [K]
    165164real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil                      ! Soil temperature [K]
    166 !   Outputs
    167 ! ---------
     165! Outputs
     166! -------
    168167real, dimension(ngrid), intent(out) :: delta_m_h2o ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]
    169 !   Locals
    170 !---------
     168! Locals
     169!-------
    171170integer                       :: ig, islope, ilay, iref ! loop index
    172171real, dimension(ngrid,nslope) :: rho                    ! density of water ice [kg/m^3]
    173172real, dimension(ngrid,nslope) :: Tice                   ! ice temperature [k]
    174 !   Code
    175 ! ------
     173! Code
     174! ----
    176175rho = 0.
    177176Tice = 0.
    … …  
    179178do ig = 1,ngrid
    180179    do islope = 1,nslope
    181         call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope))
     180        print*, ig, islope, ice_table_depth(ig,islope)
     181        call compute_Tice_pem(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope))
    182182        rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012
    183183    enddo
    … …  
    224224real,                         intent(out) :: depth_geothermal ! depth where the ice table stops [m]
    225225real, dimension(nsoilmx_PEM), intent(out) :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
    226 !   Locals
    227 !---------
     226! Locals
     227!-------
    228228real, dimension(nsoilmx_PEM) :: eta                          ! constriction
    229229real, dimension(nsoilmx_PEM) :: dz_psat                      ! first derivative of the vapor pressure at saturation
    … …  
    237237integer                      :: index_firstice               ! first index where ice appears (i.e., f > 0)
    238238real                         :: massfillabove, massfillafter ! h2O mass above and after index_geothermal
    239 !   Constant
    240 !-----------
     239! Constant
     240!---------
    241241real, parameter :: pvap2rho = 18.e-3/8.314
    242 !   Code
    243 !-------
     242! Code
     243!-----
    244244! 0. Compute constriction over the layer
    245245! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
    … …  
    364364END SUBROUTINE constriction
    365365
     366!-----------------------------------------------------------------------
     367SUBROUTINE compute_Tice_pem(nsoil, ptsoil, ptsurf, ice_depth, Tice)
     368!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     369!!!
     370!!! Purpose: Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells.
     371!!!
     372!!! Author: LL
     373!!!
     374!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     375
     376use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
     377
     378implicit none
     379
     380! Inputs
     381! ------
     382integer,                 intent(in) :: nsoil     ! Number of soil layers
     383real, dimension(nsoil),  intent(in) :: ptsoil    ! Soil temperature (K)
     384real,                    intent(in) :: ptsurf    ! Soil temperature (K)
     385real,                    intent(in) :: ice_depth ! Ice depth (m)
     386
     387! Outputs
     388! ------
     389real, intent(out) :: Tice ! Ice temperatures (K)
     390
     391! Locals
     392! ------
     393integer :: ik       ! Loop variables
     394integer :: indexice ! Index of the ice
     395
     396! Code
     397!-----
     398indexice = -1
     399if (ice_depth >= mlayer_PEM(nsoil - 1)) then
     400    Tice = ptsoil(nsoil)
     401else
     402    if(ice_depth < mlayer_PEM(0)) then
     403        indexice = 0.
     404    else     
     405        do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations
     406            if (mlayer_PEM(ik) <= ice_depth .and. mlayer_PEM(ik + 1) > ice_depth) then
     407                indexice = ik + 1
     408                exit
     409            endif
     410        enddo
     411    endif
     412    if (indexice < 0) then
     413        call abort_physic("compute_Tice_pem","subsurface ice is below the last soil layer",1)
     414    else
     415        if(indexice >= 1) then ! Linear inteprolation between soil temperature
     416            Tice = (ptsoil(indexice) - ptsoil(indexice + 1))/(mlayer_PEM(indexice - 1) - mlayer_PEM(indexice))*(ice_depth - mlayer_PEM(indexice)) + ptsoil(indexice + 1)
     417        else ! Linear inteprolation between the 1st soil temperature and the surface temperature
     418            Tice = (ptsoil(1) - ptsurf)/mlayer_PEM(0)*(ice_depth - mlayer_PEM(0)) + ptsoil(1)
     419        endif ! index ice >= 1
     420    endif !indexice < 0
     421endif ! icedepth > mlayer_PEM(nsoil - 1)
     422
     423END SUBROUTINE compute_Tice_pem
     424
    366425END MODULE ice_table_mod

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

    r3324 r3327  
    158158
    159159! Variables for h2o_ice evolution
    160 real, dimension(:,:),   allocatable :: h2o_ice              ! h2o ice in the PEM
    161 real, dimension(:,:,:), allocatable :: min_h2o_ice          ! Minima of h2o ice at each point for the PCM years [kg/m^2]
    162 real                                :: h2o_surf_ini         ! Initial surface of sublimating h2o ice
    163 real                                :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa]
    164 real                                :: global_avg_press_old ! constant: Global average pressure of initial/previous time step
    165 real                                :: global_avg_press_new ! constant: Global average pressure of current time step
    166 real, dimension(:,:),   allocatable :: zplev_new            ! Physical x Atmospheric field: mass of the atmospheric  layers in the pem at current time step [kg/m^2]
    167 real, dimension(:,:),   allocatable :: zplev_PCM            ! same but retrieved from the PCM [kg/m^2]
    168 real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
    169 real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
    170 integer                             :: stopPEM              ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu
    171 real, save                          :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
    172 real, dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
    173 integer                             :: timelen              ! # time samples
    174 real                                :: ave                  ! intermediate varibale to compute average
    175 real, dimension(:,:),   allocatable :: p                    ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
    176 real                                :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
     160real, dimension(:,:),    allocatable  :: h2o_ice              ! h2o ice in the PEM
     161real, dimension(:,:,:),  allocatable  :: min_h2o_ice          ! Minima of h2o ice at each point for the PCM years [kg/m^2]
     162real                                  :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
     163logical, dimension(:,:), allocatable  :: ini_h2oice_sublim    ! Logical array to know if h2o ice is sublimating
     164real                                  :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa]
     165real                                  :: global_avg_press_old ! constant: Global average pressure of initial/previous time step
     166real                                  :: global_avg_press_new ! constant: Global average pressure of current time step
     167real,   dimension(:,:),   allocatable :: zplev_new            ! Physical x Atmospheric field: mass of the atmospheric  layers in the pem at current time step [kg/m^2]
     168real,   dimension(:,:),   allocatable :: zplev_PCM            ! same but retrieved from the PCM [kg/m^2]
     169real,   dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
     170real,   dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
     171integer                               :: stopPEM              ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu
     172real, save                            :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
     173real,   dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
     174integer                               :: timelen              ! # time samples
     175real                                  :: ave                  ! intermediate varibale to compute average
     176real,   dimension(:,:),   allocatable :: p                    ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
     177real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
    177178
    178179! Variables for co2_ice evolution
    179 real, dimension(:,:),   allocatable :: co2_ice          ! co2 ice in the PEM
    180 real, dimension(:,:),   allocatable :: co2_ice_ini      ! Initial amount of co2 ice in the PEM
    181 real, dimension(:,:,:), allocatable :: min_co2_ice      ! Minimum of co2 ice at each point for the first year [kg/m^2]
    182 real                                :: co2_surf_ini     ! Initial surface of sublimating co2 ice
    183 real, dimension(:,:),   allocatable :: vmr_co2_PCM      ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
    184 real, dimension(:,:),   allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM
    185 real, dimension(:,:),   allocatable :: q_co2_PEM_phys   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
     180real,  dimension(:,:),   allocatable :: co2_ice           ! co2 ice in the PEM
     181real,  dimension(:,:),   allocatable :: co2_ice_ini       ! Initial amount of co2 ice in the PEM
     182real,  dimension(:,:,:), allocatable :: min_co2_ice       ! Minimum of co2 ice at each point for the first year [kg/m^2]
     183real                                 :: co2ice_ini_surf   ! Initial surface of sublimating co2 ice
     184logical, dimension(:,:), allocatable :: ini_co2ice_sublim ! Logical array to know if co2 ice is sublimating
     185real,   dimension(:,:),  allocatable :: vmr_co2_PCM       ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
     186real,   dimension(:,:),  allocatable :: vmr_co2_PEM_phys  ! Physics x Times co2 volume mixing ratio used in the PEM
     187real,   dimension(:,:),  allocatable :: q_co2_PEM_phys    ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
    186188
    187189! Variables for stratification (layering) evolution
    … …  
    556558! We save the places where h2o ice is sublimating
    557559! We compute the surface of h2o ice sublimating
    558 co2_surf_ini = 0.
    559 h2o_surf_ini = 0.
     560allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope))
     561co2ice_ini_surf = 0.
     562h2oice_ini_surf = 0.
    560563Total_surface = 0.
     564ini_co2ice_sublim = .false.
     565ini_h2oice_sublim = .false.
    561566do i = 1,ngrid
    562567    Total_surface = Total_surface + cell_area(i)
    563568    do islope = 1,nslope
    564         if (tend_co2_ice(i,islope) < 0.) co2_surf_ini = co2_surf_ini + cell_area(i)*subslope_dist(i,islope)
    565         if (tend_h2o_ice(i,islope) < 0.) h2o_surf_ini = h2o_surf_ini + cell_area(i)*subslope_dist(i,islope)
     569        if (tend_co2_ice(i,islope) < 0.) then
     570            ini_co2ice_sublim(i,islope) = .true.
     571            co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
     572        endif
     573        if (tend_h2o_ice(i,islope) < 0.) then
     574            ini_h2oice_sublim(i,islope) = .true.
     575            h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
     576        endif
    566577    enddo
    567578enddo
    568579
    569 write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2_surf_ini
    570 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2o_surf_ini
     580write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf
     581write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
    571582write(*,*) "Total surface of the planet =", Total_surface
    572583allocate(zplev_PCM(ngrid,nlayer + 1))
    … …  
    778789    do ig = 1,ngrid
    779790        do islope = 1,nslope
    780             if (co2_ice_ini(ig,islope) > 0.5 .and. co2_ice(ig,islope) < 1.e-10) then ! co2 ice disappeared, look for closest point without co2ice
     791            if (co2_ice_ini(ig,islope) > 0. .and. co2_ice(ig,islope) < 1.e-10) then ! co2 ice disappeared, look for closest point without co2ice
    781792                if (latitude_deg(ig) > 0) then
    782793                    do ig_loop = ig,ngrid
    783794                        do islope_loop = islope,iflat,-1
    784                             if (co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
     795                            if (co2_ice_ini(ig_loop,islope_loop) < 1.e-10 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
    785796                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
    786797                                bool_sublim = .true.
    … …  
    793804                    do ig_loop = ig,1,-1
    794805                        do islope_loop = islope,iflat
    795                             if(co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
     806                            if (co2_ice_ini(ig_loop,islope_loop) < 1.e-10 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
    796807                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
    797808                                bool_sublim = .true.
    … …  
    802813                    enddo
    803814                endif
    804                 co2_ice_ini(ig,islope) = 0
     815                co2_ice_ini(ig,islope) = 0.
    805816                if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
    806817                    albedo(ig,1,islope) = albedo_h2o_frost
    … …  
    954965!    II_g Checking the stopping criterion
    955966!------------------------
    956     call stopping_crit_h2o_ice(cell_area,h2o_surf_ini,h2o_ice,stopPEM,ngrid)
    957     call stopping_crit_co2(cell_area,co2_surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
     967    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
     968    call stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
    958969
    959970    year_iter = year_iter + dt_pem
    … …  
    11071118#ifndef CPP_STD
    11081119    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
    1109                   nlayer,nq,ptimestep,pday,0.,cell_area,albedodat,      &
     1120                  nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
    11101121                  inertiedat,def_slope,subslope_dist)
    1111     call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
     1122    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,      &
    11121123                  ptimestep,ztime_fin,tsurf,tsoil,inertiesoil,        &
    11131124                  albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
    … …  
    11151126#else
    11161127    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
    1117                   nlayer,nq,ptimestep,pday,time_phys,cell_area,         &
     1128                  nlayer,nq,ptimestep,pday,time_phys,cell_area,    &
    11181129                  albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
    1119     call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,  &
     1130    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,       &
    11201131                  ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
    11211132                  cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab,   &
    … …  
    11541165deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
    11551166deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
    1156 deallocate(co2_ice_ini,stratif)
     1167deallocate(co2_ice_ini,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
    11571168!----------------------------- END OUTPUT ------------------------------
    11581169

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

    r3206 r3327  
    55!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    66!!!
    7 !!! Purpose: Compute the soil thermal properties 
     7!!! Purpose: Compute the soil thermal properties
    88!!! Author:  LL, 04/2023
    99!!!
    … …  
    1919!   --------
    2020!
    21 !   author: LL, 04/2023 
     21!   author: LL, 04/2023
    2222!   -------
    23 !         
     23!
    2424!=======================================================================
    2525
    … …  
    3737logical, intent(in) :: ispureice           ! Boolean to check if ice is massive or just pore filling
    3838real,    intent(in) :: pore_filling        ! ice pore filling in each layer (1)
    39 real,    intent(in) :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2) 
     39real,    intent(in) :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2)
    4040real,    intent(out) :: ice_thermalinertia ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)
    4141
    4242!    Local Variables
    4343!    --------------
    44 REAL :: inertie_purewaterice = 2100 ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004) 
     44REAL :: inertie_purewaterice = 2100 ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004)
    4545!=======================================================================
    4646!   Beginning of the code
    … …  
    5353endif
    5454
    55 END SUBROUTINE 
     55END SUBROUTINE
    5656!=======================================================================
    5757
    … …  
    6464implicit none
    6565
    66 ! Input:
    67 integer, intent(in) ::  ngrid, nslope, nsoil_PEM      ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
    68 real, intent(in) :: p_avg_new                         ! Global average surface pressure [Pa]
    69 real, intent(in) :: tendencies_waterice(ngrid,nslope) ! Tendencies on the water ice [kg/m^2/year]
    70 real, intent(in) :: waterice(ngrid,nslope)            ! Surface Water ice [kg/m^2]
    71 real, intent(in) :: ice_depth(ngrid,nslope)           ! Ice table depth [m]
    72 real, intent(in) :: ice_thickness(ngrid,nslope)       ! Ice table thickness [m]
     66! Input:
     67integer,                       intent(in) :: ngrid, nslope, nsoil_PEM ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
     68real,                          intent(in) :: p_avg_new                ! Global average surface pressure [Pa]
     69real, dimension(ngrid,nslope), intent(in) :: tendencies_waterice      ! Tendencies on the water ice [kg/m^2/year]
     70real, dimension(ngrid,nslope), intent(in) :: waterice                 ! Surface Water ice [kg/m^2]
     71real, dimension(ngrid,nslope), intent(in) :: ice_depth                ! Ice table depth [m]
     72real, dimension(ngrid,nslope), intent(in) :: ice_thickness            ! Ice table thickness [m]
     73
    7374! Outputs:
    74 
    75 real, intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) ! Soil Thermal Inertia [J/m^2/K/s^1/2]
    76  
     75real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2]
     76
    7777! Constants:
    78 
    79  REAL ::  reg_inertie_thresold = 370.                      ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
    80  REAL ::  reg_inertie_minvalue = 50.                       ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
    81  REAL ::  reg_inertie_maxvalue = 370.                      ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
    82  REAL ::  ice_inertia                                      ! Inertia of water ice [SI]
    83  REAL ::  P610 = 610.0                                     ! current average pressure on Mars [Pa]
    84  REAL ::  C = 0.0015                                       ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [uniteless]
    85  REAL ::  K = 8.1*1e4                                      ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [tor, or 133.3Pa]
    86  REAL ::  Pa2Tor = 1./133.3                                ! Conversion from Pa to tor [Pa/tor]
    87 
     78real :: reg_inertie_thresold = 370. ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
     79real :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
     80real :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
     81real :: ice_inertia                 ! Inertia of water ice [SI]
     82real :: P610 = 610.0                ! current average pressure on Mars [Pa]
     83real :: C = 0.0015                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [uniteless]
     84real :: K = 8.1*1e4                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [tor, or 133.3Pa]
     85real :: Pa2Tor = 1./133.3           ! Conversion from Pa to tor [Pa/tor]
    8886
    8987! Local variables:
    90 
    91  INTEGER :: ig,islope,isoil,iref,iend                      ! Loop variables
    92  REAL :: regolith_inertia(ngrid,nslope)                    ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
    93  REAL :: delta                                             ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
    94  REAL :: ice_bottom_depth                                  ! Bottom depth of the subsurface ice [m]
    95  REAL :: d_part                                            ! Regolith particle size [micrometer]
    96  
    97 
    98 write(*,*) "Update soil propreties"
    99        
     88integer :: ig, islope, isoil, iref, iend          ! Loop variables
     89real, dimension(ngrid,nslope) :: regolith_inertia ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
     90real    :: delta                                  ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
     91real    :: ice_bottom_depth                       ! Bottom depth of the subsurface ice [m]
     92real    :: d_part                                 ! Regolith particle size [micrometer]
     93
     94write(*,*) "Update soil properties"
     95
    10096! 1. Modification of the regolith thermal inertia.
    101  do islope = 1,nslope
    102    regolith_inertia(:,islope) = inertiedat_PEM(:,1)
    103    do ig = 1,ngrid
    104       if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then
    105               regolith_inertia(ig,islope) = TI_regolith_avg
    106        endif
    107       if(reg_thprop_dependp) then
    108           if(TI_PEM(ig,1,islope).lt.reg_inertie_thresold) then
    109             d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Tor)**(0.6)))**(-1./(0.11*log10(P610*Pa2Tor/K))) ! compute particle size, in micrometer
    110             regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Tor)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Tor/K)))
    111             if(regolith_inertia(ig,islope).gt.reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
    112             if(regolith_inertia(ig,islope).lt.reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
    113           endif
    114       endif
    115    enddo
    116  enddo
    117 
    118 ! 2. Build new Thermal Inertia
    119 do  islope=1,nslope
    120    do ig = 1,ngrid
    121      do isoil = 1,index_breccia
    122          TI_PEM(ig,isoil,islope) = regolith_inertia(ig,islope)
    123      enddo
    124      if(regolith_inertia(ig,islope).lt.TI_breccia) then
     97do islope = 1,nslope
     98    regolith_inertia(:,islope) = inertiedat_PEM(:,1)
     99    do ig = 1,ngrid
     100        if (tendencies_waterice(ig,islope) < -1.e-5 .and. waterice(ig,islope) == 0) regolith_inertia(ig,islope) = TI_regolith_avg
     101        if (reg_thprop_dependp) then
     102            if (TI_PEM(ig,1,islope) < reg_inertie_thresold) then
     103                d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Tor)**(0.6)))**(-1./(0.11*log10(P610*Pa2Tor/K))) ! compute particle size, in micrometer
     104                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Tor)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Tor/K)))
     105                if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
     106                if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
     107            endif
     108        endif
     109    enddo
     110enddo
     111
     112! 2. Build new Thermal Inertia
     113do islope = 1,nslope
     114    do ig = 1,ngrid
     115        do isoil = 1,index_breccia
     116            TI_PEM(ig,isoil,islope) = regolith_inertia(ig,islope)
     117        enddo
     118        if (regolith_inertia(ig,islope) < TI_breccia) then
    125119!!! transition
    126              delta = depth_breccia
    127              TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
    128             (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
    129                       ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
    130             do isoil=index_breccia+2,index_bedrock
    131               TI_PEM(ig,isoil,islope) = TI_breccia
    132             enddo     
    133       else ! we keep the high ti values
    134             do isoil=index_breccia+1,index_bedrock
    135               TI_PEM(ig,isoil,islope) = TI_PEM(ig,index_breccia,islope)
    136             enddo 
    137        endif ! TI PEM and breccia comparison
    138 !! transition
    139        delta = depth_bedrock
    140        TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
    141             (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
    142             ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
    143        do isoil=index_bedrock+2,nsoil_PEM
     120            delta = depth_breccia
     121            TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/              &
     122                                                  (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
     123                                                  ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
     124            do isoil = index_breccia + 2,index_bedrock
     125                TI_PEM(ig,isoil,islope) = TI_breccia
     126            enddo
     127        else ! we keep the high ti values
     128            do isoil = index_breccia + 1,index_bedrock
     129                TI_PEM(ig,isoil,islope) = TI_PEM(ig,index_breccia,islope)
     130            enddo
     131        endif ! TI PEM and breccia comparison
     132!!! transition
     133        delta = depth_bedrock
     134        TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/              &
     135                                              (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
     136                                              ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
     137        do isoil = index_bedrock + 2,nsoil_PEM
    144138            TI_PEM(ig,isoil,islope) = TI_bedrock
    145        enddo   
    146       enddo ! ig
    147 ENDDO ! islope
     139        enddo
     140    enddo ! ig
     141enddo ! islope
    148142
    149143!  3. Build new TI in case of ice table
    150   do ig=1,ngrid
    151     do islope=1,nslope
    152      if (ice_depth(ig,islope).gt.-1e-6) then
    153 ! 3.0 Case where it is perennial ice
    154       if (ice_depth(ig,islope).lt. 1e-10) then
    155        call ice_thermal_properties(.true.,1.,0.,ice_inertia)
    156        do isoil = 1,nsoil_PEM
    157          TI_PEM(ig,isoil,islope)=ice_inertia
    158        enddo
    159       else
    160 
    161         call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
     144do ig = 1,ngrid
     145    do islope = 1,nslope
     146        if (ice_depth(ig,islope) > -1.e-6) then
     147        ! 3.0 Case where it is perennial ice
     148            if (ice_depth(ig,islope) < 1.e-10) then
     149                call ice_thermal_properties(.true.,1.,0.,ice_inertia)
     150                do isoil = 1,nsoil_PEM
     151                    TI_PEM(ig,isoil,islope) = ice_inertia
     152                enddo
     153            else
     154                call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
    162155        ! 3.1.1 find the index of the mixed layer
    163         iref=0 ! initialize iref
    164         do isoil=1,nsoil_PEM ! loop on layers to find the beginning of the ice table
    165           if (ice_depth(ig,islope).ge.layer_PEM(isoil)) then
    166             iref=isoil ! pure regolith layer up to here
    167           else
    168          ! correct iref was obtained in previous cycle
    169             exit
    170           endif
    171         enddo
     156                iref = 0 ! initialize iref
     157                do isoil = 1,nsoil_PEM ! loop on layers to find the beginning of the ice table
     158                    if (ice_depth(ig,islope) >= layer_PEM(isoil)) then
     159                        iref = isoil ! pure regolith layer up to here
     160                    else
     161                        exit ! correct iref was obtained in previous cycle
     162                    endif
     163                enddo
    172164        ! 3.1.2 find the index of the end of the ice table
    173         iend=0 ! initialize iend
    174         ice_bottom_depth = ice_depth(ig,islope)+ice_thickness(ig,islope)
    175         do isoil=1,nsoil_PEM ! loop on layers to find the end of the ice table
    176           if (ice_bottom_depth.ge.layer_PEM(isoil)) then
    177             iend=isoil ! pure regolith layer up to here
    178           else
    179          ! correct iref was obtained in previous cycle
    180             exit
    181           endif
    182         enddo
    183       ! 3.2 Build the new ti
    184         do isoil=1,iref
    185           TI_PEM(ig,isoil,islope) = TI_PEM(ig,1,islope)
    186         enddo
    187         if (iref.lt.nsoil_PEM) then
    188            if(iref.eq.iend) then
    189             ! Ice table begins and end in the same mixture Mixtures with three components
    190             if (iref.ne.0) then !
    191             ! mixed layer
    192               TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
    193                       (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
    194                       ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
    195                       ((layer_PEM(iref+1)-ice_bottom_depth)/(TI_PEM(ig,iref+1,islope)**2))))
    196             else ! first layer is already a mixed layer
    197               ! (ie: take layer(iref=0)=0)
    198               TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
    199                       (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
    200                       ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ &
    201                       ((layer_PEM(2)-ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
    202             endif ! of if (iref.ne.0)
    203            else
    204             if (iref.ne.0) then
    205             ! mixed layer
    206               TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
    207                       (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
    208                       ((layer_PEM(iref+1)-ice_depth(ig,islope))/(ice_inertia**2))))
    209             else ! first layer is already a mixed layer
    210               ! (ie: take layer(iref=0)=0)
    211               TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
    212                       (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
    213                       ((layer_PEM(1)-ice_depth(ig,islope))/(ice_inertia**2))))
    214             endif ! of if (iref.ne.0)       
    215            endif ! iref.eq.iend
    216       ! 3.3 Build the new ti
    217           do isoil=iref+2,iend
    218             TI_PEM(ig,isoil,islope)=ice_inertia
    219           enddo
    220 
    221 
    222           if(iend.lt.nsoil_PEM) then
    223              TI_PEM(ig,iend+1,islope)=sqrt((layer_PEM(iend+1)-layer_PEM(iend))/ &
    224                       (((ice_bottom_depth-layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
    225                       ((layer_PEM(iend+1)-ice_bottom_depth)/(ice_inertia**2))))
    226           endif
    227 
    228 
    229         endif ! of if (iref.lt.(nsoilmx))
    230       endif ! permanent glaciers
    231      endif !ice_depth(ig,islope) > 0.
    232 !     write(*,*) 'TI=', TI_PEM(ig,:,islope)
     165                iend = 0 ! initialize iend
     166                ice_bottom_depth = ice_depth(ig,islope) + ice_thickness(ig,islope)
     167                do isoil = 1,nsoil_PEM ! loop on layers to find the end of the ice table
     168                    if (ice_bottom_depth >= layer_PEM(isoil)) then
     169                        iend = isoil ! pure regolith layer up to here
     170                    else
     171                        exit ! correct iref was obtained in previous cycle
     172                    endif
     173                enddo
     174        ! 3.2 Build the new ti
     175                do isoil = 1,iref
     176                    TI_PEM(ig,isoil,islope) = TI_PEM(ig,1,islope)
     177                enddo
     178                if (iref < nsoil_PEM) then
     179                    if (iref == iend) then
     180                        ! Ice table begins and end in the same mixture Mixtures with three components
     181                        if (iref /= 0) then ! mixed layer
     182                            TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                             &
     183                                                         (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
     184                                                         ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) +            &
     185                                                         ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2))))
     186                        else ! first layer is already a mixed layer
     187                            ! (ie: take layer(iref=0)=0)
     188                            TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/ &
     189                                                  (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) +           &
     190                                                  ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) + &
     191                                                  ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
     192                        endif ! of if (iref /= 0)
     193                    else
     194                        if (iref /= 0) then ! mixed layer
     195                            TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                             &
     196                                                         (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
     197                                                         ((layer_PEM(iref + 1) - ice_depth(ig,islope))/(ice_inertia**2))))
     198                        else ! first layer is already a mixed layer
     199                            ! (ie: take layer(iref=0)=0)
     200                            TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                 &
     201                                                  (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + &
     202                                                  ((layer_PEM(1) - ice_depth(ig,islope))/(ice_inertia**2))))
     203                        endif ! of if (iref /= 0)
     204                    endif ! iref == iend
     205        ! 3.3 Build the new ti
     206                    do isoil = iref + 2,iend
     207                        TI_PEM(ig,isoil,islope) = ice_inertia
     208                    enddo
     209                    if (iend < nsoil_PEM) then
     210                        TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/                         &
     211                                                     (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
     212                                                     ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
     213                    endif
     214                endif ! of if (iref < nsoilmx)
     215            endif ! permanent glaciers
     216        endif ! ice_depth(ig,islope) > 0.
     217!        write(*,*) 'TI=', TI_PEM(ig,:,islope)
    233218    enddo !islope
    234    enddo !ig
     219enddo !ig
    235220
    236221END SUBROUTINE update_soil_thermalproperties

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

    r3206 r3327  
    1414!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1515
    16 SUBROUTINE stopping_crit_h2o_ice(cell_area,surf_ini,h2o_ice,stopPEM,ngrid)
     16SUBROUTINE stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
    1717
    1818use time_evol_mod, only: h2o_ice_crit
    … …  
    2626!
    2727!=======================================================================
    28 
    29 !   arguments:
    30 !   ----------
    31 !   INPUT
    32 integer,                       intent(in) :: ngrid     ! # of physical grid points
    33 real, dimension(ngrid),        intent(in) :: cell_area ! Area of the cells
    34 real, dimension(ngrid,nslope), intent(in) :: h2o_ice   ! Actual density of h2o ice
    35 real,                          intent(in) :: surf_ini  ! Initial surface of h2o ice that was sublimating
    36 !   OUTPUT
     28! Inputs
     29!-------
     30integer,                          intent(in) :: ngrid             ! # of physical grid points
     31real,    dimension(ngrid),        intent(in) :: cell_area         ! Area of the cells
     32real,    dimension(ngrid,nslope), intent(in) :: h2o_ice           ! Actual density of h2o ice
     33real,                             intent(in) :: h2oice_ini_surf   ! Initial surface of sublimating h2o ice
     34logical, dimension(ngrid,nslope), intent(in) :: ini_h2oice_sublim ! Grid points where h2o ice was initially sublimating
     35! Outputs
     36!--------
    3737integer, intent(inout) :: stopPEM ! Stopping criterion code
    38 !   local:
    39 !   ------
    40 integer :: i, islope ! Loop
    41 real    :: surf_now ! Current surface of h2o ice
     38! Locals
     39! ------
     40integer :: i, islope       ! Loop
     41real    :: h2oice_now_surf ! Current surface of h2o ice
    4242
    4343!=======================================================================
    44 ! Computation of the present surface of h2o ice
    45 surf_now = 0.
     44! Computation of the present surface of h2o ice still sublimating
     45h2oice_now_surf = 0.
    4646do i = 1,ngrid
    4747    do islope = 1,nslope
    48         if (h2o_ice(i,islope) > 0.) surf_now = surf_now + cell_area(i)*subslope_dist(i,islope)
     48        if (ini_h2oice_sublim(i,islope) .and. h2o_ice(i,islope) > 0.) h2oice_now_surf = h2oice_now_surf + cell_area(i)*subslope_dist(i,islope)
    4949    enddo
    5050enddo
    5151
    5252! Check of the criterion
    53 if (surf_now < surf_ini*(1. - h2o_ice_crit)) then
     53if (h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)) then
    5454    stopPEM = 1
    5555    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
     56    write(*,*) "h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)", h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)
     57    write(*,*) "Current surface of h2o ice sublimating =", h2oice_now_surf
     58    write(*,*) "Initial surface of h2o ice sublimating =", h2oice_ini_surf
    5959    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
    60 else if (surf_now > surf_ini*(1. + h2o_ice_crit)) then
     60else if (h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)) then
    6161    stopPEM = 1
    6262    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
     63    write(*,*) "h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)", h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)
     64    write(*,*) "Current surface of h2o ice sublimating =", h2oice_now_surf
     65    write(*,*) "Initial surface of h2o ice sublimating =", h2oice_ini_surf
    6666    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
    6767endif
    6868
    69 if (abs(surf_ini) < 1.e-5) stopPEM = 0
     69if (abs(h2oice_ini_surf) < 1.e-5) stopPEM = 0
    7070
    7171END SUBROUTINE stopping_crit_h2o_ice
    … …  
    7373!=======================================================================
    7474
    75 SUBROUTINE stopping_crit_co2(cell_area,surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
     75SUBROUTINE stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
    7676
    7777use time_evol_mod, only: co2_ice_crit, ps_criterion
    … …  
    8585!
    8686!=======================================================================
    87 
    88 !   arguments:
    89 !   ----------
    90 !   INPUT
    91 integer,                       intent(in) :: ngrid, nslope        ! # of grid physical grid points
    92 real, dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
    93 real, dimension(ngrid,nslope), intent(in) :: co2_ice              ! Actual density of h2o ice
    94 real,                          intent(in) :: surf_ini             ! Initial surface of co2 ice that was sublimating
    95 real,                          intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files
    96 real,                          intent(in) :: global_avg_press_new ! Planet average pressure from the PEM computations
    97 !   OUTPUT
     87! Inputs
     88!-------
     89integer,                          intent(in) :: ngrid, nslope        ! # of grid physical grid points
     90real,    dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
     91real,    dimension(ngrid,nslope), intent(in) :: co2_ice              ! Actual density of co2 ice
     92real,                             intent(in) :: co2ice_ini_surf      ! Initial surface of sublimatingco2 ice
     93logical, dimension(ngrid,nslope), intent(in) :: ini_co2ice_sublim    ! Grid points where co2 ice was initially sublimating
     94real,                             intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files
     95real,                             intent(in) :: global_avg_press_new ! Planet average pressure from the PEM computations
     96! Outputs
     97!--------
    9898integer, intent(inout) :: stopPEM ! Stopping criterion code
    9999
    100 !   local:
    101 !   ------
    102 integer :: i, islope ! Loop
    103 real    :: surf_now ! Current surface of co2 ice
     100! Locals
     101! ------
     102integer :: i, islope       ! Loop
     103real    :: co2ice_now_surf ! Current surface of co2 ice
    104104
    105105!=======================================================================
    106 ! Computation of the present surface of co2 ice
    107 surf_now = 0.
     106! Computation of the present surface of co2 ice still sublimating
     107co2ice_now_surf = 0.
    108108do i = 1,ngrid
    109109    do islope = 1,nslope
    110         if (co2_ice(i,islope) > 0.) surf_now = surf_now + cell_area(i)*subslope_dist(i,islope)
     110        if (ini_co2ice_sublim(i,islope) .and. co2_ice(i,islope) > 0.) co2ice_now_surf = co2ice_now_surf + cell_area(i)*subslope_dist(i,islope)
    111111    enddo
    112112enddo
    113113
    114114! Check of the criterion
    115 if (surf_now < surf_ini*(1. - co2_ice_crit)) then
     115if (co2ice_now_surf < co2ice_ini_surf*(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_crit)", surf_now < surf_ini*(1. - co2_ice_crit)
    119     write(*,*) "Current surface of co2 ice sublimating =", surf_now
    120     write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
     118    write(*,*) "co2ice_now_surf < co2ice_ini_surf*(1. - co2_ice_crit)", co2ice_now_surf < co2ice_ini_surf*(1. - co2_ice_crit)
     119    write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf
     120    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_ini_surf
    121121    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
    122 else if (surf_now > surf_ini*(1. + co2_ice_crit)) then
     122else if (co2ice_now_surf > co2ice_ini_surf*(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_crit)", surf_now > surf_ini*(1. + co2_ice_crit)
    126     write(*,*) "Current surface of co2 ice sublimating =", surf_now
    127     write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
     125    write(*,*) "co2ice_now_surf > co2ice_ini_surf*(1. + co2_ice_crit)", co2ice_now_surf > co2ice_ini_surf*(1. + co2_ice_crit)
     126    write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf
     127    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_ini_surf
    128128    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
    129129endif
    130130
    131 if (abs(surf_ini) < 1.e-5) stopPEM = 0
     131if (abs(co2ice_ini_surf) < 1.e-5) stopPEM = 0
    132132
    133133if (global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)) then
Note: See TracChangeset for help on using the changeset viewer.