Context Navigation

pem.F90 @ 2857

Last change on this file since 2857 was 2857, checked in by llange, 2 years ago
PEM Fixing bug when updating tracer mass mixing ratio (ccn_number and dust_number were previously set to 1) Tiny edits in evol_h2o_ice and criterion_ice_stop_mod_water LL
File size: 58.3 KB

Rev	Line
[2779]	1
	2	!------------------------
	3
	4	! I Initialisation
[2835]	5	! I_a READ run.def
	6	! I_b READ of start_evol.nc and starfi_evol.nc
	7	! I_c Subslope parametrisation
	8	! I_d READ GCM data and convert to the physical grid
	9	! I_e Initialisation of the PEM variable and soil
	10	! I_f Compute tendencies & Save initial situation
	11	! I_g Save initial PCM situation
	12	! I_h Read the PEMstart
	13	! I_i Compute orbit criterion
[2779]	14
	15	! II Run
[2794]	16	! II_a update pressure, ice and tracers
[2835]	17	! II_b Evolution of the ice
	18	! II_c CO2 glaciers flows
	19	! II_d Update surface and soil temperatures
	20	! II_e Update the tendencies
	21	! II_f Checking the stopping criterion
[2779]	22
	23	! III Output
[2835]	24	! III_a Update surface value for the PCM start files
[2794]	25	! III_b Write start and starfi.nc
	26	! III_c Write start_pem
[2779]	27
	28	!------------------------
	29
	30	PROGRAM pem
	31
	32	!module needed for INITIALISATION
[2842]	33	#ifndef CPP_STD
[2794]	34	use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
[2779]	35	use surfdat_h, only: tsurf, co2ice, emis,&
[2842]	36	qsurf,watercap, ini_surfdat_h, &
[2779]	37	albedodat, zmea, zstd, zsig, zgam, zthe, &
[2794]	38	hmons, summit, base,albedo_h2o_frost, &
[2842]	39	frost_albedo_threshold,emissiv
[2835]	40	use dimradmars_mod, only: totcloudfrac, albedo
[2842]	41	use dust_param_mod, only: tauscaling
	42	use tracer_mod, only: noms,igcm_h2o_ice ! tracer names
	43	#else
	44	use comsoil_h, only: nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
	45	use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe, &
	46	emissiv
	47	use tracer_h, only: noms,igcm_h2o_ice,igcm_co2_ice ! tracer names
	48	use phys_state_var_mod
	49	#endif
	50	use phyetat0_mod, only: phyetat0
	51	use phyredem, only: physdem0, physdem1
[2835]	52	use turb_mod, only: q2, wstar
[2779]	53	use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
	54	nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
	55	nf90_inquire_dimension,nf90_close
	56
	57	! For phyredem :
	58	USE control_mod, ONLY: iphysiq, day_step,nsplit_phys
[2842]	59	USE iniphysiq_mod, ONLY: iniphysiq
	60	USE logic_mod, ONLY: iflag_phys
	61	#ifndef CPP_STD
[2779]	62	use mod_phys_lmdz_para, only: is_parallel, is_sequential, &
	63	is_mpi_root, is_omp_root, &
	64	is_master
[2842]	65	use planete_h, only: aphelie, periheli, year_day, peri_day, &
	66	obliquit
	67	#else
	68	! USE comcstfi_mod, ONLY: rad,g,r,cpp,pi
	69	! USE inifis_mod, ONLY: inifis
	70	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
	71	obliquit
	72	#endif
[2779]	73	USE mod_const_mpi, ONLY: COMM_LMDZ
[2794]	74	USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, &
[2779]	75	subslope_dist,co2ice_slope, &
	76	tsurf_slope,tsoil_slope,fluxgrd_slope,&
	77	fluxrad_sky_slope,sky_slope,callsubslope,&
	78	co2iceflow, beta_slope, capcal_slope,&
	79	albedo_slope,emiss_slope,qsurf_slope,&
[2849]	80	iflat,major_slope,ini_comslope_h
[2842]	81	use time_phylmdz_mod, only: daysec,dtphys
	82	USE comconst_mod, ONLY: rad,g,r,cpp,pi
	83	USE infotrac
[2835]	84	USE geometry_mod, only: latitude_deg
[2779]	85
[2794]	86	use pemredem, only: pemdem1
[2856]	87	use co2glaciers_mod,only: co2glaciers_evol
[2794]	88	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
[2855]	89	use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
[2835]	90	TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM, & ! soil thermal inertia
	91	tsoil_PEM, mlayer_PEM,layer_PEM, & !number of subsurface layers, soil mid layer depths
	92	fluxgeo, co2_adsorbded_phys ! geothermal flux, mass of co2 in the regolith
[2794]	93	use adsorption_mod, only : regolith_co2adsorption
[2835]	94
	95	!!! For orbit parameters
	96	USE temps_mod_evol, ONLY: dt_pem, evol_orbit_pem, Max_iter_pem
	97	use orbit_param_criterion_mod, only : orbit_param_criterion
	98	use recomp_orb_param_mod, only: recomp_orb_param
	99
	100
[2779]	101	IMPLICIT NONE
	102
	103	include "dimensions.h"
	104	include "paramet.h"
[2794]	105
	106	INTEGER ngridmx
	107	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
	108
[2779]	109	include "comdissnew.h"
	110	include "comgeom.h"
	111	include "iniprint.h"
	112	! Same variable's name as in the GCM
	113
[2794]	114	INTEGER :: ngrid !Number of physical grid points
	115	INTEGER :: nlayer !Number of vertical layer
	116	INTEGER :: nq !Number of tracer
	117	INTEGER :: day_ini !First day of the simulation
	118	REAL :: pday !Physical day
	119	REAL :: time_phys !Same as GCM
	120	REAL :: ptimestep !Same as GCM
	121	REAL :: ztime_fin !Same as GCM
[2779]	122
	123	! Variable for reading start.nc
[2835]	124	character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM
[2779]	125	! variables dynamiques
	126	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
	127	REAL teta(ip1jmp1,llm) ! temperature potentielle
	128	REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
	129	REAL ps(ip1jmp1) ! pression au sol
[2794]	130	REAL, dimension(:),allocatable :: ps_phys !(ngrid) ! pression au sol
	131	REAL, dimension(:,:),allocatable :: ps_phys_timeseries !(ngrid x timelen) ! pression au sol instantannées
	132	REAL, dimension(:,:),allocatable :: ps_phys_timeseries_yr1 !(ngrid x timelen) ! pression au sol instantannées for the first year of the gcm
	133
[2779]	134	REAL masse(ip1jmp1,llm) ! masse d'air
	135	REAL phis(ip1jmp1) ! geopotentiel au sol
	136	REAL time_0
	137
	138	! Variable for reading starfi.nc
	139
[2835]	140	character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
[2779]	141	integer :: ncid, varid,status !Variable for handling opening of files
	142	integer :: phydimid, subdimid, nlayerdimid, nqdimid !Variable ID for Netcdf files
	143	integer :: lonvarid, latvarid, areavarid,sdvarid !Variable ID for Netcdf files
[2794]	144	integer :: apvarid,bpvarid !Variable ID for Netcdf files
[2779]	145
	146	! Variable for reading starfi.nc and writting restartfi.nc
	147
	148	REAL, dimension(:),allocatable :: longitude !Longitude read in FILE_NAME and written in restartfi
	149	REAL, dimension(:),allocatable :: latitude !Latitude read in FILE_NAME and written in restartfi
[2794]	150	REAL, dimension(:),allocatable :: ap !Coefficient ap read in FILE_NAME_start and written in restart
	151	REAL, dimension(:),allocatable :: bp !Coefficient bp read in FILE_NAME_start and written in restart
