Context Navigation

pemetat0.F90 @ 3595

Last change on this file since 3595 was 3571, checked in by jbclement, 4 weeks ago

PEM:

New way to manage the pressure: now the PEM manages only the average pressure and keeps the pressure deviation with the instantaneous pressure from the start to reconstruct the pressure at the end ('ps_avg = ps_start + ps_dev'). As a consequence, everything related to pressure in the PEM is modified accordingly.
Surface temperatures management is now simpler. It follows the strategy for the pressure (and soil temperature) described above.
Soil temperatures are now adapted to match the surface temperature changes occured during the PEM by modifying the soil temperature deviation at the end.
Few simplifications/optimizations: notably, the two PCM years are now read in one go in 'read_data_PCM_mod.F90' and only the needed variables are extracted.
Deletion of unused variables and unnecessary intermediate variables (memory saving and loop deletion in some cases).
Renaming of variables and subroutines to make everything clearer. In particular, the suffixes: '_avg' = average, '_start' = PCM start file, '_dev' = deviation, '_ini' or '0' = initial, '_dyn' = dynamical grid, '_timeseries' = daily average of last PCM year.
Cosmetic cleanings for readability.

JBC

File size: 25.0 KB

Rev	Line
[3076]	1	MODULE pemetat0_mod
[3019]	2
[3076]	3	implicit none
	4
	5	!=======================================================================
	6	contains
	7	!=======================================================================
	8
[3537]	9	SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness, &
[3571]	10	ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global, &
	11	watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys, &
	12	d_h2oice,d_co2ice,co2_ice,h2o_ice,m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
[3076]	13
[3297]	14	use iostart_PEM, only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
[3571]	15	use comsoil_h_PEM, only: soil_pem, layer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock
	16	use comsoil_h, only: inertiedat
[3076]	17	use adsorption_mod, only: regolith_adsorption, adsorption_pem
[3537]	18	use ice_table_mod, only: computeice_table_equilibrium, icetable_equilibrium, icetable_dynamic
[3076]	19	use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
	20	use soil_thermalproperties_mod, only: update_soil_thermalproperties
	21	use tracer_mod, only: mmol, igcm_h2o_vap ! tracer names and molar masses
	22	use abort_pem_mod, only: abort_pem
[3178]	23	use compute_soiltemp_mod, only: ini_tsoil_pem, compute_tsoil_pem
[3319]	24	use layering_mod, only: layering, array2stratif, nb_str_max, layering_algo
[3076]	25
[3019]	26	#ifndef CPP_STD
[3161]	27	use comcstfi_h, only: r, mugaz, pi
	28	use surfdat_h, only: watercaptag, watercap, perennial_co2ice
[3019]	29	#else
[3161]	30	use comcstfi_mod, only: r, mugaz, pi
[3019]	31	#endif
	32
[3076]	33	implicit none
[3019]	34
[3076]	35	include "callkeys.h"
[2835]	36
[3571]	37	character(*), intent(in) :: filename ! name of the startfi_PEM.nc
	38	integer, intent(in) :: ngrid ! # of physical grid points
	39	integer, intent(in) :: nsoil_PCM ! # of vertical grid points in the PCM
	40	integer, intent(in) :: nsoil_PEM ! # of vertical grid points in the PEM
	41	integer, intent(in) :: nslope ! # of sub-grid slopes
	42	integer, intent(in) :: timelen ! # time samples
	43	real, intent(in) :: timestep ! time step [s]
	44	real, intent(in) :: ps_avg_global ! mean average pressure on the planet [Pa]
	45	real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr1 ! surface temperature at the first year of PCM call [K]
	46	real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr2 ! surface temperature at the second year of PCM call [K]
	47	real, dimension(ngrid,timelen), intent(in) :: q_co2 ! MMR tracer co2 [kg/kg]
	48	real, dimension(ngrid,timelen), intent(in) :: q_h2o ! MMR tracer h2o [kg/kg]
	49	real, dimension(ngrid,timelen), intent(in) :: ps_timeseries ! surface pressure [Pa]
	50	real, dimension(ngrid,nslope), intent(in) :: d_h2oice ! tendencies on h2o ice
	51	real, dimension(ngrid,nslope), intent(in) :: d_co2ice ! tendencies on co2 ice
	52	real, dimension(ngrid,nslope), intent(in) :: watersurf_avg ! surface water ice density, yearly averaged (kg/m^3)
