Changeset 2895

Timestamp:

Feb 10, 2023, 8:35:22 PM (22 months ago)

Author:

llange

Message:

PEM
Soil temperature initialisation has been updated
Conf_PEM improved by adding some options to the users (thermal regolith depend on the pressure, depth of the subsurface layers, etc.)
Minor edits then (+ svn update with RV had some issues, so there are some "artefact changes" ...)
LL

Location:

trunk

Files:

• LMDZ.COMMON/libf/evolution/adsorption_mod.F90 (modified) (11 diffs)
• LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90 (modified) (2 diffs)
• LMDZ.COMMON/libf/evolution/conf_pem.F90 (modified) (3 diffs)
• LMDZ.COMMON/libf/evolution/pem.F90 (modified) (22 diffs)
• LMDZ.COMMON/libf/evolution/pemetat0.F90 (modified) (18 diffs)
• LMDZ.COMMON/libf/evolution/pemredem.F90 (modified) (1 diff)
• LMDZ.COMMON/libf/evolution/read_data_GCM.F90 (modified) (4 diffs)
• LMDZ.COMMON/libf/evolution/soil_pem.F90 (modified) (1 diff)
• LMDZ.COMMON/libf/evolution/soil_pem_ini.F90 (added)
• LMDZ.COMMON/libf/evolution/soil_settings_PEM.F (modified) (2 diffs)
• LMDZ.COMMON/libf/evolution/update_soil.F90 (modified) (4 diffs)
• LMDZ.MARS/README (modified) (1 diff)

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

@@ -89 +89 @@
 
 #ifndef CPP_STD
-      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km
+      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,index_breccia,index_bedrock
       USE comcstfi_h, only:  pi 
       use comslope_mod, only : subslope_dist,def_slope_mean
@@ -117 +117 @@
  REAL :: mu = 0.48                ! Jackosky et al. 1997
  real :: m_theta = 2.84e-7        ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
- real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
+! real :: as = 18.9e3             ! Specific area, Buhler & Piqueux 2021
+ real :: as = 9.48e4              ! Specific area, Zent 
+ real :: as_breccia = 1.7e4       ! Specific area of basalt, Zent 
  real ::  inertie_thresold = 800. ! TI > 800 means cementation
  real :: m_h2o = 18.01528E-3      ! Molecular weight of h2o (kg/mol)
  real :: m_co2 = 44.01E-3         ! Molecular weight of co2 (kg/mol)
  real :: m_noco2 = 33.37E-3       ! Molecular weight of non co2 (kg/mol)
- real ::  rho_regolith = 1500.    ! density of the regolith (2000 for buhler & piqueux 2021)
+ real ::  rho_regolith = 2000.    ! density of the regolith (Zent 1995, Buhler & piqueux 2021)
  real :: alpha_clapeyron = -6143.7! eq. (2) in Murphy & Koop 2005
  real :: beta_clapeyron = 28.9074 ! eq. (2) in Murphy & Koop 2005
 ! local variables
 #ifndef CPP_STD
- REAL :: deltam_reg_complete(ngrid,n_1km,nslope)         ! Difference in the mass per slope and soil layer (kg/m^3)
+ REAL :: deltam_reg_complete(ngrid,index_breccia,nslope)         ! Difference in the mass per slope and soil layer (kg/m^3)
 #endif
  real :: K                        ! Used to compute theta
@@ -195 +197 @@
  do ig = 1,ngrid
   do islope = 1,nslope
-    do iloop = 1,n_1km
+    do iloop = 1,index_breccia
       K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
       if(TI_PEM(ig,iloop,islope).lt.inertie_thresold)  then
@@ -203 +205 @@
         theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1+K*p_sat))**mu
       endif
-      dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*m_theta*rho_regolith
+        dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*m_theta*rho_regolith
    enddo
   enddo
@@ -214 +216 @@
    do islope = 1,nslope
     deltam_reg_slope(ig,islope) = 0.
-    do iloop = 1,n_1km
+    do iloop = 1,index_breccia
        if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
              deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope)) &
@@ -245 +247 @@
                                    tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg)
 #ifndef CPP_STD
-      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km
+      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,index_breccia,index_bedrock,index_breccia
       USE comcstfi_h, only: pi
       use comslope_mod, only : subslope_dist,def_slope_mean
@@ -271 +273 @@
  REAL :: beta =  -1541.5  ! Zent & Quinn 1995
  REAL ::  inertie_thresold = 800. ! TI > 800 means cementation
