source: trunk/LMDZ.COMMON/libf/evolution/pem.F90 @ 3595

Last change on this file since 3595 was 3591, checked in by jbclement, 12 days ago

PEM:

  • Making allocation/deallocation systematically and more efficient in the main program.
  • Some cleanings (variables deletion, more adapted type/dimension, etc).

JBC
File size: 62.7 KB
Line 
1!------------------------
2! I   Initialization
3!    I_a Read the "run.def"
4!    I_b Read the "start.nc" and "startfi.nc"
5!    I_c Subslope parametrisation
6!    I_d Read the PCM data and convert them to the physical grid
7!    I_e Initialization of the PEM variable and soil
8!    I_f Compute tendencies
9!    I_g Compute global surface pressure
10!    I_h Read the "startpem.nc"
11!    I_i Compute orbit criterion
12
13! II  Run
14!    II_a Update pressure, ice and tracers
15!    II_b Evolution of ice
16!    II_c Flow of glaciers
17!    II_d Update surface and soil temperatures
18!    II_e Outputs
19!    II_f Update the tendencies
20!    II_g Checking the stopping criterion
21
22! III Output
23!    III_a Update surface value for the PCM start files
24!    III_b Write the "restart.nc" and "restartfi.nc"
25!    III_c Write the "restartpem.nc"
26!------------------------
27
28PROGRAM pem
29
30use phyetat0_mod,               only: phyetat0
31use phyredem,                   only: physdem0, physdem1
32use netcdf,                     only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
33use turb_mod,                   only: q2, wstar
34use comslope_mod,               only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h
35use logic_mod,                  only: iflag_phys
36use mod_const_mpi,              only: COMM_LMDZ
37use infotrac
38use geometry_mod,               only: latitude_deg
39use conf_pem_mod,               only: conf_pem
40use pemredem,                   only: pemdem0, pemdem1
41use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
42                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
43use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2
44use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
45use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice
46use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
47                                      TI_PEM,               & ! Soil thermal inertia
48                                      tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
49                                      fluxgeo                 ! Geothermal flux for the PEM and PCM
50use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
51                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
52                                      co2_adsorbed_phys, h2o_adsorbed_phys          ! Mass of co2 and h2O adsorbed
53use time_evol_mod,              only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
54use orbit_param_criterion_mod,  only: orbit_param_criterion
55use recomp_orb_param_mod,       only: recomp_orb_param
56use ice_table_mod,              only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, &
57                                      ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi
58use soil_thermalproperties_mod, only: update_soil_thermalproperties
59use time_phylmdz_mod,           only: daysec, dtphys
60use abort_pem_mod,              only: abort_pem
61use soil_settings_PEM_mod,      only: soil_settings_PEM
62use compute_tend_mod,           only: compute_tend
63use info_PEM_mod,               only: info_PEM
64use get_timelen_PCM_mod,        only: get_timelen_PCM
65use pemetat0_mod,               only: pemetat0
66use read_data_PCM_mod,          only: read_data_PCM
67use recomp_tend_co2_mod,        only: recomp_tend_co2
68use compute_soiltemp_mod,       only: compute_tsoil_pem, shift_tsoil2surf
69use writediagpem_mod,           only: writediagpem, writediagsoilpem
70use co2condens_mod,             only: CO2cond_ps
71use layering_mod,               only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
72use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
73use version_info_mod,           only: print_version_info
74
75#ifndef CPP_STD
76    use comsoil_h,          only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil
77    use surfdat_h,          only: tsurf, qsurf, emis, emissiv, emisice, ini_surfdat_h,   &
78                                  albedodat, albedice, albedo_h2o_frost, albedo_h2o_cap, &
79                                  zmea, zstd, zsig, zgam, zthe, frost_albedo_threshold,  &
80                                  watercap, watercaptag, perennial_co2ice, albedo_perennialco2
81    use dimradmars_mod,     only: totcloudfrac, albedo
82    use dust_param_mod,     only: tauscaling
83    use tracer_mod,         only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
84    use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
85    use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
86    use comcstfi_h,         only: mugaz
87    use surfini_mod,        only: surfini
88    use comconst_mod,       only: pi, rad, g, r, cpp, kappa
89#else
90    use tracer_h,           only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
91    use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT,   &
92                                  PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
93                                  ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
94    use aerosol_mod,        only: iniaerosol
95    use planete_mod,        only: apoastr, periastr, year_day, peri_day, obliquit
96    use comcstfi_mod,       only: pi, rad, g, mugaz, r
97#endif
98
99#ifndef CPP_1D
100    use iniphysiq_mod,            only: iniphysiq
101    use control_mod,              only: iphysiq, day_step, nsplit_phys
102#else
103    use time_phylmdz_mod,         only: iphysiq, steps_per_sol
104    use regular_lonlat_mod,       only: init_regular_lonlat
105    use physics_distribution_mod, only: init_physics_distribution
106    use mod_grid_phy_lmdz,        only: regular_lonlat
107    use init_testphys1d_mod,      only: init_testphys1d
108    use comvert_mod,              only: ap, bp
109    use writerestart1D_mod,       only: writerestart1D
110#endif
111
112implicit none
113
114include "dimensions.h"
115include "paramet.h"
116include "comgeom.h"
117include "iniprint.h"
118include "callkeys.h"
119
120integer ngridmx
121parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
122
123! Same variable names as in the PCM
124integer, parameter :: nlayer = llm ! Number of vertical layer
125integer            :: ngrid        ! Number of physical grid points
126integer            :: nq           ! Number of tracer
127integer            :: day_ini      ! First day of the simulation
128real               :: pday         ! Physical day
129real               :: time_phys    ! Same as in PCM
130real               :: ptimestep    ! Same as in PCM
131real               :: ztime_fin    ! Same as in PCM
132
133! Variables to read "start.nc"
134character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM
135
136! Dynamic variables
137real, dimension(ip1jm,llm)          :: vcov          ! vents covariants
138real, dimension(ip1jmp1,llm)        :: ucov          ! vents covariants
139real, dimension(ip1jmp1,llm)        :: teta          ! Potential temperature
140real, dimension(:,:,:), allocatable :: q             ! champs advectes
141real, dimension(ip1jmp1)            :: ps_start_dyn  ! surface pressure in the start file (dynamic grid)
142real, dimension(:),     allocatable :: ps_start      ! surface pressure in the start file
143real, dimension(:),     allocatable :: ps_start0     ! surface pressure in the start file at the beginning
144real, dimension(:),     allocatable :: ps_avg        ! (ngrid) Averaged surface pressure
145real, dimension(:),     allocatable :: ps_dev        ! (ngrid x timelen) Surface pressure deviation
146real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) Instantaneous surface pressure
147real, dimension(ip1jmp1,llm)        :: masse         ! Air mass
148real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
149real                                :: time_0
150
151! Variables to read starfi.nc
152character(*), parameter :: startfi_name = "startfi.nc" ! Name of the file used to initialize the PEM
153character(2)            :: str2
154integer                 :: ncid, status                           ! Variable for handling opening of files
155integer                 :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files
156integer                 :: apvarid, bpvarid                       ! Variable ID for Netcdf files
157
158! Variables to read starfi.nc and write restartfi.nc
159real, dimension(:), allocatable :: longitude     ! Longitude read in startfi_name and written in restartfi
160real, dimension(:), allocatable :: latitude      ! Latitude read in startfi_name and written in restartfi
161real, dimension(:), allocatable :: cell_area     ! Cell_area read in startfi_name and written in restartfi
162real                            :: total_surface ! Total surface of the planet
163
164! Variables for h2o_ice evolution
165real, dimension(:,:),    allocatable  :: h2o_ice              ! h2o ice in the PEM
166real, dimension(:,:),    allocatable  :: d_h2oice             ! physical point x slope field: Tendency of evolution of perennial h2o ice
167real, dimension(:,:,:),  allocatable  :: min_h2o_ice          ! Minima of h2o ice at each point for the PCM years [kg/m^2]
168real                                  :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
169logical, dimension(:,:), allocatable  :: is_h2oice_sublim_ini ! Logical array to know if h2o ice is sublimating
170real                                  :: ps_avg_global_ini    ! constant: Global average pressure at initialization [Pa]
171real                                  :: ps_avg_global_old    ! constant: Global average pressure of previous time step
172real                                  :: ps_avg_global_new    ! constant: Global average pressure of current time step
173real,   dimension(:,:),   allocatable :: zplev_new            ! Grid points x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
174real,   dimension(:,:),   allocatable :: zplev_start0         ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2]
175real,   dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
176real,   dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
177integer                               :: stopPEM              ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu
178real                                  :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
179real,   dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Grid points x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
180integer                               :: timelen              ! # time samples
181real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
182
183! Variables for co2_ice evolution
184real,    dimension(:,:), allocatable :: co2_ice                ! co2 ice in the PEM
185real,    dimension(:,:), allocatable :: d_co2ice               ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
186real,    dimension(:,:), allocatable :: d_co2ice_ini           ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
187logical, dimension(:,:), allocatable :: is_co2ice_ini          ! Was there co2 ice initially in the PEM?
188real,  dimension(:,:,:), allocatable :: min_co2_ice            ! Minimum of co2 ice at each point for the first year [kg/m^2]
189real                                 :: co2ice_sublim_surf_ini ! Initial surface of sublimating co2 ice
190logical, dimension(:,:), allocatable :: is_co2ice_sublim_ini   ! Logical array to know if co2 ice is sublimating
191real,    dimension(:,:), allocatable :: vmr_co2_PCM            ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
192real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys       ! Grid points x Times co2 volume mixing ratio used in the PEM
193real,    dimension(:,:), allocatable :: q_co2_PEM_phys         ! Grid points x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
194
195! Variables for stratification (layering) evolution
196type(layering), dimension(:,:), allocatable :: stratif                     ! Layering (linked list of "stratum") for each grid point and slope
197type(ptrarray), dimension(:,:), allocatable :: current1, current2          ! Current active stratum in the layering
198logical,        dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm
199
200! Variables for slopes
201integer(kind = 1), dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow
202integer(kind = 1), dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow
203
204! Variables for surface and soil
205real, dimension(:,:),     allocatable :: tsurf_avg                        ! Grid points x Slope field: Averaged surface temperature [K]
206real, dimension(:,:),     allocatable :: tsurf_dev                        ! ngrid x Slope x Times field: Surface temperature deviation [K]
207real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Grid points x Slope field: Averaged surface temperature of first call of the PCM [K]
208real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Grid points x Soil x Slope field: Averaged Soil Temperature [K]
209real, dimension(:,:),     allocatable :: tsoil_avg_old                    ! Grid points x Soil field: Averaged Soil Temperature at the previous time step [K]
210real, dimension(:,:,:),   allocatable :: tsoil_dev                        ! Grid points x Soil x Slope field: Soil temperature deviation [K]
211real, dimension(:,:,:,:), allocatable :: tsoil_timeseries                 ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
212real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries             ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K]
213real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
214real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Grid points x Slope: Averaged  water surface density [kg/m^3]
215real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
216real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Grid points x Soil x Slopes: Averaged water soil density [kg/m^3]
217real, dimension(:),       allocatable :: delta_co2_adsorbed               ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
218real, dimension(:),       allocatable :: delta_h2o_adsorbed               ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
219real                                  :: totmassco2_adsorbed              ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
220real                                  :: totmassh2o_adsorbed              ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
221logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
222real, dimension(:,:),     allocatable :: icetable_thickness_old           ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
223real, dimension(:,:,:),   allocatable :: ice_porefilling_old              ! ngrid x nslope: Ice pore filling at the previous iteration [m]
224real, dimension(:),       allocatable :: delta_h2o_icetablesublim         ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
225real, dimension(:),       allocatable :: porefill                         ! Pore filling (output) to compute the dynamic ice table
226real                                  :: ssi_depth                        ! Ice table depth (output) to compute the dynamic ice table
227real, dimension(:,:),     allocatable :: zshift_surf                      ! Elevation shift for the surface [m]
228real, dimension(:,:),     allocatable :: zlag                             ! Newly built lag thickness [m]
229
230! Some variables for the PEM run
231real, parameter   :: year_step = 1   ! Timestep for the pem
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
236real              :: timestep        ! Timestep [s]
237character(100)    :: arg             ! To read command-line arguments program was invoked
238logical           :: timewall        ! Flag to use the time limit stopping criterion in case of a PEM job
239integer(kind = 8) :: cr              ! Number of clock ticks per second (count rate)
240integer(kind = 8) :: c1, c2          ! Counts of processor clock
241character(100)    :: chtimelimit     ! Time limit for the PEM job outputted by the SLURM command
242real              :: timelimit       ! Time limit for the PEM job in seconds
243real, parameter   :: antetime = 1200 ! Anticipation time to prevent reaching the time limit: 1200 s = 20 min by default
244integer           :: cstat, days, hours, minutes, seconds
245character(1)      :: sep
246
247#ifdef CPP_STD
248    real                                :: frost_albedo_threshold = 0.05 ! Frost albedo threeshold to convert fresh frost to old ice
249    real                                :: albedo_h2o_frost              ! Albedo of h2o frost
250    real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
251    real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
252    real, dimension(:,:),   allocatable :: tsoil_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
253    real, dimension(:),     allocatable :: emis_read_generic             ! Temporary variable to do the subslope transfert dimension when reading form generic
254    real, dimension(:,:),   allocatable :: albedo_read_generic           ! Temporary variable to do the subslope transfert dimension when reading form generic
255    real, dimension(:,:),   allocatable :: tsurf                         ! Subslope variable, only needed in the GENERIC case
256    real, dimension(:,:,:), allocatable :: qsurf                         ! Subslope variable, only needed in the GENERIC case
257    real, dimension(:,:,:), allocatable :: tsoil                         ! Subslope variable, only needed in the GENERIC case
258    real, dimension(:,:),   allocatable :: emis                          ! Subslope variable, only needed in the GENERIC case
259    real, dimension(:,:),   allocatable :: watercap                      ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
260    logical, dimension(:),  allocatable :: watercaptag                   ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
261    real, dimension(:,:,:), allocatable :: albedo                        ! Subslope variable, only needed in the GENERIC case
262    real, dimension(:,:,:), allocatable :: inertiesoil                   ! Subslope variable, only needed in the GENERIC case
263#endif
264
265#ifdef CPP_1D
266    integer            :: nsplit_phys
267    integer, parameter :: jjm_value = jjm - 1
268    integer            :: day_step
269
270    ! Dummy variables to use the subroutine 'init_testphys1d'
271    logical                             :: therestart1D, therestartfi
272    integer                             :: ndt, day0
273    real                                :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau
274    real, dimension(:),     allocatable :: zqsat
275    real, dimension(:,:,:), allocatable :: dq, dqdyn
276    real, dimension(nlayer)             :: play, w
277    real, dimension(nlayer + 1)         :: plev
278#else
279    integer, parameter                :: jjm_value = jjm
280    real, dimension(:),   allocatable :: ap ! Coefficient ap read in start_name and written in restart
281    real, dimension(:),   allocatable :: bp ! Coefficient bp read in start_name and written in restart
282    real, dimension(:,:), allocatable :: p  ! Grid points x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
283#endif
284
285! Loop variables
286integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
287
288! Elapsed time with system clock
289call system_clock(count_rate = cr)
290call system_clock(c1)
291timewall = .true.
292timelimit = 86400 ! 86400 seconds = 24 h by default
293timewall = .false.
294if (command_argument_count() > 0) then ! Get the number of command-line arguments
295    call get_command_argument(1,arg) ! Read the argument given to the program
296    select case (trim(adjustl(arg)))
297        case('version')
298            call print_version_info()
299            stop
300        case default ! This is the job id
301           ! Execute the system command
302            call execute_command_line('squeue -j '//trim(adjustl(arg))//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat)
303            if (cstat /= 0) then
304                call execute_command_line('qstat -f '//trim(adjustl(arg))//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt',cmdstat = cstat)
305                if (cstat > 0) then
306                    error stop 'pem: command execution failed!'
307                else if (cstat < 0) then
308                    error stop 'pem: command execution not supported (neither SLURM nor PBS/TORQUE is installed)!'
309                endif
310            endif
311            ! Read the output
312            open(1,file = 'tmp_cmdout.txt',status = 'old')
313            read(1,'(a)') chtimelimit
314            close(1)
315            chtimelimit = trim(adjustl(chtimelimit))
316            call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat)
317            if (cstat > 0) then
318                error stop 'pem: command execution failed!'
319            else if (cstat < 0) then
320                error stop 'pem: command execution not supported!'
321            endif
322            if (index(chtimelimit,'-') > 0) then ! 'chtimelimit' format is "D-HH:MM:SS"
323                read(chtimelimit,'(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds
324                timelimit = days*86400 + hours*3600 + minutes*60 + seconds
325            else if (index(chtimelimit,':') > 0 .and. len_trim(chtimelimit) > 5) then ! 'chtimelimit' format is "HH:MM:SS"
326                read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds
327                timelimit = hours*3600 + minutes*60 + seconds
328            else ! 'chtimelimit' format is "MM:SS"
329                read(chtimelimit,'(i2,a1,i2)') minutes, sep, seconds
330                timelimit = minutes*60 + seconds
331            endif
332    end select
333endif
334
335! Parallel variables
336#ifndef CPP_STD
337    is_sequential = .true.
338    is_parallel = .false.
339    is_mpi_root = .true.
340    is_omp_root = .true.
341    is_master = .true.
342#endif
343
344! Some constants
345day_ini = 0
346time_phys = 0.
347ngrid = ngridmx
348A = (1./m_co2 - 1./m_noco2)
349B = 1./m_noco2
350year_day = 669
351daysec = 88775.
352timestep = year_day*daysec*year_step
353
354!----------------------------- INITIALIZATION --------------------------
355!------------------------
356! I   Initialization
357!    I_a Read the "run.def"
358!------------------------
359#ifndef CPP_1D
360    dtphys = daysec/48. ! Dummy value (overwritten in phyetat0)
361    call conf_gcm(99,.true.)
362    call infotrac_init
363    nq = nqtot
364    allocate(q(ip1jmp1,llm,nqtot))
365    allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid))
366#else
367    allocate(q(1,llm,nqtot))
368    allocate(longitude(1),latitude(1),cell_area(1))
369
370    therestart1D = .false. ! Default value
371    inquire(file = 'start1D.txt',exist = therestart1D)
372    if (.not. therestart1D) then
373        write(*,*) 'There is no "start1D.txt" file!'
374        error stop 'Initialization cannot be done for the 1D PEM.'
375    endif
376    therestartfi = .false. ! Default value
377    inquire(file = 'startfi.nc',exist = therestartfi)
378    if (.not. therestartfi) then
379        write(*,*) 'There is no "startfi.nc" file!'
380        error stop 'Initialization cannot be done for the 1D PEM.'
381    endif
382
383    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,                 &
384                         time_0,ps_start_dyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
385                         play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
386    ps_start_dyn(2) = ps_start_dyn(1)
387    nsplit_phys = 1
388    day_step = steps_per_sol
389#endif
390
391call conf_pem(i_myear,n_myear)
392
393!------------------------
394! I   Initialization
395!    I_b Read of the "start.nc" and "starfi.nc"
396!------------------------
397! I_b.1 Read "start.nc"
398allocate(ps_start0(ngrid))
399#ifndef CPP_1D
400    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps_start_dyn,phis,time_0)
401
402    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps_start_dyn,ps_start0)
403
404    call iniconst ! Initialization of dynamical constants (comconst_mod)
405    call inigeom ! Initialization of geometry
406
407    allocate(ap(nlayer + 1))
408    allocate(bp(nlayer + 1))
409    status = nf90_open(start_name,NF90_NOWRITE,ncid)
410    status = nf90_inq_varid(ncid,"ap",apvarid)
411    status = nf90_get_var(ncid,apvarid,ap)
412    status = nf90_inq_varid(ncid,"bp",bpvarid)
413    status = nf90_get_var(ncid,bpvarid,bp)
414    status = nf90_close(ncid)
415
416    ! Initialization of physics constants and variables (comcstfi_h)
417    call iniphysiq(iim,jjm,llm,(jjm - 1)*iim + 2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
418#else
419    ps_start0(1) = ps_start_dyn(1)
420#endif
421
422! In the PCM, these values are given to the physic by the dynamic.
423! Here we simply read them in the "startfi.nc" file
424status = nf90_open(startfi_name,NF90_NOWRITE,ncid)
425status = nf90_inq_varid(ncid,"longitude",lonvarid)
426status = nf90_get_var(ncid,lonvarid,longitude)
427status = nf90_inq_varid(ncid,"latitude",latvarid)
428status = nf90_get_var(ncid,latvarid,latitude)
429status = nf90_inq_varid(ncid,"area",areavarid)
430status = nf90_get_var(ncid,areavarid,cell_area)
431status = nf90_inq_varid(ncid,"soildepth",sdvarid)
432status = nf90_get_var(ncid,sdvarid,mlayer)
433status = nf90_close(ncid)
434
435! I_b.2 Read the "startfi.nc"
436! First we read the initial state (starfi.nc)
437#ifndef CPP_STD
438    call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
439                  tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,          &
440                  watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
441
442    ! Remove unphysical values of surface tracer
443    where (qsurf < 0.) qsurf = 0.
444
445    call surfini(ngrid,nslope,qsurf)
446#else
447    call phys_state_var_init(nq)
448    if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
449    call initracer(ngrid,nq)
450    call iniaerosol()
451    allocate(tsurf_read_generic(ngrid))
452    allocate(qsurf_read_generic(ngrid,nq))
453    allocate(tsoil_read_generic(ngrid,nsoilmx))
454    allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
455    allocate(emis_read_generic(ngrid))
456    allocate(albedo_read_generic(ngrid,2))
457    allocate(qsurf(ngrid,nq,1))
458    allocate(tsurf(ngrid,1))
459    allocate(tsoil(ngrid,nsoilmx,1))
460    allocate(emis(ngrid,1))
461    allocate(watercap(ngrid,1))
462    allocate(watercaptag(ngrid))
463    allocate(albedo(ngrid,2,1))
464    allocate(inertiesoil(ngrid,nsoilmx,1))
465    call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
466                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,               &
467                  qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,        &
468                  rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
469    call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
470
471    nslope = 1
472    call ini_comslope_h(ngrid,1)
473
474    qsurf(:,:,1) = qsurf_read_generic
475    tsurf(:,1) = tsurf_read_generic
476    tsoil(:,:,1) = tsoil_read_generic
477    emis(:,1) = emis_read_generic
478    watercap(:,1) = 0.
479    watercaptag(:) = .false.
480    albedo(:,1,1) = albedo_read_generic(:,1)
481    albedo(:,2,1) = albedo_read_generic(:,2)
482    inertiesoil(:,:,1) = inertiedat
483
484    if (nslope == 1) then
485        def_slope(1) = 0
486        def_slope(2) = 0
487        def_slope_mean = 0
488        subslope_dist(:,1) = 1.
489    endif
490
491    ! Remove unphysical values of surface tracer
492    qsurf(:,:,1) = qsurf_read_generic
493    where (qsurf < 0.) qsurf = 0.
494
495    deallocate(tsurf_read_generic,qsurf_read_generic,qsoil_read_generic,emis_read_generic)
496#endif
497
498do nnq = 1,nqtot  ! Why not using ini_tracer ?
499    if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
500    if (noms(nnq) == "h2o_vap") then
501        igcm_h2o_vap = nnq
502        mmol(igcm_h2o_vap) = 18.
503    endif
504    if (noms(nnq) == "co2") igcm_co2 = nnq
505enddo
506
507!------------------------
508! I   Initialization
509!    I_c Subslope parametrisation
510!------------------------
511! Define some slope statistics
512iflat = 1
513do islope = 2,nslope
514    if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
515enddo
516write(*,*) 'Flat slope for islope = ',iflat
517write(*,*) 'corresponding criterium = ',def_slope_mean(iflat)
518
519!------------------------
520! I   Initialization
521!    I_d Read the PCM data and convert them to the physical grid
522!------------------------
523! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value
524call get_timelen_PCM("data_PCM_Y1.nc",timelen)
525
526allocate(vmr_co2_PCM(ngrid,timelen),q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
527allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid))
528allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
529allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope))
530allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_timeseries(ngrid,nsoilmx,nslope,timelen),tsoil_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
531allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
532
533call read_data_PCM("data_PCM_Y1.nc","data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,ps_avg,tsurf_avg_yr1,tsurf_avg, &
534                   tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2_PEM_phys,q_h2o_PEM_phys,watersurf_density_avg,watersoil_density_timeseries)
535
536! Compute the deviation from the average
537allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoilmx,nslope))
538ps_dev = ps_start0 - ps_avg
539tsurf_dev = tsurf - tsurf_avg
540tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:)
541
542!------------------------
543! I   Initialization
544!    I_e Initialization of the PEM variables and soil
545!------------------------
546call end_comsoil_h_PEM
547call ini_comsoil_h_PEM(ngrid,nslope)
548call end_adsorption_h_PEM
549call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
550call end_ice_table
551call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
552
553allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope),watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
554if (soil_pem) then
555    call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
556    tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg
557    watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
558    tsoil_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_timeseries
559    do l = nsoilmx + 1,nsoilmx_PEM
560        tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:)
561        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
562        tsoil_PEM_timeseries(:,l,:,:) = tsoil_timeseries(:,nsoilmx,:,:)
563    enddo
564    watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
565endif !soil_pem
566deallocate(tsoil_avg,watersoil_density_timeseries,tsoil_timeseries)
567
568!------------------------
569! I   Initialization
570!    I_f Compute tendencies
571!------------------------
572allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope))
573call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
574call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
575d_co2ice_ini = d_co2ice
576deallocate(min_co2_ice,min_h2o_ice)
577
578!------------------------
579! I   Initialization
580!    I_g Compute global surface pressure
581!------------------------
582total_surface = sum(cell_area)
583ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface
584ps_avg_global_new = ps_avg_global_ini
585write(*,*) "Total surface of the planet =", total_surface
586write(*,*) "Initial global averaged pressure =", ps_avg_global_ini
587
588!------------------------
589! I   Initialization
590!    I_h Read the "startpem.nc"
591!------------------------
592allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),stratif(ngrid,nslope))
593allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid))
594if (layering_algo) then
595    do islope = 1,nslope
596        do i = 1,ngrid
597            call ini_layering(stratif(i,islope))
598        enddo
599    enddo
600endif
601
602delta_h2o_icetablesublim = 0.
603call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
604              icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
605              ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice,watersurf_density_avg,         &
606              watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif)
607deallocate(tsurf_avg_yr1)
608
609! We save the places where h2o ice is sublimating
610! We compute the surface of h2o ice sublimating
611allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope),co2ice_disappeared(ngrid,nslope))
612co2ice_sublim_surf_ini = 0.
613h2oice_ini_surf = 0.
614is_co2ice_sublim_ini = .false.
615is_h2oice_sublim_ini = .false.
616is_co2ice_ini = .false.
617co2ice_disappeared = .false.
618do i = 1,ngrid
619    do islope = 1,nslope
620        if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
621        if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
622            is_co2ice_sublim_ini(i,islope) = .true.
623            co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope)
624        endif
625        if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
626            is_h2oice_sublim_ini(i,islope) = .true.
627            h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
628        endif
629    enddo
630enddo
631write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_sublim_surf_ini
632write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
633
634
635if (adsorption_pem) then
636    totmassco2_adsorbed = 0.
637    totmassh2o_adsorbed = 0.
638    do ig = 1,ngrid
639        do islope = 1,nslope
640            do l = 1,nsoilmx_PEM - 1
641                if (l == 1) then
642                   totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
643                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
644                   totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
645                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
646                else
647                   totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
648                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
649                   totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
650                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
651                endif
652            enddo
653        enddo
654    enddo
655    write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbed
656    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbed
657endif ! adsorption
658
659!------------------------
660! I   Initialization
661!    I_i Compute orbit criterion
662!------------------------
663#ifndef CPP_STD
664    call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
665#else
666    call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
667#endif
668
669n_myear_leg = Max_iter_pem
670if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
671
672!-------------------------- END INITIALIZATION -------------------------
673
674!-------------------------------- RUN ----------------------------------
675!------------------------
676! II  Run
677!    II_a Update pressure, ice and tracers
678!------------------------
679i_myear_leg = 0
680stopPEM = 0
681if (layering_algo) then
682    allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope))
683    new_str = .true.
684    new_lag1 = .true.
685    new_lag2 = .true.
686    do islope = 1,nslope
687        do ig = 1,ngrid
688            current1(ig,islope)%p => stratif(ig,islope)%top
689            current2(ig,islope)%p => stratif(ig,islope)%top
690        enddo
691    enddo
692endif
693
694do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
695! II.a.1. Compute updated global pressure
696    write(*,*) "Recomputing the new pressure..."
697    ps_avg_global_old = ps_avg_global_new
698    do i = 1,ngrid
699        do islope = 1,nslope
700            ps_avg_global_new = ps_avg_global_old - CO2cond_ps*g*cell_area(i)*d_co2ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/total_surface
701        enddo
702    enddo
703
704    if (adsorption_pem) then
705        do i = 1,ngrid
706            ps_avg_global_new = ps_avg_global_old - g*cell_area(i)*delta_co2_adsorbed(i)/total_surface
707        enddo
708    endif
709    write(*,*) 'Global average pressure old time step',ps_avg_global_old
710    write(*,*) 'Global average pressure new time step',ps_avg_global_new
711
712! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption)
713    allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen))
714    write(*,*) "Recomputing the pressure levels timeseries adapted to the old pressure..."
715    do l = 1,nlayer + 1
716        do ig = 1,ngrid
717            zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
718        enddo
719    enddo
720    write(*,*) "Recomputing the surface pressure timeseries adapted to the new pressure..."
721    do ig = 1,ngrid
722        ps_timeseries(ig,:) = ps_timeseries(ig,:)*ps_avg_global_new/ps_avg_global_old
723    enddo
724    write(*,*) "Recomputing the pressure levels timeseries adapted to the new pressure..."
725    allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
726    do l = 1,nlayer + 1
727        do ig = 1,ngrid
728            zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
729        enddo
730    enddo
731
732! II.a.3. Tracers timeseries
733    write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..."
734    allocate(vmr_co2_PEM_phys(ngrid,timelen))
735    l = 1
736    do ig = 1,ngrid
737        do t = 1,timelen
738            ! H2O
739            q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ &
740                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
741            if (q_h2o_PEM_phys(ig,t) < 0) then
742                q_h2o_PEM_phys(ig,t) = 1.e-30
743            else if (q_h2o_PEM_phys(ig,t) > 1) then
744                q_h2o_PEM_phys(ig,t) = 1.
745            endif
746            ! CO2
747            q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ &
748                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))                       &
749                                + ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))                       &
750                                -  (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/                     &
751                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
752            if (q_co2_PEM_phys(ig,t) < 0) then
753                q_co2_PEM_phys(ig,t) = 1.e-30
754            else if (q_co2_PEM_phys(ig,t) > 1) then
755                q_co2_PEM_phys(ig,t) = 1.
756            endif
757            mmean = 1./(A*q_co2_PEM_phys(ig,t) + B)
758            vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
759        enddo
760    enddo
761    deallocate(zplev_timeseries_new,zplev_timeseries_old)
762
763!------------------------
764! II  Run
765!    II_b Evolution of ice
766!------------------------
767    allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
768    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
769    call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
770    if (layering_algo) then
771        do islope = 1,nslope
772            do ig = 1,ngrid
773                call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),zshift_surf(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),zlag(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
774            enddo
775        enddo
776    else
777        zlag = 0.
778    endif
779
780!------------------------
781! II  Run
782!    II_c Flow of glaciers
783!------------------------
784    allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope))
785    if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
786                                                            ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
787    if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
788
789!------------------------
790! II  Run
791!    II_d Update surface and soil temperatures
792!------------------------
793! II_d.1 Update Tsurf
794    write(*,*) "Updating the new Tsurf"
795    do ig = 1,ngrid
796        do islope = 1,nslope
797            ! CO2 ice disappeared so we look for the closest point without CO2 ice
798            if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then
799                co2ice_disappeared(ig,islope) = .true.
800                if (latitude_deg(ig) > 0) then
801                    outer1: do ig_loop = ig,ngrid
802                        do islope_loop = islope,iflat,-1
803                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
804                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
805                                exit outer1
806                            endif
807                        enddo
808                    enddo outer1
809                else
810                    outer2: do ig_loop = ig,1,-1
811                        do islope_loop = islope,iflat
812                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
813                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
814                                exit outer2
815                            endif
816                        enddo
817                    enddo outer2
818                endif
819            else if (co2_ice(ig,islope) > 1.e-10 .and. d_co2ice(ig,islope) > 1.e-10) then ! Put tsurf as tcond CO2
820                call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_ini,ps_avg_global_new,tsurf_avg)
821            endif
822        enddo
823    enddo
824
825    if (soil_pem) then
826! II_d.2 Shifting soil temperature to surface
827        call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
828        deallocate(zshift_surf,zlag)
829
830! II_d.3 Update soil temperature
831        write(*,*)"Updating soil temperature"
832        allocate(tsoil_avg_old(ngrid,nsoilmx_PEM))
833        do islope = 1,nslope
834            tsoil_avg_old = tsoil_PEM(:,:,islope)
835            call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
836            call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
837
838            do t = 1,timelen
839                do ig = 1,ngrid
840                    do isoil = 1,nsoilmx_PEM
841            ! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
842            tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
843            ! Update of watersoil density
844                        watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/tsoil_PEM_timeseries(ig,isoil,islope,t) + alpha_clap_h2o)/tsoil_PEM_timeseries(ig,isoil,islope,t)*mmol(igcm_h2o_vap)/(mugaz*r)
845                        if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
846                    enddo
847                enddo
848            enddo
849        enddo
850        watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
851        deallocate(tsoil_avg_old)
852        write(*,*) "Update of soil temperature done"
853
854! II_d.4 Update the ice table
855        allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
856        if (icetable_equilibrium) then
857            write(*,*) "Compute ice table (equilibrium method)"
858            icetable_thickness_old = icetable_thickness
859            call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
860            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
861        else if (icetable_dynamic) then
862            write(*,*) "Compute ice table (dynamic method)"
863            ice_porefilling_old = ice_porefilling
864            allocate(porefill(nsoilmx_PEM))
865            do ig = 1,ngrid
866                do islope = 1,nslope
867                    call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps_avg(ig),(/sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2)/),ice_porefilling(ig,:,islope),porefill,ssi_depth)
868                    icetable_depth(ig,islope) = ssi_depth
869                    ice_porefilling(ig,:,islope) = porefill
870                enddo
871            enddo
872            deallocate(porefill)
873            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
874        endif
875        deallocate(icetable_thickness_old,ice_porefilling_old)
876
877! II_d.5 Update the soil thermal properties
878        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,ps_avg_global_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
879
880! II_d.6 Update the mass of the regolith adsorbed
881        if (adsorption_pem) then
882            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
883                                     tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,       &
884                                     h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
885
886            totmassco2_adsorbed = 0.
887            totmassh2o_adsorbed = 0.
888            do ig = 1,ngrid
889                do islope = 1,nslope
890                    do l = 1,nsoilmx_PEM
891                        if (l == 1) then
892                            totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
893                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
894                            totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
895                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
896                        else
897                            totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
898                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
899                            totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
900                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
901                        endif
902                    enddo
903                enddo
904            enddo
905            write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbed
906            write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbed
907        endif
908    endif !soil_pem
909
910!------------------------
911! II  Run
912!    II_e Outputs
913!------------------------
914    call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/))
915    do islope = 1,nslope
916        write(str2(1:2),'(i2.2)') islope
917        call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
918        call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
919        call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
920        call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
921        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope)))
922        call writediagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope)))
923        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
924        if (icetable_equilibrium) then
925            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
926            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope))
927        else if (icetable_dynamic) then
928            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
929        endif
930
931        if (soil_pem) then
932            call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope))
933            call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope))
934            if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
935            if (adsorption_pem) then
936                call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope))
937                call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope))
938            endif
939        endif
940    enddo
941    deallocate(flag_co2flow,flag_h2oflow)
942
943!------------------------
944! II  Run
945!    II_f Update the tendencies
946!------------------------
947    call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_old,ps_avg_global_new)
948    deallocate(vmr_co2_PEM_phys)
949
950!------------------------
951! II  Run
952!    II_g Checking the stopping criterion
953!------------------------
954    write(*,*) "Checking the stopping criteria..."
955    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
956    call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope)
957    i_myear_leg = i_myear_leg + dt
958    i_myear = i_myear + dt
959    if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
960    if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
961    call system_clock(c2)
962    if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
963    if (stopPEM > 0) then
964        select case (stopPEM)
965            case(1)
966                write(*,*) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above."
967            case(2)
968                write(*,*) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above."
969            case(3)
970                write(*,*) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
971            case(4)
972                write(*,*) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
973            case(5)
974                write(*,*) "STOPPING because maximum number of iterations is reached (possibly due to orbital parameters):", stopPEM
975            case(6)
976                write(*,*) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
977            case(7)
978                write(*,*) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM
979            case(8)
980                write(*,*) "STOPPING because the layering algorithm met an hasty end:", stopPEM
981            case default
982                write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM
983        end select
984        exit
985    else
986        write(*,*) "The PEM can continue!"
987        write(*,*) "Number of iterations of the PEM: i_myear_leg =", i_myear_leg
988        write(*,*) "Number of simulated Martian years: i_myear =", i_myear
989    endif
990enddo ! big time iteration loop of the pem
991deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed)
992deallocate(watersoil_density_PEM_avg,watersurf_density_avg,)
993deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries)
994deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim)
995deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
996deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
997if (layering_algo) then
998    do islope = 1,nslope
999        do i = 1,ngrid
1000            call del_layering(stratif(i,islope))
1001        enddo
1002    enddo
1003    deallocate(new_str,new_lag1,new_lag2,current1,current2)
1004endif
1005!------------------------------ END RUN --------------------------------
1006
1007!------------------------------- OUTPUT --------------------------------
1008!------------------------
1009! III Output
1010!    III_a Update surface value for the PCM start files
1011!------------------------
1012! III_a.1 Ice update for start file
1013watercap = 0.
1014perennial_co2ice = co2_ice
1015do ig = 1,ngrid
1016    ! H2O ice metamorphism
1017    if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
1018        h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
1019        qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
1020    endif
1021
1022    ! Is H2O ice still considered as an infinite reservoir for the PCM?
1023    if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
1024        ! There is enough ice to be considered as an infinite reservoir
1025        watercaptag(ig) = .true.
1026    else
1027        ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
1028        watercaptag(ig) = .false.
1029        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1030        h2o_ice(ig,:) = 0.
1031    endif
1032
1033    ! CO2 ice metamorphism
1034    if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
1035        perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
1036        qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
1037    endif
1038enddo
1039
1040! III.a.2. Tsurf update for start file
1041tsurf = tsurf_avg + tsurf_dev
1042deallocate(tsurf_dev)
1043
1044! III_a.3 Tsoil update for start file
1045if (soil_pem) then
1046    inertiesoil = TI_PEM(:,:nsoilmx,:)
1047    ! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value
1048    do isoil = 1,nsoilmx
1049        tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:)
1050    enddo
1051    tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev
1052#ifndef CPP_STD
1053    flux_geo = fluxgeo
1054#endif
1055endif
1056deallocate(tsurf_avg,tsoil_dev)
1057
1058! III_a.4 Pressure update for start file
1059allocate(ps_start(ngrid))
1060ps_start = ps_avg + ps_dev
1061deallocate(ps_avg,ps_dev)
1062
1063! III_a.5 Tracers update for start file
1064allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1))
1065do l = 1,nlayer + 1
1066    zplev_start0(:,l) = ap(l) + bp(l)*ps_start0
1067    zplev_new(:,l) = ap(l) + bp(l)*ps_start
1068enddo
1069
1070do nnq = 1,nqtot
1071    if (noms(nnq) /= "co2") then
1072        do l = 1,llm - 1
1073            do ig = 1,ngrid
1074                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
1075            enddo
1076            q(:,llm,nnq) = q(:,llm - 1,nnq)
1077        enddo
1078    else
1079        do l = 1,llm - 1
1080            do ig = 1,ngrid
1081                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
1082                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_start0(ig,l) - zplev_start0(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
1083            enddo
1084            q(:,llm,nnq) = q(:,llm - 1,nnq)
1085        enddo
1086    endif
1087enddo
1088deallocate(zplev_start0)
1089
1090! Conserving the tracers mass for start file
1091do nnq = 1,nqtot
1092    do ig = 1,ngrid
1093        do l = 1,llm - 1
1094            if (q(ig,l,nnq) > 1 .and. (noms(nnq) /= "dust_number") .and. (noms(nnq) /= "ccn_number") .and. (noms(nnq) /= "stormdust_number") .and. (noms(nnq) /= "topdust_number")) then
1095                extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1096                q(ig,l,nnq) = 1.
1097                q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1098                write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1099           endif
1100            if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1101        enddo
1102    enddo
1103enddo
1104deallocate(zplev_new)
1105
1106! III_a.6 Albedo update for start file
1107do ig = 1,ngrid
1108    if (latitude(ig) < 0.) then
1109        icap = 2 ! Southern hemisphere
1110    else
1111        icap = 1 ! Northern hemisphere
1112    endif
1113    do islope = 1,ngrid
1114        ! Bare ground
1115        albedo(ig,:,islope) = albedodat(ig)
1116        emis(ig,islope) = emissiv
1117
1118        ! CO2 ice/frost is treated after H20 ice/frost because it is considered dominant
1119        ! H2O ice
1120        if (h2o_ice(ig,islope) > 0.) then
1121            albedo(ig,:,islope) = albedo_h2o_cap
1122            emis(ig,islope) = 1.
1123        endif
1124        ! CO2 ice
1125        if (co2_ice(ig,islope) > 0.) then
1126            albedo(ig,:,islope) = albedo_perennialco2(icap)
1127            emis(ig,islope) = emisice(icap)
1128        endif
1129        ! H2O frost
1130        if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then
1131            albedo(ig,:,islope) = albedo_h2o_frost
1132            emis(ig,islope) = 1.
1133        endif
1134        ! CO2 frost
1135        if (qsurf(ig,igcm_co2,islope) > 0.) then
1136            albedo(ig,:,islope) = albedice(icap)
1137            emis(ig,islope) = emisice(icap)
1138        endif
1139    enddo
1140enddo
1141
1142! III_a.7 Orbital parameters update for start file
1143if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
1144
1145!------------------------
1146! III Output
1147!    III_b Write "restart.nc" and "restartfi.nc"
1148!------------------------
1149! III_b.1 Write "restart.nc"
1150ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1151pday = day_ini
1152ztime_fin = time_phys
1153
1154#ifndef CPP_1D
1155    allocate(p(ip1jmp1,nlayer + 1))
1156    ! Correction on teta due to surface pressure changes
1157    do l = 1,nlayer
1158        do i = 1,ip1jmp1
1159            teta(i,l)= teta(i,l)*(ps_start0(i)/ps_start(i))**kappa
1160        enddo
1161    enddo
1162    ! Correction on atmospheric pressure
1163    call pression(ip1jmp1,ap,bp,ps_start,p)
1164    ! Correction on the mass of atmosphere
1165    call massdair(p,masse)
1166    call dynredem0("restart.nc",day_ini,phis)
1167    call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps_start)
1168    write(*,*) "restart.nc has been written."
1169    deallocate(ap,bp,p)
1170#else
1171    call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1172    write(*,*) "restart1D.txt has been written."
1173#endif
1174deallocate(ps_start0,ps_start)
1175
1176! III_b.2 Write the "restartfi.nc"
1177#ifndef CPP_STD
1178    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1179                  nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1180                  inertiedat,def_slope,subslope_dist)
1181    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,      &
1182                  ptimestep,ztime_fin,tsurf,tsoil,inertiesoil,        &
1183                  albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1184                  wstar,watercap,perennial_co2ice)
1185#else
1186    if (allocated(noms)) deallocate(noms)
1187    deallocate(qsurf,tsurf,tsoil,emis,watercap,watercaptag,albedo,inertiesoil)
1188    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1189                  nlayer,nq,ptimestep,pday,time_phys,cell_area,    &
1190                  albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1191    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,       &
1192                  ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1193                  cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab,   &
1194                  tsea_ice,sea_ice)
1195#endif
1196write(*,*) "restartfi.nc has been written."
1197
1198!------------------------
1199! III Output
1200!    III_c Write the "restartpem.nc"
1201!------------------------
1202if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
1203call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
1204call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
1205             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,stratif)
1206write(*,*) "restartpem.nc has been written."
1207
1208call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
1209
1210write(*,*) "The PEM has run for", i_myear_leg, "Martian years."
1211write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1212write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
1213write(*,*) "PEM: so far, so good!"
1214
1215if (layering_algo) then
1216    do islope = 1,nslope
1217        do i = 1,ngrid
1218            call del_layering(stratif(i,islope))
1219        enddo
1220    enddo
1221endif
1222deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
1223deallocate(co2_ice,h2o_ice,stratif)
1224!----------------------------- END OUTPUT ------------------------------
1225
1226END PROGRAM pem
Note: See TracBrowser for help on using the repository browser.