Changeset 3130 for trunk/LMDZ.COMMON/libf/evolution

Timestamp:

Nov 20, 2023, 1:25:31 PM (12 months ago)

Author:

jbclement

Message:

PEM:
The perennial co2 ice is now taken into account with co2 frost (qsurf) to compute the tendency and to make the update + Rework of how co2 frost is converted to perennial co2 ice at the end of the PEM run + Correction of the value of 'threshold_co2_frost2perennial' to correspond to 10 m + Perennial co2 ice is now handled outside 'paleoclimate' in "phyetat0_mod.F90" of the Mars PCM + Some cleanings.

/!\ Commit for the PEM management of co2 ice before a rework of ice management in the PEM!
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• adsorption_mod.F90 (modified) (2 diffs)
• changelog.txt (modified) (1 diff)
• compute_tendencies_slope_mod.F90 (modified) (1 diff)
• constants_marspem_mod.F90 (modified) (1 diff)
• criterion_pem_stop_mod.F90 (modified) (3 diffs)
• evol_co2_ice_s_mod.F90 (modified) (2 diffs)
• evol_h2o_ice_s_mod.F90 (modified) (7 diffs)
• pem.F90 (modified) (16 diffs)
• read_data_PCM_mod.F90 (modified) (8 diffs)
• recomp_tend_co2_slope_mod.F90 (modified) (1 diff)
• soil_thermalproperties_mod.F90 (modified) (1 diff)

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

@@ -156 +156 @@
 #ifndef CPP_STD
 
-!0.1 Look at perenial ice
+!0.1 Look at perennial ice
   do ig = 1,ngrid
     do islope = 1,nslope
@@ -317 +317 @@
 #ifndef CPP_STD
 
-!0.1 Look at perenial ice
+!0.1 Look at perennial ice
   do ig = 1,ngrid
     do islope = 1,nslope

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

@@ -137 +137 @@
 Adapting the PEM soil grid to the one of the PCM
 Minor corrections when reading/initializing soil temperature, subsurface water ice
+
+== 20/11/2023 == JBC
+The perennial co2 ice is now taken into account with co2 frost (qsurf) to compute the tendency and to make the update + Rework of how co2 frost is converted to perennial co2 ice at the end of the PEM run + Correction of the value of 'threshold_co2_frost2perennial' to correspond to 10 m + Perennial co2 ice is now handled outside 'paleoclimate' in "phyetat0_mod.F90" of the Mars PCM + Some cleanings.
+
+/!\ Commit for the PEM management of co2 ice before a rework of ice management in the PEM!

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

@@ -25 +25 @@
 
 !   OUTPUT
-real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of perenial ice
+real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of perennial ice
 !=======================================================================
 

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

@@ -37 +37 @@
 real, parameter :: Lco2 =  5.71e5 ! Pilorget and Forget 2016
 
-! Conversion H2O/CO2 frost to perenial frost and vice versa
-real, parameter :: threshold_water_frost2perenial = 1000. !~ 1m
-real, parameter :: threshold_co2_frost2perenial = 3200.   !~ 2m
+! Conversion H2O/CO2 frost to perennial frost and vice versa
+real, parameter :: threshold_water_frost2perennial = 1000. !~ 1 m
+real, parameter :: threshold_co2_frost2perennial = 16000.   !~ 10 m
 
 END MODULE constants_marspem_mod

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

@@ -1 @@
-  module criterion_pem_stop_mod
-  implicit none
+MODULE criterion_pem_stop_mod
 
-  contains
+implicit none
 
+contains
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
-!!! Purpose: Criterions to check if the PEM needs to call the GCM !!!
+!!! Purpose: Criterions to check if the PEM needs to call the PCM
 !!! Author: RV & LL, 02/2023
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-
 SUBROUTINE criterion_waterice_stop(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice)
 
-  use time_evol_mod, only: water_ice_criterion
-  use comslope_mod,  only: subslope_dist,nslope
+use time_evol_mod, only: water_ice_criterion
+use comslope_mod,  only: subslope_dist, nslope
 
-      IMPLICIT NONE
+implicit none
 
 !=======================================================================
@@ -28 +27 @@
 !   arguments:
 !   ----------