- REAL ::  rho_regolith = 1500. ! density of the reoglith, buhler & piqueux 2021
+ REAL ::  rho_regolith = 2000. ! density of the regolith, buhler & piqueux 2021
  real :: m_co2 = 44.01E-3      ! Molecular weight of co2 (kg/mol)
  real :: m_noco2 = 33.37E-3    ! Molecular weight of h2o (kg/mol)
  real :: m_theta = 4.27e-7     ! Mass of co2 per m^2 absorbed 
- real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
-
+! real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
+ real :: as = 9.48e4           ! same as previous but from zent 
+ real :: as_breccia = 1.7e4       ! Specific area of basalt, Zent 1997
 ! Local          
  real :: A,B                                             ! Used to compute the mean mass above the surface
@@ -284 +287 @@
  INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)        ! Check if the h2o glacier is permanent
 #ifndef CPP_STD
- REAL :: deltam_reg_complete(ngrid,n_1km,nslope)         ! Difference in the mass per slope and soil layer (kg/m^3)
+ REAL :: deltam_reg_complete(ngrid,index_breccia,nslope)         ! Difference in the mass per slope and soil layer (kg/m^3)
 #endif
  REAL :: deltam_reg_slope(ngrid,nslope)                  ! Difference in the mass per slope  (kg/m^3)
@@ -362 +365 @@
  do ig = 1,ngrid
   do islope = 1,nslope
-    do iloop = 1,n_1km
+    do iloop = 1,index_breccia
      if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
      dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
@@ -374 +377 @@
         endif
      endif
-  enddo
+    enddo
  enddo
 enddo
@@ -384 +387 @@
    do islope = 1,nslope
     deltam_reg_slope(ig,islope) = 0.
-    do iloop = 1,n_1km
+    do iloop = 1,index_breccia
        if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
              deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) &

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

@@ -19 +19 @@
   real,save :: mu_PEM                           ! mu coefficient [SI]
   real,save :: fluxgeo                               ! Geothermal flux [W/m^2]
+  real,save :: depth_breccia                         ! Depth at which we have breccia [m]
+  real,save :: depth_bedrock                         ! Depth at which we have bedrock [m]
+  integer,save :: index_breccia                         ! last index of the depth grid before having breccia
+  integer,save :: index_bedrock                         ! last index of the depth grid before having bedrock
   LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
   real,save,allocatable,dimension(:) :: water_reservoir      ! Large reserve of water   [kg/m^2]
   real,save :: water_reservoir_nom
+  logical, save :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
+
 contains
 
@@ -61 +67 @@
     if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
     if (allocated(water_reservoir)) deallocate(water_reservoir)
-    if (allocated(water_reservoir)) deallocate(water_reservoir)
   end subroutine end_comsoil_h_PEM
 

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

@@ -15 +15 @@
   
   USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, &
-                Max_iter_pem, evol_orbit_pem,var_obl, var_ex, var_lsp
-  USE comsoil_h_pem, only: soil_pem,fluxgeo,water_reservoir_nom
+                Max_iter_pem, evol_orbit_pem, var_obl, var_ex, var_lsp
+  USE comsoil_h_pem, only: soil_pem,fluxgeo,water_reservoir_nom,depth_breccia,depth_bedrock,reg_thprop_dependp 
   USE adsorption_mod,only: adsorption_pem
   CHARACTER(len=20),parameter :: modname ='conf_pem' 
@@ -55 +55 @@
   var_obl = .true.
   CALL getin('var_obl',var_obl)
+  print*,'Does obliquity vary ?',var_obl
 
   var_ex = .true.
   CALL getin('var_ex',var_ex)
+  print*,'Does excentricity vary ?',var_ex
 
   var_lsp = .true.
   CALL getin('var_lsp',var_lsp)
+  print*,'Does Ls peri vary ?',var_lsp
+   
+  depth_breccia   = 10.
+  CALL getin('depth_breccia',depth_breccia)
+  print*,'Depth of breccia is set to',depth_breccia
+
+  depth_bedrock   = 1000.
+  CALL getin('depth_bedrock',depth_bedrock)
+  print*,'Depth of bedrock is set to',depth_bedrock
+
+  reg_thprop_dependp = .false.
+  CALL getin('reg_thprop_dependp',reg_thprop_dependp)
+  print*, 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
+
 
   if ((not(soil_pem)).and.adsorption_pem) then