[2779]	152	REAL, dimension(:),allocatable :: cell_area !Cell_area read in FILE_NAME and written in restartfi
	153	REAL :: Total_surface !Total surface of the planet
	154
	155	! Variable for h2o_ice evolution
	156
	157	REAL , dimension(:,:), allocatable :: tendencies_h2o_ice ! LON x LAT field : Tendency of evolution of perenial ice over a year
	158	REAL, dimension(:),allocatable :: tendencies_h2o_ice_phys ! physical point field : Tendency of evolution of perenial ice over a year
	159
	160	REAL , dimension(:,:), allocatable :: tendencies_co2_ice ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
	161	REAL, dimension(:),allocatable :: tendencies_co2_ice_phys ! physical point field : Tendency of evolution of perenial co2 ice over a year
	162
	163	REAL :: ini_surf ! Initial surface of sublimating water ice
[2794]	164	REAL :: ini_surf_h2o ! Initial surface of sublimating water ice
[2779]	165	REAL, dimension(:),allocatable :: initial_h2o_ice ! physical point field : Logical array indicating sublimating point
	166
	167	REAL :: ini_surf_co2 ! Initial surface of sublimating co2 ice
	168	REAL, dimension(:),allocatable :: initial_co2_ice ! physical point field : Logical array indicating sublimating point of co2 ice
	169
	170	REAL , dimension(:,:), allocatable :: min_h2o_ice_s_1 ! LON x LAT field : minimum of water ice at each point for the first year
	171	REAL , dimension(:,:), allocatable :: min_h2o_ice_s_2 ! LON x LAT field : minimum of water ice at each point for the second year
	172
	173	REAL , dimension(:,:), allocatable :: min_co2_ice_s_1 ! LON x LAT field : minimum of water ice at each point for the first year
	174	REAL , dimension(:,:), allocatable :: min_co2_ice_s_2 ! LON x LAT field : minimum of water ice at each point for the second year
	175
[2855]	176	REAL :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa]
	177	REAL :: global_ave_press_old ! constant: Global average pressure of initial/previous time step
	178	REAL :: global_ave_press_new ! constant: Global average pressure of current time step
[2779]	179
[2855]	180	REAL , dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2]
	181	REAL , dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2]
	182	REAL , dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
	183	REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
[2779]	184
	185	LOGICAL :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
	186	LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
[2835]	187	LOGICAL :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
	188	LOGICAL :: STOPPING_1_co2 ! Logical : is there still water ice to sublimate?
[2779]	189
[2855]	190	REAL, dimension(:,:,:),allocatable :: q_co2_GCM ! Lon x Lat x Time : mass mixing ratio of co2 in the first layer [kg/kg]
[2779]	191
[2855]	192	real,save :: m_co2, m_noco2, A , B, mmean ! Molar mass of co2, no co2 (Ar, ...), intermediate A, B for computations, mean molar mass of the layer [mol/kg]
	193	real ,allocatable :: vmr_co2_gcm_phys(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
	194	real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM
	195	real ,allocatable :: q_h2o_GCM_phys(:,:) ! Physics x Times h2o mass mixing ratio in the first layer from the GCM [kg/kg]
	196	real ,allocatable :: q_co2_GCM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer from the GCM [kg/kg]
	197	real ,allocatable :: q_co2_PEM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM [kg/kg]
	198	REAL, ALLOCATABLE :: ps_GCM(:,:,:) ! Lon x Lat x Times: surface pressure from the GCM [Pa]
	199	REAL, ALLOCATABLE :: ps_GCM_yr1(:,:,:) ! Lon x Lat x Times: surface pressure from the 1st year of the GCM [Pa]
	200	REAL, ALLOCATABLE :: vmr_co2_gcm(:,:,:) ! Lon x Lat x Times: co2 volumemixing ratio retrieve from the gcm [m^3/m^3]
	201	REAL, ALLOCATABLE :: q_h2o_GCM(:,:,:) ! Lon x Lat x Times: h2o volume mixing ratio retrieved from the GCM
	202	REAL ,allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM
	203	integer :: timelen ! # time samples
	204	REAL :: ave ! intermediate varibale to compute average
[2779]	205
[2855]	206	REAL, ALLOCATABLE :: p(:,:) ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
	207	REAL :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1
	208	REAL :: beta_clap_co2 = 3182.48 ! clapeyron's law for CO2
	209	REAL :: alpha_clap_co2 = 23.3494 ! Clapeyron's law for CO2
[2779]	210
[2855]	211
[2779]	212	!!!!!!!!!!!!!!!!!!!!!!!! SLOPE
[2855]	213	REAL ,allocatable :: watercap_slope(:,:) ! Physics x Nslope: watercap per slope
	214	REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1 ! LON x LAT field : minimum of co2 ice at each point for the first year [kg/m^2]
	215	REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2 ! LON x LAT field : minimum of co2 ice at each point for the second year [kg/m^2]
	216	REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_1 ! LON x LAT field : minimum of water ice at each point for the first year [kg/m^2]
	217	REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_2 ! LON x LAT field : minimum of water ice at each point for the second year [kg/m^2]
	218	REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM [kg/m^2]
	219	REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2]
[2835]	220	REAL, dimension(:,:),allocatable :: initial_co2_ice_sublim_slope ! physical point field : Logical array indicating sublimating point of co2 ice
[2855]	221	REAL, dimension(:,:),allocatable :: initial_h2o_ice_slope ! physical point field : Logical array indicating if there is water ice at initial state
	222	REAL, dimension(:,:),allocatable :: initial_co2_ice_slope ! physical point field : Logical array indicating if there is co2 ice at initial state
	223	REAL , dimension(:,:,:), allocatable :: tendencies_co2_ice_slope ! LON x LAT x nslope field : Tendency of evolution of perenial co2 ice over a year
	224	REAL , dimension(:,:,:), allocatable :: tendencies_h2o_ice_slope ! LON x LAT x slope field : Tendency of evolution of perenial water ice over a year
	225	REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
	226	REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
[2856]	227	REAL, dimension(:,:),allocatable :: tendencies_h2o_ice_phys_slope ! physical pointx slope field : Tendency of evolution of perenial co2 ice
[2779]	228	REAL , dimension(:,:), allocatable :: flag_co2flow(:,:) !(ngrid,nslope) ! To flag where there is a glacier flow
	229	REAL , dimension(:), allocatable :: flag_co2flow_mesh(:) !(ngrid) ! To flag where there is a glacier flow
	230
[2794]	231	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL
[2779]	232
[2835]	233	REAL, ALLOCATABLE :: tsurf_ave(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
[2855]	234	REAL, ALLOCATABLE :: tsurf_ave_phys(:,:) ! Physic x LAT x SLOPE field : Averaged Surface Temperature [K]
[2835]	235	REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
[2855]	236	REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature during 1st year of the GCM [K]
	237	REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) ! Physics x SLOPE field : Averaged Soil Temperature during 1st year [K]
[2835]	238	REAL, ALLOCATABLE :: TI_GCM(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Thermal Inertia [SI]
[2855]	239	REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:) ! LON X LAT x SLOPE XTULES field : Surface Temperature in timeseries [K]
	240	REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) ! Physic x SLOPE XTULES field : NOn averaged Surf Temperature [K]
[2794]	241
	242	REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
[2835]	243	REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
	244	REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
[2855]	245	REAL, ALLOCATABLE :: tsurf_ave_phys_yr1(:,:) ! Physic SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
	246	REAL, ALLOCATABLE :: inertiesoil(:,:) !Physic x Depth Thermal inertia of the mesh for restart [SI]
[2794]	247
[2855]	248	REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:) ! Physic x Depth x Slope Averaged GCM Thermal Inertia per slope [SI]
	249	REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Same but for the start