-
 !   INPUT
-  INTEGER, intent(in) :: ngrid                         ! # of grid physical grid points 
-  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
-  REAL,    intent(in) :: qsurf(ngrid,nslope)           ! physical point field : Actual density of water ice
-  REAL,    intent(in) :: ini_surf                      ! Initial surface of h2o ice that was sublimating
-  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope) ! Grid point that initialy were covered by h2o_ice
-
+integer,                       intent(in) :: ngrid           ! # of physical grid points
+real, dimension(ngrid),        intent(in) :: cell_area       ! Area of the cells
+real, dimension(ngrid,nslope), intent(in) :: qsurf           ! Actual density of water ice
+real,                          intent(in) :: ini_surf        ! Initial surface of h2o ice that was sublimating
+real, dimension(ngrid,nslope), intent(in) :: initial_h2o_ice ! Grid point that initialy were covered by h2o_ice
 !   OUTPUT
-  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
-
+logical, intent(out) :: STOPPING ! Is the criterion reached?
 !   local:
-!   -----
-  INTEGER :: i,islope   ! Loop
-  REAL :: present_surf  ! Initial/Actual surface of water ice
+!   ------
+integer :: i, islope    ! Loop
+real    :: present_surf ! Initial/Actual surface of water ice
 
 !=======================================================================
+! Initialisation to false
+STOPPING = .false.
 
-!   initialisation to false
-    STOPPING=.FALSE.
+! Computation of the present surface of water ice sublimating
+present_surf = 0.
+do i = 1,ngrid
+    do islope = 1,nslope
+        !if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
+        if (initial_h2o_ice(i,islope) > 0.5) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
+    enddo
+enddo
 
-!   computation of the present surface of water ice sublimating
-  present_surf=0.
-  do i=1,ngrid
-    do islope=1, nslope
-      if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
-         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
-      endif
-    enddo
-  enddo
-  
-!   check of the criterion
-  if(present_surf.LT.ini_surf*(1-water_ice_criterion) .OR. &
-     present_surf.GT.ini_surf*(1+water_ice_criterion)) then
-    STOPPING=.TRUE. 
-    write(*,*) "Reason of stopping : The surface of water ice sublimating reach the threshold:"
-    write(*,*) "Current surface of water ice sublimating=", present_surf
-    write(*,*) "Initial surface of water ice sublimating=", ini_surf
-    write(*,*) "Percentage of change accepted=", water_ice_criterion*100
-    write(*,*) "present_surf<ini_surf*(1-water_ice_criterion)", (present_surf.LT.ini_surf*(1-water_ice_criterion))
-  endif
+! Check of the criterion
+if (present_surf < ini_surf*(1. - water_ice_criterion) .or. present_surf > ini_surf*(1. + water_ice_criterion)) then
+    STOPPING = .true.
+    write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
+    write(*,*) "Current surface of water ice sublimating =", present_surf
+    write(*,*) "Initial surface of water ice sublimating =", ini_surf
+    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
+    write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", (present_surf < ini_surf*(1. - water_ice_criterion))
+endif
 
-  if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
-    STOPPING=.FALSE.
-  endif
+if (ini_surf < 1.e-5 .and. ini_surf > -1.e-5) STOPPING = .false.
+
 END SUBROUTINE criterion_waterice_stop
 
-! ------------------------------------------------------------------------------------------------
+! ----------------------------------------------------------------------
 
 SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_co2_ice,global_ave_press_GCM,global_ave_press_new,nslope)
 
-  use time_evol_mod, only: co2_ice_criterion,ps_criterion
-  use comslope_mod,  only: subslope_dist
+use time_evol_mod, only: co2_ice_criterion, ps_criterion
+use comslope_mod,  only: subslope_dist
 
-      IMPLICIT NONE
+implicit none
 
 !=======================================================================
 !
-!  Routine that checks if the criterion to stop the PEM is reached
+!  Routine that checks if the co2 and pressure criteria to stop the PEM are reached
 !
 !=======================================================================
@@ -92 +84 @@
 !   arguments:
 !   ----------
-
 !   INPUT
