Context Navigation

pemetat0.F90 @ 3983

Last change on this file since 3983 was 3983, checked in by jbclement, 11 days ago

PEM:

Removing completely the ice metamorphism computed by a threshold at the end of the PCM (which was commented).
Addition of a module "metamorphism" to compute the PCM frost at the PEM beginning and give it back to the PCM at the PEM end. The frost is considered as the ice given by the PCM "startfi.nc" which is above the yearly minimum. Thereby, metamorphism is performed through this operation.
Ice reservoirs representation in the PEM is modernized.

JBC

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

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

Download in other formats:

Original Format