[2794]	53	! outputs
[3571]	54	real, dimension(ngrid), intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
	55	real, dimension(ngrid), intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
	56	real, dimension(ngrid,nslope), intent(out) :: h2o_ice ! h2o ice amount [kg/m^2]
	57	real, dimension(ngrid,nslope), intent(out) :: co2_ice ! co2 ice amount [kg/m^2]
	58	type(layering), dimension(ngrid,nslope), intent(inout) :: stratif ! stratification (layerings)
	59	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! soil (mid-layer) thermal inertia in the PEM grid [SI]
	60	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM ! soil (mid-layer) temperature [K]
	61	real, dimension(ngrid,nslope), intent(inout) :: icetable_depth ! Ice table depth [m]
	62	real, dimension(ngrid,nslope), intent(inout) :: icetable_thickness ! Ice table thickness [m]
	63	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling ! Subsurface ice pore filling [m3/m3]
	64	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
	65	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
	66	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg ! surface water ice density, yearly averaged (kg/m^3)
[2794]	67	! local
[3571]	68	real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem ! soil temperature saved in the start [K]
	69	real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM ! soil thermal inertia saved in the start [SI]
	70	logical :: found, found2 ! check if variables are found in the start
	71	integer :: iloop, ig, islope, isoil ! index for loops
	72	real :: delta ! Depth of the interface regolith-breccia, breccia -bedrock [m]
	73	character(2) :: num ! intermediate string to read PEM start sloped variables
	74	logical :: startpem_file ! boolean to check if we read the startfile or not
	75	real, dimension(:,:,:,:), allocatable :: stratif_array ! Array for stratification (layerings)
[2835]	76
[3076]	77	#ifdef CPP_STD
[3149]	78	logical, dimension(ngrid) :: watercaptag
[3159]	79	watercaptag = .false.
[2985]	80	#endif
	81
[2794]	82	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	83	!!!
	84	!!! Purpose: read start_pem. Need a specific iostart_PEM
	85	!!!
[2888]	86	!!! Order: 0. Previous year of the PEM run
[3149]	87	!!! Ice initialization
[3076]	88	!!! 1. Thermal Inertia
[2888]	89	!!! 2. Soil Temperature
	90	!!! 3. Ice table
	91	!!! 4. Mass of CO2 & H2O adsorbed
[2794]	92	!!!
	93	!!! /!\ This order must be respected !
	94	!!! Author: LL
	95	!!!
	96	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2779]	97
[3170]	98	!0.1 Check if the start_PEM exist.
[2863]	99
[3143]	100	inquire(file = filename,exist = startpem_file)
[2863]	101
[3498]	102	write(,)'Is there a "startpem" file?',startpem_file
[2863]	103
[3170]	104	!0.2 Set to default values
[3537]	105	icetable_depth = -1. ! by default, no ice table
	106	icetable_thickness = -1.
	107	ice_porefilling = 0.
[2794]	108	!1. Run
	109	if (startpem_file) then
[3214]	110	write(,)'There is a "startpem.nc"!'
[3076]	111	! open pem initial state file:
	112	call open_startphy(filename)
[2835]	113
[3149]	114	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3430]	115	#ifndef CPP_STD
	116	! h2o ice
	117	h2o_ice = 0.
	118	call get_field("h2o_ice",h2o_ice,found)
	119	if (.not. found) then
	120	write(,)'Pemetat0: failed loading <h2o_ice>'
	121	write(,)'will reconstruct the values from watercaptag'
	122	write(,)'with default value ''ini_huge_h2oice'''
	123	do ig = 1,ngrid
	124	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
	125	enddo
	126	else
	127	! The variations of infinite reservoirs during the PCM years are taken into account
	128	h2o_ice = h2o_ice + watercap
	129	endif
	130
	131	! co2 ice
	132	co2_ice = perennial_co2ice
	133	#endif
	134
	135	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3297]	136	! Stratification (layerings)