-  INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points 
-  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
-  REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
-  REAL,    intent(in) :: ini_surf                      ! Initial surface of co2 ice that was sublimating
-  REAL,    intent(in) :: initial_co2_ice(ngrid,nslope) ! Grid point that initialy were covered by co2_ice
-  REAL,    intent(in) :: global_ave_press_GCM          ! Planet average pressure from the GCM start files 
-  REAL,    intent(in) :: global_ave_press_new          ! Planet average pressure from the PEM computations
-
+integer,                       intent(in) :: ngrid, nslope        ! # of grid physical grid points
+real, dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
+real, dimension(ngrid,nslope), intent(in) :: qsurf                ! Actual density of water ice
+real,                          intent(in) :: ini_surf             ! Initial surface of co2 ice that was sublimating
+real, dimension(ngrid,nslope), intent(in) :: initial_co2_ice      ! Grid point that initialy were covered by co2_ice
+real,                          intent(in) :: global_ave_press_GCM ! Planet average pressure from the PCM start files
+real,                          intent(in) :: global_ave_press_new ! Planet average pressure from the PEM computations
 !   OUTPUT
-  LOGICAL, intent(out) :: STOPPING_ice              ! Logical : is the criterion for ice reached?
-  LOGICAL, intent(out) :: STOPPING_ps               ! Logical : is the criterion for pressure reached ?
+logical, intent(out) :: STOPPING_ice ! Is the criterion for co2 ice reached?
+logical, intent(out) :: STOPPING_ps  ! Is the criterion for pressure reached?
 !   local:
-!   -----
-  INTEGER :: i,islope   ! Loop
-  REAL :: present_surf  ! Initial/Actual surface of water ice
+!   ------
+integer :: i, islope    ! Loop
+real    :: present_surf ! Initial/Actual surface of water ice
 
 !=======================================================================
+! Initialisation to false
+STOPPING_ice = .false.
+STOPPING_ps = .false.
 
-!   initialisation to false
-    STOPPING_ice=.FALSE.
-    STOPPING_ps =.FALSE.
-!   computation of the actual surface
-  present_surf=0.
-  do i=1,ngrid
-   do islope=1,nslope
-      if (initial_co2_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
-         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
-      endif
-   enddo
-  enddo
-  
-!   check of the criterion
-  if(present_surf.LT.ini_surf*(1-co2_ice_criterion) .OR. &
-     present_surf.GT.ini_surf*(1+co2_ice_criterion)) then
-    STOPPING_ice=.TRUE. 
-    write(*,*) "Reason of stopping : The surface of co2 ice sublimating reach the threshold:"
-    write(*,*) "Current surface of co2 ice sublimating=", present_surf
-    write(*,*) "Initial surface of co2 ice sublimating=", ini_surf
-    write(*,*) "Percentage of change accepted=", co2_ice_criterion*100
-    write(*,*) "present_surf<ini_surf*(1-co2_ice_criterion)", (present_surf.LT.ini_surf*(1-co2_ice_criterion))
-  endif
+! Computation of the actual surface
+present_surf = 0.
+do i = 1,ngrid
+    do islope = 1,nslope
+        if (initial_co2_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
+    enddo
+enddo
 
-  if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
-       STOPPING_ice=.FALSE.
-  endif
+! Check of the criterion
+if (present_surf < ini_surf*(1. - co2_ice_criterion) .or. present_surf > ini_surf*(1. + co2_ice_criterion)) then
+    STOPPING_ice = .true.
+    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold"
+    write(*,*) "Current surface of co2 ice sublimating =", present_surf
+    write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
+    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
+    write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", (present_surf < ini_surf*(1. - co2_ice_criterion))
+endif
 
-  if(global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion) .OR. &
-     global_ave_press_new.GT.global_ave_press_GCM*(1+ps_criterion)) then
-    STOPPING_ps=.TRUE. 
-    write(*,*) "Reason of stopping : The global pressure reach the threshold:"
-    write(*,*) "Current global pressure=", global_ave_press_new
-    write(*,*) "GCM global pressure=", global_ave_press_GCM
-    write(*,*) "Percentage of change accepted=", ps_criterion*100
-    write(*,*) "global_ave_press_new<global_ave_press_GCM*(ps_criterion)", (global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion))
-  endif
+if (abs(ini_surf) < 1.e-5) STOPPING_ice = .false.
+
+if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion) .or. global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then
+    STOPPING_ps = .true.
+    write(*,*) "Reason of stopping: the global pressure reach the threshold"
+    write(*,*) "Current global pressure =", global_ave_press_new
+    write(*,*) "PCM global pressure =", global_ave_press_GCM
+    write(*,*) "Percentage of change accepted =", ps_criterion*100.
+    write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion))
+endif
 
 END SUBROUTINE criterion_co2_stop
 