[2794]	250
[2855]	251	REAL,ALLOCATABLE :: ice_depth(:,:) ! Physic x SLope: Ice table depth [m]
	252	REAL,ALLOCATABLE :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI]
	253	REAL,ALLOCATABLE :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
	254	REAL,ALLOCATABLE :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
	255	REAL,ALLOCATABLE :: alph_locslope(:,:) ! Physic x Soil: Intermediate to compute Tsoil [1]
	256	REAL,ALLOCATABLE :: beta_locslope(:,:) ! Physic x Soil : Intermediate tocompute Tsoil [K]
	257	REAL,ALLOCATABLE :: watersurf_density_timeseries(:,:,:,:) ! Lon x Lat x Slope x Times: water surface density, time series [kg/m^3]
	258	REAL,ALLOCATABLE :: watersoil_density_timeseries(:,:,:,:,:) ! Lon x Lat x Soil x Slope x Times water soil density, time series [kg /m^3]
	259	REAL,ALLOCATABLE :: watersurf_density_phys_timeseries(:,:,:) ! Physic x Slope x Times, water surface density, time series [kg/m^3]
	260	REAL,ALLOCATABLE :: watersurf_density_phys_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
	261	REAL,ALLOCATABLE :: watersoil_density_phys_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
	262	REAL,ALLOCATABLE :: watersoil_density_phys_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
	263	REAL,ALLOCATABLE :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K]
	264	REAL, ALLOCATABLE :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
	265	REAL :: totmass_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
	266	REAL :: alpha_clap_h2o = -6143.7 ! coeffcient to compute psat, from Murphie et Kood 2005 [K]
	267	REAL :: beta_clap_h2o = 28.9074 ! coefficient to compute psat, from Murphie et Kood 2005 [1]
	268	LOGICAL :: bool_sublim ! logical to check if there is sublimation or not
[2794]	269
[2855]	270	!! Some parameters for the PEM run
	271	REAL, PARAMETER :: year_step = 1 ! timestep for the pem
	272	INTEGER :: year_iter ! number of iteration
	273	INTEGER :: year_iter_max ! maximum number of iterations before stopping
	274	REAL :: timestep ! timestep [s]
[2842]	275	#ifdef CPP_STD
	276	! INTEGER :: nsplit_phys=1
	277	! LOGICAL :: iflag_phys=.true.
[2855]	278	REAL :: frost_albedo_threshold=0.05 ! frost albedo threeshold to convert fresh frost to old ice
	279	REAL :: albedo_h2o_frost ! albedo of h2o frost
	280	REAL,ALLOCATABLE :: co2ice(:) ! Physics: co2 ice mesh averaged [kg/m^2]
[2842]	281	#endif
	282
[2779]	283	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	284
	285	! Loop variable
[2855]	286	INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil
[2842]	287	#ifndef CPP_STD
[2779]	288	! Parallel variables
	289	is_sequential=.true.
	290	is_parallel=.false.
	291	is_mpi_root=.true.
	292	is_omp_root=.true.
	293	is_master=.true.
[2842]	294	#endif
[2779]	295
	296	day_ini=0 !test
	297	time_phys=0. !test
	298
[2835]	299	! Some constants
[2794]	300
[2835]	301	ngrid=ngridmx
	302	nlayer=llm
	303
[2794]	304	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
	305	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
	306	A =(1/m_co2 - 1/m_noco2)
	307	B=1/m_noco2
	308
[2835]	309	year_day=669
	310	daysec=88775.
	311	dtphys=0
	312	timestep=year_day*daysec/year_step
[2794]	313
[2779]	314	!------------------------
	315
	316	! I Initialisation
	317	! I_a READ run.def
	318
	319	!------------------------
	320
	321	!----------------------------READ run.def ---------------------
	322	CALL conf_gcm( 99, .TRUE. )
	323
[2835]	324	call infotrac_init
	325	nq=nqtot
[2794]	326
[2779]	327	!------------------------
	328
	329	! I Initialisation
	330	! I_a READ run.def
[2835]	331	! I_b READ of start_evol.nc and starfi_evol.nc
[2779]	332
[2835]	333	!------------------------
[2779]	334
[2835]	335	!----------------------------Initialisation : READ some constant of startfi_evol.nc ---------------------
[2779]	336
[2835]	337	! In the gcm, these values are given to the physic by the dynamic.
	338	! Here we simply read them in the startfi_evol.nc file
[2779]	339	status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
	340
	341	allocate(longitude(ngrid))
	342	allocate(latitude(ngrid))
	343	allocate(cell_area(ngrid))
	344
	345	status = nf90_inq_varid(ncid, "longitude", lonvarid)
	346	status = nf90_get_var(ncid, lonvarid, longitude)
	347
	348	status = nf90_inq_varid(ncid, "latitude", latvarid)
	349	status = nf90_get_var(ncid, latvarid, latitude)
	350
	351	status = nf90_inq_varid(ncid, "area", areavarid)
	352	status = nf90_get_var(ncid, areavarid, cell_area)
	353
	354	call ini_comsoil_h(ngrid)
	355
	356	status = nf90_inq_varid(ncid, "soildepth", sdvarid)
	357	status = nf90_get_var(ncid, sdvarid, mlayer)
	358
	359	status =nf90_close(ncid)
	360
	361	!----------------------------READ start.nc ---------------------
	362
	363	allocate(q(ip1jmp1,llm,nqtot))
	364	CALL dynetat0(FILE_NAME_start,vcov,ucov, &
	365	teta,q,masse,ps,phis, time_0)
	366
[2794]	367	CALL iniconst !new
[2779]	368	CALL inigeom
[2794]	369	allocate(ap(nlayer+1))
	370	allocate(bp(nlayer+1))
	371	status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid)
	372	status = nf90_inq_varid(ncid, "ap", apvarid)
	373	status = nf90_get_var(ncid, apvarid, ap)
	374	status = nf90_inq_varid(ncid, "bp", bpvarid)
	375	status = nf90_get_var(ncid, bpvarid, bp)
	376	status =nf90_close(ncid)
[2779]	377
[2835]	378	CALL iniphysiq(iim,jjm,llm, &
[2779]	379	(jjm-1)*iim+2,comm_lmdz, &
	380	daysec,day_ini,dtphys/nsplit_phys, &
	381	rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
	382	iflag_phys)
	383
[2835]	384	!----------------------------READ startfi.nc ---------------------
[2779]	385
	386	! First we read the initial state (starfi.nc)
	387
[2835]	388	allocate(watercap_slope(ngrid,nslope))
	389	allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
[2794]	390	allocate(inertiesoil(ngrid,nsoilmx))
[2779]	391
[2842]	392	#ifndef CPP_STD
[2835]	393	CALL phyetat0 (FILE_NAME,0,0, &
[2779]	394	nsoilmx,ngrid,nlayer,nq, &
	395	day_ini,time_phys, &
	396	tsurf,tsoil,albedo,emis, &
	397	q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar, &
[2794]	398	watercap,inertiesoil,nslope,tsurf_slope, &
[2779]	399	tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
[2849]	400	subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM_start, &
[2779]	401	qsurf_slope,watercap_slope)
	402
[2835]	403	if(soil_pem) then
	404	deallocate(TI_GCM_start) !not used then
	405	endif
[2779]	406
[2835]	407	! Remove unphysical values of surface tracer
[2794]	408	DO i=1,ngrid
[2835]	409	DO nnq=1,nqtot
	410	DO islope=1,nslope
	411	if(qsurf_slope(i,nnq,islope).LT.0) then
	412	qsurf_slope(i,nnq,islope)=0.
	413	endif
[2794]	414	enddo
	415	enddo