@@ -69 +85 @@
   if ((not(soil_pem)).and.fluxgeo.gt.0.) then
        print*,'Soil is not activated but Flux Geo > 0.'
-       call abort_physic(modname,"Soil is not activated but Flux Geo > 0.'",1)
+       call abort_physic(modname,"Soil is not activated but Flux Geo > 0.",1)
+  endif
+ 
+  if ((not(soil_pem)).and.reg_thprop_dependp) then
+     print*,'Regolith properties vary with Ps only when soil is set to true'
+     call abort_physic(modname,'Regolith properties vary with Ps only when soil is set to true',1)
   endif
 
   if (evol_orbit_pem.and.year_bp_ini.eq.0.) then
-       print*,'You want to follow the file ob_ex_lsp.asc for changing orb parameters,'
-       print*,'but you did not specify from which year to start.'
-       call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini=0",1)
-  endif
-  
-   water_reservoir_nom = 1e4
+     print*,'You want to follow the file ob_ex_lsp.asc for changing orb parameters,'
+     print*,'but you did not specify from which year to start.'
+     call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini=0",1)
+  endif 
+
+  water_reservoir_nom = 1e4
   CALL getin('water_reservoir_nom',water_reservoir_nom)
+
   END SUBROUTINE conf_pem
 

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

@@ -78 +78 @@
                            co2iceflow, beta_slope, capcal_slope,&
                            albedo_slope,emiss_slope,qsurf_slope,&
-                           iflat,major_slope,ini_comslope_h
+                           iflat,major_slope,ini_comslope_h,fluxgeo_slope
       use time_phylmdz_mod, only: daysec,dtphys
       USE comconst_mod, ONLY: rad,g,r,cpp,pi
@@ -89 +89 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
       use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
-                              TI_PEM,inertiedat_PEM,         & ! soil thermal inertia          
-                              tsoil_PEM, mlayer_PEM,layer_PEM,                  & !number of subsurface layers, soil mid layer depths
-                              fluxgeo,water_reservoir                             ! geothermal flux
+                              TI_PEM,inertiedat_PEM,         &                        ! soil thermal inertia          
+                              tsoil_PEM, mlayer_PEM,layer_PEM,                  &     ! Soil temp, number of subsurface layers, soil mid layer depths
+                              fluxgeo, &                             ! geothermal flux for the PEM and GCM
+                              water_reservoir                                         ! Water ressources
       use adsorption_mod, only : regolith_adsorption,adsorption_pem, &   ! bool to check if adsorption, main subroutine
                                  ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! allocate arrays                                
@@ -159 +160 @@
 
 ! Variable for h2o_ice evolution
-      REAL :: ini_surf_h2o                                         ! Initial surface of sublimating water ice
+      REAL :: ini_surf_h2o
       REAL :: ini_surf_co2                                         ! Initial surface of sublimating co2 ice
 
@@ -258 +259 @@
      REAL :: beta_clap_h2o = 28.9074                                     ! coefficient to compute psat, from Murphie et Kood 2005 [1]
      LOGICAL :: bool_sublim                                              ! logical to check if there is sublimation or not
+      
 
 !! Some parameters for the PEM run
@@ -265 +267 @@
     REAL :: timestep                     ! timestep [s]
     REAL :: watercap_sum                 ! total mass of water cap [kg/m^2] 
-
-    REAL WC_sum
 
 #ifdef CPP_STD
@@ -482 +482 @@
      allocate(flag_co2flow_mesh(ngrid))
 
-     flag_co2flow(:,:) = 0.     
-     flag_co2flow_mesh(:) = 0.
+       flag_co2flow(:,:) = 0.     
+       flag_co2flow_mesh(:) = 0.
 
 
@@ -499 +499 @@
 
      call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
+
 
      allocate(vmr_co2_gcm(iim+1,jjm+1,timelen))
@@ -530 +531 @@
      allocate(watersoil_density_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
      allocate(watersoil_density_phys_PEM_ave(ngrid,nsoilmx_PEM,nslope))
+
      print *, "Downloading data Y1..."
 
-     call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm,vmr_co2_gcm,ps_GCM_yr1,min_co2_ice_slope_1,min_h2o_ice_slope_1,&   
+     call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm,vmr_co2_gcm,ps_GCM_yr1,min_co2_ice_slope_1,min_h2o_ice_slope_1,&      
                        nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &      
                        watersurf_density_timeseries,watersoil_density_timeseries)