-
 END MODULE

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

@@ -1 +1 @@
 MODULE evol_co2_ice_s_mod
 
-IMPLICIT NONE
+implicit none
 
-CONTAINS
+!=======================================================================
+contains
+!=======================================================================
 
-SUBROUTINE evol_co2_ice_s(qsurf,tendencies_co2_ice_phys,&
-                             iim_input,jjm_input,ngrid,cell_area,nslope)
+SUBROUTINE evol_co2_ice_s(qsurf,tendencies_co2_ice_phys,iim_input,jjm_input,ngrid,cell_area,nslope)
 
-  use time_evol_mod, only: dt_pem
+use time_evol_mod, only: dt_pem
 
-      IMPLICIT NONE
+implicit none
 
 !=======================================================================
@@ -17 +18 @@
 !
 !=======================================================================
-
 !   arguments:
 !   ----------
-
 !   INPUT
-
-  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope      ! # of grid points along longitude/latitude/ total
-  REAL, intent(in) ::  cell_area(ngrid)
+integer,                intent(in) :: iim_input, jjm_input, ngrid, nslope ! # of grid points along longitude/latitude/ total
+REAL, dimension(ngrid), intent(in) :: cell_area
 
 !   OUTPUT
-  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                   ! physical point field : Previous and actual density of water ice
-  REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
+real, dimension(ngrid,nslope), intent(inout) :: qsurf                   ! physical point field: Previous and actual density of water ice
+real, dimension(ngrid,nslope), intent(inout) :: tendencies_co2_ice_phys ! physical point field: Evolution of perennial ice over one year
 
 !   local:
 !   ----
-
-  INTEGER :: i,j,ig0,islope                                     ! loop variable
-
+integer :: i, j, ig0, islope ! loop variable
+!=======================================================================
 ! Evolution of the CO2 ice for each physical point
-  do i=1,ngrid
+do i=1,ngrid
     do islope=1,nslope
-      qsurf(i,islope)=qsurf(i,islope)+tendencies_co2_ice_phys(i,islope)*dt_pem
-      if (qsurf(i,islope).lt.0) then
-        qsurf(i,islope)=0.
-        tendencies_co2_ice_phys(i,islope)=0.
-      endif
+        qsurf(i,islope) = qsurf(i,islope) + tendencies_co2_ice_phys(i,islope)*dt_pem
+        if (qsurf(i,islope) < 0) then
+            qsurf(i,islope) = 0.
+            tendencies_co2_ice_phys(i,islope) = 0.
+        endif
     enddo
-  enddo
+enddo
 
 END SUBROUTINE evol_co2_ice_s

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

@@ -1 +1 @@
 MODULE evol_h2o_ice_s_mod
 
-IMPLICIT NONE
+implicit none
 
-CONTAINS
+!=======================================================================
+contains
+!=======================================================================
 
 SUBROUTINE evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf,tendencies_h2o_ice_phys,STOPPING)
 
-  use time_evol_mod,          only: dt_pem
-  use comslope_mod,           only: subslope_dist, def_slope_mean
-  use criterion_pem_stop_mod, only: criterion_waterice_stop
+use time_evol_mod,          only: dt_pem
+use comslope_mod,           only: subslope_dist, def_slope_mean
+use criterion_pem_stop_mod, only: criterion_waterice_stop
+
 #ifndef CPP_STD
     use comcstfi_h,   only: pi
@@ -16 +19 @@
 #endif
 
-  IMPLICIT NONE
+implicit none
 
 !=======================================================================
@@ -23 +26 @@
 !
 !=======================================================================
-
 !   arguments:
 !   ----------
-
 !   INPUT