[2835]	416	enddo
[2794]	417
[2842]	418	#else
	419	call phys_state_var_init(nq)
	420	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
	421	call initracer(ngrid,nq)
	422	call iniaerosol()
	423	call phyetat0(.true., &
	424	ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq, &
	425	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
	426	cloudfrac,totcloudfrac,hice, &
	427	rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
	428	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
	429
	430	call ini_comslope_h(ngrid,nsoilmx,nq)
	431
	432	allocate(co2ice(ngrid))
	433	co2ice(:)=qsurf(:,igcm_co2_ice)
	434	co2ice_slope(:,1)=co2ice(:)
	435	tsurf_slope(:,1)=tsurf(:)
	436
	437	if (nslope.eq.1) then
	438	def_slope(1) = 0
	439	def_slope(2) = 0
	440	def_slope_mean=0
	441	subslope_dist(:,1) = 1.
	442	endif
	443
	444	! Remove unphysical values of surface tracer
	445	DO i=1,ngrid
	446	DO nnq=1,nqtot
	447	qsurf_slope(i,nnq,1)=qsurf(i,nnq)
	448	if(qsurf(i,nnq).LT.0) then
	449	qsurf(i,nnq)=0.
	450	endif
	451	enddo
	452	enddo
	453	#endif
	454
	455	DO nnq=1,nqtot
	456	if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
	457	ENDDO
	458
[2835]	459	!------------------------
[2794]	460
[2835]	461	! I Initialisation
	462	! I_a READ run.def
	463	! I_b READ of start_evol.nc and starfi_evol.nc
	464	! I_c Subslope parametrisation
[2794]	465
[2835]	466	!------------------------
[2794]	467
[2835]	468	!----------------------------Subslope parametrisation definition ---------------------
[2794]	469
[2835]	470	! Define some slope statistics
	471	iflat=1
	472	DO islope=2,nslope
	473	IF(abs(def_slope_mean(islope)).lt. &
	474	abs(def_slope_mean(iflat))) THEN
	475	iflat = islope
	476	ENDIF
	477	ENDDO
	478
	479	PRINT*,'Flat slope for islope = ',iflat
	480	PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
	481
[2812]	482
[2835]	483	allocate(flag_co2flow(ngrid,nslope))
	484	allocate(flag_co2flow_mesh(ngrid))
[2812]	485
[2779]	486	flag_co2flow(:,:) = 0.
	487	flag_co2flow_mesh(:) = 0.
	488
[2842]	489
[2835]	490	!---------------------------- READ GCM data ---------------------
[2779]	491
[2794]	492	! I Initialisation
[2835]	493	! I_a READ run.def
	494	! I_b READ of start_evol.nc and starfi_evol.nc
	495	! I_c Subslope parametrisation
	496	! I_d READ GCM data and convert to the physical grid
[2779]	497
[2794]	498	!------------------------
	499
	500	! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the GCM run, saving only the minimum value
	501
[2835]	502	call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
[2794]	503
[2779]	504	allocate(min_h2o_ice_s_1(iim+1,jjm+1))
	505	allocate(min_co2_ice_s_1(iim+1,jjm+1))
[2794]	506	allocate(vmr_co2_gcm(iim+1,jjm+1,timelen))
	507	allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
	508	allocate(q_co2_GCM(iim+1,jjm+1,timelen))
	509	allocate(ps_GCM(iim+1,jjm+1,timelen))
	510	allocate(ps_GCM_yr1(iim+1,jjm+1,timelen))
[2779]	511	allocate(min_co2_ice_slope_1(iim+1,jjm+1,nslope))
	512	allocate(min_h2o_ice_slope_1(iim+1,jjm+1,nslope))
[2794]	513	allocate(tsurf_ave(iim+1,jjm+1,nslope))
	514	allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope))
	515	allocate(tsurf_ave_phys(ngrid,nslope))
	516	allocate(tsurf_ave_phys_yr1(ngrid,nslope))
	517	allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen))
	518	allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen))
	519	allocate(co2_ice_GCM_phys_slope(ngrid,nslope,timelen))
	520	allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen))
	521	allocate(Tsurfave_before_saved(ngrid,nslope))
[2835]	522	allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
	523	allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
	524	allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
	525	allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
	526	allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
	527	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
[2794]	528	allocate(delta_co2_adsorbded(ngrid))
[2849]	529	allocate(watersurf_density_timeseries(iim+1,jjm+1,nslope,timelen))
	530	allocate(watersoil_density_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
	531	allocate(watersurf_density_phys_timeseries(ngrid,nslope,timelen))
	532	allocate(watersurf_density_phys_ave(ngrid,nslope))
	533	allocate(watersoil_density_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
	534	allocate(watersoil_density_phys_PEM_ave(ngrid,nsoilmx_PEM,nslope))
[2835]	535	print *, "Downloading data Y1..."
[2779]	536
[2855]	537	call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm, min_h2o_ice_s_1,min_co2_ice_s_1,vmr_co2_gcm,ps_GCM_yr1,min_co2_ice_slope_1,min_h2o_ice_slope_1,&
[2849]	538	nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &
	539	watersurf_density_timeseries,watersoil_density_timeseries)
[2794]	540
	541	! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value
	542
[2835]	543	print *, "Downloading data Y1 done"
	544
[2779]	545	allocate(min_h2o_ice_s_2(iim+1,jjm+1))
	546	allocate(min_co2_ice_s_2(iim+1,jjm+1))
	547	allocate(min_co2_ice_slope_2(iim+1,jjm+1,nslope))
	548	allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
	549
[2835]	550	print *, "Downloading data Y2"
	551
[2855]	552	call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm ,min_h2o_ice_s_2,min_co2_ice_s_2,vmr_co2_gcm,ps_GCM,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
[2849]	553	nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &
	554	watersurf_density_timeseries,watersoil_density_timeseries)
[2779]	555
[2835]	556	print *, "Downloading data Y2 done"
[2779]	557
[2794]	558	! The variables in the dynamic grid are transfered to the physical grid
	559
[2779]	560	allocate(vmr_co2_gcm_phys(ngrid,timelen))
[2794]	561	allocate(vmr_co2_pem_phys(ngrid,timelen))
	562	allocate(q_h2o_GCM_phys(ngrid,timelen))
	563	allocate(q_h2o_PEM_phys(ngrid,timelen))
	564	allocate(q_co2_GCM_phys(ngrid,timelen))
	565	allocate(q_co2_PEM_phys(ngrid,timelen))
[2835]	566	allocate(ps_phys(ngrid))
	567	allocate(ps_phys_timeseries(ngrid,timelen))
	568	allocate(ps_phys_timeseries_yr1(ngrid,timelen))
[2779]	569
[2794]	570	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys)
	571	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys)
	572	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys)
[2835]	573	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
	574	call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
	575	call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
	576	CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
	577	CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
[2779]	578
[2794]	579	deallocate(vmr_co2_gcm)
	580	deallocate(q_h2o_GCM)
	581	deallocate(q_co2_GCM)
[2835]	582	deallocate(ps_GCM)
	583	deallocate(ps_GCM_yr1)
	584	deallocate(tsurf_ave)
	585	deallocate(tsurf_ave_yr1)
[2794]	586
[2835]	587	q_co2_PEM_phys(:,:)= q_co2_GCM_phys(:,:)
	588	q_h2o_PEM_phys(:,:)= q_h2o_GCM_phys(:,:)
[2794]	589
[2835]	590	!------------------------
[2794]	591
[2835]	592	! I Initialisation
	593	! I_a READ run.def
	594	! I_b READ of start_evol.nc and starfi_evol.nc
	595	! I_c Subslope parametrisation
	596	! I_d READ GCM data and convert to the physical grid
	597	! I_e Initialisation of the PEM variable and soil
