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

Last change on this file since 3742 was 3742, checked in by jbclement, 3 months ago

PEM:
Resulting modifications due to changes in the Mars PCM with r3739 and r3741.
JBC

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