-
-  INTEGER, intent(in) :: ngrid                                  ! # of grid points along longitude/latitude grid;
-  INTEGER, intent(in) :: nslope                                 ! # of subslope
-  REAL, intent(in) ::  cell_area(ngrid)                         ! Area of each mesh grid (m^2)
-  REAL, intent(in) :: delta_h2o_adsorbded(ngrid)                ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2) 
-  REAL, intent(in) :: delta_h2o_icetablesublim(ngrid)                ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2) 
+integer,                intent(in) :: ngrid                    ! # of grid points along longitude/latitude grid;
+integer,                intent(in) :: nslope                   ! # of subslope
+real, dimension(ngrid), intent(in) :: cell_area                ! Area of each mesh grid (m^2)
+real, dimension(ngrid), intent(in) :: delta_h2o_adsorbded      ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2)
+real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2)
 
 !   OUTPUT
-  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                   ! physical point field : Previous and actual density of water ice (kg/m^2)
-  REAL, intent(inout) ::  tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year (kg/m^2/year)
-  LOGICAL, INTENT(INOUT) :: STOPPING                            ! Stopping criterion
+real, dimension(ngrid,nslope), intent(inout) :: qsurf                   ! physical point field: Previous and actual density of water ice (kg/m^2)
+real, dimension(ngrid,nslope), intent(inout) :: tendencies_h2o_ice_phys ! physical point field: Evolution of perennial ice over one year (kg/m^2/year)
+logical,                       intent(inout) :: STOPPING                ! Stopping criterion
 
 !   local:
-!   ----
-
-  INTEGER :: i,j,islope                                         ! loop variable
-  REAL :: pos_tend, neg_tend, real_coefficient,negative_part    ! Variable to conserve water
-  REAL ::  new_tendencies(ngrid,nslope)                         ! Tendencies computed in order to conserve water ice on the surface, only exchange between surface are done
-  
+!   ------
+integer :: i,j,islope                                          ! loop variable
+real    :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve water
+real    ::  new_tendencies(ngrid,nslope)                       ! Tendencies computed in order to conserve water ice on the surface, only exchange between surface are done
 !=======================================================================
 
@@ -53 +51 @@
   pos_tend=0.
   neg_tend=0.
-if (ngrid.NE.1) then ! to make sure we are not in 1D 
+if (ngrid.NE.1) then ! to make sure we are not in 1D
   ! We compute the amount of water accumulating and sublimating
   do i=1,ngrid
@@ -77 +75 @@
   enddo
   ! We adapt the tendencies to conserve water and do only exchange between grid points
-   if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then ! We are sublimating more in the planet than condensing 
+   if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then ! We are sublimating more in the planet than condensing
      do i=1,ngrid
        do islope=1,nslope
@@ -97 +95 @@
        enddo
      enddo
-   elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then 
+   elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then
     write(*,*) "Reason of stopping : There is either no water ice sublimating or no water ice increasing !!"
     write(*,*) "Tendencies on ice sublimating=", neg_tend
@@ -123 +121 @@
     enddo
   enddo
-  
-  
+
+
   if(pos_tend.eq.0) then
    real_coefficient = 0.
-  else 
-   real_coefficient = negative_part/pos_tend ! We compute a coefficient by which we should remove the ice that has been added 
-                                             ! to places even if this ice was contributing to an unphysical negative amount 
+  else
+   real_coefficient = negative_part/pos_tend ! We compute a coefficient by which we should remove the ice that has been added
+                                             ! to places even if this ice was contributing to an unphysical negative amount
                                              ! of ice at other places
   endif

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

@@ -42 +42 @@
 use criterion_pem_stop_mod,       only: criterion_waterice_stop, criterion_co2_stop
 use constants_marspem_mod,        only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
-                                        m_noco2, threshold_water_frost2perenial, threshold_co2_frost2perenial
+                                        m_noco2, threshold_water_frost2perennial, threshold_co2_frost2perennial
 use evol_co2_ice_s_mod,           only: evol_co2_ice_s
 use evol_h2o_ice_s_mod,           only: evol_h2o_ice_s
@@ -77 +77 @@
                                   albedodat, zmea, zstd, zsig, zgam, zthe,     &
                                   hmons, summit, base,albedo_h2o_frost,        &