[2794]	598
[2835]	599	!------------------------
[2794]	600
[2835]	601	!---------------------------- Initialisation of the PEM soil and values ---------------------
[2794]	602
	603	call end_comsoil_h_PEM
	604	call ini_comsoil_h_PEM(ngrid,nslope)
	605
	606	allocate(ice_depth(ngrid,nslope))
[2849]	607	ice_depth(:,:) = 0.
[2794]	608	allocate(TI_GCM_phys(ngrid,nsoilmx,nslope))
	609
	610	DO islope = 1,nslope
[2835]	611	if(soil_pem) then
[2794]	612	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope))
[2835]	613	endif !soil_pem
	614	DO t=1,timelen
	615	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
	616	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
	617	enddo
[2794]	618	ENDDO
	619
[2835]	620	deallocate(co2_ice_GCM_slope)
	621	deallocate(TI_GCM)
[2849]	622	deallocate(tsurf_GCM_timeseries)
[2794]	623
[2835]	624
[2849]	625	if(soil_pem) then
[2794]	626	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
	627	DO islope = 1,nslope
[2849]	628	DO t=1,timelen
	629	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersurf_density_timeseries(:,:,islope,t),watersurf_density_phys_timeseries(:,islope,t))
	630	ENDDO
[2794]	631	DO l=1,nsoilmx
	632	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
[2835]	633	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
	634	DO t=1,timelen
	635	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
[2849]	636	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersoil_density_timeseries(:,:,l,islope,t),watersoil_density_phys_PEM_timeseries(:,l,islope,t))
[2835]	637	ENDDO
[2849]	638
[2794]	639	ENDDO
	640	DO l=nsoilmx+1,nsoilmx_PEM
	641	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
	642	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
[2849]	643	DO t=1,timelen
	644	watersoil_density_phys_PEM_timeseries(:,l,islope,t) = watersoil_density_phys_PEM_timeseries(:,nsoilmx,islope,t)
	645	ENDDO
[2794]	646	ENDDO
	647	ENDDO
[2849]	648	watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
	649	watersurf_density_phys_ave(:,:) = SUM(watersurf_density_phys_timeseries(:,:,:),3)/timelen
	650	deallocate(watersurf_density_timeseries)
	651	deallocate(watersurf_density_phys_timeseries)
	652	deallocate(watersoil_density_timeseries)
[2835]	653	deallocate(tsoil_ave_yr1)
	654	deallocate(tsoil_ave)
[2794]	655	deallocate(tsoil_GCM_timeseries)
[2835]	656	endif !soil_pem
[2794]	657
[2779]	658	!------------------------
	659
[2794]	660	! I Initialisation
[2835]	661	! I_a READ run.def
	662	! I_b READ of start_evol.nc and starfi_evol.nc
	663	! I_c Subslope parametrisation
	664	! I_d READ GCM data and convert to the physical grid
	665	! I_e Initialisation of the PEM variable and soil
	666	! I_f Compute tendencies & Save initial situation
[2779]	667
[2835]	668	!----- Compute tendencies from the PCM run
[2794]	669
[2779]	670	allocate(tendencies_h2o_ice(iim+1,jjm+1))
	671	allocate(tendencies_h2o_ice_phys(ngrid))
	672	allocate(tendencies_co2_ice(iim+1,jjm+1))
	673	allocate(tendencies_co2_ice_phys(ngrid))
	674	allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
	675	allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
[2794]	676	allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope))
[2779]	677	allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope))
	678	allocate(tendencies_h2o_ice_phys_slope(ngrid,nslope))
	679
[2835]	680	! Compute the tendencies of the evolution of ice over the years
[2779]	681
	682	call compute_tendencies(tendencies_h2o_ice,min_h2o_ice_s_1,&
	683	min_h2o_ice_s_2,iim,jjm,ngrid,tendencies_h2o_ice_phys)
	684
	685	call compute_tendencies(tendencies_co2_ice,min_co2_ice_s_1,&
	686	min_co2_ice_s_2,iim,jjm,ngrid,tendencies_co2_ice_phys)
	687
	688	call compute_tendencies_slope(tendencies_co2_ice_slope,min_co2_ice_slope_1,&
	689	min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope)
	690
[2794]	691	tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:)
	692
[2779]	693	call compute_tendencies_slope(tendencies_h2o_ice_slope,min_h2o_ice_slope_1,&
	694	min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope)
	695
[2835]	696	!------------------------
[2779]	697
[2835]	698	! I Initialisation
	699	! I_a READ run.def
	700	! I_b READ of start_evol.nc and starfi_evol.nc
	701	! I_c Subslope parametrisation
	702	! I_d READ GCM data and convert to the physical grid
	703	! I_e Initialisation of the PEM variable and soil
	704	! I_f Compute tendencies & Save initial situation
	705	! I_g Save initial PCM situation
	706
	707	!---------------------------- Save initial PCM situation ---------------------
	708
[2779]	709	allocate(initial_h2o_ice(ngrid))
	710	allocate(initial_co2_ice(ngrid))
[2794]	711	allocate(initial_co2_ice_sublim_slope(ngrid,nslope))
[2779]	712	allocate(initial_co2_ice_slope(ngrid,nslope))
	713	allocate(initial_h2o_ice_slope(ngrid,nslope))
	714
[2794]	715	! We save the places where water ice is sublimating
[2835]	716	! We compute the surface of water ice sublimating
	717	ini_surf=0.
	718	ini_surf_co2=0.
	719	ini_surf_h2o=0.
	720	Total_surface=0.
[2779]	721	do i=1,ngrid
[2835]	722	Total_surface=Total_surface+cell_area(i)
	723	if (tendencies_h2o_ice_phys(i).LT.0) then
[2779]	724	initial_h2o_ice(i)=1.
[2835]	725	ini_surf=ini_surf+cell_area(i)
	726	else
[2779]	727	initial_h2o_ice(i)=0.
[2835]	728	endif
[2779]	729	do islope=1,nslope
	730	if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
[2794]	731	initial_co2_ice_sublim_slope(i,islope)=1.
[2835]	732	ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
[2794]	733	else
	734	initial_co2_ice_sublim_slope(i,islope)=0.
	735	endif
	736	if (co2ice_slope(i,islope).GT.0) then
[2779]	737	initial_co2_ice_slope(i,islope)=1.
	738	else
	739	initial_co2_ice_slope(i,islope)=0.
	740	endif
	741	if (tendencies_h2o_ice_phys_slope(i,islope).LT.0) then
	742	initial_h2o_ice_slope(i,islope)=1.
[2835]	743	ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
[2779]	744	else
	745	initial_h2o_ice_slope(i,islope)=0.
	746	endif
	747	enddo
	748	enddo
	749
[2835]	750	print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
	751	print *, "Total initial surface of h2o ice sublimating=", ini_surf
	752	print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
	753	print *, "Total surface of the planet=", Total_surface
[2794]	754
[2835]	755	allocate(zplev_gcm(ngrid,nlayer+1))
[2779]	756
[2835]	757	DO l=1,nlayer+1
	758	DO ig=1,ngrid
	759	zplev_gcm(ig,l) = ap(l) + bp(l)*ps_phys(ig)
	760	ENDDO
	761	ENDDO
[2779]	762
[2794]	763	global_ave_press_old=0.
	764	do i=1,ngrid
	765	global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface
	766	enddo
[2779]	767
[2794]	768	global_ave_press_GCM=global_ave_press_old
[2835]	769	global_ave_press_new=global_ave_press_old
	770	print *, "Initial global average pressure=", global_ave_press_GCM
[2779]	771
	772	!------------------------
	773
