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

Last change on this file since 3440 was 3432, checked in by jbclement, 15 months ago

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