@@ -661 +663 @@
 
 !----- Compute tendencies from the PCM run
+
      allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
      allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
@@ -765 +768 @@
     do ig = 1,ngrid
       do islope = 1,nslope
-          qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
+!          qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
       enddo
     enddo
@@ -787 +790 @@
      write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
      write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
+     endif ! adsorption
      deallocate(tsoil_ave_phys_yr1)  
-    endif !soil_pem
      deallocate(tsurf_ave_phys_yr1)
      deallocate(ps_phys_timeseries_yr1)  
@@ -837 +840 @@
      enddo
        print *, 'Global average pressure old time step',global_ave_press_old
-       print *, 'Global average pressure new time step',global_ave_press_new
+
+     call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
      
      if(adsorption_pem) then
       do i=1,ngrid
         global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
-      enddo 
-      print *, 'Global average pressure old time step',global_ave_press_old
-      print *, 'Global average pressure new time step',global_ave_press_new
+       enddo 
+       print *, 'Global average pressure old time step',global_ave_press_old
+       print *, 'Global average pressure new time step',global_ave_press_new
      endif
-
-     call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
 
 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
@@ -916 +918 @@
      call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
 
-     DO islope=1, nslope
-       write(str2(1:2),'(i2.2)') islope
-       call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf_slope(:,igcm_h2o_ice,islope))
-       call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice_phys_slope(:,islope))
-     ENDDO
 
       print *, "Evolution of co2 ice"
      call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
+
+      DO islope=1, nslope
+        write(str2(1:2),'(i2.2)') islope
+        call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf_slope(:,igcm_h2o_ice,islope))
+        call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice_phys_slope(:,islope))
+      ENDDO
 
 !------------------------
@@ -939 +942 @@
                          global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
 
-     DO islope=1, nslope
-       write(str2(1:2),'(i2.2)') islope
-       call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,co2ice_slope(:,islope))
-       call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice_phys_slope(:,islope))
-     ENDDO
+      DO islope=1, nslope
+        write(str2(1:2),'(i2.2)') islope
+        call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,co2ice_slope(:,islope))
+        call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice_phys_slope(:,islope))
+      ENDDO
 
 !------------------------
@@ -1071 +1074 @@
       print *, "Update soil propreties"
 ! II_d.4 Update the soil thermal properties
-      call update_soil(ngrid,nslope,nsoilmx,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
+      call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
         ice_depth,TI_PEM)
 
@@ -1114 +1117 @@
                                    initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
 
-      year_iter=year_iter+dt_pem
+      year_iter=year_iter+dt_pem     
 
       print *, "Checking all the stopping criterion."
       if (STOPPING_water) then
         print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
-        criterion_stop=1
         criterion_stop=1
       endif
@@ -1150 +1152 @@
       endif
 
+
+
       global_ave_press_old=global_ave_press_new
 
@@ -1181 +1185 @@
 
 ! H2O ice
+
      DO ig=1,ngrid
        if(watercaptag(ig)) then
@@ -1193 +1198 @@
            endif
            watercap_sum=watercap_sum+(watercap_slope(ig,islope)-watercap_slope_saved)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-!           watercap_slope(ig,islope)=0.
+           watercap_slope(ig,islope)=0.
          enddo
            water_reservoir(ig)=water_reservoir(ig)+watercap_sum
@@ -1240 +1245 @@
 
    if(soil_pem) then
+     fluxgeo_slope(:,:) = fluxgeo
      call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
      tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
@@ -1279 +1285 @@
       ENDDO
     ENDDO
+
 
       DO ig = 1,ngrid
@@ -1346 +1353 @@
        enddo
      enddo