[2794]	774	! I Initialisation
[2835]	775	! I_a READ run.def
	776	! I_b READ of start_evol.nc and starfi_evol.nc
	777	! I_c Subslope parametrisation
	778	! I_d READ GCM data and convert to the physical grid
	779	! I_e Initialisation of the PEM variable and soil
	780	! I_f Compute tendencies & Save initial situation
	781	! I_g Save initial PCM situation
	782	! I_h Read the PEMstart
[2779]	783
[2835]	784	!---------------------------- Read the PEMstart ---------------------
[2779]	785
[2855]	786	call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
	787	tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsoil_phys_PEM_timeseries,&
[2849]	788	tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM, co2_adsorbded_phys,delta_co2_adsorbded,&
	789	watersurf_density_phys_ave,watersoil_density_phys_PEM_ave)
[2779]	790
[2835]	791	if(soil_pem) then
	792	totmass_adsorbded = 0.
	793
[2794]	794	do ig = 1,ngrid
	795	do islope =1, nslope
	796	do l = 1,nsoilmx_PEM - 1
	797	totmass_adsorbded = totmass_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
	798	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.) &
	799	cell_area(ig)
	800	enddo
	801	enddo
	802	enddo
	803
[2835]	804	write(,) "Tot mass in the regolith=", totmass_adsorbded
	805	deallocate(tsoil_ave_phys_yr1)
	806	endif !soil_pem
[2794]	807	deallocate(tsurf_ave_phys_yr1)
[2835]	808	deallocate(ps_phys_timeseries_yr1)
[2794]	809
[2835]	810	!------------------------
[2794]	811
[2835]	812	! I Initialisation
	813	! I_a READ run.def
	814	! I_b READ of start_evol.nc and starfi_evol.nc
	815	! I_c Subslope parametrisation
	816	! I_d READ GCM data and convert to the physical grid
	817	! I_e Initialisation of the PEM variable and soil
	818	! I_f Compute tendencies & Save initial situation
	819	! I_g Save initial PCM situation
	820	! I_h Read the PEMstar
	821	! I_i Compute orbit criterion
[2842]	822	#ifndef CPP_STD
[2835]	823	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
[2842]	824	#else
	825	CALL iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
	826	#endif
[2794]	827
[2835]	828	if(evol_orbit_pem) then
	829	call orbit_param_criterion(year_iter_max)
	830	else
	831	year_iter_max=Max_iter_pem
	832	endif
[2794]	833
[2835]	834	!--------------------------- END INITIALISATION ---------------------
[2794]	835
	836	!---------------------------- RUN ---------------------
	837
	838	!------------------------
	839
	840	! II Run
	841	! II_a update pressure,ice and tracers
	842
	843	!------------------------
[2835]	844	year_iter=0
	845	do while (year_iter.LT.year_iter_max)
[2794]	846
[2835]	847	! II.a.1. Compute updated global pressure
	848	print *, "Recomputing the new pressure..."
[2779]	849	do i=1,ngrid
	850	do islope=1,nslope
[2794]	851	global_ave_press_new=global_ave_press_new-gcell_area(i)tendencies_co2_ice_phys_slope(i,islope)subslope_dist(i,islope)/cos(pidef_slope_mean(islope)/180.)/Total_surface
[2857]	852	enddo
[2849]	853	enddo
	854	print *, 'Global average pressure old time step',global_ave_press_old
	855	print *, 'Global average pressure new time step',global_ave_press_new
	856
	857	do i=1,ngrid
[2835]	858	if(soil_pem) then
	859	global_ave_press_new = global_ave_press_new -gcell_area(i)delta_co2_adsorbded(i)/Total_surface
	860	endif
[2779]	861	enddo
[2849]	862	print *, 'Global average pressure old time step',global_ave_press_old
	863	print *, 'Global average pressure new time step',global_ave_press_new
[2835]	864
	865	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
[2794]	866	allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
[2835]	867	print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
	868	DO l=1,nlayer+1
	869	DO ig=1,ngrid
	870	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:)
	871	ENDDO
	872	ENDDO
[2779]	873
[2835]	874	! II.a.3. Surface pressure timeseries
	875	print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
[2794]	876	do i = 1,ngrid
	877	ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old
	878	enddo
[2779]	879
[2835]	880	! II.a.4. New pressure levels timeseries
[2794]	881	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
[2835]	882	print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
	883	do l=1,nlayer+1
	884	do ig=1,ngrid
	885	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:)
[2794]	886	enddo
[2835]	887	enddo
[2779]	888
[2835]	889	! II.a.5. Tracers timeseries
	890	print *, "Recomputing of tracer VMR timeseries for the new pressure..."
[2794]	891
[2835]	892	l=1
	893	DO ig=1,ngrid
	894	DO t = 1, timelen
[2794]	895	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))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
	896	if(q_h2o_PEM_phys(ig,t).LT.0) then
	897	q_h2o_PEM_phys(ig,t)=1E-30
	898	endif
	899	if(q_h2o_PEM_phys(ig,t).GT.1) then
	900	q_h2o_PEM_phys(ig,t)=1.
	901	endif
[2835]	902	enddo
	903	enddo
[2794]	904
[2835]	905	DO ig=1,ngrid
	906	DO t = 1, timelen
[2794]	907	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))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) &
	908	+ ( (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) - &
	909	(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t)) ) / &
	910	(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
	911	if (q_co2_PEM_phys(ig,t).LT.0) then
	912	q_co2_PEM_phys(ig,t)=1E-30
	913	elseif (q_co2_PEM_phys(ig,t).GT.1) then
	914	q_co2_PEM_phys(ig,t)=1.
	915	endif
	916	mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
	917	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
[2835]	918	ENDDO
	919	ENDDO
[2794]	920
	921	deallocate(zplev_new_timeseries)
	922	deallocate(zplev_old_timeseries)
	923
[2835]	924	! II Run
	925	! II_a update pressure, ice and tracers
	926	! II_b Evolution of the ice
[2794]	927
[2835]	928	! II.b. Evolution of the ice
	929	print *, "Evolution of h2o ice"
[2794]	930	call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
	931
[2857]	932
[2835]	933	print *, "Evolution of co2 ice"
[2779]	934	call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
	935
[2794]	936	!------------------------
	937
	938	! II Run
	939	! II_a update pressure, ice and tracers
[2835]	940	! II_b Evolution of the ice
	941	! II_c CO2 glaciers flows
[2794]	942
	943	!------------------------
	944
[2856]	945	print *, "Co2 glacier flows"
[2779]	946
	947
	948
[2856]	949	call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_phys_timeseries,&
	950	global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
	951
	952
	953
	954
	955
[2794]	956	!------------------------
[2779]	957
[2794]	958	! II Run
	959	! II_a update pressure, ice and tracers
[2835]	960	! II_b Evolution of the ice
	961	! II_c CO2 glaciers flows
	962	! II_d Update surface and soil temperatures
[2779]	963
[2794]	964	!------------------------
[2779]	965
[2835]	966	! II_d.1 Update Tsurf
[2794]	967
[2835]	968	print *, "Updating the new Tsurf"
[2794]	969	bool_sublim=0
	970	Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:)
[2835]	971	DO ig = 1,ngrid
[2794]	972	DO islope = 1,nslope
[2857]	973	if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-10) THEN !co2ice disappeared, look for closest point without co2ice
[2835]	974	if(latitude_deg(ig).gt.0) then
	975	do ig_loop=ig,ngrid
	976	DO islope_loop=islope,iflat,-1
[2857]	977	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
[2835]	978	tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
	979	bool_sublim=1
	980	exit
[2794]	981	endif
	982	enddo
[2835]	983	if (bool_sublim.eq.1) then
	984	exit
	985	endif
	986	enddo
	987	else
	988	do ig_loop=ig,1,-1
	989	DO islope_loop=islope,iflat
[2857]	990	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
[2835]	991	tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
	992	bool_sublim=1
	993	exit
