Changeset 2863


Ignore:
Timestamp:
Jan 12, 2023, 12:14:38 AM (2 years ago)
Author:
llange
Message:

PEM
H2O adsorption is now computed. Total mass of H2o adsorbded is written in the restart_PEM.
+Minor fixes and edit
LL & RV

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

Legend:

Unmodified
Added
Removed

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

    r2849 r2863  
    33  contains
    44
    5   subroutine regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, theta_h2o_adsorbded,m_h2o_adsorbed)
    6 #ifndef CPP_STD
    7  use vertical_layers_mod, ONLY: ap,bp
    8  use comsoil_h_PEM, only: n_1km
     5!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     6!!!
     7!!! Purpose: Compute CO2 and H2O adsorption, following the methods from Zent & Quinn 1995, Jackosky et al., 1997
     8!!!
     9!!! Author: LL, 01/2023
     10!!!
     11!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     12
     13  subroutine regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
     14                                m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg)
     15#ifndef CPP_STD
     16! inputs
     17 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! size dimension: physics x subslope x soil x timeseries
     18 REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) !tendancies on the glaciers [1]
     19 REAL,INTENT(IN) :: waterice(ngrid,nslope)              ! water ice at the surface [kg/m^2]
     20 REAL,INTENT(IN) :: co2ice(ngrid,nslope)                ! co2 ice at the surface [kg/m^2]
     21 REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope)      ! Soil Thermal inertia (J/K/^2/s^1/2)
     22 REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)   ! Soil temperature (K)
     23 REAL,INTENT(IN) :: ps(ngrid,timelen)                   ! Average surface pressure [Pa]
     24 REAL,INTENT(IN) :: q_co2(ngrid,timelen)                ! Mass mixing ratio of co2 in the first layer (kg/kg)
     25 REAL,INTENT(IN) :: q_h2o(ngrid,timelen)                ! Mass mixing ratio of H2o in the first layer (kg/kg)
     26
     27! outputs
     28 REAL,INTENT(INOUT) :: m_h2o_completesoil(ngrid,nsoil_PEM,nslope) ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope)     
     29 REAL,INTENT(OUT) :: delta_mh2oreg(ngrid)                         ! Difference density of h2o adsorbed (kg/m^3)
     30
     31 REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)   ! Density of co2 adsorbed (kg/m^3)
     32 REAL,INTENT(OUT) :: delta_mco2reg(ngrid)                            ! Difference density of co2 adsorbed (kg/m^3)
     33 
     34! local variables
     35 REAL :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules
     36! -------------
     37
     38! Compute H2O adsorption, then CO2 adsorption
     39
     40  call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
     41                                   theta_h2o_adsorbded,m_h2o_completesoil,delta_mh2oreg)
     42
     43
     44  call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o, &
     45                                                          tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mco2reg)
     46
     47#endif
     48   RETURN
     49  end subroutine
     50
     51!------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     52
     53  subroutine regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
     54                                   theta_h2o_adsorbded,m_h2o_completesoil,delta_mreg)
     55
     56!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     57!!!
     58!!! Purpose: Compute H2O adsorption, following the methods from Jackosky et al., 1997
     59!!!
     60!!! Author: LL, 01/2023
     61!!!
     62!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     63
     64#ifndef CPP_STD
     65      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km
     66      USE comcstfi_h, only:  pi
     67      use comslope_mod, only : subslope_dist,def_slope_mean
     68      use vertical_layers_mod, ONLY: ap,bp
    969#endif
    1070
    1171 implicit none
    12 
     72! inputs
    1373 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! size dimension
     74 REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) !tendancies on the glaciers ()
     75 REAL,INTENT(IN) :: waterice(ngrid,nslope)              ! water ice at the surface [kg/m^2]
     76 REAL,INTENT(IN) :: co2ice(ngrid,nslope)                ! co2 ice at the surface [kg/m^2]
    1477 REAL,INTENT(IN) :: ps(ngrid,timelen)                   ! surface pressure (Pa)         
    1578 REAL,INTENT(IN) :: q_co2(ngrid,timelen)                ! Mass mixing ratio of co2 in the first layer (kg/kg)
    … …  
    1881 REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)   ! Soil temperature (K)
    1982
    20 ! output
    21  REAL,INTENT(OUT) ::  m_h2o_adsorbed(ngrid,nsoil_PEM,nslope)     ! Density of h2o adsorbed (kg/m^3)
     83! outputs
     84 REAL,INTENT(INOUT) :: m_h2o_completesoil(ngrid,nsoil_PEM,nslope) ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope)     
    2285 REAL,INTENT(OUT) :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules
    23 
    24 ! constant
     86 REAL,INTENT(OUT) :: delta_mreg(ngrid)                         ! Difference density of h2o adsorbed (kg/m^3)
     87
     88! constants
    2589 REAL :: Ko =  1.57e-8            ! Jackosky et al. 1997
    2690 REAL :: e = 2573.9               ! Jackosky et al. 1997
    2791 REAL :: mu = 0.48                ! Jackosky et al. 1997
    28  REAL ::  inertie_thresold = 800. ! TI > 800 means cementation
     92 real :: m_theta = 2.84e-7        ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
     93 real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
     94 real ::  inertie_thresold = 800. ! TI > 800 means cementation
    2995 real :: m_h2o = 18.01528E-3      ! Molecular weight of h2o (kg/mol)
    3096 real :: m_co2 = 44.01E-3         ! Molecular weight of co2 (kg/mol)
    3197 real :: m_noco2 = 33.37E-3       ! Molecular weight of non co2 (kg/mol)
    32  REAL ::  rho_regolith = 2000.    ! density of the reoglith, Buhler & Piqueux 2021
     98 real ::  rho_regolith = 1500.    ! density of the regolith (2000 for buhler & piqueux 2021)
    3399 real :: alpha_clapeyron = -6143.7! eq. (2) in Murphy & Koop 2005
    34100 real :: beta_clapeyron = 28.9074 ! eq. (2) in Murphy & Koop 2005
    35  real :: mi = 2.84e-7             ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
    36  real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
    37 
    38 ! local variable
     101! local variables
     102#ifndef CPP_STD
     103 REAL :: deltam_reg_complete(ngrid,n_1km,nslope)         ! Difference in the mass per slope and soil layer (kg/m^3)
     104#endif
     105 real :: K                        ! Used to compute theta
     106 integer ig,iloop, islope,isoil,it   ! for loops
     107 INTEGER :: ispermanent_co2glaciers(ngrid,nslope)        ! Check if the co2 glacier is permanent
     108 INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)        ! Check if the h2o glacier is permanent
     109 REAL :: deltam_reg_slope(ngrid,nslope)  ! Difference density of h2o adsorbed per slope (kg/m^3)
     110 REAL :: dm_h2o_regolith_slope(ngrid,nsoil_PEM,nslope)   ! elementary h2o mass adsorded per mesh per slope
    39111 real :: A,B                                                   ! Used to compute the mean mass above the surface
    40  real :: K                                                     ! Used to compute theta
    41112 real :: p_sat                                                 ! saturated vapor pressure of ice
    42  integer ig,iloop, islope,isoil,it                             ! for loops
    43113 real,allocatable :: mass_mean(:,:)                            ! mean mass above the surface
    44114 real,allocatable :: zplev_mean(:,:)                           ! pressure above the surface
    45115 real,allocatable :: pvapor(:,:)                               ! partial pressure above the surface
    46  real, allocatable :: pvapor_avg(:)                            ! yearly average vapor pressure above the surface
     116 real, allocatable :: pvapor_avg(:)                            ! yearly
    47117
    48118! 0. Some initializations
     119#ifndef CPP_STD
    49120
    50121     allocate(mass_mean(ngrid,timelen))
    … …  
    52123     allocate(pvapor(ngrid,timelen))
    53124     allocate(pvapor_avg(ngrid))
    54 
    55 
    56 
    57       m_h2o_adsorbed(:,:,:) = 0.
    58       theta_h2o_adsorbded(:,:,:) = 0.
    59       A =(1/m_co2 - 1/m_noco2)
    60       B=1/m_noco2
    61 #ifndef CPP_STD
    62 ! 1. Compute rho surface yearly averaged
    63 
    64 !   1.1 Compute the partial pressure of vapor
     125     A =(1/m_co2 - 1/m_noco2)
     126     B=1/m_noco2
     127     theta_h2o_adsorbded(:,:,:) = 0.
     128     dm_h2o_regolith_slope(:,:,:) = 0.
     129
     130!0.1 Look at perenial ice
     131  do ig = 1,ngrid
     132    do islope = 1,nslope
     133     if((abs(tend_h2oglaciers(ig,islope)).gt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then
     134        ispermanent_h2oglaciers(ig,islope) = 1
     135     else
     136        ispermanent_h2oglaciers(ig,islope) = 0
     137     endif
     138
     139     if((abs(tend_co2glaciers(ig,islope)).gt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
     140        ispermanent_co2glaciers(ig,islope) = 1
     141     else
     142        ispermanent_co2glaciers(ig,islope) = 0
     143     endif
     144    enddo
     145   enddo
     146
     147!   0.2 Compute the partial pressure of vapor
    65148!a. the molecular mass into the column
    66149     do ig = 1,ngrid
    67150       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
    68151     enddo
     152
    69153
    70154! b. pressure level
    … …  
    74158       enddo
    75159     enddo
    76 
    77160! c. Vapor pressure
    78161     pvapor(:,:) = mass_mean(:,:)/m_h2o*q_h2o(:,:)*zplev_mean(:,:)
    79162     pvapor_avg(:) = sum(pvapor(:,:),2)/timelen
    80 
    81163     deallocate(pvapor)
    82164     deallocate(zplev_mean)
    83165     deallocate(mass_mean)
    84166
    85 ! 2. we compute the mass of co2 adsorded in each layer of the meshes 
     167! 1. we compute the mass of H2O adsorded in each layer of the meshes 
    86168
    87169 do ig = 1,ngrid
    88170  do islope = 1,nslope
    89171    do iloop = 1,n_1km
    90        K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
    91      if(TI_PEM(ig,iloop,islope).lt.inertie_thresold)  then
    92 
    93        theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu
    94        m_h2o_adsorbed(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith
    95      else
     172      K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
     173      if(TI_PEM(ig,iloop,islope).lt.inertie_thresold)  then
     174        theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu
     175      else
    96176        p_sat =exp(alpha_clapeyron/tsoil_PEM(ig,iloop,islope) +beta_clapeyron) ! we assume fixed temperature in the ice ... not really:q good but ...
    97177        theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1+K*p_sat))**mu
    98         m_h2o_adsorbed(ig,iloop,islope) =as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith
    99      endif
     178      endif
     179      dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*m_theta*rho_regolith
    100180   enddo
    101181  enddo
    102182 enddo
    103183
     184! 2. Check the exchange between the atmosphere and the regolith
     185
     186  do ig = 1,ngrid
     187   delta_mreg(ig) = 0.
     188   do islope = 1,nslope
     189    deltam_reg_slope(ig,islope) = 0.
     190    do iloop = 1,n_1km
     191       if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
     192             deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope)) &
     193                                                   *(layer_PEM(iloop+1) - layer_PEM(iloop))
     194       else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
     195             deltam_reg_complete(ig,iloop,islope) = 0.
     196       endif
     197       deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
     198    enddo
     199   delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
     200   enddo
     201  enddo
     202   m_h2o_completesoil(:,:,:) = dm_h2o_regolith_slope(:,:,:)
     203
     204
    104205 RETURN
    105206#endif
    106207  end subroutine
    107208
    108   SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_completesoil,delta_mreg)
     209!------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     210!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     211!!!
     212!!! Purpose: Compute CO2  following the methods from Zent & Quinn 1995
     213!!!
     214!!! Author: LL, 01/2023
     215!!!
     216!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     217
     218  SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,&
     219                                   tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg)
    109220#ifndef CPP_STD
    110221      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km
    111       USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
     222      USE comcstfi_h, only: pi
    112223      use comslope_mod, only : subslope_dist,def_slope_mean
    113224      use vertical_layers_mod, ONLY: ap,bp
    … …  
    115226
    116227      IMPLICIT NONE
    117 ! Input: 
     228! Inputs: 
    118229 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen             ! size dimension
    119230 REAL,INTENT(IN) :: ps(ngrid,timelen)                               ! Average surface pressure [Pa]
    … …  
    127238! Outputs:
    128239 REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)   ! Density of co2 adsorbed (kg/m^3)
    129  REAL,INTENT(INOUT) :: delta_mreg(ngrid)                            ! Difference density of co2 adsorbed (kg/m^3)
     240 REAL,INTENT(OUT) :: delta_mreg(ngrid)                              ! Difference density of co2 adsorbed (kg/m^3)
    130241 
    131242! Constants:
    … …  
    134245 REAL :: beta =  -1541.5  ! Zent & Quinn 1995
    135246 REAL ::  inertie_thresold = 800. ! TI > 800 means cementation
    136  REAL ::  rho_regolith = 2000. ! density of the reoglith, buhler & piqueux 2021
     247 REAL ::  rho_regolith = 1500. ! density of the reoglith, buhler & piqueux 2021
    137248 real :: m_co2 = 44.01E-3      ! Molecular weight of co2 (kg/mol)
    138249 real :: m_noco2 = 33.37E-3    ! Molecular weight of h2o (kg/mol)
    … …  
    144255 INTEGER :: ig,islope,iloop,it                           ! for loops
    145256 REAL :: dm_co2_regolith_slope(ngrid,nsoil_PEM,nslope)   ! elementary mass adsorded per mesh per slope
    146  INTEGER :: ispermanent_co2glaciers(ngrid,nslope)        ! Check if the glacier is permanent
    147  INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)        ! Check if the glacier is permanent
     257 INTEGER :: ispermanent_co2glaciers(ngrid,nslope)        ! Check if the co2 glacier is permanent
     258 INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)        ! Check if the h2o glacier is permanent
    148259#ifndef CPP_STD
    149260 REAL :: deltam_reg_complete(ngrid,n_1km,nslope)         ! Difference in the mass per slope and soil layer (kg/m^3)
    … …  
    152263 REAL :: m_h2o_adsorbed(ngrid,nsoil_PEM,nslope)          ! Density of CO2 adsorbed (kg/m^3)
    153264 REAL :: theta_h2o_adsorbed(ngrid,nsoil_PEM,nslope)     ! Fraction of the pores occupied by H2O molecules
     265 REAL :: delta_mh2o(ngrid)                              ! Difference density of h2o adsorbed (kg/m^3)
    154266!timelen array are allocated because heavy ...
    155267 real,allocatable :: mass_mean(:,:)                            ! mean mass above the surface
    … …  
    194306!a. the molecular mass into the column
    195307     do ig = 1,ngrid
    196        mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
     308       mass_mean(ig,:) = 1./(A*q_co2(ig,:) +B)
    197309     enddo
    198310
    … …  
    214326
    215327! 1. Compute the fraction of the pores occupied by H2O
    216  call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,theta_h2o_adsorbed, m_h2o_adsorbed)
     328
     329 call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
     330                                   theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o)
     331
     332
    217333
    218334! 2.  we compute the mass of co2 adsorded in each layer of the meshes 
    … …  
    225341                                             (alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
    226342     else
    227         if(abs(m_co2_completesoil(ig,iloop,islope)).lt.1-10) then !!! we are at first call
     343        if(abs(m_co2_completesoil(ig,iloop,islope)).lt.(1e-10)) then !!! we are at first call
    228344          dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) &
    229345                                                 /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
    … …  
    236352enddo
    237353
    238 ! 2. Check the exchange between the atmosphere and the regolith
     354! 3. Check the exchange between the atmosphere and the regolith
    239355
    240356  do ig = 1,ngrid
    … …  
    245361       if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
    246362             deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) &
    247                                                    *(layer_PEM(iloop+1) - layer_PEM(iloop)) 
     363                                                   *(layer_PEM(iloop+1) - layer_PEM(iloop))
    248364       else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
    249365             deltam_reg_complete(ig,iloop,islope) = 0.
    … …  
    254370   enddo
    255371  enddo
    256    m_co2_completesoil(:,:,:) = dm_co2_regolith_slope(:,:,:)
     372  m_co2_completesoil(:,:,:) = dm_co2_regolith_slope(:,:,:)
    257373
    258374!=======================================================================

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

    r2855 r2863  
    2020  real,parameter :: fluxgeo = 30e-3             ! Geothermal flux [W/m^2]
    2121  real, save, allocatable :: co2_adsorbded_phys(:,:,:)  ! co2 that is in the regolith [kg/m^2]
     22  real, save, allocatable :: h2o_adsorbded_phys(:,:,:)  ! h2o that is in the regolith [kg/m^2]
    2223  LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
    2324
    … …  
    3031  integer,intent(in) :: nslope ! number of slope within a mesh
    3132
    32     allocate(layer_PEM(nsoilmx_PEM)) !soil layer depths
    33     allocate(mlayer_PEM(0:nsoilmx_PEM-1)) ! soil mid-layer depths
    34     allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope)) ! soil thermal inertia
    35     allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope)) ! soil temperatures
     33    allocate(layer_PEM(nsoilmx_PEM))
     34    allocate(mlayer_PEM(0:nsoilmx_PEM-1))
     35    allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
     36    allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope))
    3637    allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1))
    3738    allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1))
    … …  
    4041    allocate(alph_PEM(ngrid,nsoilmx_PEM-1,nslope))
    4142    allocate(beta_PEM(ngrid,nsoilmx_PEM-1,nslope))
    42     allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) ! soil thermal inertia
     43    allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
    4344    allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
     45    allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
    4446  end subroutine ini_comsoil_h_PEM
    4547
    … …  
    6163    if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
    6264    if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys)
     65    if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys)
    6366  end subroutine end_comsoil_h_PEM
    6467

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

    r2859 r2863  
    11MODULE conf_pem_mod
    22
    3 IMPLICIT NONE
    43
    5 CONTAINS
    64
    7   SUBROUTINE conf_pem
    85
    9 #ifdef CPP_IOIPSL
    10   use IOIPSL, only: getin
    11 #else
    12   ! if not using IOIPSL, we still need to use (a local version of) getin
    13   use ioipsl_getincom, only: getin
    14 #endif
    15  
    16   USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, alpha_criterion, &
    17                 Max_iter_pem, evol_orbit_pem
    18   USE comsoil_h_pem, only: soil_pem
    196
    20  
    21 !PEM parameter
    22   year_bp_ini=0.
    23   CALL getin('year_bp_ini', year_bp_ini)
    247
    25   dt_pem=1
    26   CALL getin('dt_pem', dt_pem)
    278
    28   alpha_criterion=0.2
    29   CALL getin('alpha_criterion', alpha_criterion)
    309
    31   evol_orbit_pem=.false.
    32   CALL getin('evol_orbit_pem', evol_orbit_pem)
    3310
    34   Max_iter_pem=99999999
    35   CALL getin('Max_iter_pem', Max_iter_pem)
    3611
    37   soil_pem=.true.
    38   CALL getin('soil_pem', soil_pem)
    3912
    40   END SUBROUTINE conf_pem
    4113
    4214END MODULE conf_pem_mod
     15

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

    r2842 r2863  
    1111!=======================================================================
    1212!
    13 !  Routine that compute the evolution of the water ice
     13!  Routine that compute the evolution of the CO2 ice
    1414!
    1515!=======================================================================
    … …  
    3535  STOPPING=.false.
    3636
    37 ! Evolution of the water ice for each physical point
     37! Evolution of the CO2 ice for each physical point
    3838  do i=1,ngrid
    3939    do islope=1,nslope

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

    r2857 r2863  
    102102  enddo
    103103 
    104 
    105   real_coefficient=negative_part/pos_tend
    106   if(pos_tend.eq.0) real_coefficient = 0.
    107 
     104 
     105  if(pos_tend.eq.0) then
     106   real_coefficient = 0.
     107  else
     108   real_coefficient = negative_part/pos_tend
     109  endif
    108110  do i=1,ngrid
    109111    do islope=1, nslope

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

    r2859 r2863  
    8383      USE infotrac
    8484      USE geometry_mod, only: latitude_deg
    85 
    8685      use conf_pem_mod, only: conf_pem
    8786      use pemredem, only:  pemdem1
    … …  
    9190                              TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM,         & ! soil thermal inertia         
    9291                              tsoil_PEM, mlayer_PEM,layer_PEM,                  & !number of subsurface layers, soil mid layer depths
    93                               fluxgeo, co2_adsorbded_phys                         ! geothermal flux, mass of co2 in the regolith
    94       use adsorption_mod, only : regolith_co2adsorption
     92                              fluxgeo, co2_adsorbded_phys, h2o_adsorbded_phys     ! geothermal flux, mass of co2 & h2o in the regolith
     93      use adsorption_mod, only : regolith_adsorption
    9594
    9695!!! For orbit parameters
    … …  
    264263     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)                     ! Surface temperature saved from previous time step [K]
    265264     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)                         ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    266      REAL :: totmass_adsorbded                                           ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
     265     REAL, ALLOCATABLE :: delta_h2o_adsorbded(:)                         ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
     266     REAL :: totmassco2_adsorbded                                           ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
     267     REAL :: totmassh2o_adsorbded                                           ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
    267268     REAL :: alpha_clap_h2o = -6143.7                                    ! coeffcient to compute psat, from Murphie et Kood 2005 [K]
    268269     REAL :: beta_clap_h2o = 28.9074                                     ! coefficient to compute psat, from Murphie et Kood 2005 [1]
    … …  
    322323!----------------------------READ run.def ---------------------
    323324      CALL conf_gcm( 99, .TRUE. )
    324       CALL conf_pem
    325 
     325      CALL conf_pem     
     326 
    326327      call infotrac_init
    327328      nq=nqtot
    … …  
    529530     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    530531     allocate(delta_co2_adsorbded(ngrid))
     532     allocate(delta_h2o_adsorbded(ngrid))
    531533     allocate(watersurf_density_timeseries(iim+1,jjm+1,nslope,timelen))
    532534     allocate(watersoil_density_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
    … …  
    786788!---------------------------- Read the PEMstart ---------------------
    787789
    788       call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
     790      call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
    789791      tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsoil_phys_PEM_timeseries,&
    790       tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM, co2_adsorbded_phys,delta_co2_adsorbded,&
    791       watersurf_density_phys_ave,watersoil_density_phys_PEM_ave)
     792      tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,&
     793      watersurf_density_phys_ave,watersoil_density_phys_PEM_ave, &
     794      co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded)
    792795
    793796    if(soil_pem) then
    794      totmass_adsorbded = 0.
    795 
     797     totmassco2_adsorbded = 0.
     798     totmassh2o_adsorbded = 0.
    796799     do ig = 1,ngrid
    797800      do islope =1, nslope
    798801       do l = 1,nsoilmx_PEM - 1
    799           totmass_adsorbded = totmass_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     802          totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     803          subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
     804          cell_area(ig)
     805          totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
    800806          subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
    801807          cell_area(ig)
    … …  
    804810     enddo
    805811
    806      write(*,*) "Tot mass in the regolith=", totmass_adsorbded
     812     write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
     813     write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
    807814     deallocate(tsoil_ave_phys_yr1) 
    808815    endif !soil_pem
    … …  
    856863       print *, 'Global average pressure old time step',global_ave_press_old
    857864       print *, 'Global average pressure new time step',global_ave_press_new
    858 
    859      do i=1,ngrid
    860        if(soil_pem) then
    861          global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
    862        endif
    863      enddo
     865     
     866     if(soil_pem) then
     867      do i=1,ngrid
     868        global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
     869       enddo
     870     endif
    864871       print *, 'Global average pressure old time step',global_ave_press_old
    865872       print *, 'Global average pressure new time step',global_ave_press_new
    … …  
    10971104
    10981105! II_d.5 Update the mass of the regolith adsorbded
    1099    call regolith_co2adsorption(ngrid,nslope,nsoilmx_PEM,timelen,ps_phys_timeseries,tsoil_PEM,TI_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope, &
    1100                                co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), q_co2_PEM_phys,q_h2o_PEM_phys,co2_adsorbded_phys,delta_co2_adsorbded)
     1106
     1107    call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),co2ice_slope, &
     1108                             tsoil_PEM,TI_PEM,ps_phys_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
     1109                             h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
    11011110
    11021111  endif !soil_pem
    … …  
    12641273     enddo
    12651274
    1266      write(*,*) 'rapport ps',global_ave_press_new/global_ave_press_GCM
    1267 
    12681275! III_a.5 Tracer (for start)
    12691276     allocate(zplev_new(ngrid,nlayer+1))
    … …  
    13791386
    13801387         call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
    1381                       tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys)
     1388                      tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys)
    13821389
    13831390      print *, "restartfi_PEM.nc has been written"

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

    r2855 r2863  
    1    SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
     1   SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
    22                      tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
    3                       global_ave_pressure, m_co2_regolith_phys,deltam_co2_regolith_phys, watersurf_ave,watersoil_ave)
     3                      global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys, m_h2o_regolith_phys,deltam_h2o_regolith_phys)
    44   
    55   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
    66   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
    77   use comsoil_h, only:  volcapa,inertiedat
    8    use adsorption_mod, only : regolith_co2adsorption
     8   use adsorption_mod, only : regolith_adsorption
    99   USE temps_mod_evol, ONLY: year_PEM
     10
    1011
    1112  implicit none
    1213 
    1314  character(len=*), intent(in) :: filename                    ! name of the startfi_PEM.nc
    14   LOGICAL,intent(in) :: startpem_file                         ! boolean to check if we read the startfile or not
    1515  integer,intent(in) :: ngrid                                 ! # of physical grid points
    1616  integer,intent(in) :: nsoil_GCM                             ! # of vertical grid points in the GCM
    … …  
    4040  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed [kg/m^2]
    4141  real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
     42  real,intent(inout) :: m_h2o_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of h2o adsorbed [kg/m^2]
     43  real,intent(out) :: deltam_h2o_regolith_phys(ngrid)                 ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
    4244! local
    4345   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)                     ! soil temperature saved in the start [K]
    4446   real :: TI_startPEM(ngrid,nsoil_PEM,nslope)                        ! soil thermal inertia saved in the start [SI]
    4547   LOGICAL :: found                                                   ! check if variables are found in the start
     48   LOGICAL :: found2                                                  ! check if variables are found in the start
    4649   integer :: iloop,ig,islope,it,isoil                                ! index for loops
    4750   REAL :: TI_breccia = 750.                                          ! Thermal inertia of Breccia following Wood 2009 [SI]
    … …  
    5861   real :: alpha_clap_h2o = -6143.7                                   ! Intermediate coefficient to compute psat using clapeyron law, Murphie et al. 2005 [K^-1]
    5962   real :: beta_clap_h2o = 28.9074                                    ! Intermediate coefficient to compute psat using clapeyron law, Murphie et al. 2005 [1]
    60 
     63   LOGICAL :: startpem_file                                           ! boolean to check if we read the startfile or not
    6164
    6265!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    7477!!!
    7578!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     79
     80
     81!0. Check if the start_PEM exist.
     82
     83inquire(file=filename,exist =  startpem_file)
     84
     85write(*,*)'Start PEM=',startpem_file
     86
    7687
    7788!1. Run
    … …  
    247258!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    248259
    249 !4. CO2 Adsorption
    250 DO islope=1,nslope
    251   write(num,fmt='(i2.2)') islope
     260!4. CO2 & H2O Adsorption
     261  DO islope=1,nslope
     262   write(num,fmt='(i2.2)') islope
    252263   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
    253    if(.not.found) then
    254       write(*,*)'PEM settings: failed loading <m_co2_regolith_phys>'
    255       write(*,*)'will reconstruct the values of co2 adsorbded'
    256    m_co2_regolith_phys(:,:,:) = 0.
    257    call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
    258 
    259    deltam_co2_regolith_phys(:) = 0.
    260    exit
    261    else
    262    call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
    263 
    264    endif
    265 ENDDO
    266 
    267 
    268 print *,'PEMETAT0: CO2 adsorption done  '
     264    if((.not.found)) then
     265       m_co2_regolith_phys(:,:,:) = 0.
     266    endif
     267    exit
     268  ENDDO
     269
     270  DO islope=1,nslope
     271   write(num,fmt='(i2.2)') islope
     272   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
     273    if((.not.found2)) then
     274       m_h2o_regolith_phys(:,:,:) = 0.
     275    endif
     276    exit
     277  ENDDO
     278
     279    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
     280                                m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
     281
     282    if((.not.found)) then
     283      deltam_co2_regolith_phys(:) = 0.
     284    endif
     285    if((.not.found2)) then
     286       deltam_h2o_regolith_phys(:) = 0. 
     287    endif
     288print *,'PEMETAT0: CO2 & H2O adsorption done  '
    269289
    270290endif ! soil_pem
    … …  
    402422!d) Regolith adsorbed
    403423
    404    m_co2_regolith_phys(:,:,:) = 0.
    405    call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
    406    deltam_co2_regolith_phys(:) = 0.
     424
     425    m_co2_regolith_phys(:,:,:) = 0.
     426    m_h2o_regolith_phys(:,:,:) = 0.
     427
     428    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
     429                                m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
     430   
     431    deltam_co2_regolith_phys(:) = 0.
     432    deltam_h2o_regolith_phys(:) = 0. 
     433   
    407434
    408435print *,'PEMETAT0: CO2 adsorption done  '

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

    r2855 r2863  
    7272
    7373subroutine pemdem1(filename,year_iter,nsoil_PEM,ngrid,nslope, &
    74                     tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table,m_co2_regolith)
     74                    tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table, &
     75                    m_co2_regolith,m_h2o_regolith)
    7576  ! write time-dependent variable to restart file
    7677  use iostart_PEM, only : open_restartphy, close_restartphy, &
    … …  
    9495  integer,intent(IN) :: nslope
    9596  real,intent(in) :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
     97  real,intent(in) :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
    9698  integer :: islope
    9799  CHARACTER*2 :: num 
    … …  
    124126                 inertiesoil_slope_PEM(:,:,islope),year_tot)
    125127
    126   call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
     128 call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
    127129                    m_co2_regolith(:,:,islope), year_tot)
     130  call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbded in the regolith", &
     131                    m_h2o_regolith(:,:,islope), year_tot)
    128132
    129133ENDDO

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

    r2835 r2863  
    4545    do islope=1,nslope
    4646      ave=0.
    47       do t=1,timelen
    48         ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100)))**4  &
    49               -(beta/(alpha-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*global_ave_press_GCM/global_ave_press_new/100)))**4
    50       enddo
     47      if(abs(tendencies_co2_ice_phys(i,islope)).gt.1e-4) then
     48        do t=1,timelen
     49           ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4  &
     50              -(beta/(alpha-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*global_ave_press_GCM/global_ave_press_new/100.)))**4
     51         enddo
     52      endif
     53      if(ave.lt.1e-4) ave = 0.
    5154      tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)-coef*ave/timelen
    5255    enddo
    5356  enddo
    5457
    55 
    5658END SUBROUTINE recomp_tend_co2_slope

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

    r2849 r2863  
    55 USE comsoil_h_PEM, only: layer_PEM,n_1km,inertiedat_PEM
    66 USE vertical_layers_mod, ONLY: ap,bp
     7 USE comsoil_h_PEM, only: n_1km
    78 implicit none
    89! Input:
    … …  
    7172! 2. Modification of the regolith thermal inertia.
    7273
    73   do ig=1,ngrid
    74     do islope=1,nslope
    75      do iloop =1,n_1km
    76        d(ig,iloop,islope) = ((inertiedat_PEM(ig,iloop)*inertiedat_PEM(ig,iloop))/(volcapa*alpha*P610**0.6))**(-1/(0.11*log10(P610/beta)))
    77    if(TI_PEM(ig,iloop,islope).lt.inertie_thresold) then !           we are modifying the regolith properties, not ice
    78           TI_PEM(ig,iloop,islope) = sqrt(volcapa*alpha*(p_avg_new**0.6)* d(ig,iloop,islope)**(-0.11*log10(p_avg_new/beta)))
    79 
    80    endif
    81      enddo
    82    enddo
    83 enddo
     74!  do ig=1,ngrid
     75!    do islope=1,nslope
     76!     do iloop =1,n_1km
     77!       d(ig,iloop,islope) = ((inertiedat_PEM(ig,iloop)*inertiedat_PEM(ig,iloop))/(volcapa*alpha*P610**0.6))**(-1/(0.11*log10(P610/beta)))
     78!   if(TI_PEM(ig,iloop,islope).lt.inertie_thresold) then !           we are modifying the regolith properties, not ice
     79!          TI_PEM(ig,iloop,islope) = sqrt(volcapa*alpha*(p_avg_new**0.6)* d(ig,iloop,islope)**(-0.11*log10(p_avg_new/beta)))
     80!
     81!   endif
     82!     enddo
     83!   enddo
     84!enddo
    8485
    8586!  3. Build new TI for the PEM
    … …  
    8788  do ig=1,ngrid
    8889    do islope=1,nslope
    89 
     90      do iloop = 1,n_1km
     91       TI_PEM(ig,iloop,islope) = TI_PEM(ig,1,islope)
     92      enddo
    9093  if (ice_depth(ig,islope).gt. -1.e-10) then
    9194! 3.0 FIrst if permanent ice
    … …  
    108111      ! 4.2 Build the new ti
    109112      do iloop=1,iref
    110          TI_PEM(ig,iloop,islope) =TI_PEM(ig,iloop,islope)
     113         TI_PEM(ig,iloop,islope) =TI_PEM(ig,1,islope)
    111114      enddo
    112115      if (iref.lt.nsoil_PEM) then
Note: See TracChangeset for help on using the changeset viewer.