[3430]	137	nb_str_max = 1
[3319]	138	if (layering_algo) then
	139	found = inquire_dimension("nb_str_max")
	140	if (.not. found) then
	141	write(,) 'Pemetat0: failed loading <nb_str_max>!'
	142	write(,) '''nb_str_max'' is set to 3!'
	143	else
	144	nb_str_max = int(inquire_dimension_length('nb_str_max'))
	145	endif
	146	allocate(stratif_array(ngrid,nslope,nb_str_max,6))
	147	stratif_array = 0.
	148	do islope = 1,nslope
	149	write(num,'(i2.2)') islope
	150	call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
	151	call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
	152	found = found .or. found2
	153	call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
	154	found = found .or. found2
	155	call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
	156	found = found .or. found2
	157	call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
	158	found = found .or. found2
	159	call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
	160	found = found .or. found2
	161	if (.not. found) then
	162	write(,) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
	163	endif ! not found
	164	enddo ! islope
	165	call array2stratif(stratif_array,ngrid,nslope,stratif)
	166	deallocate(stratif_array)
[3297]	167	endif
	168
[3039]	169	if (soil_pem) then
[3076]	170
[3149]	171	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2794]	172	!1. Thermal Inertia
	173	! a. General case
[3149]	174	do islope = 1,nslope
	175	write(num,'(i2.2)') islope
	176	call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
	177	if (.not. found) then
	178	write(,)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
	179	write(,)'will reconstruct the values of TI_PEM'
[2779]	180
[3149]	181	do ig = 1,ngrid
	182	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
	183	!!! transition
	184	delta = depth_breccia
[3430]	185	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
	186	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
	187	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
	188	do iloop=index_breccia + 2,index_bedrock
[3149]	189	TI_PEM(ig,iloop,islope) = TI_breccia
	190	enddo
[3159]	191	else ! we keep the high TI values
[3149]	192	do iloop = index_breccia + 1,index_bedrock
	193	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
	194	enddo
	195	endif ! TI PEM and breccia comparison
	196	!! transition
	197	delta = depth_bedrock
	198	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
	199	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
	200	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
	201	do iloop = index_bedrock + 2,nsoil_PEM
	202	TI_PEM(ig,iloop,islope) = TI_bedrock
[3076]	203	enddo
	204	enddo
[3149]	205	else ! found
	206	do iloop = nsoil_PCM + 1,nsoil_PEM
[3571]	207	TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope) ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
[3149]	208	enddo
	209	endif ! not found
	210	enddo ! islope
[2779]	211
[3149]	212	write(,) 'PEMETAT0: THERMAL INERTIA done'
[2779]	213
[2794]	214	! b. Special case for inertiedat, inertiedat_PEM
[3149]	215	call get_field("inertiedat_PEM",inertiedat_PEM,found)
[3456]	216	if (.not. found) then
[3149]	217	do iloop = 1,nsoil_PCM
	218	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
	219	enddo
	220	!!! zone de transition
	221	delta = depth_breccia
	222	do ig = 1,ngrid
[3430]	223	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
[3493]	224	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
	225	(((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
[3149]	226	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
[2779]	227
[3430]	228	do iloop = index_breccia + 2,index_bedrock
[3149]	229	inertiedat_PEM(ig,iloop) = TI_breccia
	230	enddo
	231	else
[3430]	232	do iloop=index_breccia + 1,index_bedrock
[3149]	233	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
	234	enddo
	235	endif ! comparison ti breccia
	236	enddo ! ig
[2779]	237
[3149]	238	!!! zone de transition
	239	delta = depth_bedrock
	240	do ig = 1,ngrid
[3430]	241	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
	242	(((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
	243	((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
[3149]	244	enddo
[2779]	245
[3149]	246	do iloop = index_bedrock + 2, nsoil_PEM
	247	do ig = 1,ngrid
	248	inertiedat_PEM(ig,iloop) = TI_bedrock
	249	enddo
	250	enddo
	251	endif ! not found
[2779]	252
[3076]	253	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2794]	254	!2. Soil Temperature
[3571]	255	do islope = 1,nslope
[3149]	256	write(num,fmt='(i2.2)') islope
[3571]	257	call get_field("tsoil_PEM_slope"//num,tsoil_startpem(:,:,islope),found)
[3456]	258	if (.not. found) then
[3149]	259	write(,)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
	260	write(,)'will reconstruct the values of Tsoil'
[3178]	261	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
	262	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
[3149]	263	else
[3571]	264	! predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startpem is adapted firstly
	265	! for PCM year 1 and then for PCM year 2 in order to build the evolution of the profile at depth
	266	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
	267	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
	268	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope))
[3203]	269
[3532]	270	do iloop = nsoil_PCM + 1,nsoil_PEM
[3571]	271	tsoil_PEM(:,iloop,islope) = tsoil_startpem(:,iloop,islope)
[3149]	272	enddo
	273	endif !found
[2794]	274
[3532]	275	do isoil = nsoil_PCM + 1,nsoil_PEM
	276	watersoil_avg(:,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(:,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(:,isoil,islope)mmol(igcm_h2o_vap)/(mugazr)
[3149]	277	enddo
	278	enddo ! islope
	279	write(,) 'PEMETAT0: SOIL TEMP done'
[2794]	280
[3076]	281	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2794]	282	!3. Ice Table
[3493]	283	if (icetable_equilibrium) then
[3537]	284	call get_field("icetable_depth",icetable_depth,found)
[3170]	285	if (.not. found) then
[3537]	286	write(,)'PEM settings: failed loading <icetable_depth>'
[3170]	287	write(,)'will reconstruct the values of the ice table given the current state'
[3537]	288	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
[3571]	289	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
[3170]	290	do islope = 1,nslope
[3178]	291	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
[3170]	292	enddo
	293	endif
	294	write(,) 'PEMETAT0: ICE TABLE done'
[3493]	295	else if (icetable_dynamic) then
[3532]	296	call get_field("ice_porefilling",ice_porefilling,found)
[3537]	297	if (.not. found) write(,)'PEM settings: failed loading <ice_porefilling>'
	298	call get_field("icetable_depth",icetable_depth,found)
[3532]	299	if (.not. found) then
[3537]	300	write(,)'PEM settings: failed loading <icetable_depth>'
[3493]	301	write(,)'will reconstruct the values of the ice table given the current state'
[3571]	302	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
[3493]	303	do islope = 1,nslope
	304	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
	305	enddo
	306	endif
	307	write(,) 'PEMETAT0: ICE TABLE done'
[3149]	308	endif
[2794]	309
[3076]	310	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2863]	311	!4. CO2 & H2O Adsorption
[3149]	312	if (adsorption_pem) then
	313	do islope = 1,nslope
	314	write(num,fmt='(i2.2)') islope
	315	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
	316	if (.not. found) then
	317	m_co2_regolith_phys = 0.
	318	exit
	319	endif
	320	enddo
[3571]	321	do islope = 1,nslope
[3149]	322	write(num,fmt='(i2.2)') islope
	323	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
	324	if (.not. found2) then
	325	m_h2o_regolith_phys = 0.
	326	exit
	327	endif
	328	enddo
[2794]	329
[3571]	330	call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
	331	m_h2o_regolith_phys,deltam_h2o_regolith_phys,m_co2_regolith_phys,deltam_co2_regolith_phys)
[3076]	332
[3149]	333	if (.not. found) deltam_co2_regolith_phys = 0.
[3456]	334	if (.not. found2) deltam_h2o_regolith_phys = 0.
[2849]	335
[3149]	336	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
	337	endif ! adsorption_pem
	338	endif ! soil_pem
[3076]	339
[3149]	340	call close_startphy
[2779]	341
[3493]	342	else !No startpem, let's build all by hand
[2794]	343
[3297]	344	write(,)'There is no "startpem.nc"!'
	345
[3149]	346	! h2o ice
	347	h2o_ice = 0.
[3214]	348	write(,)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
	349	do ig = 1,ngrid
[3308]	350	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
[3214]	351	enddo
[2888]	352
[3149]	353	! co2 ice
[3214]	354	write(,)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
[3149]	355	co2_ice = perennial_co2ice
[2794]	356
[3430]	357	! Stratification (layerings)
	358	nb_str_max = 1
	359	if (layering_algo) write(,)'So the stratification (layerings) is initialized with 1 sub-surface stratum.'
	360
[3039]	361	if (soil_pem) then
[2794]	362
	363	!a) Thermal inertia
[3149]	364	do islope = 1,nslope
	365	do ig = 1,ngrid
	366	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
	367	!!! transition
	368	delta = depth_breccia
	369	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
	370	(((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
	371	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
	372	do iloop=index_breccia + 2,index_bedrock
	373	TI_PEM(ig,iloop,islope) = TI_breccia
	374	enddo
	375	else ! we keep the high ti values
	376	do iloop=index_breccia + 1,index_bedrock
	377	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
	378	enddo
	379	endif
	380	!! transition
	381	delta = depth_bedrock
	382	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
	383	(((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
	384	((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
	385	do iloop = index_bedrock + 2,nsoil_PEM
	386	TI_PEM(ig,iloop,islope) = TI_bedrock
	387	enddo
	388	enddo
	389	enddo
[2794]	390
[3149]	391	do iloop = 1,nsoil_PCM
	392	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
	393	enddo
	394	!!! zone de transition
	395	delta = depth_breccia
	396	do ig = 1,ngrid
	397	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
	398	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
	399	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
	400	((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
	401	do iloop = index_breccia + 2,index_bedrock
	402	inertiedat_PEM(ig,iloop) = TI_breccia
	403	enddo
	404	else
	405	do iloop = index_breccia + 1,index_bedrock
	406	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
	407	enddo
	408	endif
	409	enddo
[2794]	410
[3149]	411	!!! zone de transition
	412	delta = depth_bedrock
	413	do ig = 1,ngrid
	414	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
	415	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
	416	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
	417	enddo
[2779]	418
[3149]	419	do iloop = index_bedrock + 2,nsoil_PEM
	420	do ig = 1,ngrid
	421	inertiedat_PEM(ig,iloop) = TI_bedrock
	422	enddo
	423	enddo
[2779]	424
[3149]	425	write(,) 'PEMETAT0: TI done'
[2779]	426
[3076]	427	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	428	!b) Soil temperature
[3149]	429	do islope = 1,nslope
[3178]	430	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
	431	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
[2794]	432
[3149]	433	! First raw initialization
	434	do isoil = nsoil_PCM + 1,nsoil_PEM
	435	do ig = 1,ngrid
	436	watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)mmol(igcm_h2o_vap)/(mugazr)
	437	enddo
	438	enddo
	439	enddo !islope
	440	write(,) 'PEMETAT0: TSOIL done'
[3039]	441
[3076]	442	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3039]	443	!c) Ice table
[3493]	444	if (icetable_equilibrium) then
[3537]	445	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
[3571]	446	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
[3170]	447	do islope = 1,nslope
[3178]	448	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
[3170]	449	enddo
	450	write(,) 'PEMETAT0: Ice table done'
[3493]	451	else if (icetable_dynamic) then
[3571]	452	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
[3493]	453	do islope = 1,nslope
	454	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
	455	enddo
	456	write(,) 'PEMETAT0: Ice table done'
[3170]	457	endif
[3532]	458
[3076]	459	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2794]	460	!d) Regolith adsorbed
[3149]	461	if (adsorption_pem) then
	462	m_co2_regolith_phys = 0.
	463	m_h2o_regolith_phys = 0.
[3571]	464	call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
[3149]	465	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
	466	deltam_co2_regolith_phys = 0.
	467	deltam_h2o_regolith_phys = 0.
	468	endif
[2863]	469
[3149]	470	write(,) 'PEMETAT0: CO2 adsorption done'
	471	endif !soil_pem
[2893]	472	endif ! of if (startphy_file)
	473
[3076]	474	if (soil_pem) then ! Sanity check
	475	do ig = 1,ngrid
	476	do islope = 1,nslope
	477	do iloop = 1,nsoil_PEM
	478	if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
	479	enddo
	480	enddo
	481	enddo
	482	endif !soil_pem
[2794]	483
[3076]	484	END SUBROUTINE pemetat0
	485
	486	END MODULE pemetat0_mod

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

Download in other formats:

Original Format