-
-     DO ig=1,ngrid
-       if(watercaptag(ig)) then
-         WC_sum=0.
-         DO islope=1,nslope
-           if(qsurf_slope(ig,igcm_h2o_ice,islope).GT. (watercap_slope(ig,islope)+water_reservoir(ig))) then
-             qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig))
-           else
-             watercap_slope(ig,islope)=watercap_slope(ig,islope)+qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)/nslope)
-             qsurf_slope(ig,igcm_h2o_ice,islope)=0.
-           endif
-           WC_sum=WC_sum+watercap_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-           watercap_slope(ig,islope)=0.
-         enddo  
-           water_reservoir(ig)=water_reservoir(ig)+WC_sum
-       endif
-     enddo
-
-! H2o ice
-      DO ig = 1,ngrid
-        qsurf(ig,igcm_h2o_ice) = 0.
-        DO islope = 1,nslope
-          qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
-                                   * subslope_dist(ig,islope) /          &
-                                  cos(pi*def_slope_mean(islope)/180.)
-        ENDDO
-      ENDDO ! of DO ig=1,ngrid
-
-     DO ig=1,ngrid
-!       DO islope=1,nslope
-         if(qsurf(ig,igcm_h2o_ice).GT.500) then
-            watercaptag(ig)=.true.
-       DO islope=1,nslope
-            qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-250
-            water_reservoir(ig)=water_reservoir(ig)+250
-       ENDDO            
-         endif
-!       enddo
-     enddo
-
-     DO ig=1,ngrid
-!       DO islope=1,nslope
-         if(water_reservoir(ig).LT. 10) then
-            watercaptag(ig)=.false.
-            qsurf(ig,igcm_h2o_ice)=qsurf(ig,igcm_h2o_ice)+water_reservoir(ig)
-       DO islope=1,nslope
-            qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+water_reservoir(ig)
-       ENDDO
-         endif
-!       enddo
-     enddo
-
-      DO ig = 1,ngrid
-        watercap(ig) = 0.
-        DO islope = 1,nslope
-          watercap(ig) = watercap(ig) + watercap_slope(ig,islope) &
-                                   * subslope_dist(ig,islope) /          &
-                                  cos(pi*def_slope_mean(islope)/180.)
-        ENDDO
-      ENDDO ! of DO ig=1,ngrid
-
      
 !------------------------

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

@@ -5 +5 @@
    
    use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
-   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM,water_reservoir_nom
+   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM,water_reservoir_nom,depth_breccia,depth_bedrock,index_breccia,index_bedrock
    use comsoil_h, only:  volcapa,inertiedat 
    use adsorption_mod, only : regolith_adsorption,adsorption_pem
@@ -91 +91 @@
 write(*,*)'Is start PEM?',startpem_file
 
-startpem_file = .true.
+
 !1. Run
 
@@ -119 +119 @@
 
       do ig = 1,ngrid
-          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then 
+          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then 
 !!! transition
-             delta = 50.
-             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
-            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
-                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
-            do iloop=nsoil_GCM+2,n_1km
+             delta = depth_breccia
+             TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
+            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
+                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
+            do iloop=index_breccia+2,index_bedrock
               TI_PEM(ig,iloop,islope) = TI_breccia
             enddo      
           else ! we keep the high ti values
-            do iloop=nsoil_GCM+1,n_1km
-              TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
+            do iloop=index_breccia+1,index_bedrock
+              TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
             enddo 
           endif ! TI PEM and breccia comparison
 !! transition
-          delta = 1000.
-          TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
-               (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
-               ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
-          do iloop=n_1km+2,nsoil_PEM
+          delta = depth_bedrock
+          TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
+               (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
+               ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
+          do iloop=index_bedrock+2,nsoil_PEM
             TI_PEM(ig,iloop,islope) = TI_bedrock
           enddo   
       enddo
-   else
+   else ! found
      do iloop = nsoil_GCM+1,nsoil_PEM
        TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
@@ -158 +158 @@
  enddo 
 !!! zone de transition
-delta = 50.
+delta = depth_breccia
 do ig = 1,ngrid
-if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then 
-inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
-            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
-                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
-
- do iloop = nsoil_GCM+2,n_1km 
+if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then 
+inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
+            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
+                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
+
+ do iloop = index_breccia+2,index_bedrock 
        inertiedat_PEM(ig,iloop) = TI_breccia
   enddo
 
 else
-   do iloop=nsoil_GCM+1,n_1km
+   do iloop=index_breccia+1,index_bedrock
       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
    enddo 
@@ -177 +177 @@
 
 !!! zone de transition
-delta = 1000.
+delta = depth_bedrock
 do ig = 1,ngrid
-inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
-            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
-                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
+inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
+            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
+                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
 enddo
 
- do iloop = n_1km+2, nsoil_PEM
+ do iloop = index_bedrock+2, nsoil_PEM
    do ig = 1,ngrid
       inertiedat_PEM(ig,iloop) = TI_bedrock
@@ -201 +201 @@
       write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
       write(*,*)'will reconstruct the values of Tsoil'
-      do ig = 1,ngrid
-        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
-        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))   
-
-       do iloop=nsoil_GCM+2,n_1km
-            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
-      enddo
-        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
-        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))  
-
-       do iloop=n_1km+2,nsoil_PEM
-            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
-      enddo
-          do iloop = nsoil_GCM+1,nsoil_PEM
-             tsoil_PEM(ig,iloop,islope) =  tsoil_PEM(ig,nsoil_GCM,islope)
-          enddo
-      enddo
-!     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-!     do iloop = 1,2000000
-!      write(*,*) 'iloop,islope',islope,iloop
-!      call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-!     enddo
+!      do ig = 1,ngrid
+!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
+!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))   
+!       do iloop=index_breccia+2,index_bedrock
+!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
+!      enddo
+!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
+!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))  
+!
+!      do iloop=index_bedrock+2,nsoil_PEM
+!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
+!      enddo
+!      enddo
+
+
+     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+
     
    else