[2794]	994	endif
	995	enddo
[2835]	996	if (bool_sublim.eq.1) then
	997	exit
	998	endif
	999	enddo
	1000	endif
	1001	initial_co2_ice_slope(ig,islope)=0
[2857]	1002	if ((co2ice_slope(ig,islope).lt.1e-10).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then
[2835]	1003	albedo_slope(ig,1,islope) = albedo_h2o_frost
	1004	albedo_slope(ig,2,islope) = albedo_h2o_frost
	1005	else
	1006	albedo_slope(ig,1,islope) = albedodat(ig)
	1007	albedo_slope(ig,2,islope) = albedodat(ig)
	1008	emiss_slope(ig,islope) = emissiv
	1009	endif
[2857]	1010	elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice_phys_slope(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
[2835]	1011	ave=0.
	1012	do t=1,timelen
[2849]	1013	if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-3) then
[2855]	1014	ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
[2794]	1015	else
[2835]	1016	ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
[2794]	1017	endif
[2835]	1018	enddo
	1019	tsurf_ave_phys(ig,islope)=ave/timelen
	1020	endif
	1021	enddo
	1022	enddo
[2794]	1023
[2835]	1024	do t = 1,timelen
	1025	tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:))
	1026	enddo
	1027	! for the start
	1028	do ig = 1,ngrid
	1029	do islope = 1,nslope
	1030	tsurf_slope(ig,islope) = tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
[2794]	1031	enddo
[2835]	1032	enddo
[2794]	1033
[2849]	1034
[2835]	1035	if(soil_pem) then
[2794]	1036
[2835]	1037	! II_d.2 Update soil temperature
[2794]	1038
[2835]	1039	allocate(TI_locslope(ngrid,nsoilmx_PEM))
	1040	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
	1041	allocate(Tsurf_locslope(ngrid))
	1042	allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
	1043	allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
[2794]	1044
[2835]	1045	print *,"Updating soil temperature"
[2794]	1046
	1047	! Soil averaged
	1048	do islope = 1,nslope
	1049	TI_locslope(:,:) = TI_PEM(:,:,islope)
	1050	alph_locslope(:,:) = alph_PEM(:,:,islope)
	1051	beta_locslope(:,:) = beta_PEM(:,:,islope)
	1052	do t = 1,timelen
[2849]	1053	Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
	1054	Tsurf_locslope(:) = tsurf_phys_GCM_timeseries(:,islope,t)
[2794]	1055
[2849]	1056	call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
	1057	call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
[2794]	1058
[2849]	1059
	1060	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
	1061	do ig = 1,ngrid
	1062	do isoil = 1,nsoilmx_PEM
	1063	watersoil_density_phys_PEM_timeseries(ig,l,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
	1064	if(isnan(Tsoil_locslope(ig,isoil))) then
	1065	write(,)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil)
[2857]	1066	stop
[2849]	1067	endif
[2794]	1068	enddo
[2835]	1069	enddo
	1070	enddo
[2849]	1071	alph_PEM(:,:,islope) = alph_locslope(:,:)
	1072	beta_PEM(:,:,islope) = beta_locslope(:,:)
	1073	enddo
[2794]	1074
[2849]	1075
	1076	tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
	1077	watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
	1078	print *, "Update of soil temperature done"
	1079
	1080	deallocate(TI_locslope)
	1081	deallocate(Tsoil_locslope)
	1082	deallocate(Tsurf_locslope)
	1083	deallocate(alph_locslope)
	1084	deallocate(beta_locslope)
	1085
[2835]	1086	write(,) "Compute ice table"
[2794]	1087
[2835]	1088	! II_d.3 Update the ice table
[2849]	1089	call computeice_table(ngrid,nslope,nsoilmx_PEM,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
[2794]	1090
[2835]	1091	print *, "Update soil propreties"
	1092	! II_d.4 Update the soil thermal properties
	1093	call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
	1094	ice_depth,TI_PEM)
[2794]	1095
[2835]	1096	! II_d.5 Update the mass of the regolith adsorbded
[2794]	1097	call regolith_co2adsorption(ngrid,nslope,nsoilmx_PEM,timelen,ps_phys_timeseries,tsoil_PEM,TI_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope, &
	1098	co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), q_co2_PEM_phys,q_h2o_PEM_phys,co2_adsorbded_phys,delta_co2_adsorbded)
	1099
[2835]	1100	endif !soil_pem
[2794]	1101
	1102	!------------------------
	1103
	1104	! II Run
	1105	! II_a update pressure, ice and tracers
[2835]	1106	! II_b Evolution of the ice
	1107	! II_c CO2 glaciers flows
	1108	! II_d Update surface and soil temperatures
	1109	! II_e Update the tendencies
[2794]	1110
	1111	!------------------------
	1112
[2835]	1113	print *, "Adaptation of the new co2 tendencies to the current pressure"
[2794]	1114	call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,&
	1115	global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
	1116
[2835]	1117	!------------------------
[2779]	1118
[2835]	1119	! II Run
	1120	! II_a update pressure, ice and tracers
	1121	! II_b Evolution of the ice
	1122	! II_c CO2 glaciers flows
	1123	! II_d Update surface and soil temperatures
	1124	! II_e Update the tendencies
	1125	! II_f Checking the stopping criterion
[2779]	1126
[2835]	1127	!------------------------
[2794]	1128	call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
[2779]	1129
[2794]	1130	call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
	1131
[2835]	1132	year_iter=year_iter+dt_pem
[2794]	1133
[2835]	1134	print *, "Checking all the stopping criterion."
	1135	if (STOPPING_water) then
	1136	print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
	1137	endif
[2849]	1138	if (year_iter.ge.year_iter_max) then
	1139	print *, "STOPPING because maximum number of iterations reached"
	1140	endif
[2835]	1141	if (STOPPING_1_water) then
	1142	print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
	1143	endif
	1144	if (STOPPING_co2) then
	1145	print *, "STOPPING because surface of co2 ice sublimating is too low or global pressure changed too much, see message above", STOPPING_co2
	1146	endif
	1147	if (STOPPING_1_co2) then
	1148	print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
	1149	endif
[2794]	1150
[2857]	1151
	1152
[2779]	1153	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2) then
	1154	exit
[2835]	1155	else
	1156	print *, "We continue!"
	1157	print *, "Number of iteration of the PEM : year_iter=", year_iter
[2779]	1158	endif
	1159
	1160	global_ave_press_old=global_ave_press_new
	1161
[2835]	1162	enddo ! big time iteration loop of the pem
[2779]	1163
	1164
[2794]	1165	!---------------------------- END RUN PEM ---------------------
	1166
	1167	!---------------------------- OUTPUT ---------------------
	1168
	1169	!------------------------
	1170
	1171	! III Output
[2835]	1172	! III_a Update surface value for the PCM start files
[2794]	1173
	1174	!------------------------
	1175
[2835]	1176	! III_a.1 Ice update (for startfi)
	1177	! Co2 ice
[2779]	1178	DO ig = 1,ngrid
[2835]	1179	co2ice(ig) = 0.
	1180	DO islope = 1,nslope
	1181	co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) &
	1182	* subslope_dist(ig,islope) / &
	1183	cos(pi*def_slope_mean(islope)/180.)
	1184	ENDDO
[2842]	1185	#ifdef CPP_STD
	1186	qsurf(ig,igcm_co2_ice)=co2ice(ig)
	1187	#endif
[2779]	1188	ENDDO ! of DO ig=1,ngrid
[2835]	1189	! H2o ice
	1190	DO ig = 1,ngrid
	1191	qsurf(ig,igcm_h2o_ice) = 0.
	1192	DO islope = 1,nslope
	1193	qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
	1194	* subslope_dist(ig,islope) / &
	1195	cos(pi*def_slope_mean(islope)/180.)
	1196	ENDDO
