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

Last change on this file since 3116 was 3114, checked in by jbclement, 21 months ago

PEM:
Following r3113, addition of 'nqsoil' and 'qsoil' in the arguments of the subroutines 'phyetat0' and 'physdem1' to be able to compile.
JBC

File size: 57.4 KB
Line 
1!------------------------
2! I   Initialization
3!    I_a READ run.def
4!    I_b READ of start_evol.nc and starfi_evol.nc
5!    I_c Subslope parametrisation
6!    I_d READ PCM data and convert to the physical grid
7!    I_e Initialization of the PEM variable and soil
8!    I_f Compute tendencies & Save initial situation
9!    I_g Save initial PCM 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 the ice
16!    II_c CO2 & H2O glaciers flows
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 restart_evol.nc and restartfi_evol.nc
25!    III_c Write 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, major_slope, 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: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
42use criterion_pem_stop_mod,       only: criterion_waterice_stop, criterion_co2_stop
43use constants_marspem_mod,        only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
44                                        m_noco2, threshold_water_frost2perenial, threshold_co2_frost2perenial
45use evol_co2_ice_s_mod,           only: evol_co2_ice_s
46use evol_h2o_ice_s_mod,           only: evol_h2o_ice_s
47use comsoil_h_PEM,                only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
48                                        TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
49                                        tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
50                                        fluxgeo,                          & ! Geothermal flux for the PEM and PCM
51                                        water_reservoir                     ! Water ressources
52use adsorption_mod,               only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
53                                        ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
54                                        co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
55use time_evol_mod,                only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
56use orbit_param_criterion_mod,    only: orbit_param_criterion
57use recomp_orb_param_mod,         only: recomp_orb_param
58use ice_table_mod,                only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
59                                        ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
60use soil_thermalproperties_mod,   only: update_soil_thermalproperties
61use time_phylmdz_mod,             only: daysec, dtphys, day_end
62use abort_pem_mod,                only: abort_pem
63use soil_settings_PEM_mod,        only: soil_settings_PEM
64use compute_tendencies_slope_mod, only: compute_tendencies_slope
65use info_PEM_mod,                 only: info_PEM
66use interpolate_TIPEM_TIGCM_mod,  only: interpolate_TIPEM_TIGCM
67use nb_time_step_PCM_mod,         only: nb_time_step_PCM
68use pemetat0_mod,                 only: pemetat0
69use read_data_PCM_mod,            only: read_data_PCM
70use recomp_tend_co2_slope_mod,    only: recomp_tend_co2_slope
71use soil_pem_compute_mod,         only: soil_pem_compute
72use writediagpem_mod,             only: writediagpem
73
74#ifndef CPP_STD
75    use comsoil_h,          only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, volcapa, inertiesoil, nqsoil, qsoil
76    use surfdat_h,          only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
77                                  albedodat, zmea, zstd, zsig, zgam, zthe,     &
78                                  hmons, summit, base,albedo_h2o_frost,        &
79                                  frost_albedo_threshold, emissiv, watercaptag, perenial_co2ice, &
80                                  emisice, albedice
81    use dimradmars_mod,     only: totcloudfrac, albedo
82    use dust_param_mod,     only: tauscaling
83    use tracer_mod,         only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
84    use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
85    use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
86    use paleoclimate_mod,   only: albedo_perenialco2
87    use comcstfi_h,         only: pi, rad, g, cpp, mugaz, r
88#else
89    use tracer_h,           only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
90    use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT,   &
91                                  PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
92                                  ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
93    use aerosol_mod,        only: iniaerosol
94    use planete_mod,        only: apoastr, periastr, year_day, peri_day, obliquit
95    use comcstfi_mod,       only: pi, rad, g, cpp, mugaz, r
96#endif
97
98#ifndef CPP_1D
99    use iniphysiq_mod,            only: iniphysiq
100    use control_mod,              only: iphysiq, day_step, nsplit_phys
101#else
102    use time_phylmdz_mod,         only: iphysiq, day_step
103    use regular_lonlat_mod,       only: init_regular_lonlat
104    use physics_distribution_mod, only: init_physics_distribution
105    use mod_grid_phy_lmdz,        only: regular_lonlat
106    use init_testphys1d_mod,      only: init_testphys1d
107    use comvert_mod,              only: ap, bp
108    use writerestart1D_mod,       only: writerestart1D
109#endif
110
111implicit none
112
113include "dimensions.h"
114include "paramet.h"
115include "comgeom.h"
116include "iniprint.h"
117include "callkeys.h"
118
119integer ngridmx
120parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
121
122! Same variable names as in the PCM
123integer, parameter :: nlayer = llm ! Number of vertical layer
124integer            :: ngrid        ! Number of physical grid points
125integer            :: nq           ! Number of tracer
126integer            :: day_ini      ! First day of the simulation
127real               :: pday         ! Physical day
128real               :: time_phys    ! Same as PCM
129real               :: ptimestep    ! Same as PCM
130real               :: ztime_fin    ! Same as PCM
131
132! Variables to read start.nc
133character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM
134
135! Dynamic variables
136real, dimension(ip1jm,llm)          :: vcov          ! vents covariants
137real, dimension(ip1jmp1,llm)        :: ucov          ! vents covariants
138real, dimension(ip1jmp1,llm)        :: teta          ! temperature potentielle
139real, dimension(:,:,:), allocatable :: q             ! champs advectes
140real, dimension(ip1jmp1)            :: ps            ! pression  au sol
141real, dimension(:),     allocatable :: ps_start_GCM  ! (ngrid) pression  au sol
142real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) ! pression  au sol instantannées
143real, dimension(ip1jmp1,llm)        :: masse         ! masse d'air
144real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
145real                                :: time_0
146
147! Variables to read starfi.nc
148character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" ! Name of the file used for initialsing the PEM
149character(2)                   :: str2
150integer                        :: ncid, varid, status                      ! Variable for handling opening of files
151integer                        :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files
152integer                        :: lonvarid, latvarid, areavarid, sdvarid   ! Variable ID for Netcdf files
153integer                        :: apvarid, bpvarid                         ! Variable ID for Netcdf files
154
155! Variables to read starfi.nc and write restartfi.nc
156real, dimension(:), allocatable :: longitude     ! Longitude read in FILE_NAME and written in restartfi
157real, dimension(:), allocatable :: latitude      ! Latitude read in FILE_NAME and written in restartfi
158real, dimension(:), allocatable :: cell_area     ! Cell_area read in FILE_NAME and written in restartfi
159real                            :: Total_surface ! Total surface of the planet
160
161! Variables for h2o_ice evolution
162real                                :: ini_surf_h2o         ! Initial surface of sublimating h2o ice
163real                                :: ini_surf_co2         ! Initial surface of sublimating co2 ice
164real                                :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa]
165real                                :: global_ave_press_old ! constant: Global average pressure of initial/previous time step
166real                                :: global_ave_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_gcm            ! 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
171logical                             :: STOPPING_water       ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached?
172logical                             :: STOPPING_1_water     ! Logical: is there still water ice to sublimate?
173logical                             :: STOPPING_co2         ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached?
174logical                             :: STOPPING_pressure    ! Logical: is the criterion (% of change in the surface pressure) reached?
175integer                             :: criterion_stop       ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
176real, save                          :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
177real, dimension(:,:),   allocatable :: vmr_co2_gcm          ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
178real, dimension(:,:),   allocatable :: vmr_co2_pem_phys     ! Physics x Times  co2 volume mixing ratio used in the PEM
179real, 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]
180real, dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
181integer                             :: timelen              ! # time samples
182real                                :: ave                  ! intermediate varibale to compute average
183real, dimension(:,:),   allocatable :: p                    ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
184real                                :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
185
186! Variables for slopes
187real, dimension(:,:),   allocatable :: min_co2_ice_1          ! ngrid field: minimum of co2 ice at each point for the first year [kg/m^2]
188real, dimension(:,:),   allocatable :: min_co2_ice_2          ! ngrid field: minimum of co2 ice at each point for the second year [kg/m^2]
189real, dimension(:,:),   allocatable :: min_h2o_ice_1          ! ngrid field: minimum of water ice at each point for the first year [kg/m^2]
190real, dimension(:,:),   allocatable :: min_h2o_ice_2          ! ngrid field: minimum of water ice at each point for the second year [kg/m^2]
191real, dimension(:,:,:), allocatable :: co2_ice_GCM            ! Physics x NSLOPE x Times field: co2 ice given by the PCM  [kg/m^2]
192real, dimension(:,:),   allocatable :: initial_co2_ice_sublim ! physical point field: Logical array indicating sublimating point of co2 ice
193real, dimension(:,:),   allocatable :: initial_h2o_ice        ! physical point field: Logical array indicating if there is water ice at initial state
194real, dimension(:,:),   allocatable :: initial_co2_ice        ! physical point field: Logical array indicating if there is co2 ice at initial state
195real, dimension(:,:),   allocatable :: tendencies_co2_ice     ! physical point x slope field: Tendency of evolution of perenial co2 ice over a year
196real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the PCM
197real, dimension(:,:),   allocatable :: tendencies_h2o_ice     ! physical point x slope  field: Tendency of evolution of perenial 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_ave                          ! Physic x SLOPE field : Averaged Surface Temperature [K]
205real, dimension(:,:,:),   allocatable :: tsoil_ave                          ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
206real, dimension(:,:,:),   allocatable :: tsurf_GCM_timeseries               ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
207real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries          ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
208real, dimension(:,:,:,:), allocatable :: tsoil_GCM_timeseries               ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
209real, dimension(:,:),     allocatable :: tsurf_ave_yr1                      ! Physic x SLOPE field : Averaged Surface Temperature of first call of the PCM [K]
210real, dimension(:,:),     allocatable :: TI_locslope                        ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
211real, dimension(:,:),     allocatable :: Tsoil_locslope                     ! Physic x Soil: intermediate when computing Tsoil [K]
212real, dimension(:),       allocatable :: Tsurf_locslope                     ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
213real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries       ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
214real, dimension(:,:),     allocatable :: watersurf_density_ave              ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
215real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries   ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
216real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_ave          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
217real, dimension(:,:),     allocatable :: Tsurfave_before_saved              ! Surface temperature saved from previous time step [K]
218real, dimension(:),       allocatable :: delta_co2_adsorbded                ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
219real, dimension(:),       allocatable :: delta_h2o_adsorbded                ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
220real                                  :: totmassco2_adsorbded               ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
221real                                  :: totmassh2o_adsorbded               ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
222logical                               :: bool_sublim                        ! logical to check if there is sublimation or not
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]
233real            :: watercap_sum  ! total mass of water cap [kg/m^2]
234real            :: water_sum     ! total mass of water in the mesh [kg/m^2]
235
236#ifdef CPP_STD
237    real                                :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
238    real                                :: albedo_h2o_frost              ! albedo of h2o frost
239    real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimensiion when reading form generic
240    real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimensiion when reading form generic
241    real, dimension(:,:),   allocatable :: tsoil_read_generic            ! Temporary variable to do the subslope transfert dimensiion when reading form generic
242    real, dimension(:),     allocatable :: emis_read_generic             ! Temporary variable to do the subslope transfert dimensiion when reading form generic
243    real, dimension(:,:),   allocatable :: albedo_read_generic           ! Temporary variable to do the subslope transfert dimensiion when reading form generic
244    real, dimension(:,:),   allocatable :: tsurf                         ! Subslope variable, only needed in the GENERIC case
245    real, dimension(:,:,:), allocatable :: qsurf                         ! Subslope variable, only needed in the GENERIC case
246    real, dimension(:,:,:), allocatable :: tsoil                         ! Subslope variable, only needed in the GENERIC case
247    real, dimension(:,:),   allocatable :: emis                          ! Subslope variable, only needed in the GENERIC case
248    real, dimension(:,:),   allocatable :: watercap                      ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
249    logical, dimension(:),  allocatable :: watercaptag                   ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
250    real, dimension(:,:,:), allocatable :: albedo                        ! Subslope variable, only needed in the GENERIC case
251    real, dimension(:,:,:), allocatable :: inertiesoil                   ! Subslope variable, only needed in the GENERIC case
252#endif
253
254#ifdef CPP_1D
255    integer                           :: nsplit_phys
256    integer, parameter                :: jjm_value = jjm - 1
257    integer                           :: ierr
258
259    ! Dummy variables to use the subroutine 'init_testphys1d'
260    logical                             :: startfiles_1D, therestart1D, therestartfi
261    integer                             :: ndt, day0
262    real                                :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau
263    real, dimension(:),     allocatable :: zqsat
264    real, dimension(:,:,:), allocatable :: dq, dqdyn
265    real, dimension(nlayer)             :: play, w
266    real, dimension(nlayer + 1)         :: plev
267#else
268    integer, parameter                :: jjm_value = jjm
269    real, dimension(:), allocatable   :: ap ! Coefficient ap read in FILE_NAME_start and written in restart
270    real, dimension(:), allocatable   :: bp ! Coefficient bp read in FILE_NAME_start and written in restart
271#endif
272
273! Loop variables
274integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil, icap
275
276! Parallel variables
277#ifndef CPP_STD
278    is_sequential = .true.
279    is_parallel = .false.
280    is_mpi_root = .true.
281    is_omp_root = .true.
282    is_master = .true.
283#endif
284
285! Some constants
286day_ini = 0    ! test
287time_phys = 0. ! test
288ngrid = ngridmx
289A = (1/m_co2 - 1/m_noco2)
290B = 1/m_noco2
291year_day = 669
292daysec = 88775.
293timestep = year_day*daysec/year_step
294
295!----------------------------- INITIALIZATION --------------------------
296!------------------------
297! I   Initialization
298!    I_a READ run.def
299!------------------------
300#ifndef CPP_1D
301    dtphys = 0
302    call conf_gcm(99,.true.)
303    call infotrac_init
304    nq = nqtot
305    allocate(q(ip1jmp1,llm,nqtot))
306    allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid))
307#else
308    allocate(q(1,llm,nqtot))
309    allocate(longitude(1),latitude(1),cell_area(1))
310    call init_testphys1d(.true.,ngrid,nlayer,610.,nq,q,time_0,ps(1),ucov,vcov,teta,startfiles_1D,therestart1D, &
311                         therestartfi,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,                   &
312                         play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
313    ps(2) = ps(1)
314    nsplit_phys = 1
315#endif
316
317call conf_pem(i_myear,n_myear)
318
319!------------------------
320! I   Initialization
321!    I_b READ of start_evol.nc and starfi_evol.nc
322!------------------------
323! I_b.1 READ start_evol.nc
324allocate(ps_start_GCM(ngrid))
325#ifndef CPP_1D
326    call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0)
327
328    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
329
330    call iniconst !new
331    call inigeom
332
333    allocate(ap(nlayer + 1))
334    allocate(bp(nlayer + 1))
335    status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid)
336    status = nf90_inq_varid(ncid,"ap",apvarid)
337    status = nf90_get_var(ncid,apvarid,ap)
338    status = nf90_inq_varid(ncid,"bp",bpvarid)
339    status = nf90_get_var(ncid,bpvarid,bp)
340    status = nf90_close(ncid)
341
342    call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys, &
343                   rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
344#else
345    ps_start_GCM(1) = ps(1)
346#endif
347
348! In the PCM, these values are given to the physic by the dynamic.
349! Here we simply read them in the startfi_evol.nc file
350status = nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
351
352status = nf90_inq_varid(ncid,"longitude",lonvarid)
353status = nf90_get_var(ncid,lonvarid,longitude)
354
355status = nf90_inq_varid(ncid,"latitude",latvarid)
356status = nf90_get_var(ncid,latvarid,latitude)
357
358status = nf90_inq_varid(ncid,"area",areavarid)
359status = nf90_get_var(ncid,areavarid,cell_area)
360
361status = nf90_inq_varid(ncid,"soildepth",sdvarid)
362status = nf90_get_var(ncid,sdvarid,mlayer)
363
364status = nf90_close(ncid)
365
366! I_b.2 READ startfi_evol.nc
367! First we read the initial state (starfi.nc)
368#ifndef CPP_STD
369    call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
370                  tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,       &
371                  watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist)
372
373    ! Remove unphysical values of surface tracer
374    where (qsurf < 0.) qsurf = 0.
375
376    call surfini(ngrid,qsurf)
377#else
378    call phys_state_var_init(nq)
379    if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
380    call initracer(ngrid,nq)
381    call iniaerosol()
382    allocate(tsurf_read_generic(ngrid))
383    allocate(qsurf_read_generic(ngrid,nq))
384    allocate(tsoil_read_generic(ngrid,nsoilmx))
385    allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
386    allocate(emis_read_generic(ngrid))
387    allocate(tsurf(ngrid,1))
388    allocate(qsurf(ngrid,nq,1))
389    allocate(tsoil(ngrid,nsoilmx,1))
390    allocate(emis(ngrid,1))
391    allocate(watercap(ngrid,1))
392    allocate(watercaptag(ngrid))
393    allocate(albedo_read_generic(ngrid,2))
394    allocate(albedo(ngrid,2,1))
395    allocate(inertiesoil(ngrid,nsoilmx,1))
396    call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
397                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,            &
398                  qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,     &
399                  rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
400    call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
401
402    nslope = 1
403    call ini_comslope_h(ngrid,1)
404
405    qsurf(:,:,1) = qsurf_read_generic(:,:)
406    tsurf(:,1) = tsurf_read_generic(:)
407    tsoil(:,:,1) = tsoil_read_generic(:,:)
408    emis(:,1) = emis_read_generic(:)
409    watercap(:,1) = 0.
410    watercaptag(:) = .false.
411    albedo(:,1,1) = albedo_read_generic(:,1)
412    albedo(:,2,1) = albedo_read_generic(:,2)
413    inertiesoil(:,:,1) = inertiedat(:,:)
414
415    if (nslope == 1) then
416        def_slope(1) = 0
417        def_slope(2) = 0
418        def_slope_mean = 0
419        subslope_dist(:,1) = 1.
420    endif
421
422    ! Remove unphysical values of surface tracer
423    qsurf(:,:,1) = qsurf_read_generic(:,:)
424    where (qsurf < 0.) qsurf = 0.
425#endif
426
427do nnq = 1,nqtot  ! Why not using ini_tracer ?
428    if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
429    if (noms(nnq) == "h2o_vap") then
430        igcm_h2o_vap = nnq
431        mmol(igcm_h2o_vap)=18.
432    endif
433    if (noms(nnq) == "co2") igcm_co2 = nnq
434enddo
435r = 8.314511*1000./mugaz
436
437!------------------------
438! I   Initialization
439!    I_c Subslope parametrisation
440!------------------------
441! Define some slope statistics
442iflat = 1
443do islope = 2,nslope
444    if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
445enddo
446
447write(*,*) 'Flat slope for islope = ',iflat
448write(*,*) 'corresponding criterium = ',def_slope_mean(iflat)
449
450allocate(flag_co2flow(ngrid,nslope))
451allocate(flag_co2flow_mesh(ngrid))
452allocate(flag_h2oflow(ngrid,nslope))
453allocate(flag_h2oflow_mesh(ngrid))
454
455flag_co2flow(:,:) = 0
456flag_co2flow_mesh(:) = 0
457flag_h2oflow(:,:) = 0
458flag_h2oflow_mesh(:) = 0
459
460!------------------------
461! I   Initialization
462!    I_d READ PCM data and convert to the physical grid
463!------------------------
464! 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
465call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
466
467allocate(tsoil_ave(ngrid,nsoilmx,nslope))
468allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
469allocate(vmr_co2_gcm(ngrid,timelen))
470allocate(ps_timeseries(ngrid,timelen))
471allocate(min_co2_ice_1(ngrid,nslope))
472allocate(min_h2o_ice_1(ngrid,nslope))
473allocate(min_co2_ice_2(ngrid,nslope))
474allocate(min_h2o_ice_2(ngrid,nslope))
475allocate(tsurf_ave_yr1(ngrid,nslope))
476allocate(tsurf_ave(ngrid,nslope))
477allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
478allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
479allocate(q_co2_PEM_phys(ngrid,timelen))
480allocate(q_h2o_PEM_phys(ngrid,timelen))
481allocate(co2_ice_GCM(ngrid,nslope,timelen))
482allocate(watersurf_density_ave(ngrid,nslope))
483allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
484allocate(Tsurfave_before_saved(ngrid,nslope))
485allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
486allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
487allocate(delta_co2_adsorbded(ngrid))
488allocate(porefillingice_thickness_prev_iter(ngrid,nslope))
489allocate(delta_h2o_icetablesublim(ngrid))
490allocate(delta_h2o_adsorbded(ngrid))
491allocate(vmr_co2_pem_phys(ngrid,timelen))
492
493write(*,*) "Downloading data Y1..."
494call read_data_PCM("data_PCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
495                   tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,           &
496                   co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
497write(*,*) "Downloading data Y1 done"
498
499! 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
500write(*,*) "Downloading data Y2"
501call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
502                   tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,              &
503                   co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
504write(*,*) "Downloading data Y2 done"
505
506!------------------------
507! I   Initialization
508!    I_e Initialization of the PEM variables and soil
509!------------------------
510call end_comsoil_h_PEM
511call ini_comsoil_h_PEM(ngrid,nslope)
512call end_adsorption_h_PEM
513call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
514call end_ice_table_porefilling
515call ini_ice_table_porefilling(ngrid,nslope)
516
517if (soil_pem) then
518    call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
519    do l = 1,nsoilmx
520        tsoil_PEM(:,l,:) = tsoil_ave(:,l,:)
521        tsoil_phys_PEM_timeseries(:,l,:,:) = tsoil_GCM_timeseries(:,l,:,:)
522        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,l,:,:)
523    enddo
524    do l = nsoilmx + 1,nsoilmx_PEM
525        tsoil_PEM(:,l,:) = tsoil_ave(:,nsoilmx,:)
526        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
527    enddo
528    watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
529endif !soil_pem
530deallocate(tsoil_ave,tsoil_GCM_timeseries)
531
532!------------------------
533! I   Initialization
534!    I_f Compute tendencies & Save initial situation
535!------------------------
536allocate(tendencies_co2_ice(ngrid,nslope))
537allocate(tendencies_co2_ice_ini(ngrid,nslope))
538allocate(tendencies_h2o_ice(ngrid,nslope))
539
540! Compute the tendencies of the evolution of ice over the years
541call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice)
542tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:)
543call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice)
544
545deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2)
546
547!------------------------
548! I   Initialization
549!    I_g Save initial PCM situation
550!------------------------
551allocate(initial_co2_ice_sublim(ngrid,nslope))
552allocate(initial_co2_ice(ngrid,nslope))
553allocate(initial_h2o_ice(ngrid,nslope))
554
555! We save the places where water ice is sublimating
556! We compute the surface of water ice sublimating
557ini_surf_co2 = 0.
558ini_surf_h2o = 0.
559Total_surface = 0.
560do i = 1,ngrid
561    Total_surface = Total_surface + cell_area(i)
562    do islope = 1,nslope
563        if (tendencies_co2_ice(i,islope) < 0) then
564            initial_co2_ice_sublim(i,islope) = 1.
565            ini_surf_co2 = ini_surf_co2 + cell_area(i)*subslope_dist(i,islope)
566        else
567            initial_co2_ice_sublim(i,islope) = 0.
568        endif
569        if (qsurf(i,igcm_co2,islope) > 0) then
570            initial_co2_ice(i,islope) = 1.
571        else
572            initial_co2_ice(i,islope) = 0.
573        endif
574        if (tendencies_h2o_ice(i,islope) < 0) then
575            initial_h2o_ice(i,islope) = 1.
576            ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
577        else
578            initial_h2o_ice(i,islope) = 0.
579        endif
580    enddo
581enddo
582
583write(*,*) "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
584write(*,*) "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
585write(*,*) "Total surface of the planet=", Total_surface
586allocate(zplev_gcm(ngrid,nlayer + 1))
587
588do ig = 1,ngrid
589    zplev_gcm(ig,:) = ap(:) + bp(:)*ps_start_GCM(ig)
590enddo
591
592global_ave_press_old = 0.
593do i = 1,ngrid
594    global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface
595enddo
596
597global_ave_press_GCM = global_ave_press_old
598global_ave_press_new = global_ave_press_old
599write(*,*) "Initial global average pressure=", global_ave_press_GCM
600
601!------------------------
602! I   Initialization
603!    I_h Read the startpem.nc
604!------------------------
605call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
606              porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,          &
607              tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),                      &
608              qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave,             &
609              co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
610
611delta_h2o_icetablesublim(:) = 0.
612
613do ig = 1,ngrid
614    do islope = 1,nslope
615        qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
616    enddo
617enddo
618
619if (adsorption_pem) then
620    totmassco2_adsorbded = 0.
621    totmassh2o_adsorbded = 0.
622    do ig = 1,ngrid
623        do islope = 1,nslope
624            do l = 1,nsoilmx_PEM - 1
625                totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
626                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
627                totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
628                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
629            enddo
630        enddo
631    enddo
632
633    write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
634    write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
635endif ! adsorption
636deallocate(tsurf_ave_yr1)
637
638!------------------------
639! I   Initialization
640!    I_i Compute orbit criterion
641!------------------------
642#ifndef CPP_STD
643    call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
644#else
645    call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
646#endif
647
648if (evol_orbit_pem) then
649    call orbit_param_criterion(i_myear,year_iter_max)
650else
651    year_iter_max = Max_iter_pem
652endif
653!-------------------------- END INITIALIZATION -------------------------
654
655!-------------------------------- RUN ----------------------------------
656!------------------------
657! II  Run
658!    II_a Update pressure, ice and tracers
659!------------------------
660year_iter = 0
661
662do while (year_iter < year_iter_max .and. i_myear < n_myear)
663! II.a.1. Compute updated global pressure
664    write(*,*) "Recomputing the new pressure..."
665    do i = 1,ngrid
666        do islope = 1,nslope
667            global_ave_press_new = global_ave_press_new - g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
668        enddo
669    enddo
670
671    if (adsorption_pem) then
672        do i = 1,ngrid
673            global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
674        enddo
675    endif
676    write(*,*) 'Global average pressure old time step',global_ave_press_old
677    write(*,*) 'Global average pressure new time step',global_ave_press_new
678
679! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
680    allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
681    write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
682    do l = 1,nlayer + 1
683        do ig = 1,ngrid
684            zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
685        enddo
686    enddo
687
688! II.a.3. Surface pressure timeseries
689    write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..."
690    do ig = 1,ngrid
691        ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
692    enddo
693
694! II.a.4. New pressure levels timeseries
695    allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
696    write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
697    do l = 1,nlayer + 1
698        do ig = 1,ngrid
699            zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
700        enddo
701    enddo
702
703! II.a.5. Tracers timeseries
704    write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..."
705
706    l = 1
707    do ig = 1,ngrid
708        do t = 1,timelen
709            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))/ &
710                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
711            if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30
712            if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1.
713        enddo
714    enddo
715
716    do ig = 1,ngrid
717        do t = 1,timelen
718            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))/ &
719                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))   &
720                                + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))   &
721                                -  (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
722                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
723            if (q_co2_PEM_phys(ig,t) < 0) then
724                q_co2_PEM_phys(ig,t) = 1.e-30
725            elseif (q_co2_PEM_phys(ig,t) > 1) then
726                q_co2_PEM_phys(ig,t) = 1.
727            endif
728            mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
729            vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
730        enddo
731    enddo
732
733    deallocate(zplev_new_timeseries,zplev_old_timeseries)
734
735!------------------------
736! II  Run
737!    II_b Evolution of the ice
738!------------------------
739    write(*,*) "Evolution of h2o ice"
740    call evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,STOPPING_1_water)
741
742    write(*,*) "Evolution of co2 ice"
743    call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
744
745!------------------------
746! II  Run
747!    II_c CO2 & H2O glaciers flows
748!------------------------
749    write(*,*) "CO2 glacier flows"
750    if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
751                                               global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
752
753    write(*,*) "H2O glacier flows"
754    if (h2oglaciersflow) call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
755
756!------------------------
757! II  Run
758!    II_d Update surface and soil temperatures
759!------------------------
760! II_d.1 Update Tsurf
761    write(*,*) "Updating the new Tsurf"
762    bool_sublim = .false.
763    Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
764    do ig = 1,ngrid
765        do islope = 1,nslope
766            if (initial_co2_ice(ig,islope) > 0.5 .and. qsurf(ig,igcm_co2,islope) < 1.e-10) then !co2ice disappeared, look for closest point without co2ice
767                if (latitude_deg(ig) > 0) then
768                    do ig_loop = ig,ngrid
769                        do islope_loop = islope,iflat,-1
770                            if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
771                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
772                                bool_sublim = .true.
773                                exit
774                            endif
775                        enddo
776                        if (bool_sublim) exit
777                    enddo
778                else
779                    do ig_loop = ig,1,-1
780                        do islope_loop = islope,iflat
781                            if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
782                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
783                                bool_sublim = .true.
784                                exit
785                            endif
786                        enddo
787                        if (bool_sublim) exit
788                    enddo
789                endif
790                initial_co2_ice(ig,islope) = 0
791                if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then
792                    albedo(ig,1,islope) = albedo_h2o_frost
793                    albedo(ig,2,islope) = albedo_h2o_frost
794                else
795                    albedo(ig,1,islope) = albedodat(ig)
796                    albedo(ig,2,islope) = albedodat(ig)
797                    emis(ig,islope) = emissiv
798                endif
799            else if ((qsurf(ig,igcm_co2,islope) > 1.e-3) .and. (tendencies_co2_ice(ig,islope) > 1.e-10)) then !Put tsurf as tcond co2
800                ave = 0.
801                do t = 1,timelen
802                    if (co2_ice_GCM(ig,islope,t) > 1.e-3) then
803                        ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
804                    else
805                        ave = ave + tsurf_GCM_timeseries(ig,islope,t)
806                    endif
807                enddo
808                tsurf_ave(ig,islope) = ave/timelen
809                ! set the surface albedo to be the ice albedo
810                if (latitude_deg(ig) > 0) then
811                   icap = 1
812                else
813                   icap = 2
814                endif
815                albedo(ig,1,islope) = albedice(icap)
816                albedo(ig,2,islope) = albedice(icap)
817                emis(ig,islope) = emisice(icap)
818            endif
819        enddo
820    enddo
821
822    do t = 1,timelen
823        tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:))
824    enddo
825    ! for the start
826    do ig = 1,ngrid
827        do islope = 1,nslope
828            tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope))
829        enddo
830    enddo
831
832    if (soil_pem) then
833
834! II_d.2 Update soil temperature
835        allocate(TI_locslope(ngrid,nsoilmx_PEM))
836        allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
837        allocate(Tsurf_locslope(ngrid))
838        write(*,*)"Updating soil temperature"
839
840        ! Soil averaged
841        do islope = 1,nslope
842            TI_locslope(:,:) = TI_PEM(:,:,islope)
843            do t = 1,timelen
844                Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
845                Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
846                call soil_pem_compute(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
847                call soil_pem_compute(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
848                tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
849                do ig = 1,ngrid
850                    do isoil = 1,nsoilmx_PEM
851                        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)
852                        if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1)
853                    enddo
854                enddo
855            enddo
856        enddo
857        tsoil_PEM(:,:,:) = sum(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
858        watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
859
860        write(*,*) "Update of soil temperature done"
861
862        deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
863        write(*,*) "Compute ice table"
864
865! II_d.3 Update the ice table
866        porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:)
867        call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
868
869        call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
870
871        write(*,*) "Update soil propreties"
872
873! II_d.4 Update the soil thermal properties
874        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
875
876! II_d.5 Update the mass of the regolith adsorbded
877        if (adsorption_pem) then
878            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
879                                     qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, &
880                                     q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
881
882            totmassco2_adsorbded = 0.
883            totmassh2o_adsorbded = 0.
884            do ig = 1,ngrid
885                do islope =1, nslope
886                    do l = 1,nsoilmx_PEM - 1
887                        totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
888                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
889                        totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
890                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
891                    enddo
892                enddo
893            enddo
894            write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
895            write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
896        endif
897    endif !soil_pem
898
899!------------------------
900! II  Run
901!    II_e Outputs
902!------------------------
903    call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_ave_press_new/))
904    do islope = 1,nslope
905        write(str2(1:2),'(i2.2)') islope
906        call writediagpem(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
907        call writediagpem(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
908        call writediagpem(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
909        call writediagpem(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
910        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
911        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
912    enddo
913
914!------------------------
915! II  Run
916!    II_f Update the tendencies
917!------------------------
918    write(*,*) "Adaptation of the new co2 tendencies to the current pressure"
919    call recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),emis,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
920                               global_ave_press_GCM,global_ave_press_new)
921
922!------------------------
923! II  Run
924!    II_g Checking the stopping criterion
925!------------------------
926    call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
927
928    call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
929                            initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
930
931    year_iter = year_iter + dt_pem
932    i_myear = i_myear + dt_pem
933
934    write(*,*) "Checking all the stopping criterion."
935    if (STOPPING_water) then
936        write(*,*) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
937        criterion_stop = 1
938    endif
939    if (STOPPING_1_water) then
940        write(*,*) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
941        criterion_stop = 1
942    endif
943    if (STOPPING_co2) then
944        write(*,*) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
945        criterion_stop = 2
946    endif
947    if (STOPPING_pressure) then
948        write(*,*) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
949        criterion_stop = 3
950    endif
951    if (year_iter >= year_iter_max) then
952        write(*,*) "STOPPING because maximum number of iterations reached"
953        criterion_stop = 4
954    endif
955    if (i_myear >= n_myear) then
956        write(*,*) "STOPPING because maximum number of Martian years to be simulated reached"
957        criterion_stop = 5
958    endif
959
960    if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure)  then
961        exit
962    else
963        write(*,*) "We continue!"
964        write(*,*) "Number of iterations of the PEM: year_iter =", year_iter
965        write(*,*) "Number of simulated Martian years: i_myear =", i_myear
966    endif
967
968    global_ave_press_old = global_ave_press_new
969
970enddo  ! big time iteration loop of the pem
971!------------------------------ END RUN --------------------------------
972
973!------------------------------- OUTPUT --------------------------------
974!------------------------
975! III Output
976!    III_a Update surface value for the PCM start files
977!------------------------
978! III_a.1 Ice update (for startfi)
979
980! H2O ice
981do ig = 1,ngrid
982    if (watercaptag(ig)) then
983        watercap_sum = 0.
984        do islope = 1,nslope
985            if (qsurf(ig,igcm_h2o_ice,islope) > (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy
986                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) ! put back ancien frost
987            else
988!               2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
989                watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
990                qsurf(ig,igcm_h2o_ice,islope)=0.
991            endif
992            watercap_sum = watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
993            watercap(ig,islope) = 0.
994        enddo
995        water_reservoir(ig) = water_reservoir(ig) + watercap_sum
996    endif
997enddo
998
999do ig = 1,ngrid
1000    water_sum = 0.
1001    do islope = 1,nslope
1002        water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
1003    enddo
1004    if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
1005        if (water_sum > threshold_water_frost2perenial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
1006            watercaptag(ig) = .true.
1007            water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2perenial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
1008            do islope = 1,nslope
1009                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.)
1010            enddo
1011        endif
1012    else ! let's check that the infinite source of water disapear
1013        if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perenial) then
1014            watercaptag(ig) = .false.
1015            do islope = 1,nslope
1016                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
1017            enddo
1018            water_reservoir(ig) = 0.
1019        endif
1020    endif
1021enddo
1022
1023! CO2 ice
1024do ig = 1,ngrid
1025    do islope = 1,nslope
1026        if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then
1027            perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
1028            qsurf(ig,igcm_co2,islope)  = 0.5*qsurf(ig,igcm_co2,islope)
1029            albedo(ig,1,islope) = albedo_perenialco2
1030            albedo(ig,2,islope) = albedo_perenialco2
1031        endif
1032    enddo
1033enddo
1034
1035! III_a.2  Tsoil update (for startfi)
1036if (soil_pem) then
1037    call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
1038    tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
1039endif
1040
1041! III_a.4 Pressure (for start)
1042ps(:) = ps(:)*global_ave_press_new/global_ave_press_GCM
1043ps_start_GCM(:) = ps_start_GCM(:)*global_ave_press_new/global_ave_press_GCM
1044
1045! III_a.5 Tracer (for start)
1046allocate(zplev_new(ngrid,nlayer + 1))
1047
1048do l = 1,nlayer + 1
1049    zplev_new(:,l) = ap(l) + bp(l)*ps_start_GCM(:)
1050enddo
1051
1052do nnq = 1,nqtot
1053    if (noms(nnq) /= "co2") then
1054        do l = 1,llm - 1
1055            do ig = 1,ngrid
1056                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
1057            enddo
1058            q(:,llm,nnq) = q(:,llm - 1,nnq)
1059        enddo
1060    else
1061        do l = 1,llm - 1
1062            do ig = 1,ngrid
1063                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
1064                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
1065            enddo
1066            q(:,llm,nnq) = q(:,llm - 1,nnq)
1067        enddo
1068    endif
1069enddo
1070
1071! Conserving the tracers mass for PCM start files
1072do nnq = 1,nqtot
1073    do ig = 1,ngrid
1074        do l = 1,llm - 1
1075            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
1076                extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1077                q(ig,l,nnq) = 1.
1078                q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1079                write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1080           endif
1081            if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1082        enddo
1083    enddo
1084enddo
1085
1086if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
1087
1088!------------------------
1089! III Output
1090!    III_b Write restart_evol.nc and restartfi_evol.nc
1091!------------------------
1092! III_b.1 Write restart_evol.nc
1093ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1094pday = day_ini
1095ztime_fin = time_phys
1096
1097allocate(p(ip1jmp1,nlayer + 1))
1098#ifndef CPP_1D
1099    call pression (ip1jmp1,ap,bp,ps,p)
1100    call massdair(p,masse)
1101    call dynredem0("restart_evol.nc",day_ini,phis)
1102    call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps)
1103    write(*,*) "restart_evol.nc has been written"
1104#else
1105    call writerestart1D('restart1D_evol.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1106    write(*,*) "restart1D_evol.txt has been written"
1107#endif
1108
1109! III_b.2 Write restartfi_evol.nc
1110#ifndef CPP_STD
1111    call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
1112                  nlayer,nq,ptimestep,pday,0.,cell_area,albedodat,      &
1113                  inertiedat,def_slope,subslope_dist)
1114
1115    call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1116                  ptimestep,ztime_fin,tsurf,tsoil,inertiesoil,        &
1117                  albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1118                  wstar,watercap,perenial_co2ice)
1119#else
1120    call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
1121                  nlayer,nq,ptimestep,pday,time_phys,cell_area,         &
1122                  albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1123
1124    call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,  &
1125                  ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1126                  cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab,   &
1127                  tsea_ice,sea_ice)
1128#endif
1129write(*,*) "restartfi_evol.nc has been written"
1130
1131!------------------------
1132! III Output
1133!    III_c Write restartpem.nc
1134!------------------------
1135call pemdem0("restartpem.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
1136             float(day_ini),0.,nslope,def_slope,subslope_dist)
1137call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
1138             TI_PEM, porefillingice_depth,porefillingice_thickness,         &
1139             co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
1140write(*,*) "restartpem.nc has been written"
1141
1142call info_PEM(year_iter,criterion_stop,i_myear,n_myear)
1143
1144write(*,*) "The PEM has run for", year_iter, "Martian years."
1145write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1146write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
1147write(*,*) "LL & RV & JBC: so far, so good!"
1148
1149deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
1150deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved)
1151deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
1152deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
1153!----------------------------- END OUTPUT ------------------------------
1154
1155END PROGRAM pem
Note: See TracBrowser for help on using the repository browser.