@@ -240 +236 @@
        tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
      enddo
-   endif
+   endif !found
 
     do it = 1,timelen
@@ -263 +259 @@
 
 
-
+  call get_field("ice_table",ice_table,found)
+   if(.not.found) then
+      write(*,*)'PEM settings: failed loading <ice_table>'
+      write(*,*)'will reconstruct the values of the ice table given the current state'
      call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
-     call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
-
+     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
+     do islope = 1,nslope
+     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+     enddo
+   endif
 
 print *,'PEMETAT0: ICE TABLE  DONE'
@@ -338 +340 @@
       do ig = 1,ngrid
 
-          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then 
+          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then 
 !!! transition
-             delta = 50.
-             TI_PEM(ig,nsoil_GCM+1,islope) =sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
-            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
-                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
-
-          do iloop=nsoil_GCM+2,n_1km
+             delta = depth_breccia
+             TI_PEM(ig,index_breccia+1,islope) =sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
+            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
+                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
+
+          do iloop=index_breccia+2,index_bedrock
             TI_PEM(ig,iloop,islope) = TI_breccia
          enddo      
          else ! we keep the high ti values
-           do iloop=nsoil_GCM+1,n_1km
-                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
+           do iloop=index_breccia+1,index_bedrock
+                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
            enddo 
          endif
 
 !! transition
-             delta = 1000.
-             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
-            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
-                      ((layer_PEM(n_1km+1)-delta)/(TI_breccia**2))))
-          do iloop=n_1km+2,nsoil_PEM
+             delta = depth_bedrock
+             TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
+            (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
+                      ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
+          do iloop=index_bedrock+2,nsoil_PEM
             TI_PEM(ig,iloop,islope) = TI_bedrock
          enddo   
@@ -369 +371 @@
  enddo
 !!! zone de transition
-delta = 50.
+delta = depth_breccia 
 do ig = 1,ngrid
-      if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then 
-inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
-            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
-                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
-
-
- do iloop = nsoil_GCM+2,n_1km 
+      if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then 
+inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
+            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
+                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
+
+
+ do iloop = index_breccia+2,index_bedrock 
 
        inertiedat_PEM(ig,iloop) = TI_breccia
@@ -383 +385 @@
  enddo
 else 
-   do iloop = nsoil_GCM+1,n_1km 
-       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
+   do iloop = index_breccia+1,index_bedrock 
+       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
     enddo
 
@@ -392 +394 @@
 
 !!! zone de transition
-delta = 1000.
+delta = depth_bedrock
 do ig = 1,ngrid
-inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
-            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
-                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
+inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
+            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
+                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
 enddo
 
 
 
- do iloop = n_1km+2, nsoil_PEM
+ do iloop = index_bedrock+2, nsoil_PEM
    do ig = 1,ngrid
       inertiedat_PEM(ig,iloop) = TI_bedrock
@@ -414 +416 @@
 !b) Soil temperature 
     do islope = 1,nslope
-     do ig = 1,ngrid
-        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
-        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))
-
-       do iloop=nsoil_GCM+2,n_1km
-            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
-      enddo
-        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
-        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))
-
-       do iloop=n_1km+2,nsoil_PEM
-            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
-      enddo
+!     do ig = 1,ngrid
+!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
+!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
+!
+!       do iloop=index_breccia+2,index_bedrock
+!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
+!      enddo
+!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
+!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
+
+!       do iloop=index_bedrock+2,nsoil_PEM
+!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
+!      enddo
      