-                                  frost_albedo_threshold, emissiv, watercaptag, perenial_co2ice, &
+                                  frost_albedo_threshold, emissiv, watercaptag, perennial_co2ice, &
                                   emisice, albedice
     use dimradmars_mod,     only: totcloudfrac, albedo
@@ -84 +84 @@
     use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
     use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
-    use paleoclimate_mod,   only: albedo_perenialco2
+    use paleoclimate_mod,   only: albedo_perennialco2
     use comcstfi_h,         only: pi, rad, g, cpp, mugaz, r
 #else
@@ -161 +161 @@
 ! Variables for h2o_ice evolution
 real                                :: ini_surf_h2o         ! Initial surface of sublimating h2o ice
-real                                :: ini_surf_co2         ! Initial surface of sublimating co2 ice
 real                                :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa]
 real                                :: global_ave_press_old ! constant: Global average pressure of initial/previous time step
@@ -171 +170 @@
 logical                             :: STOPPING_water       ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached?
 logical                             :: STOPPING_1_water     ! Logical: is there still water ice to sublimate?
-logical                             :: STOPPING_co2         ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached?
 logical                             :: STOPPING_pressure    ! Logical: is the criterion (% of change in the surface pressure) reached?
 integer                             :: criterion_stop       ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
 real, save                          :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
-real, dimension(:,:),   allocatable :: vmr_co2_gcm          ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
-real, dimension(:,:),   allocatable :: vmr_co2_pem_phys     ! Physics x Times  co2 volume mixing ratio used in the PEM
-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]
 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]
 integer                             :: timelen              ! # time samples
@@ -183 +178 @@
 real, dimension(:,:),   allocatable :: p                    ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
 real                                :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
+
+! Variables for co2_ice evolution
+real                              :: ini_surf_co2         ! Initial surface of sublimating co2 ice
+logical                           :: STOPPING_co2         ! Logical: is the criterion (% of change in the surface of sublimating co2 ice) reached?
+real, dimension(:,:), allocatable :: vmr_co2_gcm          ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
+real, dimension(:,:), allocatable :: vmr_co2_pem_phys     ! Physics x Times  co2 volume mixing ratio used in the PEM
+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]
+real, dimension(:,:), allocatable :: co2frost_ini         ! Initial amount of co2 frost (at the start of the PEM run)
+real, dimension(:,:), allocatable :: perennial_co2ice_ini ! Initial amoun of perennial co2 ice (at the start of the PEM run)
+
+
 
 ! Variables for slopes
@@ -193 +199 @@
 real, dimension(:,:),   allocatable :: initial_h2o_ice        ! physical point field: Logical array indicating if there is water ice at initial state
 real, dimension(:,:),   allocatable :: initial_co2_ice        ! physical point field: Logical array indicating if there is co2 ice at initial state
-real, dimension(:,:),   allocatable :: tendencies_co2_ice     ! physical point x slope field: Tendency of evolution of perenial co2 ice over a year
-real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the PCM
-real, dimension(:,:),   allocatable :: tendencies_h2o_ice     ! physical point x slope  field: Tendency of evolution of perenial h2o ice
+real, dimension(:,:),   allocatable :: tendencies_co2_ice     ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
+real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perennial co2 ice over a year in the PCM
+real, dimension(:,:),   allocatable :: tendencies_h2o_ice     ! physical point x slope  field: Tendency of evolution of perennial h2o ice
 real, dimension(:,:),   allocatable :: flag_co2flow           ! (ngrid,nslope): Flag where there is a CO2 glacier flow
 real, dimension(:),     allocatable :: flag_co2flow_mesh      ! (ngrid)       : Flag where there is a CO2 glacier flow
@@ -383 +389 @@
     call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
                   tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,       &
-                  watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist)
+                  watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
 
     ! Remove unphysical values of surface tracer
@@ -588 +594 @@
         if (tendencies_h2o_ice(i,islope) < 0) then
             initial_h2o_ice(i,islope) = 1.
-            ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
+            ini_surf_h2o = ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
         else
             initial_h2o_ice(i,islope) = 0.
@@ -595 +601 @@
 enddo
 
