Context Navigation

pemetat0.F90 @ 3770

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

PEM:
Revision of the layering structure and algorithm:

'stratum' components are now expressed in height
deletion of the (redundant) thickness feature of 'stratum'
'stratif' in the PEM main program is renamed 'deposits'
addition of several procedures to get useful information about a stratum (major component or thickness)
all subsurface layers are now represented in the layering structure by strata with negative top elevation
simplification of the different situations arising from the input tendencies
porosity is determined for the entire stratum (and not anymore for each component) based on dominant component
sublimation of CO2 and H2O ice is now handled simultaneously (more realistic) in a stratum
linking the layering algorithm with the PEM initilization/finalization regarding PCM data and with the PEM stopping criteria
making separate cases for glaciers vs layering management
H2O sublimation flux correction is now handled with the layering when a dust lag layer layer is growing
update of 'h2o_ice_depth' and 'zdqsdif' accordingly at the PEM end for the PCM

JBC

File size: 27.1 KB

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

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

Download in other formats:

Original Format