-           do iloop = nsoil_GCM+1,nsoil_PEM
-             tsoil_PEM(ig,iloop,islope) =  tsoil_PEM(ig,nsoil_GCM,islope)
-          enddo
-       enddo
-!     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-!     do iloop = 1,2000000
-!      write(*,*) 'islope,iloop',islope,iloop
-!      call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-!     enddo
-  
+!       enddo
+
+      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+      call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+
    
     do it = 1,timelen
@@ -446 +443 @@
         enddo
      enddo
+ 
       do isoil = nsoil_GCM+1,nsoil_PEM
         do ig = 1,ngrid
@@ -451 +449 @@
         enddo
       enddo
-enddo
+enddo !islope
 print *,'PEMETAT0: TSOIL DONE  '
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !c) Ice table    
      call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
-     call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
+     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
+     do islope = 1,nslope
+     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+     enddo
+    
 
 
@@ -483 +485 @@
 endif !soil_pem
 
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    
+
+
 !. e) water reservoir
 
@@ -489 +497 @@
       do ig=1,ngrid
         if(watercaptag(ig)) then
-           water_reservoir=water_reservoir_nom
+           water_reservoir(ig)=water_reservoir_nom
         else
-           water_reservoir=0.
+           water_reservoir(ig)=0.
         endif
       enddo
@@ -498 +506 @@
 
 endif ! of if (startphy_file)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 if(soil_pem) then

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

@@ -119 +119 @@
   call put_field('water_reservoir','water_reservoir', &
            water_reservoir,year_tot)
-  call put_field("water_reservoir","water_reservoir", &
-                 water_reservoir,year_tot)
-
     if(soil_pem) then
   

trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90

@@ -77 +77 @@
 
       allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
-      allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen))
 
   print *, "Opening ", fichnom, "..."
@@ -123 +122 @@
 DO islope=1,nslope
        write(num,fmt='(i2.2)') islope
-       call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))    
+!       call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
+        watercap_slope(:,:,:,:)= 0.
 ENDDO            
      print *, "Downloading data for watercap_slope done"
@@ -187 +187 @@
 ! Compute the minimum over the year for each point
   print *, "Computing the min of h2o_ice_slope"
-  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope+watercap_slope,4)
+!  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope+watercap_slope,4)
+  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
   print *, "Computing the min of co2_ice_slope"
   min_co2_ice_slope(:,:,:)=minval(co2_ice_slope,4)
@@ -210 +211 @@
             min_co2_ice_slope(i,j,islope)  = 0.
           endif
-!          if (min_h2o_ice_slope(i,j,islope).LT.0) then
-!            min_h2o_ice_slope(i,j,islope)  = 0.
-!          endif
+          if (min_h2o_ice_slope(i,j,islope).LT.0) then
+            min_h2o_ice_slope(i,j,islope)  = 0.
+          endif
        ENDDO
     ENDDO

trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90

@@ -15 +15 @@
 !  Note: depths of layers and mid-layers, soil thermal inertia and 
 !        heat capacity are commons in comsoil_PEM.h
-!        A convergence loop is added until equilibrium
 !-----------------------------------------------------------------------
 

trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F

@@ -4 +4 @@
 
 !      use netcdf
-      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
-      use comsoil_h, only: inertiedat,layer,mlayer, volcapa
+      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, 
+     &   depth_breccia,depth_bedrock,index_breccia,index_bedrock  
+      use comsoil_h, only: inertiedat,layer,mlayer, volcapa                       
+
+
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
@@ -86 +89 @@
       enddo
 
+
+! 3. Index for subsurface layering
+! ------------------
+       write(*,*) 'depth=',depth_breccia,depth_bedrock
+      index_breccia= 1
+      do iloop = 1,nsoil_PEM-1
+        if (depth_breccia.ge.layer_PEM(iloop)) then
+          index_breccia=iloop 
+        else
+         exit
+        endif
+      enddo
+
+      index_bedrock= 1
+      do iloop = 1,nsoil_PEM-1
+        if (depth_bedrock.ge.layer_PEM(iloop)) then
+          index_bedrock=iloop 
+        else
+         exit
+        endif
+      enddo
+
+
       end subroutine soil_settings_PEM

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

@@ -1 @@
-   SUBROUTINE update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,TI_PEM)
+   SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,TI_PEM)
 #ifndef CPP_STD
  USE comsoil_h, only:  inertiedat, volcapa