-write(*,*) "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
-write(*,*) "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
+write(*,*) "Total initial surface of co2 ice sublimating (slope) =", ini_surf_co2
+write(*,*) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o
 write(*,*) "Total surface of the planet=", Total_surface
 allocate(zplev_gcm(ngrid,nlayer + 1))
@@ -625 +631 @@
 delta_h2o_icetablesublim(:) = 0.
 
+! We save the initial values for the co2 frost and perennial ice
+allocate(co2frost_ini(ngrid,nslope),perennial_co2ice_ini(ngrid,nslope))
+co2frost_ini = qsurf(:,igcm_co2,:)
+perennial_co2ice_ini = perennial_co2ice
+
 do ig = 1,ngrid
     do islope = 1,nslope
         qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
+        qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) + perennial_co2ice(ig,islope)
     enddo
 enddo
@@ -1000 +1012 @@
 !               2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
                 watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
-                qsurf(ig,igcm_h2o_ice,islope)=0.
+                qsurf(ig,igcm_h2o_ice,islope) = 0.
             endif
-            watercap_sum = watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+            watercap_sum = watercap_sum + watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
             watercap(ig,islope) = 0.
         enddo
@@ -1015 +1027 @@
     enddo
     if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
-        if (water_sum > threshold_water_frost2perenial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
+        if (water_sum > threshold_water_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
             watercaptag(ig) = .true.
-            water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2perenial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
+            water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
             do islope = 1,nslope
-                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.)
+                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.)
             enddo
         endif
     else ! let's check that the infinite source of water disapear
-        if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perenial) then
+        if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perennial) then
             watercaptag(ig) = .false.
             do islope = 1,nslope
@@ -1036 +1048 @@
 do ig = 1,ngrid
     do islope = 1,nslope
-        if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then
-            perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
-            qsurf(ig,igcm_co2,islope)  = 0.5*qsurf(ig,igcm_co2,islope)
-            albedo(ig,1,islope) = albedo_perenialco2
-            albedo(ig,2,islope) = albedo_perenialco2
+        if (qsurf(ig,igcm_co2,islope) == threshold_co2_frost2perennial) then
+        ! If co2 ice is equal to the threshold, then everything is transformed into perennial ice
+            perennial_co2ice(ig,islope) = threshold_co2_frost2perennial
+            qsurf(ig,igcm_co2,islope) = 0.
+        else if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perennial) then
+        ! If co2 ice is superior to the threshold, then co2 frost is equal to the threshold (max possible value)
+        ! and the leftover is transformed into perennial ice
+            perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope) - threshold_co2_frost2perennial
+            qsurf(ig,igcm_co2,islope) = threshold_co2_frost2perennial
+        else
+        ! If co2 ice is inferior to the threshold, then we compare with the initial state
+            if (qsurf(ig,igcm_co2,islope) > perennial_co2ice_ini(ig,islope)) then
+            ! If co2 ice higher than the initial perennial ice, then the change is affected only to the frost
+                perennial_co2ice(ig,islope) = perennial_co2ice_ini(ig,islope)
+                qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) - perennial_co2ice_ini(ig,islope)
+            else
+            ! If co2 ice is lower than the initial perennial ice, then there is no frost anymore
+                perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope)
+                qsurf(ig,igcm_co2,islope) = 0.
+            endif
         endif
     enddo
@@ -1124 +1151 @@
                   nlayer,nq,ptimestep,pday,0.,cell_area,albedodat,      &
                   inertiedat,def_slope,subslope_dist)
-
     call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
                   ptimestep,ztime_fin,tsurf,tsoil,inertiesoil,        &
                   albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
-                  wstar,watercap,perenial_co2ice)
+                  wstar,watercap,perennial_co2ice)
 #else
     call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
                   nlayer,nq,ptimestep,pday,time_phys,cell_area,         &
                   albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
-
     call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,  &
                   ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
@@ -1162 +1187 @@
 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
+deallocate(co2frost_ini,perennial_co2ice_ini)
 !----------------------------- END OUTPUT ------------------------------
 

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

