Context Navigation

pemetat0.F90 @ 3896

Last change on this file since 3896 was 3850, checked in by jbclement, 7 months ago

PEM:

Handling the case where the new date is positive (present-day climate) in "recomp_orb_param_mod.F90", especially because of rounding errors.
Correction in "visu_evol_layering.py" to get the proper extent of the stratication + replacing 'imshow' by 'pcolormesh' to better process the data.
Addition of a graphic to plot the changes in obliquity directly onto the evolution of the stratification over time for better analysis. Obliquity curve is colored according to stratification gain/loss.
Correction of numbering in "lib_launchPEM.sh" when requesting a relaunch.
Update of "PEMrun.job" due to new parsing of arguments introduced in r3836/r3837.

JBC

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

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

Download in other formats:

Original Format