[2779]	1197	ENDDO ! of DO ig=1,ngrid
	1198
[2849]	1199	! III_a.2 Tsoil update (for startfi)
[2779]	1200
[2835]	1201	if(soil_pem) then
	1202	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
	1203	tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
	1204	else
	1205	TI_GCM_phys(:,:,:)=TI_GCM_start(:,:,:)
	1206	endif !soil_pem
[2779]	1207
[2794]	1208
[2842]	1209	#ifndef CPP_STD
[2794]	1210	DO ig = 1,ngrid
[2835]	1211	DO iloop = 1,nsoilmx
[2794]	1212	tsoil(ig,iloop) = 0.
	1213	inertiesoil(ig,iloop) = 0.
[2835]	1214	DO islope = 1,nslope
	1215	tsoil(ig,iloop) = tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
[2794]	1216	* subslope_dist(ig,islope)
[2835]	1217	inertiesoil(ig,iloop) = inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
[2794]	1218	* subslope_dist(ig,islope)
	1219	ENDDO
[2835]	1220	ENDDO
[2794]	1221	ENDDO ! of DO ig=1,ngrid
	1222
[2835]	1223	! III_a.3 Surface optical properties (for startfi)
[2794]	1224
	1225	DO ig = 1,ngrid
	1226	DO l = 1,2
	1227	albedo(ig,l) =0.
	1228	DO islope = 1,nslope
	1229	albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) &
	1230	*subslope_dist(ig,islope)
	1231	ENDDO
	1232	ENDDO
[2835]	1233	ENDDO
[2794]	1234
	1235	DO ig = 1,ngrid
[2835]	1236	emis(ig) =0.
	1237	DO islope = 1,nslope
	1238	emis(ig)= emis(ig)+emiss_slope(ig,islope) &
[2794]	1239	*subslope_dist(ig,islope)
[2835]	1240	ENDDO
	1241	ENDDO
[2849]	1242
	1243
	1244	DO ig = 1,ngrid
	1245	tsurf(ig) = 0.
	1246	DO islope = 1,nslope
	1247	tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)tsurf_slope(ig,islope))*4 &
	1248	* subslope_dist(ig,islope)
	1249	ENDDO
	1250	tsurf(ig) = tsurf(ig)**(0.25)/emis(ig)
	1251	ENDDO ! of DO ig=1,ngrid
	1252
[2842]	1253	#endif
[2794]	1254
[2835]	1255	! III_a.4 Pressure (for start)
	1256	do i=1,ip1jmp1
	1257	ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM
	1258	enddo
[2794]	1259
[2835]	1260	do i = 1,ngrid
	1261	ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_GCM
	1262	enddo
	1263
[2857]	1264	write(,) 'rapport ps',global_ave_press_new/global_ave_press_GCM
	1265
[2835]	1266	! III_a.5 Tracer (for start)
	1267	allocate(zplev_new(ngrid,nlayer+1))
	1268
	1269	do l=1,nlayer+1
	1270	do ig=1,ngrid
	1271	zplev_new(ig,l) = ap(l) + bp(l)*ps_phys(ig)
	1272	enddo
	1273	enddo
	1274
	1275	DO nnq=1,nqtot
	1276	if (noms(nnq).NE."co2") then
	1277	DO l=1,llm-1
	1278	DO ig=1,ngrid
	1279	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))
	1280	ENDDO
	1281	q(:,llm,nnq)=q(:,llm-1,nnq)
	1282	ENDDO
	1283	else
	1284	DO l=1,llm-1
	1285	DO ig=1,ngrid
	1286	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)) &
	1287	+ ( (zplev_new(ig,l)-zplev_new(ig,l+1)) - &
	1288	(zplev_gcm(ig,l)-zplev_gcm(ig,l+1)) ) / &
	1289	(zplev_new(ig,l)-zplev_new(ig,l+1))
	1290	ENDDO
	1291	q(:,llm,nnq)=q(:,llm-1,nnq)
	1292	ENDDO
	1293	endif
	1294	ENDDO
	1295
	1296	! Conserving the tracers's mass for GCM start files
	1297	DO nnq=1,nqtot
	1298	DO ig=1,ngrid
	1299	DO l=1,llm-1
[2857]	1300	if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") ) then
[2835]	1301	extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
	1302	q(ig,l,nnq)=1.
	1303	q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
[2857]	1304	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq).NE."dust_number",noms(nnq).NE."ccn_number"
[2835]	1305	endif
	1306	if(q(ig,l,nnq).LT.0) then
	1307	q(ig,l,nnq)=1E-30
	1308	endif
	1309	ENDDO
	1310	enddo
	1311	enddo
	1312
[2779]	1313	!------------------------
[2835]	1314	if(evol_orbit_pem) then
	1315	call recomp_orb_param(year_iter)
	1316	endif
[2779]	1317
	1318	! III Output
[2835]	1319	! III_a Update surface value for the PCM start files
[2794]	1320	! III_b Write start and starfi.nc
[2779]	1321
	1322	!------------------------
	1323
[2794]	1324	! III_b.1 WRITE restart.nc
[2779]	1325
	1326	ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys
	1327	pday=day_ini
	1328	ztime_fin=0.
	1329
	1330	allocate(p(ip1jmp1,nlayer+1))
[2835]	1331	CALL pression (ip1jmp1,ap,bp,ps,p)
	1332	CALL massdair(p,masse)
[2779]	1333
[2794]	1334	CALL dynredem0("restart_evol.nc", day_ini, phis)
	1335
	1336	CALL dynredem1("restart_evol.nc", &
[2779]	1337	time_0,vcov,ucov,teta,q,masse,ps)
[2835]	1338	print *, "restart_evol.nc has been written"
[2779]	1339
[2794]	1340	! III_b.2 WRITE restartfi.nc
[2842]	1341	#ifndef CPP_STD
[2835]	1342	call physdem0("restartfi_evol.nc",longitude,latitude, &
[2794]	1343	nsoilmx,ngrid,nlayer,nq, &
	1344	ptimestep,pday,0.,cell_area, &
	1345	albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
	1346	hmons,summit,base,nslope,def_slope, &
	1347	subslope_dist)
[2779]	1348
[2835]	1349	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
[2779]	1350	ptimestep,ztime_fin, &
	1351	tsurf,tsoil,co2ice,albedo,emis, &
	1352	q2,qsurf,tauscaling,totcloudfrac,wstar, &
[2794]	1353	watercap,inertiesoil,nslope,co2ice_slope, &
[2779]	1354	tsurf_slope,tsoil_slope, albedo_slope, &
[2794]	1355	emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys)
[2842]	1356	#else
	1357	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
	1358	ptimestep,pday,time_phys,cell_area, &
	1359	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[2779]	1360
[2842]	1361	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
	1362	ptimestep,ztime_fin, &
	1363	tsurf,tsoil,emis,q2,qsurf, &
	1364	cloudfrac,totcloudfrac,hice, &
	1365	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
	1366	#endif
	1367
[2835]	1368	print *, "restartfi_evol.nc has been written"
[2794]	1369	!------------------------
[2779]	1370
[2794]	1371	! III Output
[2835]	1372	! III_a Update surface value for the PCM start files
[2794]	1373	! III_b Write start and starfi.nc
	1374	! III_c Write start_pem
[2779]	1375
[2794]	1376	!------------------------
[2779]	1377
[2835]	1378	call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
[2794]	1379	tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys)
[2779]	1380
[2835]	1381	print *, "restartfi_PEM.nc has been written"
	1382	print *, "The PEM had run for ", year_iter, " years."
[2794]	1383	print *, "LL & RV : So far so good"
	1384
[2779]	1385	END PROGRAM pem

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format