@@ -66 +66 @@
 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn             ! h2o ice per slope of the concatenated file [kg/m^2]
 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: perennial_co2ice
 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm               ! CO2 volume mixing ratio in the first layer [mol/m^3]
 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                    ! Surface Pressure [Pa]
@@ -106 +107 @@
     call get_var3("watercap",watercap(:,:,1,:))
     write(*,*) "Data for watercap downloaded!"
+
+    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
+    write(*,*) "Data for perennial_co2ice downloaded!"
 #endif
 
@@ -142 +146 @@
     enddo
     write(*,*) "Data for watercap downloaded!"
+
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
+    enddo
+    write(*,*) "Data for perennial_co2ice downloaded!"
 #endif
 
@@ -179 +189 @@
 write(*,*) "Computing the min of co2_ice..."
 where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
-min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn,4)
+min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn + perennial_co2ice,4)
 
 ! Compute averages over the year for each point
@@ -204 +214 @@
             endif
             if (q_h2o_dyn(i,j,t) < 0) then
-                q_h2o_dyn(i,j,t) = 1.e-30
+                q_h2o_dyn(i,j,t) = 1.e-10
             else if (q_h2o_dyn(i,j,t) > 1) then
                 q_h2o_dyn(i,j,t) = 1.
@@ -282 +292 @@
 !=======================================================================
 
-SUBROUTINE get_var1(filename,v)
-
-implicit none
-
-character(len = *), intent(in)  :: filename
+SUBROUTINE get_var1(var,v)
+
+implicit none
+
+character(len = *), intent(in)  :: var
 real, dimension(:), intent(out) :: v
 
-call error_msg(NF90_INQ_VARID(fID,filename,vID),"inq",filename)
-call error_msg(NF90_GET_VAR(fID,vID,v),"get",filename)
+call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
+call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
 
 END SUBROUTINE get_var1
@@ -296 +306 @@
 !=======================================================================
 
-SUBROUTINE get_var3(filename,v) ! on U grid
-
-implicit none
-
-character(len = *),     intent(in)  :: filename
+SUBROUTINE get_var3(var,v) ! on U grid
+
+implicit none
+
+character(len = *),     intent(in)  :: var
 real, dimension(:,:,:), intent(out) :: v
 
-call error_msg(NF90_INQ_VARID(fID,filename,vID),"inq",filename)
-call error_msg(NF90_GET_VAR(fID,vID,v),"get",filename)
+call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
+call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
 
 END SUBROUTINE get_var3
@@ -310 +320 @@
 !=======================================================================
 
-SUBROUTINE get_var4(filename,v)
-
-implicit none
-
-character(len = *),       intent(in)  :: filename
+SUBROUTINE get_var4(var,v)
+
+implicit none
+
+character(len = *),       intent(in)  :: var
 real, dimension(:,:,:,:), intent(out) :: v
 
-call error_msg(NF90_INQ_VARID(fID,filename,vID),"inq",filename)
-call error_msg(NF90_GET_VAR(fID,vID,v),"get",filename)
+call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
+call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
 
 END SUBROUTINE get_var4

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

@@ -29 +29 @@
 real,                           intent(in) :: global_ave_press_GCM        ! global averaged pressure at previous timestep
 real,                           intent(in) :: global_ave_press_new        ! global averaged pressure at current timestep
-real, dimension(ngrid,nslope),  intent(in) :: tendencies_co2_ice_phys_ini ! physical point field : Evolution of perenial ice over one year
+real, dimension(ngrid,nslope),  intent(in) :: tendencies_co2_ice_phys_ini ! physical point field : Evolution of perennial ice over one year
 real, dimension(ngrid,nslope),  intent(in) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope(kg/m^2)
 real, dimension(ngrid,nslope),  intent(in) :: emissivity_slope            ! Emissivity per mesh and sub-grid slope(1)
 !   OUTPUT
-real, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
+real, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perennial ice over one year
 
 !   local:

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

@@ -147 +147 @@
     do islope=1,nslope
      if (ice_depth(ig,islope).gt.-1e-6) then
-! 3.0 Case where it is perenial ice
+! 3.0 Case where it is perennial ice
       if (ice_depth(ig,islope).lt. 1e-10) then
        call ice_thermal_properties(.true.,1.,0.,ice_inertia)