Changeset 3498 for trunk


Ignore:
Timestamp:
Nov 7, 2024, 2:48:08 PM (2 weeks ago)
Author:
jbclement
Message:

PEM:

  • Correction of the variable name for the ice table depth in "pemetat0.F90". So it is now got as intended from the "startpem.nc" file;
  • Renaming of the tendencies in the PEM with the prefix 'd_' instead of 'tend_';
  • Modification of the PEM time step type from integer to real. As a consequence, all time variables are now of real type. This change adds the possibility to consider fractions of year as time step.

JBC

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

Legend:

Unmodified
Added
Removed

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

    r3493 r3498  
    2727      enddo
    2828!      return
    29       END
     29      END SUBROUTINE tridag
    3030!  (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31&#5(+.

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

    r3456 r3498  
    4040
    4141!=======================================================================
    42 SUBROUTINE regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
     42SUBROUTINE regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
    4343                               m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg)
    4444
    … …  
    4747! Inputs
    4848integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM, timelen  ! size dimension: physics x subslope x soil x timeseries
    49 real,    dimension(ngrid,nslope),           intent(in) :: tend_h2oglaciers, tend_co2glaciers ! tendancies on the glaciers [1]
     49real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! tendancies on the glaciers [1]
    5050real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! water ice at the surface [kg/m^2]
    5151real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! co2 ice at the surface [kg/m^2]
    … …  
    6666! -------------
    6767! Compute H2O adsorption, then CO2 adsorption
    68 call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
     68call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
    6969                            theta_h2o_adsorbded,m_h2o_completesoil,delta_mh2oreg)
    70 call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o, &
     70call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o, &
    7171                            tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mco2reg)
    7272
    … …  
    7474
    7575!=======================================================================
    76 SUBROUTINE regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
     76SUBROUTINE regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
    7777                                  theta_h2o_adsorbded,m_h2o_completesoil,delta_mreg)
    7878
    … …  
    100100! Inputs
    101101integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM,timelen   ! Size dimension
    102 real,    dimension(ngrid,nslope),           intent(in) :: tend_h2oglaciers, tend_co2glaciers ! Tendencies on the glaciers ()
     102real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! Tendencies on the glaciers ()
    103103real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
    104104real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! CO2 ice at the surface [kg/m^2]
    … …  
    151151    do ig = 1,ngrid
    152152        do islope = 1,nslope
    153             if (abs(tend_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
    154             if (abs(tend_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
     153            if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
     154            if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
    155155        enddo
    156156    enddo
    … …  
    224224!!!
    225225!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    226 SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg)
     226SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg)
    227227
    228228use comsoil_h_PEM,         only: layer_PEM, index_breccia, index_breccia
    … …  
    244244real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Mean Soil Temperature [K]
    245245real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Mean Thermal Inertia [USI]
    246 real,    dimension(ngrid,nslope),           intent(in) :: tend_h2oglaciers, tend_co2glaciers ! Tendencies on the glaciers ()
     246real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! Tendencies on the glaciers ()
    247247real,    dimension(ngrid,timelen),          intent(in) :: q_co2, q_h2o                       ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
    248248real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
    … …  
    292292    do ig = 1,ngrid
    293293        do islope = 1,nslope
    294             if (abs(tend_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
    295             if (abs(tend_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
     294            if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
     295            if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
    296296        enddo
    297297    enddo
    … …  
    317317
    318318    ! 1. Compute the fraction of the pores occupied by H2O
    319     call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
     319    call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
    320320                                theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o)
    321321

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

    r3495 r3498  
    472472    - The execution command line in the job scripts that should be modified by the user according to the set-up is now given as an argument at the beginning to be more identifiable and adaptable;
    473473    - Making the job scripts more robust to detect a successful end.
     474
     475== 07/11/2024 == JBC
     476- Correction of the variable name for the ice table depth in "pemetat0.F90". So it is now got as intended from the "startpem.nc" file;
     477- Renaming of the tendencies in the PEM with the prefix 'd_' instead of 'tend_';
     478- Modification of the PEM time step type from integer to real. As a consequence, all time variables are now of real type. This change adds the possibility to consider fractions of year as time step.

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

    r3367 r3498  
    77!=======================================================================
    88
    9 SUBROUTINE compute_tend(ngrid,nslope,min_ice,tendencies_ice)
     9SUBROUTINE compute_tend(ngrid,nslope,min_ice,d_ice)
    1010
    1111implicit none
    … …  
    2525
    2626!   OUTPUT
    27 real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! Difference between the minima = evolution of perennial ice
     27real, dimension(ngrid,nslope), intent(out) :: d_ice ! Difference between the minima = evolution of perennial ice
    2828!=======================================================================
    2929! We compute the difference
    30 tendencies_ice = min_ice(:,:,2) - min_ice(:,:,1)
     30d_ice = min_ice(:,:,2) - min_ice(:,:,1)
    3131
    3232! If the difference is too small, then there is no evolution
    33 where (abs(tendencies_ice) < 1.e-10) tendencies_ice = 0.
     33where (abs(d_ice) < 1.e-10) d_ice = 0.
    3434
    3535! If the minimum over the last year is 0, then we have no perennial ice
    36 where (abs(min_ice(:,:,2)) < 1.e-10) tendencies_ice = 0.
     36where (abs(min_ice(:,:,2)) < 1.e-10) d_ice = 0.
    3737
    3838END SUBROUTINE compute_tend

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

    r3490 r3498  
    2323#endif
    2424
    25 use time_evol_mod,         only: year_bp_ini, dt_pem, h2o_ice_crit, co2_ice_crit, ps_criterion, Max_iter_pem, &
    26                           evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years, ecritpem
     25use time_evol_mod,         only: year_bp_ini, dt, h2o_ice_crit, co2_ice_crit, ps_criterion, Max_iter_pem, &
     26                                 evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years, ecritpem
    2727use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_huge_h2oice, depth_breccia, depth_bedrock, reg_thprop_dependp
    2828use adsorption_mod,        only: adsorption_pem
    … …  
    3434implicit none
    3535
    36 integer, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
     36real, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
    3737
    3838character(20), parameter :: modname ='conf_pem'
    … …  
    6464year_bp_ini = 0.
    6565call getin('year_earth_bp_ini',year_earth_bp_ini)
    66 year_bp_ini = int(year_earth_bp_ini/convert_years)
     66year_bp_ini = year_earth_bp_ini/convert_years
    6767
    6868var_obl = .true.
    … …  
    9191call getin('ps_criterion',ps_criterion)
    9292
    93 dt_pem = 1
    94 call getin('dt_pem', dt_pem)
     93dt = 1.
     94call getin('dt',dt)
    9595
    9696!#---------- Subsurface parameters ----------#
    … …  
    128128    call abort_physic(modname,"Ice table must be used when soil_pem = T",1)
    129129endif
     130if (icetable_equilibrium .or. icetable_dynamic) then
     131    write(*,*) 'Ice table is asked to be computed both by the equilibrium and dynamic method.'
     132    write(*,*) 'The dynamic method is then chosen.'
     133endif
    130134
    131135if ((.not. soil_pem) .and. adsorption_pem) then
    … …  
    144148endif
    145149
    146 if (evol_orbit_pem .and. year_bp_ini == 0.) then
    147     write(*,*) 'You want to follow the file obl_ecc_lsp.asc for changing orb parameters,'
     150if (evol_orbit_pem .and. abs(year_bp_ini) < 1.e-10) then
     151    write(*,*) 'You want to follow the file "obl_ecc_lsp.asc" for changing orbital parameters,'
    148152    write(*,*) 'but you did not specify from which year to start.'
    149     call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini=0",1)
     153    call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini = 0",1)
    150154endif
    151155

  • trunk/LMDZ.COMMON/libf/evolution/deftank/README

    r3495 r3498  
    77      (ii)  nPCM_ini -> the number of initial PCM runs (at least 2);
    88      (iii) nPCM -> the number of PCM runs between each PEM run (usually 2);
    9       (iv)  mode -> the launching mode (0 = "processing scripts"; any other values = "submitting jobs"). The former option is usually used to process the script on a local machine while the latter is used to submit jobs on supercomputer.
     9      (iv)  mode -> the launching mode (0 = "processing scripts"; any other values = "submitting jobs"). The former option is usually used to process the script on a local machine while the latter is used to submit jobs on a supercomputer.
    1010  The script can take an argument:
    1111      - If there is no argument, then the script initiates a PEM simulation from scratch.

  • trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh

    r3495 r3498  
    1414###############################################
    1515# Set the number of years to be simulated, either Martian or Earth years:
    16 n_mars_years=100
    17 #n_earth_years=300
     16n_mars_years=100.
     17#n_earth_years=300.
    1818
    1919# Set the number of initial PCM runs (>= 2):
    … …  
    2323nPCM=2
    2424
    25 # Set the launching mode (0 = "processing scripts"; any other values = "submitting jobs"). The former option is usually used to process the script on a local machine while the latter is used to submit jobs on supercomputer:
     25# Set the launching mode (0 = "processing scripts"; any other values = "submitting jobs"). The former option is usually used to process the script on a local machine while the latter is used to submit jobs on a supercomputer:
    2626mode=1
    2727########################################################################

  • trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh

    r3495 r3498  
    5252
    5353    if [ -v n_mars_years ] && [ ! -z "$n_mars_years" ]; then
    54         if [ $n_mars_years -lt 1 ]; then
     54        if [ $n_mars_years -lt 0. ]; then
    5555            echo "Error: the value of 'n_mars_years' must be >0!"
    5656            errlaunch
    5757        fi
    5858    elif [ -v n_earth_years ] && [ ! -z "$n_earth_years" ]; then
    59         if [ $n_earth_years -lt 1 ]; then
     59        if [ $n_earth_years -lt 0. ]; then
    6060            echo "Error: the value of 'n_earth_years' must be >0!"
    6161            errlaunch
    … …  
    132132        echo "Number of years to be simulated: $n_myear Martian years."
    133133    elif [ -v n_earth_years ]; then
    134         n_myear=$(echo "($n_earth_years/$convert_years + 0.999999)/1" | bc) # Ceiling of n_earth_years/convert_years
     134        n_myear=$(echo "$n_earth_years/$convert_years" | bc -l)
    135135        echo "Number of years to be simulated: $n_earth_years Earth years = $n_myear Martian years."
    136136    fi
    … …  
    141141    echo "This is a chained simulation for PEM and PCM runs in $dir on $machine by $user."
    142142    convertyears
    143     i_myear=0
     143    i_myear=0.
    144144    iPEM=1
    145145    iPCM=1
    … …  
    176176                errlaunch
    177177            fi
    178         else # Mode: launching jobs
     178        else # Mode: submitting jobs
    179179            cp PCMrun.job PCMrun${iPCM}.job
    180180            sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPCM}/" PCMrun${iPCM}.job
    … …  
    187187        fi
    188188        ((iPCM++))
    189         ((i_myear++))
     189        i_myear = $(echo "$i_myear + 1." | bc - l)
    190190        ((ii++))
    191191    else
    … …  
    201201                    errlaunch
    202202                fi
    203             else # Mode: launching jobs
     203            else # Mode: submitting jobs
    204204                cp PCMrun.job PCMrun${iPCM}.job
    205205                sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPCM}/" PCMrun${iPCM}.job
    … …  
    209209            fi
    210210            ((iPCM++))
    211             ((i_myear++))
     211            i_myear = $(echo "$i_myear + 1." | bc - l)
    212212        else
    213213            endlaunch
    … …  
    226226                errlaunch
    227227            fi
    228         else # Mode: launching jobs
     228        else # Mode: submitting jobs
    229229            sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPEM}/" PEMrun.job
    230230            jobID=$(eval "$submit_job PEMrun.job")
    … …  
    255255                errlaunch
    256256            fi
    257         else # Mode: launching jobs
     257        else # Mode: submitting jobs
    258258            sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPEM}/" PEMrun.job
    259259            jobID=$(eval "$submit_dependjob=afterok:${jobID} PEMrun.job")

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

    r3339 r3498  
    3939# ps_criterion = 0.15
    4040
    41 # Time step length of the PEM in Martian years? Default = 1
    42 # dt_pem=1
     41# Time step length of the PEM in Martian years? Default = 1.
     42# dt=1
    4343
    4444#---------- Subsurface parameters ----------#

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

    r3368 r3498  
    77!=======================================================================
    88
    9 SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
     9SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
    1010
    11 use time_evol_mod, only: dt_pem
     11use time_evol_mod, only: dt
    1212
    1313implicit none
    … …  
    2525!   OUTPUT
    2626real, dimension(ngrid,nslope), intent(inout) :: co2_ice      ! Previous and actual density of CO2 ice
    27 real, dimension(ngrid,nslope), intent(inout) :: tend_co2_ice ! Evolution of perennial ice over one year
     27real, dimension(ngrid,nslope), intent(inout) :: d_co2ice ! Evolution of perennial ice over one year
    2828
    2929!   local:
    … …  
    3535
    3636co2_ice_old = co2_ice
    37 co2_ice = co2_ice + tend_co2_ice*dt_pem
     37co2_ice = co2_ice + d_co2ice*dt
    3838where (co2_ice < 0.)
    3939    co2_ice = 0.
    40     tend_co2_ice = -co2_ice_old/dt_pem
     40    d_co2ice = -co2_ice_old/dt
    4141end where
    4242
    … …  
    4545!=======================================================================
    4646
    47 SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
     47SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
    4848
    49 use time_evol_mod,          only: dt_pem
     49use time_evol_mod,          only: dt
    5050use comslope_mod,           only: subslope_dist, def_slope_mean
    5151
    … …  
    7373!   OUTPUT
    7474real, dimension(ngrid,nslope), intent(inout) :: h2o_ice      ! Previous and actual density of h2o ice (kg/m^2)
    75 real, dimension(ngrid,nslope), intent(inout) :: tend_h2o_ice ! Evolution of perennial ice over one year (kg/m^2/year)
     75real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Evolution of perennial ice over one year (kg/m^2/year)
    7676integer,                       intent(inout) :: stopPEM      ! Stopping criterion code
    7777
    … …  
    8080integer                       :: i, islope                                           ! Loop variables
    8181real                          :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve h2o
    82 real, dimension(ngrid,nslope) :: new_tendencies                                      ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
     82real, dimension(ngrid,nslope) :: new_tend                                      ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
    8383!=======================================================================
    8484if (ngrid /= 1) then ! Not in 1D
    … …  
    9999        do islope = 1,nslope
    100100            if (h2o_ice(i,islope) > 0.) then
    101                 if (tend_h2o_ice(i,islope) > 0.) then
    102                     pos_tend = pos_tend + tend_h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
     101                if (d_h2oice(i,islope) > 0.) then
     102                    pos_tend = pos_tend + d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
    103103                else
    104                     neg_tend = neg_tend - tend_h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
     104                    neg_tend = neg_tend - d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
    105105                endif
    106106            endif
    … …  
    108108    enddo
    109109
    110     new_tendencies = 0.
     110    new_tend = 0.
    111111    if (abs(pos_tend) < 1.e-10 .or. abs(neg_tend) < 1.e-10) then
    112112        write(*,*) "Reason of stopping: there is no sublimating or condensing h2o ice!"
    … …  
    118118        ! We adapt the tendencies to conserve h2o and do only exchange between grid points
    119119        if (neg_tend > pos_tend .and. pos_tend > 0.) then ! More sublimation on the planet than condensation
    120             where (tend_h2o_ice < 0.) ! We lower the sublimating rate by a coefficient
    121                 new_tendencies = tend_h2o_ice*pos_tend/neg_tend
     120            where (d_h2oice < 0.) ! We lower the sublimating rate by a coefficient
     121                new_tend = d_h2oice*pos_tend/neg_tend
    122122            elsewhere ! We don't change the condensing rate
    123                 new_tendencies = tend_h2o_ice
     123                new_tend = d_h2oice
    124124            end where
    125125        else if (neg_tend < pos_tend .and. neg_tend > 0.) then ! More condensation on the planet than sublimation
    126             where (tend_h2o_ice < 0.) ! We don't change the sublimating rate
    127                 new_tendencies = tend_h2o_ice
     126            where (d_h2oice < 0.) ! We don't change the sublimating rate
     127                new_tend = d_h2oice
    128128            elsewhere ! We lower the condensing rate by a coefficient
    129                 new_tendencies = tend_h2o_ice*neg_tend/pos_tend
     129                new_tend = d_h2oice*neg_tend/pos_tend
    130130            end where
    131131        endif
    … …  
    133133
    134134    ! Evolution of the h2o ice for each physical point
    135     h2o_ice = h2o_ice + new_tendencies*dt_pem
     135    h2o_ice = h2o_ice + new_tend*dt
    136136
    137137    ! We compute the amount of h2o that is sublimated in excess
    … …  
    142142                negative_part = negative_part - h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
    143143                h2o_ice(i,islope) = 0.
    144                 tend_h2o_ice(i,islope) = 0.
     144                d_h2oice(i,islope) = 0.
    145145            endif
    146146        enddo
    … …  
    157157    do islope = 1,nslope
    158158        do i = 1,ngrid
    159              if (new_tendencies(i,islope) > 0.) h2o_ice(i,islope) = h2o_ice(i,islope) - new_tendencies(i,islope)*real_coefficient*dt_pem*cos(def_slope_mean(islope)*pi/180.)
     159             if (new_tend(i,islope) > 0.) h2o_ice(i,islope) = h2o_ice(i,islope) - new_tend(i,islope)*real_coefficient*dt*cos(def_slope_mean(islope)*pi/180.)
    160160        enddo
    161161     enddo
    162162else ! ngrid == 1, i.e. in 1D
    163     h2o_ice = h2o_ice + tend_h2o_ice*dt_pem
     163    h2o_ice = h2o_ice + d_h2oice*dt
    164164endif
    165165

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

    r3493 r3498  
    1010!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1111
    12 logical, save                           :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium
    13 logical, save                           :: icetable_dynamic     ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method
    14 real, allocatable, dimension(:,:)       :: icetable_depth       ! ngrid x nslope: Depth of the ice table [m]
    15 real, allocatable, dimension(:,:)       :: icetable_thickness   ! ngrid x nslope: Thickness of the ice table [m]
    16 real, allocatable, dimension(:,:,:)     :: ice_porefilling      ! the amout of porefilling in each layer in each grid [m^3/m^3]
     12logical                             :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium
     13logical                             :: icetable_dynamic     ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method
     14real, allocatable, dimension(:,:)   :: icetable_depth       ! ngrid x nslope: Depth of the ice table [m]
     15real, allocatable, dimension(:,:)   :: icetable_thickness   ! ngrid x nslope: Thickness of the ice table [m]
     16real, allocatable, dimension(:,:,:) :: ice_porefilling      ! the amout of porefilling in each layer in each grid [m^3/m^3]
    1717
    1818!-----------------------------------------------------------------------

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

    r3387 r3498  
    99!=======================================================================
    1010
    11 SUBROUTINE info_PEM(year_iter,stopPEM,i_myear,n_myear)
     11SUBROUTINE info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
    1212
    1313!=======================================================================
    … …  
    2424
    2525!----- Arguments
    26 integer, intent(in) :: year_iter, stopPEM ! # of year and reason to stop
    27 integer, intent(in) :: i_myear, n_myear          ! Current simulated Martian year and maximum number of Martian years to be simulated
     26integer, intent(in) :: stopPEM          ! Reason to stop
     27real,    intent(in) :: i_myear_leg      ! # of years
     28real,    intent(in) :: i_myear, n_myear ! Current simulated Martian year and maximum number of Martian years to be simulated
    2829
    2930!----- Local variables
    … …  
    3536inquire(file = 'info_PEM.txt',exist = ok)
    3637if (ok) then
    37     write(ich1,'(i0)') i_myear
    38     write(ich2,'(i0)') n_myear
     38    write(ich1,'(f0.4)') i_myear
     39    write(ich2,'(f0.4)') n_myear
    3940    write(fch,'(f0.4)') convert_years ! 4 digits to the right of the decimal point to respect the precision of Martian year in "launch_pem.sh"
    4041    write(ich3,'(i0)') iPCM
    … …  
    5152    ! Martian date, Number of Martians years done by the PEM run, Number of Martians years done by the chainded simulation, Code of the stopping criterion
    5253    ! The conversion ratio from Planetary years to Earth years is given in the header of the file
    53     write(1,*) year_bp_ini + i_myear, year_iter, i_myear, stopPEM
     54    write(1,*) year_bp_ini + i_myear, i_myear_leg, i_myear, stopPEM
    5455    close(1)
    5556else

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

    r3430 r3498  
    1818
    1919! Physical parameters
    20 real, parameter :: tend_dust = 5.78e-2          ! Tendency of dust [kg.m-2.y-1]
     20real, parameter :: d_dust = 5.78e-2             ! Tendency of dust [kg.m-2.y-1]
    2121real, parameter :: dry_lag_porosity = 0.5       ! Porosity of dust lag
    2222real, parameter :: dry_regolith_porosity = 0.45 ! Porosity of regolith
    … …  
    452452!=======================================================================
    453453! Layering algorithm
    454 SUBROUTINE make_layering(this,tend_co2ice,tend_h2oice,tend_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)
     454SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)
    455455
    456456implicit none
    … …  
    461461logical,                intent(inout) :: new_str, new_lag1, new_lag2
    462462integer,                intent(inout) :: stopPEM
    463 real,     intent(in) :: tend_co2ice, tend_h2oice, tend_dust
    464 
    465 !---- Local variables
    466 real                   :: htend_co2ice, htend_h2oice, htend_dust
     463real,     intent(in) :: d_co2ice, d_h2oice, d_dust
     464
     465!---- Local variables
     466real                   :: dh_co2ice, dh_h2oice, dh_dust
    467467real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
    468468real                   :: thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
    … …  
    471471
    472472!---- Code
    473 htend_co2ice = tend_co2ice/rho_co2ice
    474 htend_h2oice = tend_h2oice/rho_h2oice
    475 htend_dust = tend_dust/rho_dust
    476 
    477 if (htend_dust < 0.) then ! Dust lifting only
    478     if (abs(htend_co2ice) > eps .or. abs(htend_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
     473dh_co2ice = d_co2ice/rho_co2ice
     474dh_h2oice = d_h2oice/rho_h2oice
     475dh_dust = d_dust/rho_dust
     476
     477if (dh_dust < 0.) then ! Dust lifting only
     478    if (abs(dh_co2ice) > eps .or. abs(dh_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
    479479    write(*,'(a)') ' Stratification -> Dust lifting'
    480     h2lift_tot = abs(htend_dust)
     480    h2lift_tot = abs(dh_dust)
    481481    do while (h2lift_tot > 0. .and. associated(current1))
    482482        ! Is the considered stratum a dust lag?
    … …  
    508508
    509509!------ Dust sedimentation only
    510     if (abs(htend_co2ice) < eps .and. abs(htend_h2oice) < eps) then
     510    if (abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then
    511511        write(*,'(a)') ' Stratification -> Dust sedimentation'
    512512        ! New stratum at the layering top by sedimentation of dust
    513         thickness = htend_dust/(1. - dry_regolith_porosity)
     513        thickness = dh_dust/(1. - dry_regolith_porosity)
    514514        if (thickness > 0.) then ! Only if the stratum is building up
    515515             if (new_str) then
    516                  call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,htend_dust/thickness,dry_regolith_porosity)
     516                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity)
    517517                 new_str = .false.
    518518             else
    … …  
    523523
    524524!------ Condensation of CO2 ice + H2O ice
    525     else if ((htend_co2ice >= 0. .and. htend_h2oice > 0.) .or. (htend_co2ice > 0. .and. htend_h2oice >= 0.)) then
     525    else if ((dh_co2ice >= 0. .and. dh_h2oice > 0.) .or. (dh_co2ice > 0. .and. dh_h2oice >= 0.)) then
    526526        write(*,'(a)') ' Stratification -> CO2 and H2O ice condensation'
    527527        ! New stratum at the layering top by condensation of CO2 and H2O ice
    528         thickness = htend_co2ice/(1. - co2ice_porosity) + htend_h2oice/(1. - h2oice_porosity) + htend_dust
     528        thickness = dh_co2ice/(1. - co2ice_porosity) + dh_h2oice/(1. - h2oice_porosity) + dh_dust
    529529        if (thickness > 0.) then ! Only if the stratum is building up
    530530             if (new_str) then
    531                  call add_stratum(this,thickness,this%top%top_elevation + thickness,htend_co2ice/thickness,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_co2ice/thickness + htend_h2oice/thickness + htend_dust/thickness))
     531                 call add_stratum(this,thickness,this%top%top_elevation + thickness,dh_co2ice/thickness,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_co2ice/thickness + dh_h2oice/thickness + dh_dust/thickness))
    532532                 new_str = .false.
    533533             else
    … …  
    538538
    539539!------ Sublimation of CO2 ice + Condensation of H2O ice
    540     else if (htend_co2ice < 0. .and. htend_h2oice > 0.) then
     540    else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then
    541541        write(*,'(a)') ' Stratification -> Sublimation of CO2 ice + Condensation of H2O ice'
    542542        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
    543         h2subl_tot = abs(htend_co2ice)
     543        h2subl_tot = abs(dh_co2ice)
    544544        do while (h2subl_tot > 0. .and. associated(current1))
    545545            h_co2ice_old = current1%co2ice_volfrac*current1%thickness
    … …  
    585585        endif
    586586        ! New stratum at the layering top by condensation of H2O ice
    587         thickness = htend_h2oice/(1. - h2oice_porosity) + htend_dust
     587        thickness = dh_h2oice/(1. - h2oice_porosity) + dh_dust
    588588        if (thickness > 0.) then ! Only if the stratum is building up
    589589             if (new_str) then
    590                  call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_h2oice/thickness + htend_dust/thickness))
     590                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness))
    591591                 new_str = .false.
    592592             else
    … …  
    597597
    598598!------ Sublimation of CO2 ice + H2O ice
    599     else if ((htend_co2ice <= 0. .and. htend_h2oice < 0.) .or. (htend_co2ice < 0. .and. htend_h2oice <= 0.)) then
     599    else if ((dh_co2ice <= 0. .and. dh_h2oice < 0.) .or. (dh_co2ice < 0. .and. dh_h2oice <= 0.)) then
    600600        write(*,'(a)') ' Stratification -> Sublimation of CO2 and H2O ice'
    601601        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
    602         h2subl_tot = abs(htend_co2ice)
     602        h2subl_tot = abs(dh_co2ice)
    603603        do while (h2subl_tot > 0. .and. associated(current1))
    604604            h_co2ice_old = current1%co2ice_volfrac*current1%thickness
    … …  
    644644        endif
    645645        ! H2O ice sublimation in the considered stratum + New stratum for dust lag
    646         h2subl_tot = abs(htend_h2oice)
     646        h2subl_tot = abs(dh_h2oice)
    647647        do while (h2subl_tot > 0. .and. associated(current2))
    648648            h_co2ice_old = current2%co2ice_volfrac*current2%thickness
    … …  
    689689
    690690!------ Condensation of CO2 ice + Sublimation of H2O ice
    691     else if (htend_co2ice > 0. .and. htend_h2oice < 0.) then
     691    else if (dh_co2ice > 0. .and. dh_h2oice < 0.) then
    692692        error stop 'Impossible situation: CO2 ice condensation + H2O ice sublimation!'
    693693    endif

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

    r3399 r3498  
    1515!=======================================================================
    1616
    17 SUBROUTINE orbit_param_criterion(i_myear,year_iter_max)
     17SUBROUTINE orbit_param_criterion(i_myear,n_myear_leg)
    1818
    1919#ifdef CPP_IOIPSL
    … …  
    3838! Input Variables
    3939!--------------------------------------------------------
    40 integer, intent(in) :: i_myear ! Martian year of the simulation
     40real, intent(in) :: i_myear ! Martian year of the simulation
    4141
    4242!--------------------------------------------------------
    4343! Output Variables
    4444!--------------------------------------------------------
    45 integer, intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much
     45real, intent(out) :: n_myear_leg ! Maximum number of iteration of the PEM before orb changes too much
    4646
    4747!--------------------------------------------------------
    4848! Local variables
    4949!--------------------------------------------------------
    50 integer :: Year                                           ! Year of the simulation
     50real    :: Year                                           ! Year of the simulation
    5151real    :: Lsp                                            ! Ls perihelion
    5252real    :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
    5353real    :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
    5454real    :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
    55 integer :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
     55real    :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
    5656integer :: ilask                                          ! Loop variable
    5757real    :: xa, xb, ya, yb                                 ! Variables for interpolation
    … …  
    113113! If we do not want some orb parameter to change, they should not be a stopping criterion,
    114114! So the number of iteration corresponding is set to maximum
    115 max_obl_iter = 1000000
    116 max_ecc_iter = 1000000
    117 max_lsp_iter = 1000000
     115max_obl_iter = 1000000.
     116max_ecc_iter = 1000000.
     117max_lsp_iter = 1000000.
    118118
    119119! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters
    … …  
    136136        yb = obllask(ilask - 1)
    137137        if (yb < min_obl) then ! The minimum accepted is overtaken
    138             max_obl_iter = floor((min_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
     138            max_obl_iter = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
    139139            found_obl = .true.
    140140        else if (max_obl < yb) then ! The maximum accepted is overtaken
    141             max_obl_iter = floor((max_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
     141            max_obl_iter = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
    142142            found_obl = .true.
    143143        endif
    … …  
    149149        yb = ecclask(ilask - 1)
    150150        if (yb < min_ecc) then ! The minimum accepted is overtaken
    151             max_ecc_iter = floor((min_ecc - ya)*(xb - xa)/(yb - ya) + xa) - Year
     151            max_ecc_iter = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
    152152            found_ecc = .true.
    153153        else if (max_ecc < yb) then ! The maximum accepted is overtaken
    154             max_ecc_iter = floor((max_ecc - ya)*(xb - xa)/(yb - ya) + xa) - Year
     154            max_ecc_iter = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
    155155            found_ecc = .true.
    156156        endif
    … …  
    187187        endif
    188188        if (yb < min_lsp) then ! The minimum accepted is overtaken
    189             max_lsp_iter = floor((min_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
     189            max_lsp_iter = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
    190190            found_lsp = .true.
    191191        else if (max_lsp < yb) then ! The maximum accepted is overtaken
    192             max_lsp_iter = floor((max_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
     192            max_lsp_iter = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
    193193            found_lsp = .true.
    194194        endif
    … …  
    206206write(*,*) 'Max. number of iterations for the Lsp parameter  =', max_lsp_iter
    207207
    208 year_iter_max = min(max_obl_iter,max_ecc_iter,max_lsp_iter)
    209 year_iter_max = max(year_iter_max,1)
    210 write(*,*) 'So the max. number of iterations for the orbital criteria =', year_iter_max
     208n_myear_leg = min(max_obl_iter,max_ecc_iter,max_lsp_iter)
     209if (n_myear_leg < 1.) n_myear_leg = 1. ! n_myear_leg should be at least equal to 1
     210write(*,*) 'So the max. number of iterations for the orbital criteria =', n_myear_leg
    211211
    212212END SUBROUTINE orbit_param_criterion

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

    r3493 r3498  
    5151                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
    5252                                      co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
    53 use time_evol_mod,              only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
     53use time_evol_mod,              only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
    5454use orbit_param_criterion_mod,  only: orbit_param_criterion
    5555use recomp_orb_param_mod,       only: recomp_orb_param
    … …  
    7070use writediagpem_mod,           only: writediagpem, writediagsoilpem
    7171use co2condens_mod,             only: CO2cond_ps
    72 use layering_mod,               only: tend_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
     72use layering_mod,               only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
    7373
    7474#ifndef CPP_STD
    … …  
    195195
    196196! Variables for slopes
    197 real, dimension(:,:),   allocatable :: tend_co2_ice       ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
    198 real, dimension(:,:),   allocatable :: tend_co2_ice_ini   ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
    199 real, dimension(:,:),   allocatable :: tend_h2o_ice       ! physical point x slope field: Tendency of evolution of perennial h2o ice
    200 real, dimension(:,:),   allocatable :: flag_co2flow       ! (ngrid,nslope): Flag where there is a CO2 glacier flow
    201 real, dimension(:),     allocatable :: flag_co2flow_mesh  ! (ngrid)       : Flag where there is a CO2 glacier flow
    202 real, dimension(:,:),   allocatable :: flag_h2oflow       ! (ngrid,nslope): Flag where there is a H2O glacier flow
    203 real, dimension(:),     allocatable :: flag_h2oflow_mesh  ! (ngrid)       : Flag where there is a H2O glacier flow
     197real, dimension(:,:),   allocatable :: d_co2ice          ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
     198real, dimension(:,:),   allocatable :: d_co2ice_ini      ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
     199real, dimension(:,:),   allocatable :: d_h2oice          ! physical point x slope field: Tendency of evolution of perennial h2o ice
     200real, dimension(:,:),   allocatable :: flag_co2flow      ! (ngrid,nslope): Flag where there is a CO2 glacier flow
     201real, dimension(:),     allocatable :: flag_co2flow_mesh ! (ngrid)       : Flag where there is a CO2 glacier flow
     202real, dimension(:,:),   allocatable :: flag_h2oflow      ! (ngrid,nslope): Flag where there is a H2O glacier flow
     203real, dimension(:),     allocatable :: flag_h2oflow_mesh ! (ngrid)       : Flag where there is a H2O glacier flow
    204204
    205205! Variables for surface and soil
    … …  
    230230! Some variables for the PEM run
    231231real, parameter :: year_step = 1   ! Timestep for the pem
    232 integer         :: year_iter       ! Number of iteration
    233 integer         :: year_iter_max   ! Maximum number of iterations before stopping
    234 integer         :: i_myear         ! Global number of Martian years of the chained simulations
    235 integer         :: n_myear         ! Maximum number of Martian years of the chained simulations
     232real            :: i_myear_leg       ! Number of iteration
     233real            :: n_myear_leg     ! Maximum number of iterations before stopping
     234real            :: i_myear         ! Global number of Martian years of the chained simulations
     235real            :: n_myear         ! Maximum number of Martian years of the chained simulations
    236236real            :: timestep        ! Timestep [s]
    237237character(20)   :: job_id          ! Job id provided as argument passed on the command line when the program was invoked
    … …  
    604604!    I_f Compute tendencies
    605605!------------------------
    606 allocate(tend_co2_ice(ngrid,nslope),tend_h2o_ice(ngrid,nslope))
    607 allocate(tend_co2_ice_ini(ngrid,nslope))
     606allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope))
     607allocate(d_co2ice_ini(ngrid,nslope))
    608608
    609609! Compute the tendencies of the evolution of ice over the years
    610 call compute_tend(ngrid,nslope,min_co2_ice,tend_co2_ice)
    611 call compute_tend(ngrid,nslope,min_h2o_ice,tend_h2o_ice)
    612 tend_co2_ice_ini = tend_co2_ice
     610call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
     611call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
     612d_co2ice_ini = d_co2ice
    613613deallocate(min_co2_ice,min_h2o_ice)
    614614
    … …  
    646646call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
    647647              icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
    648               ps_timeseries,tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,               &
    649               global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,          &
     648              ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice,                       &
     649              global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,         &
    650650              delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
    651651
    … …  
    662662    do islope = 1,nslope
    663663        if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
    664         if (tend_co2_ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
     664        if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
    665665            ini_co2ice_sublim(i,islope) = .true.
    666666            co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
    667667        endif
    668         if (tend_h2o_ice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
     668        if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
    669669            ini_h2oice_sublim(i,islope) = .true.
    670670            h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
    … …  
    712712#endif
    713713
    714 year_iter_max = Max_iter_pem
    715 if (evol_orbit_pem) call orbit_param_criterion(i_myear,year_iter_max)
     714n_myear_leg = Max_iter_pem
     715if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
    716716
    717717!-------------------------- END INITIALIZATION -------------------------
    … …  
    722722!    II_a Update pressure, ice and tracers
    723723!------------------------
    724 year_iter = 0
     724i_myear_leg = 0
    725725stopPEM = 0
    726726if (layering_algo) then
    … …  
    737737endif
    738738
    739 do while (year_iter < year_iter_max .and. i_myear < n_myear)
     739do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
    740740! II.a.1. Compute updated global pressure
    741741    write(*,*) "Recomputing the new pressure..."
    742742    do i = 1,ngrid
    743743        do islope = 1,nslope
    744             global_avg_press_new = global_avg_press_new - CO2cond_ps*g*cell_area(i)*tend_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
     744            global_avg_press_new = global_avg_press_new - CO2cond_ps*g*cell_area(i)*d_co2ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
    745745        enddo
    746746    enddo
    … …  
    817817!    II_b Evolution of ice
    818818!------------------------
    819     call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
    820     call evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
     819    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
     820    call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
    821821    if (layering_algo) then
    822822        do islope = 1,nslope
    823823            do ig = 1,ngrid
    824                 call make_layering(stratif(ig,islope),tend_co2_ice(ig,islope),tend_h2o_ice(ig,islope),tend_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
     824                call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
    825825            enddo
    826826        enddo
    … …  
    832832!------------------------
    833833    if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
    834                                             global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
     834                                                            global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
    835835    if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
    836836
    … …  
    878878                    emis(ig,islope) = emissiv
    879879                endif
    880             else if ((co2_ice(ig,islope) > 1.e-3) .and. (tend_co2_ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
     880            else if ((co2_ice(ig,islope) > 1.e-3) .and. (d_co2ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
    881881                ave = 0.
    882882                do t = 1,timelen
    … …  
    945945
    946946! II_d.4 Update the soil thermal properties
    947         call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM)
     947        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM)
    948948
    949949! II_d.5 Update the mass of the regolith adsorbed
    950950        if (adsorption_pem) then
    951             call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice,                   &
     951            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,                   &
    952952                                     h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
    953953                                     h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
    … …  
    986986        call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
    987987        call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
    988         call writediagpem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope))
    989         call writediagpem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope))
     988        call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
     989        call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
    990990        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
    991991        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
    … …  
    10121012!    II_f Update the tendencies
    10131013!------------------------
    1014     call recomp_tend_co2_slope(ngrid,nslope,timelen,tend_co2_ice,tend_co2_ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries, &
     1014    call recomp_tend_co2_slope(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries, &
    10151015                               global_avg_press_PCM,global_avg_press_new)
    10161016
    … …  
    10231023    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
    10241024    call stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
    1025     year_iter = year_iter + dt_pem
    1026     i_myear = i_myear + dt_pem
    1027     if (stopPEM <= 0 .and. year_iter >= year_iter_max) stopPEM = 5
     1025    i_myear_leg = i_myear_leg + dt
     1026    i_myear = i_myear + dt
     1027    if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
    10281028    if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
    10291029    call system_clock(c2)
    … …  
    10531053    else
    10541054        write(*,*) "The PEM can continue!"
    1055         write(*,*) "Number of iterations of the PEM: year_iter =", year_iter
     1055        write(*,*) "Number of iterations of the PEM: i_myear_leg =", i_myear_leg
    10561056        write(*,*) "Number of simulated Martian years: i_myear =", i_myear
    10571057    endif
    … …  
    11511151enddo
    11521152
    1153 if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
     1153if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
    11541154
    11551155!------------------------
    … …  
    12131213write(*,*) "restartpem.nc has been written"
    12141214
    1215 call info_PEM(year_iter,stopPEM,i_myear,n_myear)
    1216 
    1217 write(*,*) "The PEM has run for", year_iter, "Martian years."
     1215call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
     1216
     1217write(*,*) "The PEM has run for", i_myear_leg, "Martian years."
    12181218write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
    12191219write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
    1220 write(*,*) "LL & RV & JBC: so far, so good!"
     1220write(*,*) "PEM: so far, so good!"
    12211221
    12221222if (layering_algo) then

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

    r3493 r3498  
    88
    99SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness,      &
    10                     ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, &
     10                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice, &
    1111                    global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                         &
    1212                    m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
    … …  
    4949real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
    5050real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
    51 real, dimension(ngrid,nslope),  intent(in) :: tend_h2o_ice        ! tendencies on h2o ice
    52 real, dimension(ngrid,nslope),  intent(in) :: tend_co2_ice        ! tendencies on co2 ice
     51real, dimension(ngrid,nslope),  intent(in) :: d_h2oice        ! tendencies on h2o ice
     52real, dimension(ngrid,nslope),  intent(in) :: d_co2ice        ! tendencies on co2 ice
    5353real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg       ! surface water ice density, yearly averaged (kg/m^3)
    5454! outputs
    … …  
    105105inquire(file = filename,exist = startpem_file)
    106106
    107 write(*,*)'Is start PEM?',startpem_file
     107write(*,*)'Is there a "startpem" file?',startpem_file
    108108
    109109!0.2 Set to default values
    … …  
    315315                write(*,*)'will reconstruct the values of the ice table given the current state'
    316316                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
    317                 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     317                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
    318318                do islope = 1,nslope
    319319                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    … …  
    327327                write(*,*)'will reconstruct the values of the ice table given the current state'
    328328                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
    329                 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     329                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
    330330                do islope = 1,nslope
    331331                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    … …  
    360360            enddo
    361361
    362             call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
     362            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
    363363                                        m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
    364364
    … …  
    504504        if (icetable_equilibrium) then
    505505            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
    506             call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     506            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
    507507            do islope = 1,nslope
    508508                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    … …  
    512512            ice_porefilling = 0.
    513513            ice_table_depth = 0.
    514             call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     514            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
    515515            do islope = 1,nslope
    516516                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    … …  
    524524            m_co2_regolith_phys = 0.
    525525            m_h2o_regolith_phys = 0.
    526             call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
     526            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
    527527                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
    528528            deltam_co2_regolith_phys = 0.

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

    r3493 r3498  
    7373
    7474character(*),                            intent(in) :: filename
    75 integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope, i_myear
     75integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope
     76real,                                    intent(in) :: i_myear
    7677real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_slope_PEM       ! under mesh bining according to slope
    7778real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: inertiesoil_slope_PEM ! under mesh bining according to slope

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

    r3206 r3498  
    1313!=======================================================================
    1414
    15 SUBROUTINE recomp_orb_param(i_myear,year_iter)
     15SUBROUTINE recomp_orb_param(i_myear,i_myear_leg)
    1616
    1717use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp
    … …  
    3030! Input Variables
    3131!--------------------------------------------------------
    32 integer, intent(in) :: i_myear   ! Number of simulated Martian years
    33 integer, intent(in) :: year_iter ! Number of iterations of the PEM
     32real, intent(in) :: i_myear   ! Number of simulated Martian years
     33real, intent(in) :: i_myear_leg ! Number of iterations of the PEM
    3434
    3535!--------------------------------------------------------
    … …  
    4040! Local variables
    4141!--------------------------------------------------------
    42 integer :: Year           ! Year of the simulation
     42real    :: Year           ! Year of the simulation
    4343real    :: Lsp            ! Ls perihelion
    4444integer :: ilask          ! Loop variable
    … …  
    6161Lsp = lsperi*180./pi         ! We convert in degree
    6262
    63 write(*,*) 'recomp_orb_param, Old year =', Year - year_iter
     63write(*,*) 'recomp_orb_param, Old year =', Year - i_myear_leg
    6464write(*,*) 'recomp_orb_param, Old obl. =', obliquit
    6565write(*,*) 'recomp_orb_param, Old ecc. =', e_elips

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

    r3381 r3498  
    77!=======================================================================
    88
    9 SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope,emissivity_slope, &
     9SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_phys_ini,co2ice_slope,emissivity_slope, &
    1010                                 vmr_co2_PCM,vmr_co2_pem,ps_PCM_2,global_avg_press_PCM,global_avg_press_new)
    1111
    … …  
    2929real,                           intent(in) :: global_avg_press_PCM        ! global averaged pressure at previous timestep
    3030real,                           intent(in) :: global_avg_press_new        ! global averaged pressure at current timestep
    31 real, dimension(ngrid,nslope),  intent(in) :: tendencies_co2_ice_phys_ini ! physical point field: Evolution of perennial ice over one year
     31real, dimension(ngrid,nslope),  intent(in) :: d_co2ice_phys_ini ! physical point field: Evolution of perennial ice over one year
    3232real, dimension(ngrid,nslope),  intent(in) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope (kg/m^2)
    3333real, dimension(ngrid,nslope),  intent(in) :: emissivity_slope            ! Emissivity per mesh and sub-grid slope(1)
    3434!   OUTPUT
    35 real, dimension(ngrid,nslope), intent(inout) ::  tendencies_co2_ice_phys ! physical point field: Evolution of perennial ice over one year
     35real, dimension(ngrid,nslope), intent(inout) ::  d_co2ice_phys ! physical point field: Evolution of perennial ice over one year
    3636
    3737!   local:
    … …  
    4747        coef = sols_per_my*sec_per_sol*emissivity_slope(i,islope)*sigmaB/Lco2
    4848        ave = 0.
    49         if (co2ice_slope(i,islope) > 1.e-4 .and. abs(tendencies_co2_ice_phys(i,islope)) > 1.e-5) then
     49        if (co2ice_slope(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
    5050            do t=1,timelen
    5151                ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM_2(i,t)/100.)))**4 &
    … …  
    5353            enddo
    5454            if (ave < 1e-4) ave = 0.
    55             tendencies_co2_ice_phys(i,islope) = tendencies_co2_ice_phys_ini(i,islope) - coef*ave/timelen
     55            d_co2ice_phys(i,islope) = d_co2ice_phys_ini(i,islope) - coef*ave/timelen
    5656        endif
    5757    enddo

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

    r3214 r3498  
    33implicit none
    44
    5 integer :: year_bp_ini    ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
    6 integer :: dt_pem         ! Time step used by the PEM in Planetary years
     5real    :: year_bp_ini    ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
     6real    :: dt             ! Time step used by the PEM in Planetary years
    77real    :: convert_years  ! Conversion ratio from Planetary years to Earth years
    88real    :: h2o_ice_crit   ! Percentage of change of the surface of h2o ice sublimating before stopping the PEM
    99real    :: co2_ice_crit   ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
    1010real    :: ps_criterion   ! Percentage of change of averaged surface pressure before stopping the PEM
    11 integer :: Max_iter_pem   ! Maximal number of iteration when converging to a steady state, read in evol.def
     11integer :: Max_iter_pem   ! Maximal number of iterations when converging to a steady state, read in evol.def
    1212logical :: evol_orbit_pem ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
    1313logical :: var_obl        ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity

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

    r3345 r3498  
    5151use mod_phys_lmdz_para, only: is_parallel, is_mpi_root, is_master, gather
    5252use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat, nbp_lev, grid_type, unstructured
    53 use time_evol_mod,      only: ecritpem, dt_pem
     53use time_evol_mod,      only: ecritpem, dt
    5454
    5555implicit none
    … …  
    255255! to writediagpem
    256256!------------------------------------------------------------------------
    257 if (nom == firstnom) zitau = zitau + dt_pem
     257if (nom == firstnom) zitau = zitau + 1
    258258
    259259!--------------------------------------------------------
    … …  
    600600use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat
    601601use mod_grid_phy_lmdz,  only: grid_type, unstructured
    602 use time_evol_mod,      only: ecritpem, dt_pem
     602use time_evol_mod,      only: ecritpem, dt
    603603use iniwritesoil_mod,   only: iniwritesoil
    604604
    … …  
    747747if (name.eq.firstname) then
    748748  ! if we run across 'firstname', then it is a new time step
    749   zitau = zitau + dt_pem
     749  zitau = zitau + dt
    750750endif
    751751
Note: See TracChangeset for help on using the changeset viewer.