- USE comsoil_h_PEM, only: layer_PEM,n_1km,inertiedat_PEM
+ USE comsoil_h_PEM, only: layer_PEM,n_1km,inertiedat_PEM,depth_breccia,depth_bedrock,index_breccia,index_bedrock,reg_thprop_dependp
  USE vertical_layers_mod, ONLY: ap,bp
  USE comsoil_h_PEM, only: n_1km
  implicit none
 ! Input: 
- INTEGER,INTENT(IN) ::  ngrid, nslope,nsoil_GCM, nsoil_PEM
+ INTEGER,INTENT(IN) ::  ngrid, nslope, nsoil_PEM
  REAL,INTENT(IN) :: p_avg_new
  REAL,INTENT(IN) :: tendencies_waterice(ngrid,nslope)
@@ -41 +41 @@
 
 
-
  do islope = 1,nslope
    regolith_inertia(:,islope) = inertiedat_PEM(:,1)
    do ig = 1,ngrid
-
       if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then
               regolith_inertia(ig,islope) = inertie_averaged
        endif
-        write(*,*) 'ig,islope',ig,islope,inertie_thresold,TI_PEM(ig,1,islope)
-       if(TI_PEM(ig,1,islope).lt.inertie_thresold) then
-!          regolith_inertia(ig,islope) = regolith_inertia(ig,islope)*(p_avg_new/P610)**0.3
-       endif 
-       TI_breccia_new = TI_breccia !*(p_avg_new/P610)**0.3
+      if(reg_thprop_dependp) then
+          if(TI_PEM(ig,1,islope).lt.inertie_thresold) then
+            regolith_inertia(ig,islope) = regolith_inertia(ig,islope)*(p_avg_new/P610)**0.3
+          endif 
+         TI_breccia_new = TI_breccia*(p_avg_new/P610)**0.3
+      else 
+         TI_breccia_new = TI_breccia
+      endif
    enddo
  enddo
@@ -62 +63 @@
 do  islope=1,nslope
    do ig = 1,ngrid
-     do iloop = 1,nsoil_GCM
+     do iloop = 1,index_breccia
          TI_PEM(ig,iloop,islope) = regolith_inertia(ig,islope)
      enddo
      if(regolith_inertia(ig,islope).lt.TI_breccia_new) then 
 !!! transition
-             delta = 50.
-             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
-            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
-                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia_new**2))))
-            do iloop=nsoil_GCM+2,n_1km
+             delta = depth_breccia
+             TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
+            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
+                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia_new**2))))
+            do iloop=index_breccia+2,index_bedrock
               TI_PEM(ig,iloop,islope) = TI_breccia_new
             enddo      
       else ! we keep the high ti values
-            do iloop=nsoil_GCM+1,n_1km
-              TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
+            do iloop=index_breccia+1,index_bedrock
+              TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
             enddo 
        endif ! TI PEM and breccia comparison
 !! transition
-       delta = 1000.
-       TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
-            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
-            ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
-       do iloop=n_1km+2,nsoil_PEM
+       delta = depth_bedrock
+       TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
+            (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
+            ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
+       do iloop=index_bedrock+2,nsoil_PEM
             TI_PEM(ig,iloop,islope) = TI_bedrock
        enddo   
@@ -94 +95 @@
   do ig=1,ngrid
     do islope=1,nslope
-      do iloop = 1,n_1km
-       TI_PEM(ig,iloop,islope) = TI_PEM(ig,1,islope)
-      enddo
+!      do iloop = 1,index_bedrock
+!       TI_PEM(ig,iloop,islope) = TI_PEM(ig,1,islope)
+!      enddo
   if (ice_depth(ig,islope).gt. -1.e-10) then
 ! 3.0 FIrst if permanent ice
    if (ice_depth(ig,islope).lt. 1e-10) then
       do iloop = 1,nsoil_PEM
-       TI_PEM(ig,iloop,islope)=max(ice_inertia,inertiedat_PEM(ig,iloop))
+       TI_PEM(ig,iloop,islope)=ice_inertia
       enddo
      else 

trunk/LMDZ.MARS/README

@@ -3874 +3874 @@
 If false, it is set to constant (to the value of the tab_cntrl in the start)
 
+== 10/02/2023 == LL
+PEM
+Soil temperature initialisation has been updated
+Conf_PEM improved by adding some options to the users (thermal regolith depend on the pressure, depth of the subsurface layers, etc.)
+Minor edits then (+ svn update with RV had some issues, so there are some "artefact changes" ...)
+