Context Navigation

source: branches/ORCHIDEE_3_CMIP6/ORCHIDEE/src_stomate/stomate_soil_carbon_discretization.f90 @ 8367

Last change on this file since 8367 was 6731, checked in by bertrand.guenet, 4 years ago
bug correction: missing variables in stomate_Cforcing
Property svn:executable set to ``*
File size: 188.7 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_soil_carbon_discretization
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see
8	! ORCHIDEE/ORCHIDEE_CeCILL.LIC
9	!
10	!>\BRIEF Calculate permafrost soil carbon dynamics following POPCRAN by Dmitry Khvorstyanov
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! SVN :
17	!! $HeadURL:
18	!svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-MICT/ORCHIDEE/src_stomate/stomate_soilcarbon.f90
19	!$
20	!! $Date: 2013-10-14 15:38:24 +0200 (Mon, 14 Oct 2013) $
21	!! $Revision: 1536 $
22	!! \n
23	!_
24	!================================================================================================================================
25
26	MODULE stomate_soil_carbon_discretization
27
28	! modules used:
29	USE ioipsl_para
30	USE constantes_soil_var
31	USE constantes_soil
32	USE constantes_var
33	USE pft_parameters
34	USE vertical_soil
35	USE stomate_data
36	USE grid
37	USE mod_orchidee_para
38	USE xios_orchidee
39
40	IMPLICIT NONE
41	PRIVATE
42	PUBLIC stomate_soil_carbon_discretization_deep_somcycle, stomate_soil_carbon_discretization_clear, &
43	stomate_soil_carbon_discretization_microactem, calc_vert_int_som
44
45	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: zf_soil !! depths of full levels (m)
46	!$OMP THREADPRIVATE(zf_soil)
47	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: zi_soil !! depths of intermediate levels (m)
48	!$OMP THREADPRIVATE(zi_soil)
49	REAL(r_std), SAVE :: mu_soil
50	!$OMP THREADPRIVATE(mu_soil)
51	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: alphaO2_soil
52	!$OMP THREADPRIVATE(alphaO2_soil)
53	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: betaO2_soil
54	!$OMP THREADPRIVATE(betaO2_soil)
55	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: alphaCH4_soil
56	!$OMP THREADPRIVATE(alphaCH4_soil)
57	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: betaCH4_soil
58	!$OMP THREADPRIVATE(betaCH4_soil)
59
60	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: heights_snow !! total thickness of snow levels (m)
61	!$OMP THREADPRIVATE(heights_snow)
62	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: zf_snow !! depths of full levels (m)
63	!$OMP THREADPRIVATE(zf_snow)
64	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: zi_snow !! depths of intermediate levels (m)
65	!$OMP THREADPRIVATE(zi_snow)
66	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: zf_snow_nopftdim !! depths of full levels (m)
67	!$OMP THREADPRIVATE(zf_snow_nopftdim)
68	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: zi_snow_nopftdim !! depths of intermediate levels (m)
69	!$OMP THREADPRIVATE(zi_snow_nopftdim)
70
71	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zf_coeff_snow
72	!$OMP THREADPRIVATE(zf_coeff_snow)
73	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zi_coeff_snow
74	!$OMP THREADPRIVATE(zi_coeff_snow)
75	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mu_snow
76	!$OMP THREADPRIVATE(mu_snow)
77	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: alphaO2_snow
78	!$OMP THREADPRIVATE(alphaO2_snow)
79	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: betaO2_snow
80	!$OMP THREADPRIVATE(betaO2_snow)
81	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: alphaCH4_snow
82	!$OMP THREADPRIVATE(alphaCH4_snow)
83	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: betaCH4_snow
84	!$OMP THREADPRIVATE(betaCH4_snow)
85
86	real(r_std), allocatable, save, dimension(:,:,:,:) :: deepSOM_pftmean !! Deep soil organic matter profiles, mean over all PFTs
87	!$OMP THREADPRIVATE(deepSOM_pftmean)
88
89	INTEGER(i_std), SAVE :: yr_len = 360
90	!$OMP THREADPRIVATE(yr_len)
91	!! Arrays related to cryoturbation processes
92	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: diff_k !! Diffusion constant (m^2/s)
93	!$OMP THREADPRIVATE(diff_k)
94	REAL(r_std), DIMENSION(:,:), ALLOCATABLE, SAVE :: xe_a
95	!$OMP THREADPRIVATE(xe_a)
96	REAL(r_std), DIMENSION(:,:), ALLOCATABLE, SAVE :: xe_s
97	!$OMP THREADPRIVATE(xe_s)
98	REAL(r_std), DIMENSION(:,:), ALLOCATABLE, SAVE :: xe_p
99	!$OMP THREADPRIVATE(xe_p)
100	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: xc_cryoturb
101	!$OMP THREADPRIVATE(xc_cryoturb)
102	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: xd_cryoturb
103	!$OMP THREADPRIVATE(xd_cryoturb)
104	REAL(r_std), DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: alpha_a
105	!$OMP THREADPRIVATE(alpha_a)
106	REAL(r_std), DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: alpha_s
107	!$OMP THREADPRIVATE(alpha_s)
108	REAL(r_std), DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: alpha_p
109	!$OMP THREADPRIVATE(alpha_p)
110	REAL(r_std), DIMENSION(:,:), ALLOCATABLE, SAVE :: mu_soil_rev
111	!$OMP THREADPRIVATE(mu_soil_rev)
112
113	REAL(r_std), DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: beta_a
114	!$OMP THREADPRIVATE(beta_a)
115	REAL(r_std), DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: beta_s
116	!$OMP THREADPRIVATE(beta_s)
117	REAL(r_std), DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: beta_p
118	!$OMP THREADPRIVATE(beta_p)
119	LOGICAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: cryoturb_location
120	!$OMP THREADPRIVATE(cryoturb_location)
121	LOGICAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: bioturb_location
122	!$OMP THREADPRIVATE(bioturb_location)
123	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: airvol_soil
124	!$OMP THREADPRIVATE(airvol_soil)
125	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: totporO2_soil !! total oxygen porosity in the soil
126	!$OMP THREADPRIVATE(totporO2_soil)
127	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: totporCH4_soil !! total methane porosity in the soil
128	!$OMP THREADPRIVATE(totporCH4_soil)
129	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: conduct_soil
130	!$OMP THREADPRIVATE(conduct_soil)
131	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: diffO2_soil !! oxygen diffusivity in the soil (m**2/s)
132	!$OMP THREADPRIVATE(diffO2_soil)
133	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: diffCH4_soil !! methane diffusivity in the soil (m**2/s)
134	!$OMP THREADPRIVATE(diffCH4_soil)
135	REAL(r_std), DIMENSION(:,:,:),ALLOCATABLE, SAVE :: airvol_snow
136	!$OMP THREADPRIVATE(airvol_snow)
137	REAL(r_std), DIMENSION(:,:,:),ALLOCATABLE, SAVE :: totporO2_snow !! total oxygen porosity in the snow
138	!$OMP THREADPRIVATE(totporO2_snow)
139	REAL(r_std), DIMENSION(:,:,:),ALLOCATABLE, SAVE :: totporCH4_snow !! total methane porosity in the snow
140	!$OMP THREADPRIVATE(totporCH4_snow)
141	REAL(r_std), DIMENSION(:,:,:),ALLOCATABLE, SAVE :: conduct_snow
142	!$OMP THREADPRIVATE(conduct_snow)
143	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: diffCH4_snow !! methane diffusivity in the snow (m**2/s)
144	!$OMP THREADPRIVATE(diffCH4_snow)
145	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: diffO2_snow !! oxygen diffusivity in the snow (m**2/s)
146	!$OMP THREADPRIVATE(diffO2_snow)
147	REAL(r_std), DIMENSION(:,:), ALLOCATABLE, SAVE :: altmax_lastyear !! active layer thickness
148	!$OMP THREADPRIVATE(altmax_lastyear)
149	REAL(r_std), DIMENSION(:,:), ALLOCATABLE, SAVE :: alt
150	!$OMP THREADPRIVATE(alt)
151	INTEGER(i_std), DIMENSION(:,:), ALLOCATABLE, SAVE :: alt_ind !! active layer thickness
152	!$OMP THREADPRIVATE(alt_ind)
153	INTEGER(i_std), DIMENSION(:,:),ALLOCATABLE, SAVE :: altmax_ind !! Maximum over the year active layer thickness
154	!$OMP THREADPRIVATE(altmax_ind)
155	INTEGER(i_std), DIMENSION(:,:),ALLOCATABLE, SAVE :: altmax_ind_lastyear
156	!$OMP THREADPRIVATE(altmax_ind_lastyear)
157	REAL(r_std), DIMENSION(:,:),ALLOCATABLE, SAVE :: z_root !! Rooting depth
158	!$OMP THREADPRIVATE(z_root)
159	INTEGER(i_std), DIMENSION(:,:),ALLOCATABLE, SAVE :: rootlev !! The deepest model level within the rooting depth
160	!$OMP THREADPRIVATE(rootlev)
161	REAL(r_std),DIMENSION(:,:),ALLOCATABLE, SAVE :: lalo_global !! Geogr. coordinates (latitude,longitude) (degrees)
162	!$OMP THREADPRIVATE(lalo_global)
163	LOGICAL,DIMENSION(:,:),ALLOCATABLE, SAVE :: veget_mask_2d !! whether there is vegetation
164	!$OMP THREADPRIVATE(veget_mask_2d)
165
166	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: fc !! flux fractions within carbon pools
167	!$OMP THREADPRIVATE(fc)
168	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: fr !! fraction of decomposed carbon that goes into the atmosphere
169	!$OMP THREADPRIVATE(fr)
170	CONTAINS
171
172	!!
173	!================================================================================================================================
174	!! SUBROUTINE : stomate_soil_carbon_discretization_deep_somcycle
175	!!
176	!>\BRIEF Recalculate vegetation cover and LAI
177	!!
178	!!\n DESCRIPTION :
179	!!
180	!! RECENT CHANGE(S) : None
181	!!
182	!! MAIN OUTPUT VARIABLE(S): None
183	!!
184	!! REFERENCE(S) : None
185	!!
186	!! FLOWCHART :
187	!_
188	!================================================================================================================================
189
190	SUBROUTINE stomate_soil_carbon_discretization_deep_somcycle(kjpindex, index, itau, time_step, lalo, clay, &
191	tsurf, tprof, hslong_in, snow, heat_Zimov, pb, &
192	sfluxCH4_deep, sfluxCO2_deep, &
193	deepSOM_a, deepSOM_s, deepSOM_p, O2_soil, CH4_soil, O2_snow, CH4_snow, &
194	depth_organic_soil, som_in, veget_max, &
195	rprof, altmax, som, som_surf, resp_hetero_soil, fbact, CN_target, fixed_cryoturbation_depth, &
196	snowdz,snowrho, n_mineralisation)
197
198	!! 0. Variable and parameter declaration
199
200	!! 0.1 Input variables
201	INTEGER(i_std), INTENT(in) :: kjpindex
202	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
203	INTEGER(i_std), intent(in) :: itau !! time step number
204	REAL(r_std),DIMENSION(kjpindex,2),INTENT(in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
205	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! surface pressure [pa]
206	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: clay !! clay content
207	INTEGER(i_std),DIMENSION(kjpindex),INTENT(in) :: index !! Indeces of the points on the map
208	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: snow !! Snow mass [Kg/m^2]
209	REAL(r_std), DIMENSION(kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth [m]
210	REAL(r_std), DIMENSION(kjpindex,nsnow), INTENT(in) :: snowrho !! snow density (Kg/m^3)
211	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm),INTENT (in) :: tprof !! deep temperature profile
212	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm),INTENT (in) :: hslong_in !! deep long term soil humidity profile
213	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements),INTENT(in) :: som_in !! carbon going into carbon pools [gC/(m**2 of ground)/day]
214	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: veget_max !! Maximum vegetation fraction
215	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: rprof !! rooting depth (m)
216	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: tsurf !! skin temperature [K]
217	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: fbact !! turnover constant (day)
218	REAL(r_std), DIMENSION(kjpindex,nvm,ncarb), INTENT(in) :: CN_target !! C to N ratio of SOM flux from one pool to another (gN m-2 dt-1)
219	!! 0.2 Output variables
220	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: sfluxCH4_deep !! total CH4 flux [g CH4 / m**2 / s]
221	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: sfluxCO2_deep !! total CO2 flux [g C / m**2 / s]
222	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: resp_hetero_soil !! soil heterotrophic respiration (first in gC/day/m**2 of ground )
223	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT (out) :: heat_Zimov !! Heating associated with decomposition [W/m**3 soil]
224	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements), INTENT (out) :: som !! vertically-integrated (diagnostic) soil carbon pool: active, slow,
225	!! or passive, (gC/(m**2 of ground))
226	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements), INTENT (out) :: som_surf !! vertically-integrated (diagnostic) soil carbon pool: active, slow,
227	!! or passive, (gC/(m**2 of ground))
228
229	!! 0.3 Modified variables
230	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deepSOM_a !! Active soil carbon (g/m**3)
231	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deepSOM_s !! Slow soil carbon (g/m**3)
232	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deepSOM_p !! Passive soil carbon (g/m**3)
233	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: O2_snow !! oxygen in the snow (g O2/m**3 air)
234	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: O2_soil !! oxygen in the soil (g O2/m**3 air)
235	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: CH4_snow !! methane in the snow (g CH4/m**3 air)
236	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: CH4_soil !! methane in the soil (g CH4/m**3 air)
237	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: altmax !! active layer thickness (m)
238	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(inout) :: fixed_cryoturbation_depth !! depth to hold cryoturbation to for fixed runs
239	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: n_mineralisation !! net nitrogen mineralisation of decomposing SOM
240	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: depth_organic_soil !! depth to organic soil
241
242	!! 0.4 Local variables
243
244	REAL(r_std), DIMENSION(kjpindex) :: overburden
245	REAL(r_std), DIMENSION(kjpindex,nvm) :: fluxCH4,febul
246	REAL(r_std), DIMENSION(kjpindex,nvm) :: sfluxCH4
247	REAL(r_std), DIMENSION(kjpindex,nvm) :: flupmt
248	REAL(r_std), DIMENSION(kjpindex,nvm) :: MT !! depth-integrated methane consumed in methanotrophy
249	REAL(r_std), DIMENSION(kjpindex,nvm) :: MG !! depth-integrated methane released in methanogenesis
250	REAL(r_std), DIMENSION(kjpindex,nvm) :: CH4i !! depth-integrated methane
251	REAL(r_std), DIMENSION(kjpindex,nvm) :: CH4ii !! depth-integrated initial methane
252	REAL(r_std), DIMENSION(kjpindex,nvm) :: dC1i !! depth-integrated oxic decomposition carbon
253	REAL(r_std), DIMENSION(kjpindex,nvm) :: dCi !! depth-integrated soil carbon
254
255	REAL(r_std), DIMENSION(kjpindex,nvm) :: Tref !! Ref. temperature for growing season caluculation (C)
256	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: deltaCH4g, deltaCH4
257	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements) :: deltaSOM1_a,deltaSOM1_s,deltaSOM1_p,deltaSOM2,deltaSOM3
258	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: CH4ini_soil
259	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: hslong !! deep long term soil humidity profile
260	INTEGER(i_std) :: ip, il, itz, iz
261	REAL(r_std), SAVE, DIMENSION(3) :: lhc !! specific heat of soil organic matter oxidation (J/kg carbon)
262	!$OMP THREADPRIVATE(lhc)
263	REAL(r_std), SAVE :: O2m !! oxygen concentration [g/m3] below which there is anoxy
264	!$OMP THREADPRIVATE(O2m)
265	LOGICAL, SAVE :: ok_methane !! Is Methanogenesis and -trophy taken into account?
266	!$OMP THREADPRIVATE(ok_methane)
267	LOGICAL, SAVE :: ok_cryoturb !! cryoturbate the carbon?
268	!$OMP THREADPRIVATE(ok_cryoturb)
269	REAL(r_std), SAVE :: cryoturbation_diff_k_in !! input time constant of cryoturbation (m^2/y)
270	!$OMP THREADPRIVATE(cryoturbation_diff_k_in)
271	REAL(r_std), SAVE :: bioturbation_diff_k_in !! input time constant of bioturbation (m^2/y)
272	!$OMP THREADPRIVATE(bioturbation_diff_k_in)
273	REAL(r_std), SAVE :: tau_CH4troph !! time constant of methanetrophy (s)
274	!$OMP THREADPRIVATE(tau_CH4troph)
275	REAL(r_std), SAVE :: fbactratio !! time constant of methanogenesis (ratio to that of oxic)
276	!$OMP THREADPRIVATE(fbactratio)
277	LOGICAL, SAVE :: firstcall = .TRUE. !! first call?
278	!$OMP THREADPRIVATE(firstcall)
279	REAL(r_std), SAVE, DIMENSION(2) :: lhCH4 !! specific heat of methane transformation (J/kg) (/ 3.1E6, 9.4E6 /)
280	!$OMP THREADPRIVATE(lhCH4)
281	LOGICAL, SAVE :: oxlim !! O2 limitation taken into account
282	!$OMP THREADPRIVATE(oxlim)
283	REAL(r_std), PARAMETER :: refdep = 0.20_r_std !! Depth to compute reference temperature for the growing season (m). WH2000 use 0.50
284	REAL(r_std), PARAMETER :: Tgr = 5. !! Temperature when plant growing starts and this becomes constant
285	INTEGER(i_std) :: month,year,dayno !! current time parameters
286	REAL(r_std) :: scnd
287	REAL(r_std) :: organic_layer_thickness
288	INTEGER(i_std) :: ier, iv, m, jv, iele
289	CHARACTER(80) :: yedoma_map_filename
290	REAL(r_std) :: yedoma_depth, yedoma_cinit_act, yedoma_cinit_slo, yedoma_cinit_pas
291	LOGICAL :: reset_yedoma_carbon
292	LOGICAL, SAVE :: MG_useallCpools = .true. !! Do we allow all three C pools to feed methanogenesis?
293	!$OMP THREADPRIVATE(MG_useallCpools)
294	CHARACTER(LEN=10) :: part_str !! string suffix indicating an index
295	REAL(r_std), SAVE :: max_shum_value = 1.0 !! maximum saturation degree on the thermal axes
296	!$OMP THREADPRIVATE(max_shum_value)
297	REAL(r_std), DIMENSION(kjpindex) :: alt_pftmean, altmax_pftmean, tsurf_pftmean
298	REAL(r_std), DIMENSION (kjpindex,nvm) :: veget_max_bg
299
300	IF (printlev>=3) WRITE(,) 'Entering stomate_soil_carbon_discretization_deep_somcycle'
301
302	!! 0. first call
303	IF ( firstcall ) THEN
304
305	overburden(:)=1.
306	!
307	!Config Key = organic_layer_thickness
308	!Config Desc = The thickness of organic layer
309	!Config Def = 0.0
310	!Config If = OK_SOIL_CARBON_DISCRETIZATION
311	!Config Help = This parameters allows the user to prescibe the organic
312	!Config layer thickness
313	!Config Units = [-]
314	!
315	organic_layer_thickness = 0.
316	CALL getin_p('organic_layer_thickness', organic_layer_thickness)
317	depth_organic_soil(:) = overburden(:)*organic_layer_thickness
318
319	!Config Key = OK_METHANE
320	!Config Desc = Is Methanogenesis and -trophy taken into account?
321	!Config Def = n
322	!Config If = OK_SOIL_CARBON_DISCRETIZATION
323	!Config Help =
324	!Config
325	!Config Units = [flag]
326	!
327	ok_methane = .FALSE.
328	CALL getin_p('OK_METHANE',ok_methane)
329	!
330	!Config Key = HEAT_CO2_ACT
331	!Config Desc = specific heat of soil organic matter oxidation for active carbon (J/kg carbon)
332	!Config Def = 40.0E6
333	!Config If = OK_SOIL_CARBON_DISCRETIZATION
334	!Config Help =
335	!Config
336	!Config Units = [J/Kg]
337	!
338	lhc(iactive) = 40.0e6
339	CALL getin_p('HEAT_CO2_ACT',lhc(iactive))
340	!
341	!Config Key = HEAT_CO2_SLO
342	!Config Desc = specific heat of soil organic matter oxidation for slow
343	!Config carbon pool (J/kg carbon)
344	!Config Def = 30.0E6
345	!Config If = OK_SOIL_CARBON_DISCRETIZATION
346	!Config Help =
347	!Config
348	!Config Units = [J/Kg]
349	!
350	lhc(islow) = 30.0E6
351	CALL getin_p('HEAT_CO2_SLO',lhc(islow))
352	!
353	!Config Key = HEAT_CO2_PAS
354	!Config Desc = specific heat of soil organic matter oxidation for
355	!Config passive carbon pool (J/kg carbon)
356	!Config Def = 10.0E6
357	!Config If = OK_SOIL_CARBON_DISCRETIZATION
358	!Config Help =
359	!Config
360	!Config Units = [J/Kg]
361	!
362	lhc(ipassive) = 10.0e6
363	CALL getin_p('HEAT_CO2_PAS',lhc(ipassive))
364	!
365	!Config Key = TAU_CH4_TROPH
366	!Config Desc = time constant of methanetrophy
367	!Config Def = 432000
368	!Config If = OK_SOIL_CARBON_DISCRETIZATION
369	!Config Help =
370	!Config
371	!Config Units = [s]
372	!
373	tau_CH4troph = 432000
374	CALL getin_p('TAU_CH4_TROPH',tau_CH4troph)
375	!
376	!Config Key = TAU_CH4_GEN_RATIO
377	!Config Desc = time constant of methanogenesis (ratio to that of oxic)
378	!Config Def = 9.0
379	!Config If = OK_SOIL_CARBON_DISCRETIZATION
380	!Config Help =
381	!Config
382	!Config Units = [-]
383	!
384	fbactratio = 9.0
385	CALL getin_p('TAU_CH4_GEN_RATIO',fbactratio)
386	!
387	!Config Key = O2_SEUIL_MGEN
388	!Config Desc = oxygen concentration below which there is anoxy
389	!Config Def = 3.0
390	!Config If = OK_SOIL_CARBON_DISCRETIZATION
391	!Config Help =
392	!Config
393	!Config Units = [g/m3]
394	!
395	O2m = 3.0
396	CALL getin_p('O2_SEUIL_MGEN',O2m)
397	!
398	!Config Key = HEAT_CH4_GEN
399	!Config Desc = specific heat of methanogenesis
400	!Config Def = 0
401	!Config If = OK_SOIL_CARBON_DISCRETIZATION
402	!Config Help =
403	!Config
404	!Config Units = [J/kgC]
405	!
406	lhCH4(1) = 0
407	CALL getin_p('HEAT_CH4_GEN',lhCH4(1))
408	!
409	!Config Key = HEAT_CH4_TROPH
410	!Config Desc = specific heat of methanotrophy
411	!Config Def = 0
412	!Config If = OK_SOIL_CARBON_DISCRETIZATION
413	!Config Help =
414	!Config
415	!Config Units = [J/kgC]
416	!
417	lhCH4(2) = 0
418	CALL getin_p('HEAT_CH4_TROPH',lhCH4(2))
419	!
420	!Config Key = O2_LIMIT
421	!Config Desc = O2 limitation taken into account
422	!Config Def = n
423	!Config If = OK_SOIL_CARBON_DISCRETIZATION
424	!Config Help =
425	!Config
426	!Config Units = [flag]
427	!
428	oxlim=.FALSE.
429	CALL getin_p('O2_LIMIT',oxlim)
430	!
431	!Config Key = cryoturbate
432	!Config Desc = Do we allow for cyoturbation?
433	!Config Def = y
434	!Config If = OK_SOIL_CARBON_DISCRETIZATION
435	!Config Help =
436	!Config
437	!Config Units = [flag]
438	!
439	ok_cryoturb=.TRUE.
440	CALL getin_p('cryoturbate',ok_cryoturb)
441	!
442	!Config Key = cryoturbation_diff_k_in
443	!Config Desc = diffusion constant for cryoturbation
444	!Config Def = 0.001
445	!Config If = OK_SOIL_CARBON_DISCRETIZATION
446	!Config Help =
447	!Config
448	!Config Units = [m2/year]
449	!
450	cryoturbation_diff_k_in = .001
451	CALL getin_p('cryoturbation_diff_k',cryoturbation_diff_k_in)
452	!
453	!Config Key = bioturbation_diff_k_in
454	!Config Desc = diffusion constant for bioturbation
455	!Config Def = 0.0
456	!Config If = OK_SOIL_CARBON_DISCRETIZATION
457	!Config Help =
458	!Config
459	!Config Units = [m2/year]
460	!
461	bioturbation_diff_k_in = 0.0001
462	CALL getin_p('bioturbation_diff_k',bioturbation_diff_k_in)
463	!
464	!Config Key = MG_useallCpools
465	!Config Desc = Do we allow all three C pools to feed methanogenesis?
466	!Config Def = y
467	!Config If = OK_SOIL_CARBON_DISCRETIZATION
468	!Config Help =
469	!Config
470	!Config Units = [flag]
471	!
472	MG_useallCpools = .TRUE.
473	CALL getin_p('MG_useallCpools', MG_useallCpools)
474	!
475	!Config Key = max_shum_value
476	!Config Desc = maximum saturation degree on the thermal axes
477	!Config Def = 1
478	!Config If = OK_SOIL_CARBON_DISCRETIZATION
479	!Config Help =
480	!Config
481	!Config Units = [-]
482	!
483	max_shum_value=1.0
484	CALL getin_p('max_shum_value',max_shum_value)
485	hslong(:,:,:) = MAX(MIN(hslong_in(:,:,:),max_shum_value),zero)
486	!
487
488	!! Arrays allocations
489
490	ALLOCATE (veget_mask_2d(kjpindex,nvm),stat=ier)
491	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
492	'Pb in alloc for veget_mask_2d','','')
493
494	ALLOCATE(lalo_global(kjpindex,2),stat=ier)
495	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
496	'Pb in alloc for lalo_global','','')
497
498	ALLOCATE (alt(kjpindex,nvm),stat=ier)
499	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
500	'Pb in alloc for alt','','')
501
502	ALLOCATE (altmax_lastyear(kjpindex,nvm),stat=ier)
503	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
504	'Pb in alloc for altmax_lastyear','','')
505
506	ALLOCATE (alt_ind(kjpindex,nvm),stat=ier)
507	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
508	'Pb in alloc for alt_ind','','')
509
510	ALLOCATE (altmax_ind(kjpindex,nvm),stat=ier)
511	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
512	'Pb in alloc for altmax_ind','','')
513
514	ALLOCATE (altmax_ind_lastyear(kjpindex,nvm),stat=ier)
515	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
516	'Pb in alloc for altmax_ind_lastyear','','')
517
518	ALLOCATE (z_root(kjpindex,nvm),stat=ier)
519	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
520	'Pb in alloc for z_root','','')
521
522	ALLOCATE (rootlev(kjpindex,nvm),stat=ier)
523	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
524	'Pb in alloc for rootlev','','')
525
526	ALLOCATE (heights_snow(kjpindex,nvm),stat=ier)
527	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
528	'Pb in alloc for heights_snow','','')
529
530	ALLOCATE (zf_soil(0:ngrnd),stat=ier)
531	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
532	'Pb in alloc for zf_soil','','')
533
534	ALLOCATE (zi_soil(ngrnd),stat=ier)
535	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
536	'Pb in alloc for zi_soil','','')
537
538	ALLOCATE (zf_snow(kjpindex,0:nsnow,nvm),stat=ier)
539	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
540	'Pb in alloc for zf_snow','','')
541
542	ALLOCATE (zi_snow(kjpindex,nsnow,nvm),stat=ier)
543	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
544	'Pb in alloc for zi_snow','','')
545
546	ALLOCATE (zf_snow_nopftdim(kjpindex,0:nsnow),stat=ier)
547	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
548	'Pb in alloc for zf_snow_nopftdim','','')
549
550	ALLOCATE (zi_snow_nopftdim(kjpindex,nsnow),stat=ier)
551	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
552	'Pb in alloc for zi_snow_nopftdim','','')
553
554	ALLOCATE (airvol_soil(kjpindex,ngrnd,nvm),stat=ier)
555	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
556	'Pb in alloc for airvol_soil','','')
557
558	ALLOCATE (totporO2_soil(kjpindex,ngrnd,nvm),stat=ier)
559	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
560	'Pb in alloc for totporO2_soil','','')
561
562	ALLOCATE (totporCH4_soil(kjpindex,ngrnd,nvm),stat=ier)
563	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
564	'Pb in alloc for totporCH4_soil','','')
565
566	ALLOCATE (conduct_soil(kjpindex,ngrnd,nvm),stat=ier)
567	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
568	'Pb in alloc for conduct_soil','','')
569
570	ALLOCATE (diffO2_soil(kjpindex,ngrnd,nvm),stat=ier)
571	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
572	'Pb in alloc for diffO2_soil','','')
573
574	ALLOCATE (diffCH4_soil(kjpindex,ngrnd,nvm),stat=ier)
575	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
576	'Pb in alloc for diffCH4_soil','','')
577
578	ALLOCATE (airvol_snow(kjpindex,nsnow,nvm),stat=ier)
579	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
580	'Pb in alloc for airvol_snow','','')
581
582	ALLOCATE (totporO2_snow(kjpindex,nsnow,nvm),stat=ier)
583	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
584	'Pb in alloc for totporO2_snow','','')
585
586	ALLOCATE (totporCH4_snow(kjpindex,nsnow,nvm),stat=ier)
587	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
588	'Pb in alloc for totporCH4_snow','','')
589
590	ALLOCATE (conduct_snow(kjpindex,nsnow,nvm),stat=ier)
591	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
592	'Pb in alloc for conduct_snow','','')
593
594	ALLOCATE (diffO2_snow(kjpindex,nsnow,nvm),stat=ier)
595	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
596	'Pb in alloc for diffO2_snow','','')
597
598	ALLOCATE (diffCH4_snow(kjpindex,nsnow,nvm),stat=ier)
599	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
600	'Pb in alloc for diffCH4_snow','','')
601
602	ALLOCATE (deepSOM_pftmean(kjpindex,ngrnd,ncarb,nelements),stat=ier)
603	IF (ier /= 0) CALL ipslerr_p(3,'stomate_soil_carbon_discretization_deep_somcycle', &
604	'Pb in alloc for deepSOM_pftmean','','')
605
606	!! assign values for arrays
607	yr_len = NINT(one_year)
608
609	veget_max_bg(:,2:nvm) = veget_max(:,2:nvm)
610	veget_max_bg(:,1) = MAX((un - SUM(veget_max(:,2:nvm), 2)), zero)
611	!! veget_mask_2d(:,:) = veget_max_bg .GT. EPSILON(zero)
612	!! WHERE( ALL((.NOT. veget_mask_2d(:,:)), dim=2) )
613	!! veget_mask_2d(:,1) = .TRUE.
614	!! END WHERE
615	veget_mask_2d(:,:) = .TRUE.
616
617	lalo_global(:,:) = lalo(:,:)
618	alt(:,:) = 0
619	altmax_lastyear(:,:) = 0
620	alt_ind(:,:) = 0
621	altmax_ind(:,:) = 0
622	altmax_ind_lastyear(:,:) = 0
623	z_root(:,:) = 0
624	rootlev(:,:) = 0
625	febul(:,:) = 0
626	flupmt(:,:) = 0
627	! make sure gas concentrations where not defined by veget_mask are equal
628	!to initial conditions
629	DO iv = 1, ngrnd
630	WHERE ( .NOT. veget_mask_2d(:,:) )
631	O2_soil(:,iv,:) = O2_init_conc
632	CH4_soil(:,iv,:) = CH4_init_conc
633	END WHERE
634	END DO
635	DO iv = 1, nsnow
636	WHERE ( .NOT. veget_mask_2d(:,:) )
637	O2_snow(:,iv,:) = O2_surf
638	CH4_snow(:,iv,:) = CH4_surf
639	END WHERE
640	END DO
641
642	heights_snow(:,:) = zero
643	zf_soil(:) = zero
644	zi_soil(:) = zero
645	zf_snow(:,:,:) = zero
646	zi_snow(:,:,:) = zero
647	zf_snow_nopftdim(:,:) = zero
648	zi_snow_nopftdim(:,:) = zero
649	airvol_soil(:,:,:) = zero
650	totporO2_soil(:,:,:) = zero
651	totporCH4_soil(:,:,:) = zero
652	conduct_soil(:,:,:) = zero
653	diffO2_soil(:,:,:) = zero
654	diffCH4_soil(:,:,:) = zero
655	airvol_snow(:,:,:) = zero
656	totporO2_snow(:,:,:) = zero
657	totporCH4_snow(:,:,:) = zero
658	conduct_snow(:,:,:) = zero
659	diffO2_snow(:,:,:) = zero
660	diffCH4_snow(:,:,:) = zero
661
662	! get snow and soil levels
663	DO iv = 1, nvm
664	heights_snow(:,iv) = SUM(snowdz(:,1:nsnow), 2)
665	ENDDO
666	! Calculating intermediate and full depths for snow
667	call snowlevels (kjpindex, snowdz, zi_snow, zf_snow, veget_max_bg)
668
669	! here we need to put the shallow and deep soil levels together to make the complete soil levels.
670	! This requires pulling in the indices from thermosoil and deepsoil_freeze.
671	zi_soil(:) = znt(:)
672	zf_soil(1:ngrnd) = zlt(:)
673	zf_soil(0) = 0.
674
675	! allocate arrays for gas diffusion !
676	! get diffusion coefficients: heat capacity,
677	! conductivity, and oxygen diffusivity
678
679	CALL get_gasdiff (kjpindex,hslong,tprof,snow,airvol_snow, &
680	totporO2_snow,totporCH4_snow,diffO2_snow,diffCH4_snow, &
681	airvol_soil,totporO2_soil,totporCH4_soil,diffO2_soil,diffCH4_soil, snowrho)
682
683	!
684	! initialize soil temperature calculation
685	!
686	CALL soil_gasdiff_main (kjpindex,time_step,index,'initialize', &
687	pb,tsurf,tprof,diffO2_snow,diffCH4_snow, &
688	totporO2_snow,totporCH4_snow,O2_snow,CH4_snow,diffO2_soil,diffCH4_soil, &
689	totporO2_soil,totporCH4_soil,O2_soil,CH4_soil, zi_snow, zf_snow)
690
691	!
692	! calculate the coefficients
693	!
694	CALL soil_gasdiff_main (kjpindex,time_step,index,'coefficients', &
695	pb,tsurf,tprof,diffO2_snow,diffCH4_snow, &
696	totporO2_snow,totporCH4_snow,O2_snow,CH4_snow,diffO2_soil,diffCH4_soil, &
697	totporO2_soil,totporCH4_soil,O2_soil,CH4_soil, zi_snow, zf_snow)
698
699
700
701	CALL itau2ymds(itau, time_step, year, month, dayno, scnd)
702	dayno = (month-1)*30 + dayno
703	CALL altcalc (kjpindex, time_step, dayno, scnd, tprof, zi_soil, alt, alt_ind, altmax, altmax_ind, &
704	altmax_lastyear, altmax_ind_lastyear)
705
706	IF (printlev>=3 ) THEN
707	WRITE(,) 'stomate_soil_carbon_discretization_deep_somcycle: finished firstcall calcs'
708	ENDIF
709
710	! reset
711	!
712	!Config Key = reset_yedoma_carbon
713	!Config Desc = Do we reset carbon concentrations for yedoma region?
714	!Config Def = n
715	!Config If = OK_SOIL_CARBON_DISCRETIZATION
716	!Config Help =
717	!Config
718	!Config Units = [flag]
719	!
720	reset_yedoma_carbon = .false.
721	CALL getin_p('reset_yedoma_carbon',reset_yedoma_carbon)
722
723	IF (reset_yedoma_carbon) THEN
724	yedoma_map_filename = 'NONE'
725	yedoma_depth = zero
726	yedoma_cinit_act = zero
727	yedoma_cinit_slo = zero
728	yedoma_cinit_pas = zero
729	!
730	!Config Key = yedoma_map_filename
731	!Config Desc = The filename for yedoma map
732	!Config Def = yedoma_map.nc
733	!Config If = OK_SOIL_CARBON_DISCRETIZATION
734	!Config Help =
735	!Config
736	!Config Units = []
737	!
738	CALL getin_p('yedoma_map_filename', yedoma_map_filename)
739	!
740	!Config Key = yedoma_depth
741	!Config Desc = The depth for soil carbon in yedoma
742	!Config Def = 20
743	!Config If = OK_SOIL_CARBON_DISCRETIZATION
744	!Config Help =
745	!Config
746	!Config Units = [m]
747	!
748	CALL getin_p('yedoma_depth', yedoma_depth)
749	!
750	!Config Key = deepC_a_init
751	!Config Desc = Carbon concentration for active soil C pool in yedoma
752	!Config Def = 1790.1
753	!Config If = OK_SOIL_CARBON_DISCRETIZATION
754	!Config Help =
755	!Config
756	!Config Units = []
757	!
758	CALL getin_p('deepC_a_init', yedoma_cinit_act)
759	!
760	!Config Key = deepC_s_init
761	!Config Desc = Carbon concentration for slow soil C pool in yedoma
762	!Config Def = 14360.8
763	!Config If = OK_SOIL_CARBON_DISCRETIZATION
764	!Config Help =
765	!Config
766	!Config Units = []
767	!
768	CALL getin_p('deepC_s_init', yedoma_cinit_slo)
769	!
770	!Config Key = deepC_p_init
771	!Config Desc = Carbon concentration for passive soil C pool in yedoma
772	!Config Def = 1436
773	!Config If = OK_SOIL_CARBON_DISCRETIZATION
774	!Config Help =
775	!Config
776	!Config Units = []
777	!
778	CALL getin_p('deepC_p_init', yedoma_cinit_pas)
779	! intialize the yedoma carbon stocks
780	CALL initialize_yedoma_carbonstocks(kjpindex, lalo, deepSOM_a, deepSOM_s, deepSOM_p, &
781	yedoma_map_filename, yedoma_depth, yedoma_cinit_act,yedoma_cinit_slo, yedoma_cinit_pas, altmax_ind)
782	ENDIF
783
784
785	ENDIF ! firstcall
786
787	! Prepare values for arrays
788	veget_max_bg(:,2:nvm) = veget_max(:,2:nvm)
789	veget_max_bg(:,1) = MAX((un - SUM(veget_max(:,2:nvm), 2)), zero)
790
791	IF ( ANY(rootlev(:,:) .GT. ngrnd) ) THEN
792	WRITE(,) 'problems with rootlev:', rootlev
793	STOP
794	ENDIF
795
796	DO iv = 1, nvm
797	heights_snow(:,iv) = SUM(snowdz(:,1:nsnow), 2)
798	ENDDO
799	!
800	! define initial CH4 value (before the time step)
801	CH4ini_soil(:,:,:) = CH4_soil(:,:,:)
802
803	! apply maximum soil wetness criteria to prevent soils from turning to wetlands where they aren't supposed to
804	hslong(:,:,:) = MAX(MIN(hslong_in(:,:,:),max_shum_value),zero)
805
806
807	! update the gas profiles
808	!
809	CALL soil_gasdiff_main (kjpindex, time_step, index, 'diffuse', &
810	pb,tsurf,tprof,diffO2_snow,diffCH4_snow, &
811	totporO2_snow,totporCH4_snow,O2_snow,CH4_snow,diffO2_soil,diffCH4_soil, &
812	totporO2_soil,totporCH4_soil,O2_soil,CH4_soil, zi_snow, zf_snow)
813
814	! get new snow levels and interpolate gases on these levels
815	!
816	CALL snow_interpol (kjpindex,O2_snow, CH4_snow, zi_snow, zf_snow, veget_max_bg, snowdz)
817
818	! Compute active layer thickness
819	CALL itau2ymds(itau, time_step, year, month, dayno, scnd)
820	dayno = (month-1)*30 + dayno
821
822	CALL altcalc (kjpindex, time_step, dayno, scnd, tprof, zi_soil, alt, alt_ind, altmax, altmax_ind, &
823	altmax_lastyear, altmax_ind_lastyear)
824
825	! list pft-mean alt and altmax for debugging purposes
826	IF (printlev>=3) THEN
827	alt_pftmean(:) = 0.
828	altmax_pftmean(:) = 0.
829	tsurf_pftmean(:) = 0.
830	DO iv = 1, nvm
831	WHERE ( veget_mask_2d(:,iv) )
832	alt_pftmean(:) = alt_pftmean(:) + alt(:,iv)*veget_max_bg(:,iv)
833	altmax_pftmean(:) = altmax_pftmean(:) + altmax(:,iv)*veget_max_bg(:,iv)
834	tsurf_pftmean(:) = tsurf_pftmean(:) + tprof(:,1,iv)*veget_max_bg(:,iv)
835	END WHERE
836	END DO
837	END IF
838
839	! Make sure the rooting depth is within the active layer
840
841	!need to sort out the rooting depth, by each STOMATE PFT
842	WHERE ( altmax_lastyear(:,:) .LT. z_root_max .and. veget_mask_2d(:,:) )
843	z_root(:,:) = altmax_lastyear(:,:)
844	rootlev(:,:) = altmax_ind_lastyear(:,:)
845	ELSEWHERE ( veget_mask_2d(:,:) )
846	z_root(:,:) = z_root_max
847	rootlev(:,:) = altmax_ind_lastyear(:,:)
848	ENDWHERE
849
850	!
851	! Carbon input into the soil
852	!
853	CALL sominput(kjpindex,time_step,itau*time_step,tprof,tsurf,hslong,dayno,z_root,altmax_lastyear, &
854	deepSOM_a, deepSOM_s, deepSOM_p, som_in, veget_max_bg, rprof)
855	!
856
857	CALL permafrost_decomp (kjpindex, time_step, tprof, airvol_soil, &
858	oxlim, tau_CH4troph, ok_methane, fbactratio, O2m, &
859	totporO2_soil, totporCH4_soil, hslong, clay, &
860	deepSOM_a, deepSOM_s, deepSOM_p, &
861	deltaCH4g, deltaCH4, deltaSOM1_a, deltaSOM1_s, deltaSOM1_p, deltaSOM2, &
862	deltaSOM3, O2_soil, CH4_soil, fbact, CN_target, MG_useallCpools)
863
864	! calculate coefficients for cryoturbation calculation
865	IF (ok_cryoturb) CALL cryoturbate(kjpindex, time_step, dayno, altmax_ind_lastyear, deepSOM_a, deepSOM_s, deepSOM_p, &
866	'coefficients', cryoturbation_diff_k_in/(one_dayone_year),bioturbation_diff_k_in/(one_dayone_year), &
867	altmax_lastyear, fixed_cryoturbation_depth)
868
869	IF (ok_cryoturb) CALL cryoturbate(kjpindex, time_step, dayno, altmax_ind_lastyear, deepSOM_a, deepSOM_s, deepSOM_p, &
870	'diffuse', cryoturbation_diff_k_in/(one_dayone_year), bioturbation_diff_k_in/(one_dayone_year), &
871	altmax_lastyear, fixed_cryoturbation_depth)
872
873	DO ip = 1, kjpindex
874	DO iv = 1, nvm
875	IF ( veget_mask_2d(ip,iv) ) THEN
876	! oxic decomposition
877	heat_Zimov(ip,:,iv) = lhc(iactive)1.E-3deltaSOM1_a(ip,:,iv,icarbon) + &
878	lhc(islow)1.E-3deltaSOM1_s(ip,:,iv,icarbon) + &
879	lhc(ipassive)1.E-3deltaSOM1_p(ip,:,iv,icarbon)
880	!
881	! methanogenesis
882	heat_Zimov(ip,:,iv) = heat_Zimov(ip,:,iv) + lhCH4(1)1.E-3deltaSOM2(ip,:,iv,icarbon)
883	!
884	! methanotrophy
885	heat_Zimov(ip,:,iv) = heat_Zimov(ip,:,iv) + lhCH4(2)1.E-3deltaCH4(ip,:,iv) * &
886	totporCH4_soil(ip,:,iv)
887	!
888	heat_Zimov(ip,:,iv) = heat_Zimov(ip,:,iv)/time_step
889
890	!
891	fluxCH4(ip,iv) = zero
892	ELSE
893	heat_Zimov(ip,:,iv) = zero
894	fluxCH4(ip,iv) = zero
895	ENDIF
896	ENDDO
897	ENDDO
898
899	IF ( .NOT. firstcall) THEN
900	!
901	! Plant-mediated CH4 transport
902	!
903	CALL traMplan(CH4_soil,O2_soil,kjpindex,time_step,totporCH4_soil,totporO2_soil,z_root, &
904	rootlev,Tgr,Tref,hslong,flupmt, &
905	refdep, zi_soil, tprof)
906	! flupmt=zero
907	!
908	! CH4 ebullition
909	!
910
911	CALL ebullition (kjpindex,time_step,tprof,totporCH4_soil,hslong,CH4_soil,febul)
912
913	!
914	ENDIF
915
916	!
917	MT(:,:)=zero
918	MG(:,:)=zero
919	CH4i(:,:)=zero
920	CH4ii(:,:)=zero
921	dC1i(:,:)=zero
922	dCi(:,:)=zero
923	!
924	DO ip = 1, kjpindex
925	DO iv = 1, nvm
926	IF ( veget_mask_2d(ip,iv) ) THEN
927	DO il=1,ngrnd
928	MT(ip,iv) = MT(ip,iv) + deltaCH4(ip,il,iv)totporCH4_soil(ip,il,iv) &
929	( zf_soil(il) - zf_soil(il-1) )
930	MG(ip,iv) = MG(ip,iv) + deltaCH4g(ip,il,iv)totporCH4_soil(ip,il,iv) &
931	( zf_soil(il) - zf_soil(il-1) )
932	CH4i(ip,iv) = CH4i(ip,iv) + CH4_soil(ip,il,iv)totporCH4_soil(ip,il,iv) &
933	(zf_soil(il)-zf_soil(il-1))
934	CH4ii(ip,iv) = CH4ii(ip,iv) + &
935	CH4ini_soil(ip,il,iv)totporCH4_soil(ip,il,iv) &
936	(zf_soil(il)-zf_soil(il-1))
937	dC1i(ip,iv) = dC1i(ip,iv) + &
938	(deltaSOM1_a(ip,il,iv,icarbon)+deltaSOM1_s(ip,il,iv,icarbon)+deltaSOM1_p(ip,il,iv,icarbon)) * &
939	( zf_soil(il) - zf_soil(il-1) )
940	n_mineralisation(ip,iv)=n_mineralisation(ip,iv)+ &
941	(deltaSOM1_a(ip,il,iv,initrogen)+deltaSOM1_s(ip,il,iv,initrogen)+deltaSOM1_p(ip,il,iv,initrogen))* &
942	( zf_soil(il) - zf_soil(il-1) )
943	!MERGE: This dCi variables might not be necessary
944	dCi(ip,iv) = dCi(ip,iv) + &
945	(deepSOM_a(ip,il,iv,icarbon) + deepSOM_s(ip,il,iv,icarbon) + deepSOM_p(ip,il,iv,icarbon)) * &
946	( zf_soil(il) - zf_soil(il-1) )
947	END DO
948	ENDIF
949	ENDDO
950	ENDDO
951
952	!
953	!
954
955	DO ip = 1, kjpindex
956	! Total CH4 flux
957	sfluxCH4_deep(ip) = SUM(veget_max_bg(ip,:)*( CH4ii(ip,:)-CH4i(ip,:)+MG(ip,:)-MT(ip,:) ))/time_step
958	! TotalCO2 flux
959	sfluxCO2_deep(ip) = SUM(veget_max_bg(ip,:)( dC1i(ip,:) + MT(ip,:)(12./16.) ) )/time_step
960	END DO
961
962	resp_hetero_soil(:,:) = ( dC1i(:,:) + MT(:,:)(12./16.) ) one_day/time_step
963	sfluxCH4(:,:) = ( CH4ii(:,:)-CH4i(:,:)+MG(:,:)-MT(:,:) ) *one_day/time_step
964
965
966	! calculate the coefficients for the next timestep:
967	!
968	! get diffusion coefficients: heat capacity,
969	! conductivity, and oxygen diffusivity
970	!
971	CALL get_gasdiff (kjpindex,hslong,tprof,snow,airvol_snow, &
972	totporO2_snow,totporCH4_snow,diffO2_snow,diffCH4_snow, &
973	airvol_soil,totporO2_soil,totporCH4_soil,diffO2_soil,diffCH4_soil, snowrho)
974
975	!
976	! calculate the coefficients for the next time step
977	!
978	CALL soil_gasdiff_main (kjpindex,time_step,index,'coefficients', &
979	pb,tsurf,tprof,diffO2_snow,diffCH4_snow, &
980	totporO2_snow,totporCH4_snow,O2_snow,CH4_snow,diffO2_soil,diffCH4_soil, &
981	totporO2_soil,totporCH4_soil,O2_soil,CH4_soil, zi_snow, zf_snow)
982
983	call calc_vert_int_som(kjpindex, deepSOM_a, deepSOM_s, deepSOM_p, som, som_surf, zf_soil)
984	IF (printlev>=3) WRITE(,) 'after calc_vert_int_som'
985
986	! XIOS history output
987	IF ( .NOT. soilc_isspinup ) THEN
988
989	!CALL xios_orchidee_send_field ('fluxCH4', sfluxCH4)
990	!CALL xios_orchidee_send_field ('febul', (febul*one_day))
991	!CALL xios_orchidee_send_field ('flupmt', (flupmt*one_day))
992	CALL xios_orchidee_send_field ( 'alt', alt )
993	CALL xios_orchidee_send_field ( 'altmax', altmax)
994	!CALL xios_orchidee_send_field ( 'sfluxCH4_deep', sfluxCH4_deep)
995	!CALL xios_orchidee_send_field ( 'sfluxCO2_deep', sfluxCO2_deep)
996	CALL xios_orchidee_send_field ( 'pb', pb)
997
998	!CALL xios_orchidee_send_field ( 'O2_soil', O2_soil)
999	!CALL xios_orchidee_send_field ( 'CH4_soil', CH4_soil)
1000	!CALL xios_orchidee_send_field ('O2_snow', O2_snow)
1001	!CALL xios_orchidee_send_field ( 'CH4_snow', CH4_snow)
1002
1003	!CALL xios_orchidee_send_field ( 'deltaCH4g', deltaCH4g)
1004	!CALL xios_orchidee_send_field ( 'deltaCH4', deltaCH4)
1005
1006	!CALL xios_orchidee_send_field ( 'heat_Zimov', heat_Zimov)
1007
1008	!CALL xios_orchidee_send_field ( 'totporO2_soil', totporO2_soil)
1009	!CALL xios_orchidee_send_field ( 'diffO2_soil', diffO2_soil)
1010	!CALL xios_orchidee_send_field ( 'alphaO2_soil', alphaO2_soil)
1011	!CALL xios_orchidee_send_field ( 'betaO2_soil', betaO2_soil)
1012
1013	!CALL xios_orchidee_send_field ( 'totporCH4_soil', totporCH4_soil)
1014	!CALL xios_orchidee_send_field ( 'diffCH4_soil', diffCH4_soil)
1015	!CALL xios_orchidee_send_field ('alphaCH4_soil', alphaCH4_soil)
1016	!CALL xios_orchidee_send_field ( 'betaCH4_soil', betaCH4_soil)
1017
1018	ENDIF
1019
1020	IF (printlev>=3) WRITE(,) 'cdk: leaving stomate_soil_carbon_discretization_deep_somcycle'
1021
1022	IF ( firstcall ) firstcall = .FALSE.
1023
1024
1025	END SUBROUTINE stomate_soil_carbon_discretization_deep_somcycle
1026
1027	!!
1028	!================================================================================================================================
1029	!! SUBROUTINE : altcalc
1030	!!
1031	!>\BRIEF This routine calculate active layer thickness
1032	!!
1033	!! DESCRIPTION :
1034	!!
1035	!! RECENT CHANGE(S) : None
1036	!!
1037	!! MAIN OUTPUT VARIABLE(S) : alt
1038	!!
1039	!! REFERENCE(S) : None
1040	!!
1041	!! FLOWCHART11 : None
1042	!! \n
1043	!_
1044	!================================================================================================================================
1045	SUBROUTINE altcalc (kjpindex,time_step,dayno,scnd, temp, zprof, alt, alt_ind, altmax, altmax_ind, &
1046	altmax_lastyear, altmax_ind_lastyear)
1047
1048	!! 0. Variable and parameter declaration
1049
1050	!! 0.1 Input variables
1051
1052	INTEGER(i_std), INTENT(in) :: kjpindex
1053	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
1054	INTEGER(i_std), INTENT(in) :: dayno !! number of the day in the current year
1055	REAL(r_std), INTENT(in) :: scnd !! model time & time step
1056	REAL(r_std), DIMENSION(kjpindex,ngrnd, nvm), INTENT(in) :: temp !! soil temperature
1057	REAL(r_std), DIMENSION(ngrnd), INTENT(in) :: zprof !! soil depths (m)
1058
1059	!! 0.2 Output variables
1060
1061	REAL(r_std), DIMENSION(kjpindex, nvm), INTENT(out) :: alt !! active layer thickness
1062	INTEGER, DIMENSION(kjpindex, nvm), INTENT(out) :: alt_ind !! active layer index
1063
1064	!! 0.3 Modified variables
1065
1066	REAL(r_std), DIMENSION(kjpindex, nvm),INTENT(inout) :: altmax_lastyear !! Maximum active-layer thickness
1067	REAL(r_std), DIMENSION(kjpindex, nvm),INTENT(inout) :: altmax !! Maximum active-layer thickness
1068	INTEGER(i_std), DIMENSION(kjpindex, nvm),INTENT(inout) :: altmax_ind !! Maximum over the year active-layer index
1069	INTEGER(i_std), DIMENSION(kjpindex, nvm),INTENT(inout) :: altmax_ind_lastyear !! Maximum over the year active-layer index
1070
1071	!! 0.4 Local variables
1072
1073	INTEGER :: ix,iz,il,iv !! grid indices
1074	LOGICAL, SAVE :: firstcall = .TRUE.
1075	!$OMP THREADPRIVATE(firstcall)
1076	INTEGER, save :: tcounter
1077	INTEGER(i_std), SAVE :: id, id2
1078	!$OMP THREADPRIVATE(id)
1079	!$OMP THREADPRIVATE(id2)
1080	LOGICAL, SAVE :: check = .FALSE.
1081	!$OMP THREADPRIVATE(check)
1082	LOGICAL, SAVE :: newaltcalc = .FALSE.
1083	!$OMP THREADPRIVATE(newaltcalc)
1084	LOGICAL, DIMENSION(kjpindex,nvm) :: inalt, bottomlevelthawed
1085	CHARACTER(LEN=16) :: buf
1086	INTEGER :: lev
1087
1088
1089	IF ( firstcall ) THEN
1090
1091	! calculate altmax_ind from altmax
1092	altmax_ind(:,:) = 0
1093	DO ix = 1, kjpindex
1094	DO iv = 1, nvm
1095	IF ( veget_mask_2d(ix,iv) ) THEN
1096	DO il=1,ngrnd
1097	IF ( altmax(ix,iv) .GE. zprof(il) ) THEN
1098	altmax_ind(ix,iv) = altmax_ind(ix,iv) + 1
1099	END IF
1100	END DO
1101	END IF
1102	END DO
1103	END DO
1104	altmax_lastyear(:,:) = altmax(:,:)
1105	altmax_ind_lastyear(:,:) = altmax_ind(:,:)
1106	firstcall = .FALSE.
1107
1108	!Config Key = newaltcalc
1109	!Config Desc = calculate alt ?
1110	!Config Def = n
1111	!Config If = OK_SOIL_CARBON_DISCRETIZATION
1112	!Config Help =
1113	!Config Unit = [flag]
1114	CALL getin_p('newaltcalc', newaltcalc)
1115
1116	ELSE
1117	! all other timesteps
1118	IF ( .NOT. newaltcalc ) THEN
1119	DO ix = 1, kjpindex
1120	DO iv = 1, nvm
1121	IF ( veget_mask_2d(ix,iv) ) THEN
1122	iz = 1
1123	DO WHILE( temp(ix,iz,iv) > ZeroCelsius .AND. iz < ngrnd )
1124	iz = iz + 1
1125	END DO
1126	IF( iz == 1 ) THEN
1127	! it means that all is frozen
1128	alt(ix,iv) = zero
1129	ELSE
1130	alt(ix,iv) = zprof(iz-1)
1131	END IF
1132	alt_ind(ix,iv) = iz-1
1133	END IF
1134	END DO
1135	END DO
1136	ELSE
1137	! initialize for pfts that don't exist
1138	alt(:,:) = zprof(ngrnd)
1139	bottomlevelthawed(:,:) = .FALSE.
1140	! start from bottom and work up instead
1141	WHERE (temp(:,ngrnd,:) > ZeroCelsius )
1142	bottomlevelthawed(:,:) = .TRUE.
1143	alt(:,:) = zprof(ngrnd)
1144	alt_ind(:,:) = ngrnd
1145	END WHERE
1146	inalt(:,:) = .FALSE.
1147	DO iz = 1, ngrnd - 1
1148	lev = ngrnd - iz
1149	WHERE ( temp(:,lev,:) > ZeroCelsius .AND. .NOT. inalt(:,:) .AND. .NOT. bottomlevelthawed(:,:) )
1150	inalt(:,:) = .TRUE.
1151	alt(:,:) = zprof(lev)
1152	alt_ind(:,:) = lev
1153	ELSEWHERE ( temp(:,lev,:) <= ZeroCelsius .AND. inalt(:,:) .AND. .NOT. bottomlevelthawed(:,:) )
1154	inalt(:,:) = .FALSE.
1155	END WHERE
1156	END DO
1157	WHERE ( .NOT. inalt .AND. .NOT. bottomlevelthawed(:,:) )
1158	alt(:,:) = zero
1159	alt_ind(:,:) = 0
1160	END WHERE
1161	ENDIF
1162
1163	! debug
1164	IF ( check ) THEN
1165	IF (ANY(alt(:,:) .GT. zprof(ngrnd))) THEN
1166	WRITE(,) 'error: alt greater than soil depth.'
1167	ENDIF
1168	ENDIF
1169
1170	! Maximum over the year active layer thickness
1171	WHERE ( ( alt(:,:) .GT. altmax(:,:) ) .AND. veget_mask_2d(:,:) )
1172	altmax(:,:) = alt(:,:)
1173	altmax_ind(:,:) = alt_ind(:,:)
1174	ENDWHERE
1175
1176	IF ( .NOT. soilc_isspinup ) THEN
1177	! do it on the second timestep, that way when we are writing restart files it is not done before that!
1178	! now we are doing daily permafrost calcs, so just run it on the second day.
1179	IF ( ( dayno .EQ. 2) ) THEN
1180	! Reinitialize ALT_max
1181	altmax_lastyear(:,:) = altmax(:,:)
1182	altmax_ind_lastyear(:,:) = altmax_ind(:,:)
1183	altmax(:,:) = alt(:,:)
1184	altmax_ind(:,:) = alt_ind(:,:)
1185	END IF
1186	ELSE
1187
1188	! for spinup, best to set altmax_lastyear to altmax, and not boter to reset since every year is the same,
1189	! and if you try to do so, it doesn't work properly -- 06 may 2010
1190	altmax_lastyear(:,:) = altmax(:,:)
1191	altmax_ind_lastyear(:,:) = altmax_ind(:,:)
1192	END IF
1193	END IF
1194
1195	IF (printlev>=3) WRITE(,) 'leaving altcalc'
1196	END SUBROUTINE altcalc
1197
1198	!!
1199	!================================================================================================================================
1200	!! SUBROUTINE : soil_gasdiff_main
1201	!!
1202	!>\BRIEF This routine calculate oxygen and methane in the snow/soil medium
1203	!!
1204	!! DESCRIPTION :
1205	!!
1206	!! RECENT CHANGE(S) : None
1207	!!
1208	!! MAIN OUTPUT VARIABLE(S) :
1209	!!
1210	!! REFERENCE(S) : None
1211	!!
1212	!! FLOWCHART11 : None
1213	!! \n
1214	!_
1215	!================================================================================================================================
1216	SUBROUTINE soil_gasdiff_main( kjpindex,time_step,index,cur_action, &
1217	psol,tsurf,tprof,diffO2_snow,diffCH4_snow, &
1218	totporO2_snow,totporCH4_snow,O2_snow,CH4_snow,diffO2_soil,diffCH4_soil, &
1219	totporO2_soil,totporCH4_soil,O2_soil,CH4_soil, zi_snow, zf_snow)
1220
1221	!! 0. Variable and parameter declaration
1222
1223	!! 0.1 Input variables
1224
1225	INTEGER(i_std), INTENT(in) :: kjpindex !! number of grid points
1226	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
1227	CHARACTER(LEN=*), INTENT(in) :: cur_action !! what to do
1228	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: psol !! surface pressure (Pa)
1229	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: tsurf !! Surface temperature (K)
1230	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: tprof !! Soil temperature (K)
1231	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: diffO2_snow !! oxygen diffusivity (m**2/s)
1232	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: diffCH4_snow !! methane diffusivity (m**2/s)
1233	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: totporO2_snow !! total O2 porosity (Tans, 1998)
1234	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: totporCH4_snow !! total CH4 porosity (Tans, 1998)
1235	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: diffO2_soil !! oxygen diffusivity (m**2/s)
1236	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: diffCH4_soil !! methane diffusivity (m**2/s)
1237	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: totporO2_soil !! total O2 porosity (Tans, 1998)
1238	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: totporCH4_soil !! total CH4 porosity (Tans, 1998)
1239	INTEGER(i_std),DIMENSION(kjpindex),INTENT(in) :: index !! Indeces of permafrost points on the map
1240
1241	!! 0.2 Output variables
1242
1243	!! 0.3 Modified variables
1244
1245	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: O2_snow !! oxygen (g O2/m**3 air)
1246	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: CH4_snow !! methane (g CH4/m**3 air)
1247	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: O2_soil !! oxygen (g O2/m**3 air)
1248	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: CH4_soil !! methane (g CH4/m**3 air)
1249	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm), intent(inout):: zf_snow !! depths of full levels (m)
1250	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), intent(inout) :: zi_snow !! depths of intermediate levels (m)
1251
1252	!! 0.4 local variables
1253
1254	CHARACTER(LEN=50), SAVE :: last_action = 'not called'
1255	!$OMP THREADPRIVATE(last_action)
1256
1257
1258	! 1. ensure that we do not repeat actions
1259	!
1260	IF ( TRIM(cur_action) .EQ. TRIM(last_action) ) THEN
1261	!
1262	CALL ipslerr_p(3, 'soil_gasdiff_main','CANNOT TAKE THE SAME ACTION TWICE: ',cur_action, last_action)
1263	!
1264	ENDIF
1265	!
1266	! 2. decide what to do
1267	!
1268	IF ( TRIM(cur_action) .EQ. 'initialize' ) THEN
1269	!
1270	! 2.1 initialize
1271	!
1272	IF ( TRIM(last_action) .NE. 'not called' ) THEN
1273	!
1274	CALL ipslerr_p(3, 'soil_gasdiff_main','SOIL MODEL CANNOT BE INITIALIZED TWICE.', '', '')
1275	!
1276	ENDIF
1277	!
1278	CALL soil_gasdiff_alloc( kjpindex )
1279	!
1280	ELSEIF ( TRIM(cur_action) .EQ. 'diffuse' ) THEN
1281	!
1282	! 2.2 calculate soil temperatures
1283	!
1284	CALL soil_gasdiff_diff( kjpindex,time_step,index,psol,tsurf, O2_snow, CH4_snow, O2_soil, CH4_soil)
1285	!
1286	ELSEIF ( TRIM(cur_action) .EQ. 'coefficients' ) THEN
1287	!
1288	! 2.3 calculate coefficients (heat flux and apparent surface heat capacity)
1289	!
1290	CALL soil_gasdiff_coeff( kjpindex,time_step,tprof,O2_snow,CH4_snow, &
1291	diffO2_snow,diffCH4_snow,totporO2_snow,totporCH4_snow,O2_soil,CH4_soil, &
1292	diffO2_soil,diffCH4_soil,totporO2_soil,totporCH4_soil, zi_snow, zf_snow)
1293	!
1294	ELSE
1295	!
1296	! 2.4 do not know this action
1297	!
1298	CALL ipslerr_p(3,'soil_gasdiff_main', 'This action does not exists and must be implemented', &
1299	'or its wrong:',cur_action)
1300	!
1301	ENDIF
1302	!
1303	! 2.5 keep last action in mind
1304	!
1305	last_action = TRIM(cur_action)
1306
1307	IF (printlev>=3) WRITE(,) 'leaving soil_gasdiff_main'
1308	END SUBROUTINE soil_gasdiff_main
1309
1310	!!
1311	!================================================================================================================================
1312	!! SUBROUTINE : soil_gasdiff_alloc
1313	!!
1314	!>\BRIEF This routine allocate arrays related to oxygen and methane in the snow/soil medium
1315	!!
1316	!! DESCRIPTION :
1317	!!
1318	!! RECENT CHANGE(S) : None
1319	!!
1320	!! MAIN OUTPUT VARIABLE(S) :
1321	!!
1322	!! REFERENCE(S) : None
1323	!!
1324	!! FLOWCHART11 : None
1325	!! \n
1326	!_
1327	!================================================================================================================================
1328	SUBROUTINE soil_gasdiff_alloc( kjpindex )
1329
1330	!! 0. Variable and parameter declaration
1331
1332	!! 0.1 Input variables
1333
1334	INTEGER(i_std), INTENT(in) :: kjpindex
1335
1336	!! 0.2 Output variables
1337
1338	!! 0.3 Modified variables
1339
1340	!! 0.4 local variables
1341
1342	INTEGER(i_std) :: ier
1343
1344	! Allocate the variables that need to be saved after soil_gasdiff_coeff
1345
1346	ALLOCATE (alphaO2_soil(kjpindex,ngrnd,nvm),stat=ier)
1347	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for alphaO2_soil','','')
1348
1349	ALLOCATE (betaO2_soil(kjpindex,ngrnd,nvm),stat=ier)
1350	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for betaO2_soil','','')
1351
1352	ALLOCATE (alphaCH4_soil(kjpindex,ngrnd,nvm),stat=ier)
1353	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for alphaCH4_soil','','')
1354
1355	ALLOCATE (betaCH4_soil(kjpindex,ngrnd,nvm),stat=ier)
1356	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for betaCH4_soil','','')
1357
1358	ALLOCATE (alphaO2_snow(kjpindex,nsnow,nvm),stat=ier)
1359	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for alphaO2_snow','','')
1360
1361	ALLOCATE (betaO2_snow(kjpindex,nsnow,nvm),stat=ier)
1362	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for betaO2_snow','','')
1363
1364	ALLOCATE (alphaCH4_snow(kjpindex,nsnow,nvm),stat=ier)
1365	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for alphaCH4_snow','','')
1366
1367	ALLOCATE (betaCH4_snow(kjpindex,nsnow,nvm),stat=ier)
1368	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for betaCH4_snow','','')
1369
1370	ALLOCATE (zf_coeff_snow(kjpindex,0:nsnow,nvm),stat=ier)
1371	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for zf_coeff_snow','','')
1372
1373	ALLOCATE (zi_coeff_snow(kjpindex,nsnow,nvm),stat=ier)
1374	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for zi_coeff_snow','','')
1375
1376	ALLOCATE (mu_snow(kjpindex,nvm),stat=ier)
1377	IF (ier /= 0) CALL ipslerr_p(3,'soil_gasdiff_alloc', 'Pb in alloc for mu_snow','','')
1378
1379
1380	alphaO2_soil(:,:,:) = zero
1381	betaO2_soil(:,:,:) = zero
1382	alphaCH4_soil(:,:,:) = zero
1383	betaCH4_soil(:,:,:) = zero
1384	alphaO2_snow(:,:,:) = zero
1385	betaO2_snow(:,:,:) = zero
1386	alphaCH4_snow(:,:,:) = zero
1387	betaCH4_snow(:,:,:) = zero
1388	zf_coeff_snow(:,:,:) = zero
1389	zi_coeff_snow(:,:,:) = zero
1390	mu_snow(:,:) = zero
1391
1392	END SUBROUTINE soil_gasdiff_alloc
1393
1394	!!
1395	!================================================================================================================================
1396	!! SUBROUTINE : soil_gasdiff_coeff
1397	!!
1398	!>\BRIEF This routine calculate coeff related to gas diffuvisity
1399	!!
1400	!! DESCRIPTION :
1401	!!
1402	!! RECENT CHANGE(S) : None
1403	!!
1404	!! MAIN OUTPUT VARIABLE(S) :
1405	!!
1406	!! REFERENCE(S) : None
1407	!!
1408	!! FLOWCHART11 : None
1409	!! \n
1410	!_
1411	!================================================================================================================================
1412
1413	SUBROUTINE soil_gasdiff_coeff( kjpindex,time_step,tprof,O2_snow,CH4_snow, &
1414	diffO2_snow,diffCH4_snow,totporO2_snow,totporCH4_snow,O2_soil,CH4_soil, &
1415	diffO2_soil,diffCH4_soil,totporO2_soil,totporCH4_soil, zi_snow, zf_snow)
1416
1417
1418	!! 0. Variable and parameter declaration
1419
1420	!! 0.1 Input variables
1421
1422	INTEGER(i_std), INTENT(in) :: kjpindex !! number of grid points
1423	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
1424	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: tprof !! Soil temperature (K)
1425	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: diffO2_snow !! oxygen diffusivity (m**2/s)
1426	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: diffCH4_snow !! methane diffusivity (m**2/s)
1427	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: totporO2_snow !! total O2 porosity (Tans, 1998)
1428	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: totporCH4_snow !! total CH4 porosity (Tans, 1998)
1429	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: diffO2_soil !! oxygen diffusivity (m**2/s)
1430	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: diffCH4_soil !! methane diffusivity (m**2/s)
1431	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: totporO2_soil !! total O2 porosity (Tans, 1998)
1432	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: totporCH4_soil !! total CH4 porosity (Tans, 1998)
1433	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: O2_snow !! oxygen (g O2/m**3 air)
1434	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: CH4_snow !! methane (g CH4/m**3 air)
1435	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: O2_soil !! oxygen (g O2/m**3 air)
1436	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: CH4_soil !! methane (g CH4/m**3 air)
1437	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm), INTENT(in) :: zf_snow !! depths of full levels (m)
1438	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: zi_snow !! depths of intermediate levels (m)
1439
1440	!! 0.2 Output variables
1441
1442	!! 0.3 Modified variables
1443
1444	!! 0.4 local variables
1445
1446	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: xcO2_snow,xdO2_snow
1447	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: xcCH4_snow,xdCH4_snow
1448	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: xcO2_soil,xdO2_soil
1449	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: xcCH4_soil,xdCH4_soil
1450	INTEGER(i_std) :: il
1451	REAL(r_std), DIMENSION(kjpindex,nvm) :: xeO2,xeCH4
1452	LOGICAL, DIMENSION(kjpindex,nvm) :: snow_height_mask_2d
1453	LOGICAL, SAVE :: firstcall = .true.
1454	!$OMP THREADPRIVATE(firstcall)
1455
1456	! loop over materials (soil, snow), beginning at the bottom
1457	!
1458	! 1. define useful variables linked to geometry and physical properties
1459	!
1460	! 1.1 normal levels
1461	!
1462	! default value if inexistent
1463	xcO2_snow(:,:,:) = 0
1464	xdO2_snow(:,:,:) = 0
1465	xcCH4_snow(:,:,:) = 0
1466	xdCH4_snow(:,:,:) = 0
1467	xcO2_soil(:,:,:) = 0
1468	xdO2_soil(:,:,:) = 0
1469	xcCH4_soil(:,:,:) = 0
1470	xdCH4_soil(:,:,:) = 0
1471	xeO2 = 0
1472	xeCH4 = 0
1473	!
1474	snow_height_mask_2d(:,:) = ( heights_snow(:,:) .GT. hmin_tcalc )
1475	!
1476	DO il = 1,nsnow-1
1477	!
1478	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1479	!
1480	xcO2_snow(:,il,:) = ( zf_snow(:,il,:) - zf_snow(:,il-1,:) ) * &
1481	totporO2_snow(:,il,:) / time_step
1482	xcCH4_snow(:,il,:) = ( zf_snow(:,il,:) - zf_snow(:,il-1,:) ) * &
1483	totporCH4_snow(:,il,:) / time_step
1484	!
1485	xdO2_snow(:,il,:) = diffO2_snow(:,il,:) / &
1486	(zi_snow(:,il+1,:)-zi_snow(:,il,:))
1487	xdCH4_snow(:,il,:) = diffCH4_snow(:,il,:) / &
1488	(zi_snow(:,il+1,:)-zi_snow(:,il,:))
1489	!
1490	ENDWHERE
1491	END DO
1492	!
1493	DO il = 1,ngrnd-1
1494	!
1495	WHERE ( veget_mask_2d(:,:) )
1496	!
1497	xcO2_soil(:,il,:) = ( zf_soil(il) - zf_soil(il-1) ) * &
1498	totporO2_soil(:,il,:) / time_step
1499	xcCH4_soil(:,il,:) = ( zf_soil(il) - zf_soil(il-1) ) * &
1500	totporCH4_soil(:,il,:) / time_step
1501	!
1502	xdO2_soil(:,il,:) = diffO2_soil(:,il,:) / &
1503	(zi_soil(il+1)-zi_soil(il))
1504	xdCH4_soil(:,il,:) = diffCH4_soil(:,il,:) / &
1505	(zi_soil(il+1)-zi_soil(il))
1506	!
1507	ENDWHERE
1508	!
1509	ENDDO
1510	!
1511	! 1.2 for the lower boundary, define a similar geometric variable.
1512	!
1513	!snow
1514	!
1515	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1516	xcO2_snow(:,nsnow,:) = ( zf_snow(:,nsnow,:) - &
1517	zf_snow(:,nsnow-1,:) ) * &
1518	totporO2_snow(:,nsnow,:) / time_step
1519	xdO2_snow(:,nsnow,:) = diffO2_snow(:,nsnow,:) / &
1520	( zi_soil(1) + &
1521	zf_snow(:,nsnow,:) - zi_snow(:,nsnow,:) )
1522	xcCH4_snow(:,nsnow,:) = ( zf_snow(:,nsnow,:) - &
1523	zf_snow(:,nsnow-1,:) ) * &
1524	totporCH4_snow(:,nsnow,:) / time_step
1525	xdCH4_snow(:,nsnow,:) = diffCH4_snow(:,nsnow,:) / &
1526	( zi_soil(1) + &
1527	zf_snow(:,nsnow,:) - zi_snow(:,nsnow,:) )
1528	ENDWHERE
1529	!
1530	! soil
1531	!
1532	WHERE ( veget_mask_2d(:,:) ) ! removed heights_soil logic
1533	xcO2_soil(:,ngrnd,:) = &
1534	( zf_soil(ngrnd) - zf_soil(ngrnd-1) ) * &
1535	totporO2_soil(:,ngrnd,:) / time_step
1536	xdO2_soil(:,ngrnd,:) = diffO2_soil(:,ngrnd,:) / &
1537	( zf_soil(ngrnd) - zi_soil(ngrnd) )
1538	xcCH4_soil(:,ngrnd,:) = &
1539	( zf_soil(ngrnd) - zf_soil(ngrnd-1) ) * &
1540	totporCH4_soil(:,ngrnd,:) / time_step
1541	xdCH4_soil(:,ngrnd,:) = diffCH4_soil(:,ngrnd,:) / &
1542	( zf_soil(ngrnd) - zi_soil(ngrnd) )
1543	ENDWHERE
1544	!
1545	! 1.3 extrapolation factor from first levels to surface
1546	!
1547	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1548	mu_snow(:,:) = zi_snow(:,1,:) / ( zi_snow(:,2,:) - zi_snow(:,1,:) )
1549	ELSEWHERE ( veget_mask_2d(:,:) )
1550	mu_snow(:,:) = .5 ! any value
1551	ENDWHERE
1552	!
1553	mu_soil = zi_soil(1) / ( zi_soil(2) - zi_soil(1) )
1554	!
1555	! 2. bottom level: treatment depends on lower boundary condition
1556	!
1557	! soil
1558	!
1559	WHERE ( veget_mask_2d(:,:) ) ! removed heights_soil logic
1560	!
1561	xeO2(:,:) = xcO2_soil(:,ngrnd,:) + xdO2_soil(:,ngrnd-1,:)
1562	xeCH4(:,:) = xcCH4_soil(:,ngrnd,:) + xdCH4_soil(:,ngrnd-1,:)
1563	!
1564	alphaO2_soil(:,ngrnd-1,:) = xdO2_soil(:,ngrnd-1,:) / xeO2(:,:)
1565	alphaCH4_soil(:,ngrnd-1,:) = xdCH4_soil(:,ngrnd-1,:) &
1566	/ xeCH4(:,:)
1567	!
1568	betaO2_soil(:,ngrnd-1,:) = &
1569	(xcO2_soil(:,ngrnd,:)*O2_soil(:,ngrnd,:))/xeO2(:,:)
1570	betaCH4_soil(:,ngrnd-1,:) = &
1571	(xcCH4_soil(:,ngrnd,:)*CH4_soil(:,ngrnd,:))/xeCH4(:,:)
1572	!
1573	ENDWHERE
1574	!
1575	!snow
1576	!
1577	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1578	!
1579	! dernier niveau
1580	!
1581	xeO2(:,:) = xcO2_soil(:,1,:) + &
1582	(1.-alphaO2_soil(:,1,:))*xdO2_soil(:,1,:) + &
1583	xdO2_snow(:,nsnow,:)
1584	xeCH4(:,:) = xcCH4_soil(:,1,:) + &
1585	(1.-alphaCH4_soil(:,1,:))*xdCH4_soil(:,1,:) + &
1586	xdCH4_snow(:,nsnow,:)
1587	!
1588	alphaO2_snow(:,nsnow,:) = xdO2_snow(:,nsnow,:)/xeO2(:,:)
1589	alphaCH4_snow(:,nsnow,:) = xdCH4_snow(:,nsnow,:) &
1590	/xeCH4(:,:)
1591	!
1592	betaO2_snow(:,nsnow,:) = &
1593	( xcO2_soil(:,1,:)*O2_soil(:,1,:) + &
1594	xdO2_soil(:,1,:)*betaO2_soil(:,1,:) ) &
1595	/ xeO2(:,:)
1596	betaCH4_snow(:,nsnow,:) = &
1597	( xcCH4_soil(:,1,:)*CH4_soil(:,1,:) + &
1598	xdCH4_soil(:,1,:)*betaCH4_soil(:,1,:) ) &
1599	/ xeCH4(:,:)
1600	!
1601	! avant-dernier niveau
1602	!
1603	xeO2(:,:) = xcO2_snow(:,nsnow,:) + &
1604	(1.-alphaO2_snow(:,nsnow,:))*xdO2_snow(:,nsnow,:) + &
1605	xdO2_snow(:,nsnow-1,:)
1606	xeCH4(:,:) = xcCH4_snow(:,nsnow,:) + &
1607	(1.-alphaCH4_snow(:,nsnow,:))*xdCH4_snow(:,nsnow,:) &
1608	+ xdCH4_snow(:,nsnow-1,:)
1609	!
1610	alphaO2_snow(:,nsnow-1,:) = &
1611	xdO2_snow(:,nsnow-1,:) / xeO2(:,:)
1612	alphaCH4_snow(:,nsnow-1,:) = &
1613	xdCH4_snow(:,nsnow-1,:) / xeCH4(:,:)
1614	!
1615	betaO2_snow(:,nsnow-1,:) = &
1616	( xcO2_snow(:,nsnow,:)*O2_snow(:,nsnow,:) + &
1617	xdO2_snow(:,nsnow,:)*betaO2_snow(:,nsnow,:) ) &
1618	/ xeO2(:,:)
1619	betaCH4_snow(:,nsnow-1,:) = &
1620	( xcCH4_snow(:,nsnow,:)*CH4_snow(:,nsnow,:) + &
1621	xdCH4_snow(:,nsnow,:)*betaCH4_snow(:,nsnow,:) ) &
1622	/ xeCH4(:,:)
1623	!
1624	ELSEWHERE ( veget_mask_2d(:,:) )
1625	!
1626	alphaO2_snow(:,nsnow,:) = 1.
1627	alphaCH4_snow(:,nsnow,:) = 1.
1628	betaO2_snow(:,nsnow,:) = zero
1629	betaCH4_snow(:,nsnow,:) = zero
1630	!
1631	alphaO2_snow(:,nsnow-1,:) = 1.
1632	alphaCH4_snow(:,nsnow-1,:) = 1.
1633	betaO2_snow(:,nsnow-1,:) = zero
1634	betaCH4_snow(:,nsnow-1,:) = zero
1635	!
1636	ENDWHERE
1637	!
1638
1639	!
1640	! 3. the other levels
1641	!
1642	DO il = nsnow-2,1,-1 !snow
1643	!
1644	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1645	!
1646	xeO2(:,:) = xcO2_snow(:,il+1,:) + &
1647	(1.-alphaO2_snow(:,il+1,:))*xdO2_snow(:,il+1,:) + xdO2_snow(:,il,:)
1648	xeCH4(:,:) = xcCH4_snow(:,il+1,:) + &
1649	(1.-alphaCH4_snow(:,il+1,:))*xdCH4_snow(:,il+1,:) + &
1650	xdCH4_snow(:,il,:)
1651	!
1652	alphaO2_snow(:,il,:) = xdO2_snow(:,il,:) / xeO2(:,:)
1653	alphaCH4_snow(:,il,:) = xdCH4_snow(:,il,:) / xeCH4(:,:)
1654	!
1655	betaO2_snow(:,il,:) = &
1656	( xcO2_snow(:,il+1,:)*O2_snow(:,il+1,:) + &
1657	xdO2_snow(:,il+1,:)*betaO2_snow(:,il+1,:) ) / xeO2(:,:)
1658	betaCH4_snow(:,il,:) = &
1659	( xcCH4_snow(:,il+1,:)*CH4_snow(:,il+1,:) + &
1660	xdCH4_snow(:,il+1,:)*betaCH4_snow(:,il+1,:) ) / xeCH4(:,:)
1661	!
1662	ELSEWHERE ( veget_mask_2d(:,:) )
1663	!
1664	alphaO2_snow(:,il,:) = 1.
1665	alphaCH4_snow(:,il,:) = 1.
1666	!
1667	betaO2_snow(:,il,:) = zero
1668	betaCH4_snow(:,il,:) = zero
1669	!
1670	ENDWHERE
1671	!
1672	ENDDO
1673	!
1674	DO il = ngrnd-2,1,-1 !soil
1675	!
1676	WHERE ( veget_mask_2d(:,:) ) !removed heights_soil logic
1677	!
1678	xeO2(:,:) = xcO2_soil(:,il+1,:) + &
1679	(1.-alphaO2_soil(:,il+1,:))*xdO2_soil(:,il+1,:) + xdO2_soil(:,il,:)
1680	xeCH4(:,:) = xcCH4_soil(:,il+1,:) + &
1681	(1.-alphaCH4_soil(:,il+1,:))*xdCH4_soil(:,il+1,:) + &
1682	xdCH4_soil(:,il,:)
1683	!
1684	alphaO2_soil(:,il,:) = xdO2_soil(:,il,:) / xeO2(:,:)
1685	alphaCH4_soil(:,il,:) = xdCH4_soil(:,il,:) / xeCH4(:,:)
1686	!
1687	betaO2_soil(:,il,:) = &
1688	( xcO2_soil(:,il+1,:)*O2_soil(:,il+1,:) + &
1689	xdO2_soil(:,il+1,:)*betaO2_soil(:,il+1,:) ) / xeO2(:,:)
1690	betaCH4_soil(:,il,:) = &
1691	( xcCH4_soil(:,il+1,:)*CH4_soil(:,il+1,:) + &
1692	xdCH4_soil(:,il+1,:)*betaCH4_soil(:,il+1,:) ) / xeCH4(:,:)
1693	!
1694	ENDWHERE
1695	!
1696	ENDDO
1697	!
1698	! 4. store thickness of the different levels for all soil types (for security)
1699	!
1700	zf_coeff_snow(:,:,:) = zf_snow(:,:,:)
1701	zi_coeff_snow(:,:,:) = zi_snow(:,:,:)
1702
1703	!--hist out for keeping track of these
1704	IF (firstcall) THEN
1705	firstcall = .false.
1706	ELSE
1707	ENDIF
1708
1709	END SUBROUTINE soil_gasdiff_coeff
1710
1711	!!
1712	!================================================================================================================================
1713	!! SUBROUTINE : soil_gasdiff_diff
1714	!!
1715	!>\BRIEF This routine update oxygen and methane in the snow and soil
1716	!!
1717	!! DESCRIPTION :
1718	!!
1719	!! RECENT CHANGE(S) : None
1720	!!
1721	!! MAIN OUTPUT VARIABLE(S) :
1722	!!
1723	!! REFERENCE(S) : None
1724	!!
1725	!! FLOWCHART11 : None
1726	!! \n
1727	!_
1728	!================================================================================================================================
1729
1730	SUBROUTINE soil_gasdiff_diff( kjpindex,time_step,index,pb,tsurf, O2_snow, CH4_snow, O2_soil, CH4_soil)
1731
1732	!! 0. Variable and parameter declaration
1733
1734	!! 0.1 Input variables
1735
1736	INTEGER(i_std), INTENT(in) :: kjpindex !! number of grid points
1737	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
1738	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure
1739	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: tsurf !! Surface temperature
1740	INTEGER(i_std),DIMENSION(kjpindex),INTENT(in) :: index !! Indeces of the points on the map
1741	!! 0.2 Output variables
1742
1743	!! 0.3 Modified variables
1744
1745	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: O2_snow !! oxygen (g O2/m**3 air)
1746	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: CH4_snow !! methane (g CH4/m**3 air)
1747	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: O2_soil !! oxygen (g O2/m**3 air)
1748	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: CH4_soil !! methane (g CH4/m**3 air)
1749
1750	!! 0.4 local variables
1751
1752	INTEGER(i_std) :: it, ip, il, iv
1753	LOGICAL, DIMENSION(kjpindex,nvm) :: snowtop
1754	REAL(r_std), DIMENSION(kjpindex,nvm) :: O2sa, CH4sa
1755
1756	!
1757	! 1.1 Determine which is the first existing soil type.
1758	!
1759	snowtop(:,:) = .FALSE.
1760	!
1761	!ignore snow for now...
1762	WHERE ( heights_snow(:,:) .GT. hmin_tcalc )
1763	snowtop(:,:) = .TRUE.
1764	ENDWHERE
1765	!
1766	! 2.gas diffusion
1767	!
1768	! 2.1 top level
1769	!
1770	! 2.1.1 non-existing
1771	!
1772	DO iv = 1, nvm
1773	O2sa(:,iv) = pb(:)/(RRtsurf(:)) O2_surf * wO2
1774	CH4sa(:,iv) = pb(:)/(RRtsurf(:)) CH4_surf * wCH4
1775	ENDDO
1776	!
1777	WHERE ( (.NOT. snowtop(:,:)) .AND. veget_mask_2d(:,:) ) ! it equals 1 (snow) but there is no snow...
1778	!
1779	O2_snow(:,1,:) = O2sa(:,:)
1780	CH4_snow(:,1,:) = CH4sa(:,:)
1781	!
1782	O2_soil(:,1,:) = ( O2sa(:,:) + mu_soil*betaO2_soil(:,1,:) ) / &
1783	( 1. + mu_soil*(1.-alphaO2_soil(:,1,:)) )
1784	CH4_soil(:,1,:) = ( CH4sa(:,:) + mu_soil*betaCH4_soil(:,1,:) ) / &
1785	( 1. + mu_soil*(1.-alphaCH4_soil(:,1,:)) )
1786	!
1787	ENDWHERE
1788	!
1789	! 2.1.2 first existing soil type
1790	!
1791	WHERE ( snowtop(:,:) .AND. veget_mask_2d(:,:) )
1792	!
1793	O2_snow(:,1,:) = ( O2sa(:,:) + mu_snow(:,:)*betaO2_snow(:,1,:) ) / &
1794	( 1. + mu_snow(:,:)*(1.-alphaO2_snow(:,1,:)) )
1795	CH4_snow(:,1,:) = ( CH4sa(:,:) + mu_snow(:,:)*betaCH4_snow(:,1,:) ) / &
1796	( 1. + mu_snow(:,:)*(1.-alphaCH4_snow(:,1,:)) )
1797	!
1798	O2_soil(:,1,:) = &
1799	alphaO2_snow(:,nsnow,:) * O2_snow(:,nsnow,:) + &
1800	betaO2_snow(:,nsnow,:)
1801	CH4_soil(:,1,:) = &
1802	alphaCH4_snow(:,nsnow,:) * CH4_snow(:,nsnow,:) + &
1803	betaCH4_snow(:,nsnow,:)
1804	! debug: need to check for weird numbers here!
1805	ENDWHERE
1806	!
1807	! 2.2 other levels
1808	!
1809	DO il = 2, nsnow
1810
1811	WHERE ( veget_mask_2d(:,:) )
1812	!
1813	O2_snow(:,il,:) = &
1814	alphaO2_snow(:,il-1,:) * O2_snow(:,il-1,:) + &
1815	betaO2_snow(:,il-1,:)
1816	CH4_snow(:,il,:) = &
1817	alphaCH4_snow(:,il-1,:) * CH4_snow(:,il-1,:) + &
1818	betaCH4_snow(:,il-1,:)
1819	END WHERE
1820	ENDDO
1821	DO il = 2, ngrnd
1822
1823	WHERE ( veget_mask_2d(:,:) )
1824	!
1825	O2_soil(:,il,:) = &
1826	alphaO2_soil(:,il-1,:) * O2_soil(:,il-1,:) + &
1827	betaO2_soil(:,il-1,:)
1828	CH4_soil(:,il,:) = &
1829	alphaCH4_soil(:,il-1,:) * CH4_soil(:,il-1,:) + &
1830	betaCH4_soil(:,il-1,:)
1831	END WHERE
1832	ENDDO
1833
1834	END SUBROUTINE soil_gasdiff_diff
1835
1836	!!
1837	!================================================================================================================================
1838	!! SUBROUTINE : get_gasdiff
1839	!!
1840	!>\BRIEF This routine update oxygen and methane in the snow and soil
1841	!!
1842	!! DESCRIPTION :
1843	!!
1844	!! RECENT CHANGE(S) : None
1845	!!
1846	!! MAIN OUTPUT VARIABLE(S) :
1847	!!
1848	!! REFERENCE(S) : None
1849	!!
1850	!! FLOWCHART11 : None
1851	!! \n
1852	!_
1853	!================================================================================================================================
1854	SUBROUTINE get_gasdiff (kjpindex,hslong,tprof,snow,airvol_snow, &
1855	totporO2_snow,totporCH4_snow,diffO2_snow,diffCH4_snow, &
1856	airvol_soil,totporO2_soil,totporCH4_soil,diffO2_soil,diffCH4_soil, snowrho)
1857
1858	!! 0. Variable and parameter declaration
1859
1860	!! 0.1 Input variables
1861
1862	INTEGER(i_std), INTENT(in) :: kjpindex !! number of grid points
1863	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: hslong !! deep long term soil humidity profile
1864	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: tprof !! Soil temperature (K)
1865	REAL(r_std), DIMENSION(kjpindex,nsnow), INTENT(in) :: snowrho !! snow density (Kg/m^3)
1866	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: snow !! Snow mass [Kg/m^2]
1867	!! 0.2 Output variables
1868
1869	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(out) :: airvol_soil
1870	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(out) :: totporO2_soil !! total O2 porosity (Tans, 1998)
1871	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(out) :: totporCH4_soil !! total CH4 porosity
1872	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(out) :: diffO2_soil !! oxygen diffusivity (m**2/s)
1873	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(out) :: diffCH4_soil !! methane diffusivity (m**2/s)
1874	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: airvol_snow
1875	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: totporO2_snow !! total O2 porosity (Tans, 1998)
1876	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: totporCH4_snow !! total CH4 porosity (Tans, 1998)
1877	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: diffO2_snow !! oxygen diffusivity (m**2/s)
1878	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: diffCH4_snow !! methane diffusivity (m**2/s)
1879
1880	!! 0.3 Modified variables
1881
1882	!! 0.4 local variables
1883
1884	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: density_snow
1885	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: porosity_snow
1886	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: tortuosity_snow
1887	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: density_soil
1888	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: porosity_soil
1889	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: tortuosity_soil
1890	INTEGER(i_std) :: it,ip, il, iv
1891	REAL(r_std) :: x, rho_iw
1892	REAL(r_std) :: csat, fng
1893	REAL(r_std), SAVE :: cond_fact
1894	!$OMP THREADPRIVATE(cond_fact)
1895	LOGICAL, SAVE :: pr_fois=.TRUE.
1896	!$OMP THREADPRIVATE(pr_fois)
1897
1898	IF (pr_fois) THEN
1899	cond_fact=1.
1900	CALL getin_p('COND_FACT',cond_fact)
1901	WRITE(,) 'COND_FACT=',cond_fact
1902	pr_fois=.FALSE.
1903	ENDIF
1904
1905	!
1906	! 1. Three-layers snow model with snow density resolved at each snow layer
1907	!
1908	DO iv = 1, nvm
1909	density_snow(:,:,iv) = snowrho(:,:)
1910	ENDDO
1911	porosity_snow(:,:,:) = (1. - density_snow(:,:,:)/rho_ice )
1912	tortuosity_snow(:,:,:) = porosity_snow(:,:,:)**(1./3.) ! based on Sommerfeld et al., GBC, 1996
1913	diffO2_snow(:,:,:) = diffO2_air * porosity_snow(:,:,:) * tortuosity_snow(:,:,:)
1914	diffCH4_snow(:,:,:) = diffCH4_air * porosity_snow(:,:,:) * tortuosity_snow(:,:,:)
1915	airvol_snow(:,:,:) = MAX(porosity_snow(:,:,:),avm)
1916	totporO2_snow(:,:,:) = airvol_snow(:,:,:)
1917	totporCH4_snow(:,:,:) = airvol_snow(:,:,:)
1918	!
1919	! 2. soil: depends on temperature and soil humidity
1920	!
1921	DO ip = 1, kjpindex
1922	!
1923	DO iv = 1, nvm
1924	!
1925	IF ( veget_mask_2d(ip,iv) ) THEN
1926	!
1927	DO il = 1, ngrnd
1928	!
1929	! 2.1 soil dry density, porosity, and dry heat capacity
1930	!
1931	porosity_soil(ip,il,iv) = tetasat
1932	!
1933	!
1934	! 2.2 heat capacity and density as a function of
1935	! ice and water content
1936	! removed these as we are calculating thermal evolution in the sechiba subroutines
1937
1938	!
1939	! 2.3 oxygen diffusivity: soil can get waterlogged,
1940	! therefore take soil humidity into account
1941	!
1942	tortuosity_soil(ip,il,iv) = 2./3. ! Hillel, 1980
1943	airvol_soil(ip,il,iv) = porosity_soil(ip,il,iv)*(1.-hslong(ip,il,iv))
1944	totporO2_soil(ip,il,iv) = airvol_soil(ip,il,iv) + porosity_soil(ip,il,iv)BunsenO2hslong(ip,il,iv)
1945	totporCH4_soil(ip,il,iv) = airvol_soil(ip,il,iv) + porosity_soil(ip,il,iv)BunsenCH4hslong(ip,il,iv)
1946	diffO2_soil(ip,il,iv) = (diffO2_air*airvol_soil(ip,il,iv) + &
1947	diffO2_wBunsenO2hslong(ip,il,iv)porosity_soil(ip,il,iv))tortuosity_soil(ip,il,iv)
1948	diffCH4_soil(ip,il,iv) = (diffCH4_air*airvol_soil(ip,il,iv) + &
1949	diffCH4_wBunsenCH4hslong(ip,il,iv)porosity_soil(ip,il,iv))tortuosity_soil(ip,il,iv)
1950	!
1951	END DO
1952	ELSE
1953	tortuosity_soil(ip,:,iv) = EPSILON(0.)
1954	airvol_soil(ip,:,iv) = EPSILON(0.)
1955	totporO2_soil(ip,:,iv) = EPSILON(0.)
1956	totporCH4_soil(ip,:,iv) = EPSILON(0.)
1957	diffO2_soil(ip,:,iv) = EPSILON(0.)
1958	diffCH4_soil(ip,:,iv) = EPSILON(0.)
1959	END IF
1960	ENDDO
1961	ENDDO
1962
1963	END SUBROUTINE get_gasdiff
1964
1965	!!
1966	!================================================================================================================================
1967	!! SUBROUTINE : traMplan
1968	!!
1969	!>\BRIEF This routine calculates plant-mediated transport of methane
1970	!!
1971	!! DESCRIPTION :
1972	!!
1973	!! RECENT CHANGE(S) : None
1974	!!
1975	!! MAIN OUTPUT VARIABLE(S) :
1976	!!
1977	!! REFERENCE(S) : None
1978	!!
1979	!! FLOWCHART11 : None
1980	!! \n
1981	!_
1982	!================================================================================================================================
1983	SUBROUTINE traMplan(CH4,O2,kjpindex,time_step,totporCH4,totporO2,z_root,rootlev,Tgr,Tref,hslong,flupmt, &
1984	refdep, zi_soil, tprof)
1985
1986	!! 0. Variable and parameter declaration
1987
1988	!! 0.1 Input variables
1989
1990	INTEGER(i_std), INTENT(in) :: kjpindex
1991	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
1992	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: totporO2 !! total oxygen porosity
1993	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: totporCH4 !! total methane porosity
1994	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: tprof !! soil temperature (K)
1995	INTEGER(i_std),DIMENSION(kjpindex,nvm),INTENT(in) :: rootlev !! the deepest model level within the rooting depth
1996	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: z_root !! the rooting depth
1997	REAL(r_std), INTENT(in) :: Tgr !! Temperature at which plants begin to grow (C)
1998	REAL(r_std), DIMENSION(ngrnd), INTENT(in) :: zi_soil !! depths at intermediate levels
1999	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: hslong !! deep soil humidity
2000
2001	!! 0.2 Output variables
2002
2003	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: flupmt !! plant-mediated methane flux (g m-2 s-1)
2004
2005	!! 0.3 Modified variables
2006
2007	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(inout) :: Tref !! Ref. temperature for growing season caluculation (C)
2008	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: O2
2009	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: CH4
2010
2011	!! 0.4 local variables
2012	REAL(r_std), DIMENSION(kjpindex,nvm) :: CH4atm !! CH4 atm concentration
2013	REAL(r_std), DIMENSION(kjpindex,nvm) :: dCH4 !! delta CH4 per m3 air
2014	REAL(r_std), DIMENSION(kjpindex,nvm) :: dO2 !! O2 change
2015	REAL(r_std), DIMENSION(kjpindex,nvm) :: fgrow !! Plant growing state (maturity index)
2016	REAL(r_std) :: froot !! vertical distribution of roots
2017	REAL(r_std) :: Tmat !! Temperature at which plants reach maturity (C)
2018	REAL(r_std), PARAMETER :: La_min = zero
2019	REAL(r_std), PARAMETER :: La = 4.
2020	REAL(r_std), PARAMETER :: La_max = La_min + La
2021	REAL(r_std), PARAMETER :: Tveg = 10 !! Vegetation type control on the plant-mediated transport, Adjustable parameter,
2022	!! but we start from 10 following Walter et al (2001) tundra value
2023	REAL(r_std), PARAMETER :: Pox = 0.5 !! fraction of methane oxydized near the roots
2024	LOGICAL, SAVE :: firstcall=.TRUE.
2025	!$OMP THREADPRIVATE(firstcall)
2026	INTEGER(i_std) :: il,ip, iv
2027	LOGICAL, SAVE :: check = .FALSE.
2028	!$OMP THREADPRIVATE(check)
2029	REAL(r_std), INTENT(in) :: refdep !! Depth to compute reference temperature for the growing season (m)
2030	INTEGER(i_std), SAVE :: reflev = 0 !! Level closest to reference depth refdep
2031	!$OMP THREADPRIVATE(reflev)
2032
2033
2034	IF (firstcall) THEN
2035	firstcall = .FALSE.
2036
2037	! Find the level closest to refdep
2038	DO il=1,ngrnd
2039	IF (zi_soil(il) .GT. refdep .AND. reflev.EQ.0) reflev = il-1
2040	ENDDO
2041	IF (reflev.EQ.0) reflev = ngrnd
2042
2043
2044	IF (check) THEN
2045	OPEN (28,file='pmt.dat',status='unknown')
2046	OPEN (29,file='pmtf.dat',status='unknown')
2047	ENDIF
2048	ENDIF
2049
2050	! Update seasonal reference temperature trace record
2051	WHERE ( veget_mask_2d(:,:) )
2052	Tref(:,:) = tprof(:,reflev,:) - ZeroCelsius
2053	END WHERE
2054
2055	Tmat = Tgr + 10._r_std
2056	flupmt(:,:) = zero
2057	CH4atm(:,:) = zero
2058
2059
2060	! Plant growing state (maturity index)
2061	WHERE (Tref(:,:).LE.Tgr .AND. veget_mask_2d(:,:) )
2062	fgrow(:,:) = La_min
2063	ELSEWHERE (Tref(:,:).GE.Tmat .AND. veget_mask_2d(:,:) )
2064	fgrow(:,:) = La_max
2065	ELSEWHERE ( veget_mask_2d(:,:))
2066	fgrow(:,:) = La_min + La * (1 - ((Tmat - Tref(:,:))/(Tmat - Tgr))**2)
2067	ENDWHERE
2068
2069	DO ip=1,kjpindex
2070	DO iv = 1, nvm
2071	IF ( (z_root(ip,iv) .GT. 0.) .AND. veget_mask_2d(ip,iv) ) THEN ! added this to prevent pmt calcs when soil frozen
2072	DO il=1,rootlev(ip,iv)
2073	! vertical distribution of roots
2074	froot = MAX( 2 * (z_root(ip,iv) - REAL( zi_soil(il) )) / z_root(ip,iv), zero)
2075	! Methane removal from a given depth. We assume that the methane
2076	! in air pores is always in equilibrium with that dissolved
2077	! in water-filled pores. If soil humidity is low,
2078	! with root water as well
2079	! We assume that PMT is proportional to soil humidity
2080	dCH4(ip,iv) = 0.01_r_std * Tveg * froot * fgrow(ip,iv) * hslong(ip,il,iv) * (CH4(ip,il,iv) - CH4atm(ip,iv))
2081	! No transport if soil concentration is less than atmospheric
2082	IF (dCH4(ip,iv).LT.CH4atm(ip,iv)) dCH4(ip,iv) = zero
2083	! Strange thing in WH 2001: 0.01Tvegfroot*fgrow > 1
2084	! at Tveg=15, froot&fgrow=max, i.e. more CH4 is taken than available
2085	! So need to impose a limitation:
2086	IF (dCH4(ip,iv).GT.CH4(ip,il,iv)) dCH4(ip,iv) = CH4(ip,il,iv)
2087	! Methane concentration is decreased within the root layer:
2088
2089	CH4(ip,il,iv) = CH4(ip,il,iv) - dCH4(ip,iv)
2090	! O2 concentration is decreased in reaction with
2091	! dCH4Poxtime_step
2092	dO2(ip,iv) = dCH4(ip,iv)Pox wO2/wCH4 * totporCH4(ip,il,iv)/totporO2(ip,il,iv)
2093	IF ( dO2(ip,iv).LT.O2(ip,il,iv) ) O2(ip,il,iv) = O2(ip,il,iv) - dO2(ip,iv)
2094
2095	! CO2 concentration is increased by dCH4(:)*Pox
2096
2097	! Integration
2098	flupmt(ip,iv) = flupmt(ip,iv) + dCH4(ip,iv)totporCH4(ip,il,iv)/time_step (1 - Pox) * &
2099	( zf_soil(il) - zf_soil(il-1) )
2100	ENDDO
2101	END IF
2102	ENDDO
2103	ENDDO
2104
2105	IF (check) THEN
2106	WRITE(29,*) flupmt(:,:)
2107	CALL flush(28)
2108	CALL flush(29)
2109	END IF
2110
2111	END SUBROUTINE traMplan
2112
2113	!!
2114	!================================================================================================================================
2115	!! SUBROUTINE : ebullition
2116	!!
2117	!>\BRIEF This routine calculates CH4 ebullition
2118	!!
2119	!! DESCRIPTION :
2120	!!
2121	!! RECENT CHANGE(S) : None
2122	!!
2123	!! MAIN OUTPUT VARIABLE(S) :
2124	!!
2125	!! REFERENCE(S) : None
2126	!!
2127	!! FLOWCHART11 : None
2128	!! \n
2129	!_
2130	!================================================================================================================================
2131	SUBROUTINE ebullition (kjpindex,time_step,tprof,totporCH4_soil,hslong,Ch4_soil,febul)
2132
2133	!! 0. Variable and parameter declaration
2134
2135	!! 0.1 Input variables
2136
2137	INTEGER(i_std), INTENT(in) :: kjpindex
2138	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
2139	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm),INTENT(in) :: tprof !! soil temperature (K)
2140	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: totporCH4_soil !! total methane porosity
2141	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: hslong !! deep soil humidity
2142
2143	!! 0.2 Output variables
2144
2145	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: febul !! CH4 ebullition
2146
2147	!! 0.3 Modified variables
2148
2149	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: Ch4_soil !! methane
2150
2151	!! 0.4 Local variables
2152	REAL(r_std) :: dCH4, CH4d
2153	INTEGER(i_std) :: ip, il, iv
2154	REAL(r_std) :: dz
2155	REAL(r_std), PARAMETER :: tortuosity=2./3.
2156	REAL(r_std), PARAMETER :: wsize=0.01
2157	REAL(r_std), PARAMETER :: CH4wm = 12.E-3 !! CH4 concentration threshold for ebullition (8-16 mg/m3 in Walter&Heimann 2000)
2158	REAL(r_std) :: hum
2159
2160	febul(:,:)=zero
2161
2162	DO ip=1,kjpindex
2163	DO iv = 1, nvm
2164	IF ( veget_mask_2d(ip,iv) ) THEN
2165	IF (hslong(ip,1,iv).GT.ebuthr) THEN
2166	DO il = ngrnd, 1, -1
2167	CH4d = Ch4_soil(ip,il,iv) - CH4wm/BunsenCH4
2168	IF (CH4d .GT. EPSILON(0.)) THEN
2169	IF (il.GT.1) THEN
2170	dz = zi_soil(il) - zi_soil(il-1)
2171	hum = ( hslong(ip,il,iv) + hslong(ip,il-1,iv) ) / 2
2172	ELSE
2173	dz = zi_soil(1)
2174	hum = hslong(ip,1,iv)
2175	ENDIF
2176
2177	dCH4 = hum*( dz/wsize/tortuosity ) CH4d
2178	dCH4 = CH4d
2179
2180	Ch4_soil(ip,il,iv) = Ch4_soil(ip,il,iv) - dCH4
2181
2182
2183	febul(ip,iv) = febul(ip,iv) + dCH4 * totporCH4_soil(ip,il,iv) * &
2184	( zf_soil(il) - zf_soil(il-1) ) / time_step
2185
2186	ENDIF
2187	ENDDO
2188	ENDIF
2189	END IF
2190	ENDDO
2191	ENDDO
2192	END SUBROUTINE ebullition
2193
2194	!!
2195	!================================================================================================================================
2196	!! FUNCTION : stomate_soil_carbon_discretization_microactem
2197	!!
2198	!>\BRIEF This function calculates parameters describing bacterial activity (time constant tau[s]) as a function of temperature
2199	!!
2200	!! DESCRIPTION :
2201	!!
2202	!! RECENT CHANGE(S) : None
2203	!!
2204	!! MAIN OUTPUT VARIABLE(S) :
2205	!!
2206	!! REFERENCE(S) : None
2207	!!
2208	!! FLOWCHART11 : None
2209	!! \n
2210	!_
2211	!================================================================================================================================
2212	FUNCTION stomate_soil_carbon_discretization_microactem &
2213	( temp, frozen_respiration_func, moist_in, i_ind, j_ind, k_ind, zi_soil ) RESULT ( fbact )
2214
2215	!! 0. Variable and parameter declaration
2216
2217	!! 0.1 Input variables
2218
2219	INTEGER(i_std), INTENT(in) :: i_ind
2220	INTEGER(i_std), INTENT(in) :: j_ind
2221	INTEGER(i_std), INTENT(in) :: k_ind
2222	INTEGER(i_std), INTENT(in) :: frozen_respiration_func
2223	REAL, DIMENSION(i_ind, j_ind, k_ind), INTENT(in) :: moist_in
2224	REAL, DIMENSION(i_ind, j_ind, k_ind), INTENT(in) :: temp
2225	REAL, DIMENSION(j_ind), INTENT(in) :: zi_soil
2226
2227	!! 0.2 Output variables
2228
2229	!! 0.3 Modified variables
2230
2231	!! 0.4 Local variables
2232
2233	REAL, DIMENSION(i_ind, j_ind, k_ind) :: fbact
2234	REAL, DIMENSION(i_ind, j_ind, k_ind) :: tempfunc_result
2235	REAL, DIMENSION(i_ind, j_ind, k_ind) :: temp_kelvin
2236	INTEGER(i_std), PARAMETER :: ntconfun = 7
2237	REAL(r_std), DIMENSION(ntconfun) :: tconfun
2238	REAL(r_std), DIMENSION(ntconfun) :: tauconfun
2239	INTEGER :: itz
2240	INTEGER :: ii, ij, ik
2241	REAL, DIMENSION(i_ind, j_ind, k_ind) :: moistfunc_result
2242	REAL(r_std), parameter :: q10 = 2.0
2243	logical, parameter :: limit_decomp_moisture = .true.
2244
2245	temp_kelvin(:,:,:) = temp(:,:,:) + ZeroCelsius
2246	SELECT CASE(frozen_respiration_func)
2247
2248	CASE(0) ! this is the standard ORCHIDEE state
2249
2250	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2251	tempfunc_result(:,:,:) = MIN( 1._r_std, tempfunc_result(:,:,:) )
2252
2253	CASE(1) ! cutoff respiration when T < -1C
2254	WHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius ) ! normal as above
2255	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2256	ELSEWHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius - 1. ) ! linear dropoff to zero
2257	tempfunc_result(:,:,:) = (temp_kelvin(:,:,:) - (ZeroCelsius - 1.)) * &
2258	EXP( log(q10) * ( ZeroCelsius - (ZeroCelsius+30.) ) / 10. )
2259	ELSEWHERE ! zero
2260	tempfunc_result(:,:,:) = EPSILON(0.)
2261	endwhere
2262
2263	tempfunc_result(:,:,:) = MAX(MIN( 1._r_std, tempfunc_result(:,:,:) ), EPSILON(0.))
2264
2265	CASE(2) ! cutoff respiration when T < -3C
2266	WHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius ) ! normal as above
2267	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2268	ELSEWHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius - 3. ) ! linear dropoff to zero
2269	tempfunc_result(:,:,:) = ((temp_kelvin(:,:,:) - (ZeroCelsius - 3.))/3.) * &
2270	EXP( log(q10) * ( ZeroCelsius - (ZeroCelsius+30.) ) / 10. )
2271	ELSEWHERE ! zero
2272	tempfunc_result(:,:,:) = EPSILON(0.)
2273	endwhere
2274
2275	CASE(3) ! q10 = 100 when below zero
2276	WHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius ) ! normal as above
2277	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2278	ELSEWHERE
2279	tempfunc_result(:,:,:) = EXP( log(100.) * ( temp_kelvin(:,:,:) - (ZeroCelsius) ) / 10. ) * &
2280	EXP( log(q10) * ( -30. ) / 10. )
2281	endwhere
2282
2283	CASE(4) ! q10 = 1000 when below zero
2284	WHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius ) ! normal as above
2285	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2286	ELSEWHERE
2287	tempfunc_result(:,:,:) = EXP( log(1000.) * ( temp_kelvin(:,:,:) - (ZeroCelsius) ) / 10. ) * &
2288	EXP( log(q10) * ( -30. ) / 10. )
2289	endwhere
2290
2291	CASE DEFAULT
2292	CALL ipslerr_p(3,'stomate_soil_carbon_discretization_microactem', &
2293	'frozen_respiration_func can only be 0, 1, 2, 3 or 4','','')
2294	END SELECT
2295	tempfunc_result(:,:,:) = MAX(MIN( 1._r_std, tempfunc_result(:,:,:) ), EPSILON(0.))
2296
2297	!---- stomate residence times: -----!
2298	! residence times in carbon pools (days)
2299	!carbon_tau(iactive) = .149 * one_year !!!!???? 1.5 years
2300	!carbon_tau(islow) = 5.48 * one_year !!!!???? 25 years
2301	!carbon_tau(ipassive) = 241. * one_year !!!!???? 1000 years
2302	!-----------------------------------!
2303	IF ( limit_decomp_moisture ) THEN
2304	! stomate moisture control function
2305	moistfunc_result(:,:,:) = -1.1 * moist_in(:,:,:) * moist_in(:,:,:) + 2.4 * moist_in(:,:,:) - 0.29
2306	moistfunc_result(:,:,:) = max( 0.25_r_std, min( 1._r_std, moistfunc_result(:,:,:) ) )
2307	ELSE
2308	moistfunc_result(:,:,:) = 1._r_std
2309	ENDIF
2310
2311	DO ij = 1, ngrnd
2312	fbact(:,ij,:) = stomate_tau/(moistfunc_result(:,ij,:) * tempfunc_result(:,ij,:)) / EXP(-zi_soil(ij)/depth_modifier)
2313	ENDDO
2314
2315	! chaoyue: We tentatively increase the turnover of soil C in croplands,
2316	! as shown here to decrease its tau -- the residence time.
2317	DO ik = 1,nvm
2318	fbact(:,:,ik) = fbact(:,:,ik)/ decomp_factor(ik)
2319	ENDDO
2320
2321	END FUNCTION stomate_soil_carbon_discretization_microactem
2322
2323
2324	!!
2325	!================================================================================================================================
2326	!! SUBROUTINE : snowlevels
2327	!!
2328	!>\BRIEF This routine calculates depths of full levels and intermediate
2329	!! levels related to snow pack
2330	!!
2331	!! DESCRIPTION :
2332	!!
2333	!! RECENT CHANGE(S) : None
2334	!!
2335	!! MAIN OUTPUT VARIABLE(S) :
2336	!!
2337	!! REFERENCE(S) : None
2338	!!
2339	!! FLOWCHART11 : None
2340	!! \n
2341	!_
2342	!================================================================================================================================
2343
2344	SUBROUTINE snowlevels( kjpindex, snowdz, zi_snow, zf_snow, veget_max )
2345
2346	!! 0. Variable and parameter declaration
2347
2348	!! 0.1 Input variables
2349
2350	INTEGER(i_std), INTENT(in) :: kjpindex
2351	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: veget_max !! maximum vegetation fraction
2352	REAL(r_std), DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz !! snow depth [m]
2353
2354	!! 0.2 Output variables
2355
2356	!! 0.3 Modified variables
2357
2358	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm), INTENT(inout) :: zf_snow !! depths of full levels (m)
2359	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: zi_snow !! depths of intermediate levels (m)
2360
2361	!! 0.4 Local variables
2362
2363	REAL(r_std), DIMENSION(kjpindex,nvm) :: z_alpha !! parameter of the geometric series
2364	INTEGER(i_std) :: il,it, ix, iv
2365	INTEGER(i_std) :: it_beg,it_end
2366	INTEGER(i_std), PARAMETER :: niter = 10
2367	REAL(r_std), DIMENSION(kjpindex) :: dxmin
2368	INTEGER(i_std), DIMENSION(kjpindex) :: imin
2369	INTEGER(i_std) :: i,j
2370	REAL(r_std), DIMENSION(kjpindex,nvm) :: xi, xf
2371	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: snowdz_pft
2372
2373	snowdz_pft(:,:,:) = 0.0
2374	DO il = 1,nsnow
2375	DO iv = 1, nvm
2376	WHERE ( veget_mask_2d(:,iv) )
2377	snowdz_pft(:,il,iv) = snowdz(:,il)
2378	ENDWHERE
2379	ENDDO
2380	ENDDO
2381	!
2382	! calculate snow discretisation
2383	!
2384	WHERE ( veget_mask_2d(:,:) )
2385	zf_snow(:,0,:) = 0.
2386	END WHERE
2387	!
2388	DO il = 1, nsnow
2389	IF ( il .EQ. 1 ) THEN
2390	WHERE ( veget_mask_2d(:,:) )
2391
2392	zi_snow(:,il,:) = snowdz_pft(:,1,:) / 2.
2393
2394	zf_snow(:,il,:) = snowdz_pft(:,1,:)
2395
2396	END WHERE
2397	ENDIF
2398
2399	IF ( il .GT. 1 ) THEN
2400	WHERE ( veget_mask_2d(:,:) )
2401
2402	zi_snow(:,il,:) = zf_snow(:,il-1,:) + snowdz_pft(:,il,:) / 2
2403
2404	zf_snow(:,il,:) = SUM(snowdz_pft(:,1:il,:),2)
2405
2406	END WHERE
2407	ENDIF
2408
2409	ENDDO
2410
2411	DO ix = 1, kjpindex
2412	DO il = 1, nsnow
2413	zi_snow_nopftdim(ix,il) = SUM(zi_snow(ix,il,:)*veget_max(ix,:))
2414	zf_snow_nopftdim(ix,il) = SUM(zf_snow(ix,il,:)*veget_max(ix,:))
2415	END DO
2416	END DO
2417
2418	END SUBROUTINE snowlevels
2419
2420	!!
2421	!================================================================================================================================
2422	!! SUBROUTINE : snow_interpol
2423	!!
2424	!>\BRIEF This routine interpolates oxygen and methane into snow layers
2425	!!
2426	!! DESCRIPTION :
2427	!!
2428	!! RECENT CHANGE(S) : None
2429	!!
2430	!! MAIN OUTPUT VARIABLE(S) :
2431	!!
2432	!! REFERENCE(S) : None
2433	!!
2434	!! FLOWCHART11 : None
2435	!! \n
2436	!_
2437	!================================================================================================================================
2438
2439	SUBROUTINE snow_interpol (kjpindex,snowO2, snowCH4, zi_snow, zf_snow, veget_max, snowdz)
2440
2441	!! 0. Variable and parameter declaration
2442
2443	!! 0.1 Input variables
2444
2445	INTEGER(i_std), INTENT(in) :: kjpindex
2446	REAL(r_std), DIMENSION(kjpindex,nsnow), INTENT(in) :: snowdz !! snow depth at each layer [m]
2447	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: veget_max !! maximum vegetation fraction
2448
2449	!! 0.2 Output variables
2450
2451	!! 0.3 Modified variables
2452
2453	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: snowO2 !! snow oxygen (g O2/m**3 air)
2454	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: snowCH4 !! snow methane (g CH4/m**3 air), needed just for num. scheme
2455	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm), INTENT(inout) :: zf_snow !! depths at full levels
2456	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: zi_snow !! depths at intermediate levels
2457
2458	!! 0.4 Local variables
2459	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: isnow !! index of first old layer that is deeper
2460	INTEGER(i_std), DIMENSION(kjpindex,nsnow,nvm) :: i1,i2 !! indices of the layers used for the inter- or extrapolation
2461	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: snowO2o !! initial snow oxygen (g O2/m**3 air)
2462	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: snowCH4o !! initial snow methane (g CH4/m**3 air)
2463	REAL(r_std), DIMENSION(kjpindex,nvm) :: dzio !! initial distance between two levels
2464	INTEGER(i_std) :: il, it, ip, ill, iv !! indices
2465	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm) :: zfo !! initial depths at full levels
2466	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: zio !! initial depths at intermediate levels
2467
2468
2469
2470	! 1. save old discretisation and temperatures
2471
2472	zio(:,:,:) = zi_snow(:,:,:)
2473
2474	zfo(:,:,:) = zf_snow(:,:,:)
2475
2476	snowO2o(:,:,:) = snowO2(:,:,:)
2477	snowCH4o(:,:,:) = snowCH4(:,:,:)
2478
2479	! 2. new discretisation
2480
2481	CALL snowlevels( kjpindex, snowdz, zi_snow, zf_snow, veget_max)
2482
2483	! 3. for each new intermediate layer, look for the first old intermediate
2484	! layer that is deeper
2485
2486	DO il = 1, nsnow
2487
2488	isnow(:,il,:) = -1
2489
2490	DO ill = nsnow,1,-1
2491
2492	WHERE ( zio(:,ill,:) .GT. zi_snow(:,il,:) .AND. veget_mask_2d(:,:) )
2493
2494	isnow(:,il,:) = ill
2495
2496	ENDWHERE
2497
2498	ENDDO
2499
2500	ENDDO
2501
2502	! 4. determine which levels to take for the inter- or extrapolation
2503
2504
2505	DO ip = 1, kjpindex
2506	DO iv = 1, nvm
2507	IF ( veget_mask_2d(ip,iv) ) THEN
2508	DO il = 1, nsnow
2509	!
2510	IF ( isnow(ip,il,iv) .EQ. 1 ) THEN
2511	!
2512	! 4.1 first old layer is below new layer:
2513	! extrapolation from layers 1 and 2
2514	!
2515	i1(ip,il,iv) = 1
2516	i2(ip,il,iv) = 2
2517	!
2518	ELSEIF ( isnow(ip,il,iv) .EQ. -1 ) THEN
2519	!
2520	! 4.2 new layer is below last old layer:
2521	! extrapolation from layers nsnow-1 and nsnow
2522	!
2523	i1(ip,il,iv) = nsnow-1
2524	i2(ip,il,iv) = nsnow
2525	!
2526	ELSE
2527	!
2528	! 4.3 new layer is between two old layers: interpolation
2529	!
2530	i1(ip,il,iv) = isnow(ip,il,iv)-1
2531	i2(ip,il,iv) = isnow(ip,il,iv)
2532	!
2533	ENDIF
2534
2535	ENDDO
2536	ENDIF
2537	ENDDO
2538	ENDDO
2539
2540	! 5. inter- or extrapolate
2541
2542	DO ip = 1, kjpindex
2543	DO iv = 1, nvm
2544	IF ( veget_mask_2d(ip,iv) ) THEN
2545	DO il = 1, nsnow
2546	dzio(ip,iv) = zio(ip,i2(ip,il,iv),iv) - zio(ip,i1(ip,il,iv),iv)
2547
2548	IF ( dzio(ip,iv) .GT. min_stomate ) THEN
2549
2550	snowO2(ip,il,iv) = snowO2o(ip,i1(ip,il,iv),iv) + &
2551	( zi_snow(ip,il,iv) - zio(ip,i1(ip,il,iv),iv) ) / dzio(ip,iv) * &
2552	( snowO2o(ip,i2(ip,il,iv),iv) - snowO2o(ip,i1(ip,il,iv),iv) )
2553	snowCH4(ip,il,iv) = snowCH4o(ip,i1(ip,il,iv),iv) + &
2554	( zi_snow(ip,il,iv) - zio(ip,i1(ip,il,iv),iv) ) / dzio(ip,iv) * &
2555	( snowCH4o(ip,i2(ip,il,iv),iv) - snowCH4o(ip,i1(ip,il,iv),iv) )
2556
2557	ELSE
2558
2559	snowO2(ip,il,iv) = snowO2o(ip,i1(ip,il,iv),iv)
2560	snowCH4(ip,il,iv) = snowCH4o(ip,i1(ip,il,iv),iv)
2561
2562	ENDIF
2563
2564	ENDDO
2565	ENDIF
2566	ENDDO
2567
2568	ENDDO
2569	END SUBROUTINE snow_interpol
2570
2571	!!
2572	!================================================================================================================================
2573	!! SUBROUTINE : stomate_soil_carbon_discretization_clear
2574	!!
2575	!>\BRIEF
2576	!!
2577	!! DESCRIPTION :
2578	!!
2579	!! RECENT CHANGE(S) : None
2580	!!
2581	!! MAIN OUTPUT VARIABLE(S) :
2582	!!
2583	!! REFERENCE(S) : None
2584	!!
2585	!! FLOWCHART11 : None
2586	!! \n
2587	!_
2588	!================================================================================================================================
2589	SUBROUTINE stomate_soil_carbon_discretization_clear()
2590	IF (ALLOCATED(veget_mask_2d)) DEALLOCATE(veget_mask_2d)
2591	IF (ALLOCATED(heights_snow)) DEALLOCATE(heights_snow)
2592	IF (ALLOCATED(zf_soil)) DEALLOCATE(zf_soil)
2593	IF (ALLOCATED(zi_soil)) DEALLOCATE(zi_soil)
2594	IF (ALLOCATED(zf_snow)) DEALLOCATE(zf_snow)
2595	IF (ALLOCATED(zi_snow)) DEALLOCATE(zi_snow)
2596	IF (ALLOCATED(alphaO2_soil )) DEALLOCATE(alphaO2_soil )
2597	IF (ALLOCATED(betaO2_soil )) DEALLOCATE(betaO2_soil )
2598	IF (ALLOCATED(alphaCH4_soil )) DEALLOCATE(alphaCH4_soil )
2599	IF (ALLOCATED(betaCH4_soil )) DEALLOCATE(betaCH4_soil )
2600	IF (ALLOCATED(alphaO2_snow )) DEALLOCATE(alphaO2_snow )
2601	IF (ALLOCATED(betaO2_snow )) DEALLOCATE(betaO2_snow )
2602	IF (ALLOCATED(alphaCH4_snow )) DEALLOCATE(alphaCH4_snow )
2603	IF (ALLOCATED(betaCH4_snow )) DEALLOCATE(betaCH4_snow )
2604	IF (ALLOCATED(zf_coeff_snow )) DEALLOCATE(zf_coeff_snow )
2605	IF (ALLOCATED(zi_coeff_snow )) DEALLOCATE(zi_coeff_snow )
2606	IF (ALLOCATED(mu_snow )) DEALLOCATE(mu_snow )
2607	IF (ALLOCATED(deepSOM_pftmean )) DEALLOCATE(deepSOM_pftmean )
2608
2609	END SUBROUTINE stomate_soil_carbon_discretization_clear
2610
2611	!!
2612	!================================================================================================================================
2613	!! SUBROUTINE : initialize_yedoma_carbonstocks
2614	!!
2615	!>\BRIEF This routine intialize soil carbon in yedoma region
2616	!!
2617	!! DESCRIPTION :
2618	!!
2619	!! RECENT CHANGE(S) : None
2620	!!
2621	!! MAIN OUTPUT VARIABLE(S) :
2622	!!
2623	!! REFERENCE(S) : None
2624	!!
2625	!! FLOWCHART11 : None
2626	!! \n
2627	!_
2628	!================================================================================================================================
2629
2630	SUBROUTINE initialize_yedoma_carbonstocks(kjpindex, lalo, soilsom_a, soilsom_s, soilsom_p, &
2631	yedoma_map_filename, yedoma_depth, yedoma_cinit_act, yedoma_cinit_slo, yedoma_cinit_pas, altmax_ind)
2632
2633	!! 0. Variable and parameter declaration
2634
2635	!! 0.1 Input variables
2636
2637	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
2638	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo !! geographic lat/lon
2639	CHARACTER(LEN=80), INTENT (in) :: yedoma_map_filename !! yedoma map
2640	REAL(r_std), INTENT(in) :: yedoma_depth !! depth of yedoma carbon stock
2641	REAL(r_std), INTENT(in) :: yedoma_cinit_act !! initial active soil C concentration
2642	REAL(r_std), INTENT(in) :: yedoma_cinit_slo !! initial slow soil C concentration
2643	REAL(r_std), INTENT(in) :: yedoma_cinit_pas !! initial passive soil C concentration
2644	INTEGER(i_std), DIMENSION(kjpindex,nvm),INTENT(in) :: altmax_ind !! Maximum over the year active-layer index
2645
2646	!! 0.2 Output variables
2647
2648	!! 0.3 Modified variables
2649
2650	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: soilsom_a !! active soil C concentration
2651	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: soilsom_s !! slow soil C concentration
2652	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: soilsom_p !! passive soil C concentration
2653
2654	!! 0.4 Local variables
2655	REAL(r_std), DIMENSION(kjpindex) :: yedoma
2656	INTEGER(i_std) :: il, ils, ip, ix, iy, imin, jmin, ier, iv
2657	REAL(r_std) :: dlon, dlonmin, dlat, dlatmin
2658	INTEGER(i_std) :: iml, jml, lml, tml, fid
2659	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: xx,yy, yedoma_file
2660	REAL(r_std),ALLOCATABLE,DIMENSION(:) :: x,y
2661	REAL(r_std) :: lev(1), date, dt
2662	INTEGER(i_std) :: itau(1)
2663	INTEGER(i_std) :: yedoma_depth_index, iz
2664	REAL(r_std), DIMENSION(kjpindex) :: NC_yedoma_act !! NC ratio of the active pool for yedoma
2665	REAL(r_std), DIMENSION(kjpindex) :: NC_yedoma_slo !! NC ratio of the slow pool for yedoma
2666	REAL(r_std), DIMENSION(kjpindex) :: NC_yedoma_pas !! NC ratio of the passive pool for yedoma
2667	! plus bas, on prend la temperature lue dans un fichier climato si celui-ci existe
2668
2669	IF ( yedoma_map_filename .EQ. "NONE" ) THEN
2670	yedoma(:) = zero
2671	ELSE IF ( yedoma_map_filename .EQ. "EVERYWHERE" ) THEN
2672	yedoma(:) = 1.
2673	ELSE
2674	CALL flininfo(yedoma_map_filename,iml, jml, lml, tml, fid)
2675
2676	ALLOCATE (yy(iml,jml),stat=ier)
2677	IF (ier /= 0) CALL ipslerr_p(3,'initialize_yedoma_carbonstocks', 'Pb in alloc for yy','','')
2678
2679	ALLOCATE (xx(iml,jml),stat=ier)
2680	IF (ier /= 0) CALL ipslerr_p(3,'initialize_yedoma_carbonstocks', 'Pb in alloc for xx','','')
2681
2682	ALLOCATE (x(iml),stat=ier)
2683	IF (ier /= 0) CALL ipslerr_p(3,'initialize_yedoma_carbonstocks', 'Pb in alloc for x','','')
2684
2685	ALLOCATE (y(jml),stat=ier)
2686	IF (ier /= 0) CALL ipslerr_p(3,'initialize_yedoma_carbonstocks', 'Pb in alloc for y','','')
2687
2688	ALLOCATE (yedoma_file(iml,jml),stat=ier)
2689	IF (ier /= 0) CALL ipslerr_p(3,'initialize_yedoma_carbonstocks', 'Pb in alloc for yedoma_file','','')
2690
2691	CALL flinopen (yedoma_map_filename, .FALSE., iml, jml, lml, &
2692	xx, yy, lev, tml, itau, date, dt, fid)
2693	CALL flinget (fid, 'yedoma', iml, jml, lml, tml, &
2694	1, 1, yedoma_file)
2695	CALL flinclo (fid)
2696	! On suppose que le fichier est regulier.
2697	! Si ce n'est pas le cas, tant pis. Les temperatures seront mal
2698	! initialisees et puis voila. De toute maniere, il faut avoir
2699	! l'esprit mal tourne pour avoir l'idee de faire un fichier de
2700	! climatologie avec une grille non reguliere.
2701	x(:) = xx(:,1)
2702	y(:) = yy(1,:)
2703	! prendre la valeur la plus proche
2704	DO ip = 1, kjpindex
2705	dlonmin = HUGE(1.)
2706	DO ix = 1,iml
2707	dlon = MIN( ABS(lalo(ip,2)-x(ix)), ABS(lalo(ip,2)+360.-x(ix)), ABS(lalo(ip,2)-360.-x(ix)) )
2708	IF ( dlon .LT. dlonmin ) THEN
2709	imin = ix
2710	dlonmin = dlon
2711	ENDIF
2712	ENDDO
2713	dlatmin = HUGE(1.)
2714	DO iy = 1,jml
2715	dlat = ABS(lalo(ip,1)-y(iy))
2716	IF ( dlat .LT. dlatmin ) THEN
2717	jmin = iy
2718	dlatmin = dlat
2719	ENDIF
2720	ENDDO
2721	yedoma(ip) = yedoma_file(imin,jmin)
2722	ENDDO
2723	DEALLOCATE (yy)
2724	DEALLOCATE (xx)
2725	DEALLOCATE (x)
2726	DEALLOCATE (y)
2727	DEALLOCATE (yedoma_file)
2728	ENDIF
2729
2730	yedoma_depth_index = 0
2731	DO iz = 1, ngrnd
2732	IF (znt(iz) .LE. yedoma_depth ) yedoma_depth_index = yedoma_depth_index + 1
2733	END DO
2734	WRITE(,) 'yedoma_depth_index ', yedoma_depth_index, ' at depth ', yedoma_depth
2735
2736	IF ( yedoma_depth_index .GT. 0) THEN
2737	DO ix = 1, kjpindex
2738	DO iv = 2, nvm !!! no yedoma carbon for PFT zero.
2739	IF ( veget_mask_2d(ix,iv) ) THEN
2740	DO iz = 1, yedoma_depth_index
2741	IF (yedoma(ix) .GT. 0.) THEN
2742	IF ( iz .GE. altmax_ind(ix,iv) ) THEN !!! only put yedoma carbon at base of and below the active layer
2743
2744	!MERGE: provisory fix for NC ratio should be done in a cleaner way later
2745	NC_yedoma_act(:)=0.1
2746	NC_yedoma_slo(:)=0.1
2747	NC_yedoma_pas(:)=0.1
2748	CALL getin_p('NC_yedoma_act',NC_yedoma_act)
2749	CALL getin_p('NC_yedoma_slo',NC_yedoma_slo)
2750	CALL getin_p('NC_yedoma_pas',NC_yedoma_pas)
2751
2752	soilsom_a(ix, iz,iv,icarbon) = yedoma_cinit_act
2753	soilsom_s(ix, iz,iv,icarbon) = yedoma_cinit_slo
2754	soilsom_p(ix, iz,iv,icarbon) = yedoma_cinit_pas
2755	soilsom_a(ix, iz,iv,initrogen) = yedoma_cinit_act * NC_yedoma_act(ix)
2756	soilsom_s(ix, iz,iv,initrogen) = yedoma_cinit_slo * NC_yedoma_slo(ix)
2757	soilsom_p(ix, iz,iv,initrogen) = yedoma_cinit_pas * NC_yedoma_pas(ix)
2758	ELSE
2759	soilsom_a(ix, iz,iv,:) = zero
2760	soilsom_s(ix, iz,iv,:) = zero
2761	soilsom_p(ix, iz,iv,:) = zero
2762	ENDIF
2763	ELSE
2764	soilsom_a(ix, iz,iv,:) = zero
2765	soilsom_s(ix, iz,iv,:) = zero
2766	soilsom_p(ix, iz,iv,:) = zero
2767	END IF
2768	END DO
2769	ENDIF
2770	ENDDO
2771	ENDDO
2772	ENDIF
2773
2774	END SUBROUTINE initialize_yedoma_carbonstocks
2775	!!
2776	!================================================================================================================================
2777	!! SUBROUTINE : sominput
2778	!!
2779	!>\BRIEF This routine calculate carbon input to the soil
2780	!!
2781	!! DESCRIPTION :
2782	!!
2783	!! RECENT CHANGE(S) : None
2784	!!
2785	!! MAIN OUTPUT VARIABLE(S) :
2786	!!
2787	!! REFERENCE(S) : None
2788	!!
2789	!! FLOWCHART11 : None
2790	!! \n
2791	!_
2792	!================================================================================================================================
2793	SUBROUTINE sominput(kjpindex,time_step,time,tprof,tsurf,hslong, dayno,z_root,altmax, &
2794	soilsom_a, soilsom_s, soilsom_p, som_input, veget_max, rprof)
2795
2796	!! 0. Variable and parameter declaration
2797
2798	!! 0.1 Input variables
2799
2800	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
2801	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
2802	REAL(r_std), INTENT(in) :: time
2803	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: tprof !! Soil temperature (K)
2804	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: tsurf !! Surface temperature (K)
2805	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: hslong !! deep soil humidity
2806	INTEGER(i_std), INTENT(in) :: dayno !! current day of year
2807	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: z_root !! the rooting depth
2808	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: altmax !! Maximum over the year active-layer thickness
2809	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in) :: veget_max !! Maximum fraction of vegetation type
2810	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements), INTENT(in) :: som_input !! quantity of organic matter going into the SOM pools from litter decomposition (gC/(m**2 of ground)/day)
2811
2812	!! 0.2 Output variables
2813
2814
2815	!! 0.3 Modified variables
2816
2817	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: soilsom_a !! active soil organic matter
2818	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: soilsom_s !! slow soil organic matter
2819	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: soilsom_p !! passive soil organic matter
2820
2821	!! 0.4 Local variables
2822
2823	REAL(r_std) :: dt !! Time step \f$(dt_sechiba one_day^{-1})$\f
2824	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements) :: dsom_litter !! depth-integrated carbon input due to litter decomposition
2825	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements) :: som_input_finite
2826	REAL(r_std), DIMENSION(kjpindex,nvm) :: intdep !! integral depth of carbon deposition
2827	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements) :: sominp_correction
2828	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements) :: som_input_TS
2829	LOGICAL, SAVE :: firstcall = .TRUE.
2830	!$OMP THREADPRIVATE(firstcall)
2831	REAL(r_std), DIMENSION(kjpindex,nvm) :: z_lit !! litter input e-folding depth
2832	INTEGER :: il,ic,iv,ix,iele
2833	INTEGER :: ipts, ivm
2834	LOGICAL, SAVE :: check = .FALSE.
2835	!$OMP THREADPRIVATE(check)
2836	REAL(r_std), PARAMETER :: dgyrst = 96.
2837	INTEGER(i_std), SAVE :: id, id2, id3, id4
2838	!$OMP THREADPRIVATE(id)
2839	!$OMP THREADPRIVATE(id2)
2840	!$OMP THREADPRIVATE(id3)
2841	!$OMP THREADPRIVATE(id4)
2842	CHARACTER(LEN=16) :: buf
2843	INTEGER :: recn
2844	LOGICAL, SAVE :: correct_carboninput_vertprof = .TRUE.
2845	!$OMP THREADPRIVATE(correct_carboninput_vertprof)
2846	LOGICAL, SAVE :: new_carbinput_intdepzlit = .FALSE.
2847	!$OMP THREADPRIVATE(new_carbinput_intdepzlit)
2848	REAL(r_std), DIMENSION(ngrnd) :: z_thickness
2849	REAL(r_std), DIMENSION(ngrnd) :: root_prof
2850	REAL(r_std), SAVE :: minaltmax = 0.1
2851	!$OMP THREADPRIVATE(minaltmax)
2852	REAL(r_std), SAVE :: maxaltmax = 2.
2853	!$OMP THREADPRIVATE(maxaltmax)
2854	REAL(r_std), SAVE :: finerootdepthratio = 0.5 !! the ratio of fine root to overall root e-folding depth (for C inputs)
2855	!$OMP THREADPRIVATE(finerootdepthratio)
2856	REAL(r_std), SAVE :: altrootratio = 0.5 !! the maximum ratio of fine root depth to active layer thickness (for C inputs)
2857	!$OMP THREADPRIVATE(altrootratio)
2858	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: rprof !! root depth (m)
2859	INTEGER, save :: tcounter
2860	!$OMP THREADPRIVATE(tcounter)
2861	INTEGER(i_std) :: imbc, igrnd !! Indices (unitless)
2862	INTEGER(i_std) :: inbpools, icarb !! Indices (unitless)
2863	REAL(r_std), DIMENSION(kjpindex,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
2864	!! mass balance check for this routine
2865	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2866	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
2867	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2868	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: pool_start !! Start and end pool of this routine
2869	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2870	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: pool_end !! Start and end pool of this routine
2871	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2872	REAL(r_std), DIMENSION(ngrnd) :: som_profile !! vertical profile (0-1, unitless)
2873	INTEGER(r_std) :: best_layer !! number of the layer of which the depth
2874	!! best matches a specific target depth
2875	!! (0-ngrnd, unitless)
2876	REAL(r_std), DIMENSION(kjpindex,ncarb,ngrnd,nvm,nelements) :: dsom_litter_z !! depth_dependent organic matter input due to litter
2877	!_ ================================================================================================================================
2878
2879
2880	IF (firstcall) THEN
2881
2882	!Config Key = new_carbinput_intdepzlit
2883	!Config Desc = ???
2884	!Config Def = n
2885	!Config If = OK_SOIL_CARBON_DISCRETIZATION
2886	!Config Help =
2887	!Config Units = [flag]
2888	CALL getin_p('new_carbinput_intdepzlit', new_carbinput_intdepzlit)
2889
2890	!
2891	!Config Key = correct_carboninput_vertprof
2892	!Config Desc = ???
2893	!Config Def = n
2894	!Config If = OK_SOIL_CARBON_DISCRETIZATION
2895	!Config Help =
2896	!Config Units = [flag]
2897	CALL getin_p('correct_carboninput_vertprof', correct_carboninput_vertprof)
2898
2899
2900	! Diagnostic output init
2901
2902	IF (check) THEN
2903	tcounter = 1
2904	WRITE(buf,'(I3)') yr_len
2905	id2 = 0
2906	CALL fliocrfd ('alt.nc', (/'geo ','veg ','time'/), (/kjpindex, nvm, -1/), id, id2, 'REPLACE')
2907	CALL fliodefv (id,'time',(/ 3 /),units='seconds since 0000-01-01 00:00:00',v_t=flio_r8)
2908	CALL flioputa (id,'time','title','time')
2909	CALL flioputa (id,'time','calendar',TRIM(buf)//'d')
2910	CALL fliodefv (id,'alt',(/ 1,2,3 /),units='m',v_t=flio_r8)
2911
2912	CALL fliocrfd ('som_litterinput.nc', (/'geo ','carb','veg ','time'/), (/kjpindex,ncarb,nvm,-1/), id3, id4, 'REPLACE')
2913	CALL fliodefv (id3,'time',(/ 4 /),units='seconds since 0000-01-01 00:00:00',v_t=flio_r8)
2914	CALL flioputa (id3,'time','title','time')
2915	CALL flioputa (id3,'time','calendar',TRIM(buf)//'d')
2916	CALL fliodefv (id3,'dsom_litter',(/ 1,2,3,4 /),units='g C / ts',v_t=flio_r8)
2917	CALL fliodefv (id3,'som_input_TS',(/ 1,2,3,4 /),units='g C / ts',v_t=flio_r8)
2918	ENDIF ! check
2919
2920	firstcall = .FALSE.
2921	!
2922	ENDIF ! firstcall
2923
2924	! Calculate thickness of the soil layers for later use. Could be put into
2925	! firstcall loop but it will then needs to be defined as SAVE
2926	DO il = 1, ngrnd
2927	z_thickness(il) = zf_soil(il) - zf_soil(il-1)
2928	END DO
2929
2930	! 1. Litter input and decomposition
2931	!
2932	! add up the soil carbon from all veg pools, and change units from (gC/(m**2 of ground)/day) to gC/m^2 per timestep
2933	som_input_TS(:,:,:,:) = som_input(:,:,:,:)*time_step/one_day
2934
2935	!
2936	! 2. Carbon input e-folding depth. We distribute with e-depth = min(z_root,intdep)
2937	! and integral depth = min(altmax,z_org)
2938	! e-folding depth cannot be greater than integral depth
2939	IF ( .NOT. new_carbinput_intdepzlit ) THEN
2940	!WRITE(numout,*) 'z_root, ', z_root(:,:)
2941	! change to make intdep equal to z_root alone
2942	z_lit(:,:) = z_root(:,:)
2943	intdep(:,:) = z_root(:,:)
2944	ELSE
2945	!change to separate e-folding depths for roots from total depth over which to integrate
2946	! z_lit is the e-folding depth
2947	z_lit(:,:) = MIN(rprof(:,:)finerootdepthratio, altmax(:,:)altrootratio)
2948	! intdep is the maximum depth of integration;
2949	intdep(:,:) = MIN(altmax(:,:), maxaltmax)
2950	ENDIF
2951
2952	!
2953	! 3. Carbon input.
2954	!
2955	dsom_litter_z(:,:,:,:,:) = zero
2956	dsom_litter(:,:,:,:)=zero
2957	DO ipts = 1,kjpindex
2958	DO ivm = 1,nvm
2959
2960	! The original code was using zi_soil. If we use zi_soil here we introduce
2961	! inconsistencies in the mass balance because we will distribute the
2962	! carbon over zi_soil but then integrate it over zf_soil. Given the top
2963	! layers are rather thin and the number of layers seems arbitrary anyway,
2964	! it seems acceptable to move from zi_soil to zf_soil.
2965	! We always need to integrate to the exact depth of one of the layers,
2966	! otherwise we will add a lot of complexity to close the mass balance. If
2967	! intdep(:,:) .LT. zi_soil(2), we set the depth to an exact layer. A solution
2968	! was added in case intdep(:,:) .GT. zi_soil(2). In that case the nearest
2969	! exact layer was used to distribute the carbon and nitrogen.
2970	IF ( intdep(ipts,ivm) .LT. zf_soil(2) ) THEN
2971	! Litter is decomposed somehow (?) even when alt == 0. To avoid carbon loss,
2972	! we distribute this carbon within the first 2 soil layers when alt == 0
2973	! Adding EPSILON(zero) avoids problems with the next IF-statement where
2974	! zf_soil is compared against intdep.
2975	intdep(ipts,ivm) = zf_soil(2) + EPSILON(0.)
2976	ELSE
2977	! Find the layer that best represents the depth of interest
2978	best_layer = MINLOC(ABS(zf_soil(:) - intdep(ipts,ivm)),dim=1)
2979	intdep(ipts,ivm) = zf_soil(best_layer) + EPSILON(0.)
2980	ENDIF
2981	IF ( z_lit(ipts,ivm) .LT. zf_soil(2) ) THEN
2982	! Litter is decomposed somehow (?) even when alt == 0. To avoid carbon loss,
2983	! we distribute this carbon within the first 2 soil layers when alt == 0
2984	z_lit(ipts,ivm) = zf_soil(2)
2985	ELSE
2986	! Find the layer that best represents the depth of interest
2987	best_layer = MINLOC(ABS(zf_soil(:) - z_lit(ipts,ivm)),dim=1)
2988	z_lit(ipts,ivm) = zf_soil(best_layer)
2989	ENDIF
2990
2991	! Initialize
2992	som_profile(:) = zero
2993
2994	! Rewritten the calculation of the redistribution factor. Irrespective
2995	! of whether zi_soil or zf_soil is used, the original code appeared
2996	! bugged to me. It did not account for the first layer. It divided by
2997	! som_input_TS by z_lit which implies it redistributed the input
2998	! twice (the division is a uniform redistribution, the EXP an exponential
2999	! redistribution). Division by z_lit results in a change in the units
3000	! this conversion is now taken care of later in the code.
3001	DO il = 1, ngrnd
3002	! Calculate the som profile
3003	IF ( zf_soil(il) .LT. intdep(ipts,ivm) .AND. veget_mask_2d(ipts,ivm) ) THEN
3004	! The first part of the equation, i.e., 1/(1-EXP(-intdep/z_lit)
3005	! is a scaling factor to ensure that the profile down to
3006	! intdep will add up to 1 (note that z_lit determines the shape of
3007	! exponential function that will be used to redistribute the som
3008	! over the vertical layers).
3009	som_profile(il) = un / ( un - EXP( -intdep(ipts,ivm) / z_lit(ipts,ivm) ) ) * &
3010	( EXP(-zf_soil(il-1)/z_lit(ipts,ivm)) - &
3011	EXP( -zf_soil(il)/z_lit(ipts,ivm) ) )
3012	ELSE
3013	! This layer is too far down and is no longer accounted for
3014	som_profile(il) = zero
3015	END IF
3016	END DO
3017
3018	! Check for errors. If none, overcome possible precision issues
3019	IF (veget_mask_2d(ipts,ivm)) THEN
3020	IF (ABS(SUM(som_profile(:))-un).GT.min_stomate) THEN
3021	! The above calculation is less precise than expeceted
3022	!IF (err_act.GT.1) THEN
3023	WRITE(numout,*) 'ipts, ivm, ', ipts,ivm
3024	WRITE(numout,*) 'zf_soil, ', zf_soil(:)
3025	WRITE(numout,*) 'intdep, z_lit,', intdep(ipts,ivm),z_lit(ipts,ivm)
3026	WRITE(numout,*) 'som_profile, ', som_profile(:)
3027	WRITE(numout,*) 'Sum of som_profile, ', SUM(som_profile(:))
3028	CALL ipslerr_p(3,'som_input','The sum of the som profile differs from 1',&
3029	'mass balance will not be closed','')
3030	!ENDIF
3031	ELSE
3032	! The above calculation should result in a som_profile of 0.9999999 which is
3033	! pretty good but not good enough to keep the mass balance closed. Put the
3034	! residual in the top layer (just because that layer should always be present)
3035	som_profile(1) = un - SUM(som_profile(2:ngrnd))
3036	END IF
3037	END IF
3038
3039	! Use the profile to redistribute the som_input
3040	DO ic = 1, ncarb
3041	DO iele = 1, nelements
3042	DO il = 1, ngrnd
3043	! Divide by the tickness of the layer, not the tickness of the whole
3044	! profile (which is the case when using z_lit). The unit of som_input_Ts
3045	! is gC/m^2 per timestep. We want dsom_litter_z to be in gC/m^3 per
3046	! timestep so divide by the layer depth over which the som input is to
3047	! be distributed.
3048	dsom_litter_z(ipts,ic,il,ivm,iele) = som_input_TS(ipts,ic,ivm,iele) * &
3049	som_profile(il)/ z_thickness(il)
3050	! The original code already used zf_soil. This may have resulted in an
3051	! inconsitency. Given that zi_soil has now been changed to zf_soil. The
3052	! original code is consistent with the new approach that uses zf_soil
3053	dsom_litter(ipts,ic,ivm,iele) = dsom_litter(ipts,ic,ivm,iele) + &
3054	dsom_litter_z(ipts,ic,il,ivm,iele) * z_thickness(il)
3055	END DO
3056	END DO
3057	END DO
3058
3059	END DO
3060	END DO
3061
3062	! Update the active, slow and passive soilsom pools
3063	DO il = 1, ngrnd
3064	DO iele = 1, nelements
3065	WHERE ( veget_mask_2d(:,:) )
3066	soilsom_a(:,il,:,iele) = soilsom_a(:,il,:,iele) + dsom_litter_z(:,iactive,il,:,iele)
3067	soilsom_s(:,il,:,iele) = soilsom_s(:,il,:,iele) + dsom_litter_z(:,islow,il,:,iele)
3068	soilsom_p(:,il,:,iele) = soilsom_p(:,il,:,iele) + dsom_litter_z(:,ipassive,il,:,iele)
3069	END WHERE
3070	ENDDO
3071	END DO
3072
3073
3074	! Diagnostic output
3075	IF (check) THEN
3076	recn = NINT(time/time_step)
3077	tcounter = tcounter + 1
3078	WRITE(,) 'sominput check: output to .nc number',recn
3079	WRITE(,) 'time',time
3080	WRITE(,) 'time_step',time_step
3081
3082	CALL flioputv (id,'time', time, (/ tcounter /) )
3083	CALL flioputv (id,'alt', altmax(:,:), start = (/ 1, 1, tcounter /), count = (/ kjpindex, nvm, 1 /) )
3084	CALL fliosync(id)
3085
3086	CALL flioputv (id3,'time', time, (/ tcounter /) )
3087	CALL flioputv (id3,'som_input_TS', som_input_TS(:,:,:,:), start = (/ 1,1, 1, 1, tcounter /), &
3088	count = (/ kjpindex, ncarb, nvm, nelements, 1 /) )
3089	CALL flioputv (id3,'dsom_litter', dsom_litter(:,:,:,:), start = (/ 1,1, 1, 1, tcounter /), &
3090	count = (/ kjpindex, ncarb, nvm,nelements, 1 /) )
3091	CALL fliosync(id3)
3092	ENDIF
3093
3094
3095	END SUBROUTINE sominput
3096
3097	!!
3098	!================================================================================================================================
3099	!! SUBROUTINE : cryoturbate
3100	!!
3101	!>\BRIEF This routine calculates cryoturbation process
3102	!!
3103	!! DESCRIPTION :
3104	!!
3105	!! RECENT CHANGE(S) : None
3106	!!
3107	!! MAIN OUTPUT VARIABLE(S) :
3108	!!
3109	!! REFERENCE(S) : None
3110	!!
3111	!! FLOWCHART11 : None
3112	!! \n
3113	!_
3114	!================================================================================================================================
3115
3116	SUBROUTINE cryoturbate(kjpindex, time_step, dayno, altmax_ind, deepSOM_a, deepSOM_s, deepSOM_p, &
3117	action, diff_k_const, bio_diff_k_const, altmax_lastyear, fixed_cryoturbation_depth)
3118
3119	!! 0. Variable and parameter declaration
3120
3121	!! 0.1 Input variables
3122
3123	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
3124	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
3125	INTEGER(i_std), INTENT(in) :: dayno !! number of the day in the current year
3126	INTEGER(i_std), DIMENSION(kjpindex,nvm),INTENT(in) :: altmax_ind !! Maximum over the year active-layer index
3127	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: altmax_lastyear !! Maximum over the year active-layer thickness
3128	CHARACTER(LEN=*), INTENT(in) :: action !! what to do
3129	REAL(r_std), INTENT(in) :: diff_k_const
3130	REAL(r_std), INTENT(in) :: bio_diff_k_const
3131
3132	!! 0.2 Output variables
3133
3134	!! 0.3 Modified variables
3135	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deepSOM_a !! soil soil organic matter (g/m**3) active
3136	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deepSOM_s !! soil soil organic matter (g/m**3) slow
3137	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deepSOM_p !! soil soil organic matter (g/m**3) passive
3138	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(inout) :: fixed_cryoturbation_depth !! depth to hold cryoturbation to for fixed runs
3139
3140	!! 0.4 Local variables
3141	LOGICAL, SAVE :: firstcall = .TRUE.
3142	!$OMP THREADPRIVATE(firstcall)
3143	LOGICAL, SAVE :: use_new_cryoturbation
3144	!$OMP THREADPRIVATE(use_new_cryoturbation)
3145	INTEGER, SAVE :: cryoturbation_method
3146	!$OMP THREADPRIVATE(cryoturbation_method)
3147	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: altSOM_a_old !! soil organic matter (g/m**2) active integrated over active layer before cryoturbation
3148	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: altSOM_s_old !! soil organic matter (g/m**2) slow integrated over active layer before cryoturbation
3149	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: altSOM_p_old !! soil organic matter (g/m**2) passive integrated over active layer before cryoturbation
3150	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: altSOM_a
3151	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: altSOM_s
3152	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: altSOM_p
3153	INTEGER(i_std), PARAMETER :: n_totakefrom = 3 !! how many surface layers to subtract from in mass balance
3154	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: surfSOM_totake_a !! active soil organic matter to subtract from surface layers to maintain mass balance (g/m**3)
3155	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: surfSOM_totake_s !! slow soil organic matter to subtract from surface layers to maintain mass balance (g/m**3)
3156	REAL(r_std), DIMENSION(kjpindex,nvm,nelements) :: surfSOM_totake_p !! passive soil organic matter to subtract from surface layers to maintain mass balance (g/m**3)
3157	REAL(r_std), DIMENSION(kjpindex,nvm) :: error_a
3158	REAL(r_std), DIMENSION(kjpindex,nvm) :: error_s
3159	REAL(r_std), DIMENSION(kjpindex,nvm) :: error_p
3160	INTEGER(i_std) :: ip, il, ier, iv, iele
3161	CHARACTER(LEN=20), SAVE :: last_action = 'not called'
3162	!$OMP THREADPRIVATE(last_action)
3163	INTEGER(i_std) :: cryoturb_date
3164	REAL(r_std), SAVE :: max_cryoturb_alt
3165	!$OMP THREADPRIVATE(max_cryoturb_alt)
3166	REAL(r_std), SAVE :: min_cryoturb_alt
3167	!$OMP THREADPRIVATE(min_cryoturb_alt)
3168	REAL(r_std), SAVE :: bioturbation_depth
3169	!$OMP THREADPRIVATE(bioturbation_depth)
3170	LOGICAL, SAVE :: reset_fixed_cryoturbation_depth = .FALSE.
3171	!$OMP THREADPRIVATE(reset_fixed_cryoturbation_depth)
3172	LOGICAL, SAVE :: use_fixed_cryoturbation_depth = .FALSE.
3173	!$OMP THREADPRIVATE(use_fixed_cryoturbation_depth)
3174	REAL(r_std), DIMENSION(kjpindex,nvm) :: cryoturbation_depth
3175
3176
3177	! 1. ensure that we do not repeat actions
3178	!
3179	IF ( action .EQ. last_action ) THEN
3180	!
3181	WRITE(,) 'CANNOT TAKE THE SAME ACTION TWICE: ',TRIM(action)
3182	STOP
3183	!
3184	ENDIF
3185
3186	IF (firstcall) THEN
3187
3188	! 2. faire les trucs du debut
3189
3190	! 2.1 allocation des variables
3191	ALLOCATE (xe_a(kjpindex,nvm),stat=ier)
3192	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for xe_a','','')
3193
3194	ALLOCATE (xe_s(kjpindex,nvm),stat=ier)
3195	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for xe_s','','')
3196
3197	ALLOCATE (xe_p(kjpindex,nvm),stat=ier)
3198	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for xe_p','','')
3199
3200	ALLOCATE (xc_cryoturb(kjpindex,ngrnd,nvm),stat=ier)
3201	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for xc_cryoturb','','')
3202
3203	ALLOCATE (xd_cryoturb(kjpindex,ngrnd,nvm),stat=ier)
3204	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for xd_cryoturb','','')
3205
3206	ALLOCATE (alpha_a(kjpindex,ngrnd,nvm,nelements),stat=ier)
3207	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for alpha_a','','')
3208	alpha_a(:,:,:,:)=0.
3209
3210	ALLOCATE (alpha_s(kjpindex,ngrnd,nvm,nelements),stat=ier)
3211	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for alpha_s','','')
3212	alpha_s(:,:,:,:)=0.
3213
3214	ALLOCATE (alpha_p(kjpindex,ngrnd,nvm,nelements),stat=ier)
3215	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for alpha_p','','')
3216	alpha_p(:,:,:,:)=0.
3217
3218	ALLOCATE (mu_soil_rev(kjpindex,nvm),stat=ier)
3219	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for mu_soil_rev','','')
3220	mu_soil_rev(:,:)=0.
3221
3222	ALLOCATE (beta_a(kjpindex,ngrnd,nvm,nelements),stat=ier)
3223	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for beta_a','','')
3224	beta_a(:,:,:,:)=0.
3225
3226	ALLOCATE (beta_s(kjpindex,ngrnd,nvm,nelements),stat=ier)
3227	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for beta_s','','')
3228	beta_s(:,:,:,:)=0.
3229
3230	ALLOCATE (beta_p(kjpindex,ngrnd,nvm,nelements),stat=ier)
3231	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for beta_p','','')
3232	beta_p(:,:,:,:)=0.
3233
3234	ALLOCATE (diff_k(kjpindex,ngrnd,nvm),stat=ier)
3235	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for diff_k','','')
3236
3237	ALLOCATE (cryoturb_location(kjpindex,nvm),stat=ier)
3238	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for cryoturb_location','','')
3239
3240	ALLOCATE (bioturb_location(kjpindex,nvm),stat=ier)
3241	IF (ier /= 0) CALL ipslerr_p(3,'cryoturbate', 'Pb in alloc for bioturb_location','','')
3242
3243	cryoturb_location(:,:) = .false.
3244	use_new_cryoturbation = .false.
3245	!
3246	!Config Key = use_new_cryoturbation
3247	!Config Desc =
3248	!Config Def = n
3249	!Config If = OK_SOIL_CARBON_DISCRETIZATION
3250	!Config Help =
3251	!Config Units = [flag]
3252	CALL getin_p('use_new_cryoturbation', use_new_cryoturbation)
3253	!
3254	!Config Key = cryoturbation_method
3255	!Config Desc =
3256	!Config Def = 1
3257	!Config If = OK_SOIL_CARBON_DISCRETIZATION
3258	!Config Help =
3259	!Config Units = []
3260	cryoturbation_method = 4
3261	CALL getin_p('cryoturbation_method', cryoturbation_method)
3262	!
3263	!Config Key = max_cryoturb_alt
3264	!Config Desc =
3265	!Config Def = 1
3266	!Config If = OK_SOIL_CARBON_DISCRETIZATION
3267	!Config Help =
3268	!Config Units = []
3269	max_cryoturb_alt = 3.
3270	CALL getin_p('max_cryoturb_alt',max_cryoturb_alt)
3271	!
3272	!Config Key = min_cryoturb_alt
3273	!Config Desc =
3274	!Config Def = 1
3275	!Config If = OK_SOIL_CARBON_DISCRETIZATION
3276	!Config Help =
3277	!Config Units = []
3278	min_cryoturb_alt = 0.01
3279	CALL getin_p('min_cryoturb_alt',min_cryoturb_alt)
3280	!
3281	!Config Key = reset_fixed_cryoturbation_depth
3282	!Config Desc =
3283	!Config Def = n
3284	!Config If = OK_SOIL_CARBON_DISCRETIZATION
3285	!Config Help =
3286	!Config Units = [flag]
3287	CALL getin_p('reset_fixed_cryoturbation_depth',reset_fixed_cryoturbation_depth)
3288	IF (reset_fixed_cryoturbation_depth) THEN
3289	fixed_cryoturbation_depth = altmax_lastyear
3290	ENDIF
3291	!
3292	!Config Key = use_fixed_cryoturbation_depth
3293	!Config Desc =
3294	!Config Def = n
3295	!Config If = OK_SOIL_CARBON_DISCRETIZATION
3296	!Config Help =
3297	!Config Units = [flag]
3298	CALL getin_p('use_fixed_cryoturbation_depth',use_fixed_cryoturbation_depth)
3299	bioturb_location(:,:) = .false.
3300	!
3301	!Config Key = bioturbation_depth
3302	!Config Desc = maximum bioturbation depth
3303	!Config Def = 2
3304	!Config If = OK_SOIL_CARBON_DISCRETIZATION
3305	!Config Help =
3306	!Config Units = m
3307	bioturbation_depth = 2.
3308	CALL getin_p('bioturbation_depth',bioturbation_depth)
3309
3310	firstcall = .FALSE.
3311	ENDIF
3312
3313	IF ( action .EQ. 'diffuse' ) THEN
3314	! 1. calculate the total soil organic matter in the active layer
3315	altSOM_a_old(:,:,:) = zero
3316	altSOM_s_old(:,:,:) = zero
3317	altSOM_p_old(:,:,:) = zero
3318	altSOM_a(:,:,:) = zero
3319	altSOM_s(:,:,:) = zero
3320	altSOM_p(:,:,:) = zero
3321
3322	DO ip = 1, kjpindex
3323	DO iv = 1, nvm
3324	DO iele = 1, nelements
3325	IF ( cryoturb_location(ip,iv) .OR. bioturb_location(ip,iv) )THEN
3326	! 1. calculate the total soil organic matter
3327	DO il = 1, ngrnd
3328	altSOM_a_old(ip,iv,iele) = altSOM_a_old(ip,iv,iele) + deepSOM_a(ip,il,iv,iele)*(zf_soil(il)-zf_soil(il-1))
3329	altSOM_s_old(ip,iv,iele) = altSOM_s_old(ip,iv,iele) + deepSOM_s(ip,il,iv,iele)*(zf_soil(il)-zf_soil(il-1))
3330	altSOM_p_old(ip,iv,iele) = altSOM_p_old(ip,iv,iele) + deepSOM_p(ip,il,iv,iele)*(zf_soil(il)-zf_soil(il-1))
3331	ENDDO
3332
3333	! 2. diffuse the soil organic matter
3334	deepSOM_a(ip,1,iv,iele) = (deepSOM_a(ip,1,iv,iele)+mu_soil_rev(ip,iv)*beta_a(ip,1,iv,iele)) / &
3335	(1.+mu_soil_rev(ip,iv)*(1.-alpha_a(ip,1,iv,iele)))
3336	deepSOM_s(ip,1,iv,iele) = (deepSOM_s(ip,1,iv,iele)+mu_soil_rev(ip,iv)*beta_s(ip,1,iv,iele)) / &
3337	(1.+mu_soil_rev(ip,iv)*(1.-alpha_s(ip,1,iv,iele)))
3338	deepSOM_p(ip,1,iv,iele) = (deepSOM_p(ip,1,iv,iele)+mu_soil_rev(ip,iv)*beta_p(ip,1,iv,iele)) / &
3339	(1.+mu_soil_rev(ip,iv)*(1.-alpha_p(ip,1,iv,iele)))
3340
3341	DO il = 2, ngrnd
3342	deepSOM_a(ip,il,iv,iele) = alpha_a(ip,il-1,iv,iele)*deepSOM_a(ip,il-1,iv,iele) + beta_a(ip,il-1,iv,iele)
3343	deepSOM_s(ip,il,iv,iele) = alpha_s(ip,il-1,iv,iele)*deepSOM_s(ip,il-1,iv,iele) + beta_s(ip,il-1,iv,iele)
3344	deepSOM_p(ip,il,iv,iele) = alpha_p(ip,il-1,iv,iele)*deepSOM_p(ip,il-1,iv,iele) + beta_p(ip,il-1,iv,iele)
3345	ENDDO
3346
3347	! 3. recalculate the total soil organic matter
3348	DO il = 1, ngrnd
3349	altSOM_a(ip,iv,iele) = altSOM_a(ip,iv,iele) + deepSOM_a(ip,il,iv,iele)*(zf_soil(il)-zf_soil(il-1))
3350	altSOM_s(ip,iv,iele) = altSOM_s(ip,iv,iele) + deepSOM_s(ip,il,iv,iele)*(zf_soil(il)-zf_soil(il-1))
3351	altSOM_p(ip,iv,iele) = altSOM_p(ip,iv,iele) + deepSOM_p(ip,il,iv,iele)*(zf_soil(il)-zf_soil(il-1))
3352	ENDDO
3353
3354
3355	IF ( altSOM_a_old(ip,iv,iele) > min_stomate .AND. &
3356	(ABS(altSOM_a(ip,iv,iele)-altSOM_a_old(ip,iv,iele))/altSOM_a_old(ip,iv,iele).GT. min_stomate) ) THEN
3357	WRITE (numout,*) 'DZ warn: cryoturbate: total C or N not conserved','iele=',iele, 'ip=',ip,'iv=',iv, &
3358	'A,diff=',altSOM_a(ip,iv,iele),altSOM_a_old(ip,iv,iele),altSOM_a(ip,iv,iele)-altSOM_a_old(ip,iv,iele), &
3359	(altSOM_a(ip,iv,iele)-altSOM_a_old(ip,iv,iele))/altSOM_a_old(ip,iv,iele)
3360	CALL ipslerr_p (3,'cryoturbate','','','')
3361	ENDIF
3362
3363	IF ( altSOM_s_old(ip,iv,iele) > min_stomate .AND. &
3364	(ABS(altSOM_s(ip,iv,iele)-altSOM_s_old(ip,iv,iele))/altSOM_s_old(ip,iv,iele).GT. min_stomate) ) THEN
3365	WRITE (numout,*) 'DZ warn: cryoturbate: total C or N not conserved','iele=',iele,'ip=',ip,'iv=',iv, &
3366	'S,diff=',altSOM_s(ip,iv,iele),altSOM_s_old(ip,iv,iele),altSOM_s(ip,iv,iele)-altSOM_s_old(ip,iv,iele), &
3367	(altSOM_s(ip,iv,iele)-altSOM_s_old(ip,iv,iele))/altSOM_s_old(ip,iv,iele)
3368	CALL ipslerr_p (3,'cryoturbate','','','')
3369	ENDIF
3370
3371	IF ( altSOM_p_old(ip,iv,iele) > min_stomate .AND. &
3372	(ABS(altSOM_p(ip,iv,iele)-altSOM_p_old(ip,iv,iele))/altSOM_p_old(ip,iv,iele).GT. min_stomate) ) THEN
3373	WRITE (numout,*) 'DZ warn: cryoturbate: total C or N not conserved','iele=',iele, 'ip=',ip,'iv=',iv, &
3374	'P,diff=',altSOM_p(ip,iv,iele),altSOM_p_old(ip,iv,iele),altSOM_p(ip,iv,iele)-altSOM_p_old(ip,iv,iele), &
3375	(altSOM_p(ip,iv,iele)-altSOM_p_old(ip,iv,iele))/altSOM_p_old(ip,iv,iele)
3376	CALL ipslerr_p (3,'cryoturbate','','','')
3377	ENDIF
3378
3379	! 4. subtract the organic matter in the top layer(s) so that the total organic matter content of the active layer is conserved.
3380	! for now remove this correction term...
3381	! surfC_totake_a(ip,iv) = (altC_a(ip,iv)-altC_a_old(ip,iv))/(zf_soil(altmax_ind(ip,iv))-zf_soil(0))
3382	! surfC_totake_s(ip,iv) = (altC_s(ip,iv)-altC_s_old(ip,iv))/(zf_soil(altmax_ind(ip,iv))-zf_soil(0))
3383	! surfC_totake_p(ip,iv) = (altC_p(ip,iv)-altC_p_old(ip,iv))/(zf_soil(altmax_ind(ip,iv))-zf_soil(0))
3384	! deepC_a(ip,1:altmax_ind(ip,iv),iv) = deepC_a(ip,1:altmax_ind(ip,iv),iv) - surfC_totake_a(ip,iv)
3385	! deepC_s(ip,1:altmax_ind(ip,iv),iv) = deepC_s(ip,1:altmax_ind(ip,iv),iv) - surfC_totake_s(ip,iv)
3386	! deepC_p(ip,1:altmax_ind(ip,iv),iv) = deepC_p(ip,1:altmax_ind(ip,iv),iv) - surfC_totake_p(ip,iv)
3387	!
3388	! ! if negative values appear, we don't subtract the delta-C from top layers
3389	! IF (ANY(deepC_a(ip,1:altmax_ind(ip,iv),iv) .LT. zero) ) THEN
3390	! deepC_a(ip,1:altmax_ind(ip,iv),iv)=deepC_a(ip,1:altmax_ind(ip,iv),iv)+surfC_totake_a(ip,iv)
3391	! IF (altC_a(ip,iv) .GT. zero) THEN
3392	! deepC_a(ip,:,iv)=deepC_a(ip,:,iv)*altC_a_old(ip,iv)/altC_a(ip,iv)
3393	! ENDIF
3394	! ENDIF
3395	! IF (ANY(deepC_s(ip,1:altmax_ind(ip,iv),iv) .LT. zero) ) THEN
3396	! deepC_s(ip,1:altmax_ind(ip,iv),iv)=deepC_s(ip,1:altmax_ind(ip,iv),iv)+surfC_totake_s(ip,iv)
3397	! IF (altC_s(ip,iv) .GT. zero) THEN
3398	! deepC_s(ip,:,iv)=deepC_s(ip,:,iv)*altC_s_old(ip,iv)/altC_s(ip,iv)
3399	! ENDIF
3400	! ENDIF
3401	! IF (ANY(deepC_p(ip,1:altmax_ind(ip,iv),iv) .LT. zero) ) THEN
3402	! deepC_p(ip,1:altmax_ind(ip,iv),iv)=deepC_p(ip,1:altmax_ind(ip,iv),iv)+surfC_totake_p(ip,iv)
3403	! IF (altC_p(ip,iv) .GT. zero) THEN
3404	! deepC_p(ip,:,iv)=deepC_p(ip,:,iv)*altC_p_old(ip,iv)/altC_p(ip,iv)
3405	! ENDIF
3406	! ENDIF
3407
3408	! Consistency check. Potentially add to STRICT_CHECK flag
3409	IF ( ANY(deepSOM_a(ip,:,iv,iele) .LT. zero) ) THEN
3410	WRITE (numout,*) 'cryoturbate: deepSOM_a<0','iele=',iele, &
3411	'ip=',ip,'iv=',iv,'deepSOM_a=',deepSOM_a(ip,:,iv,iele)
3412	CALL ipslerr_p (3,'cryoturbate','','','')
3413	ENDIF
3414	IF ( ANY(deepSOM_s(ip,:,iv,iele) .LT. zero) ) THEN
3415	WRITE (numout,*) 'cryoturbate: deepSOM_s<0','iele=',iele, &
3416	'ip=',ip,'iv=',iv,'deepSOM_s=',deepSOM_s(ip,:,iv,iele)
3417	CALL ipslerr_p (3,'cryoturbate','','','')
3418	ENDIF
3419	IF ( ANY(deepSOM_p(ip,:,iv,iele) .LT. zero) ) THEN
3420	WRITE (numout,*) 'cryoturbate: deepSOM_p<0','iele=',iele, &
3421	'ip=',ip,'iv=',iv,'deepSOM_p=',deepSOM_p(ip,:,iv,iele)
3422	CALL ipslerr_p (3,'cryoturbate','','','')
3423	ENDIF
3424
3425	ENDIF
3426	ENDDO !End loop over nelements
3427	ENDDO
3428	ENDDO
3429
3430
3431	ELSEIF ( action .EQ. 'coefficients' ) THEN
3432	IF (firstcall) THEN
3433	WRITE(,) 'error: initilaizations have to happen before coefficients calculated. we stop.'
3434	STOP
3435	ENDIF
3436
3437	cryoturb_location(:,:) = ( altmax_lastyear(:,:) .LT. max_cryoturb_alt ) &
3438	!In the former vertical discretization scheme the first level was at 0.016 cm; now it's only 0.00048 so we set an equivalent threshold directly as a fixed depth of 1 cm,
3439	.AND. ( altmax_lastyear(:,:) .GE. min_cryoturb_alt ) .AND. veget_mask_2d(:,:)
3440	IF (use_fixed_cryoturbation_depth) THEN
3441	cryoturbation_depth(:,:) = fixed_cryoturbation_depth(:,:)
3442	ELSE
3443	cryoturbation_depth(:,:) = altmax_lastyear(:,:)
3444	ENDIF
3445
3446	bioturb_location(:,:) = ( ( altmax_lastyear(:,:) .GE. max_cryoturb_alt ) .AND. veget_mask_2d(:,:) )
3447
3448	DO ip = 1, kjpindex
3449	DO iv = 1,nvm
3450	IF ( cryoturb_location(ip,iv) ) THEN
3451	!
3452	IF (use_new_cryoturbation) THEN
3453	SELECT CASE(cryoturbation_method)
3454	CASE(1)
3455	!
3456	DO il = 1, ngrnd ! linear dropoff to zero between alt and 2*alt
3457	IF ( zi_soil(il) .LE. cryoturbation_depth(ip,iv) ) THEN
3458	diff_k(ip,il,iv) = diff_k_const
3459	ELSE
3460	diff_k(ip,il,iv) = diff_k_const*(un-MAX(MIN((zi_soil(il)/cryoturbation_depth(ip,iv))-un,un),zero))
3461	ENDIF
3462	END DO
3463	!
3464	CASE(2)
3465	!
3466	DO il = 1, ngrnd ! exponential dropoff with e-folding distace = alt, below the active layer
3467	IF ( zi_soil(il) .LE. cryoturbation_depth(ip,iv) ) THEN
3468	diff_k(ip,il,iv) = diff_k_const
3469	ELSE
3470	diff_k(ip,il,iv) = diff_k_const*(EXP(-MAX((zi_soil(il)/cryoturbation_depth(ip,iv)-un),zero)))
3471	ENDIF
3472	END DO
3473	!
3474	CASE(3)
3475	!
3476	! exponential dropoff with e-folding distace = alt, starting at surface
3477	diff_k(ip,:,iv) = diff_k_const*(EXP(-(zi_soil(:)/cryoturbation_depth(ip,iv))))
3478	!
3479	CASE(4)
3480	!
3481	DO il = 1, ngrnd ! linear dropoff to zero between alt and 3*alt
3482	IF ( zi_soil(il) .LE. cryoturbation_depth(ip,iv) ) THEN
3483	diff_k(ip,il,iv) = diff_k_const
3484	ELSE
3485	diff_k(ip,il,iv) = diff_k_const*(un-MAX(MIN((zi_soil(il)-cryoturbation_depth(ip,iv))/ &
3486	(2.*cryoturbation_depth(ip,iv)),un),zero))
3487	ENDIF
3488	IF ( zf_soil(il) .GT. max_cryoturb_alt ) THEN
3489	diff_k(ip,il,iv) = zero
3490	ENDIF
3491	END DO
3492	!
3493	IF (printlev>=3) WRITE(,) 'cryoturb method 4: ip, iv, diff_k(ip,:,iv): ', ip, iv, diff_k(ip,:,iv)
3494	CASE(5)
3495	!
3496	DO il = 1, ngrnd ! linear dropoff to zero between alt and 3m
3497	IF ( zi_soil(il) .LE. cryoturbation_depth(ip,iv) ) THEN
3498	diff_k(ip,il,iv) = diff_k_const
3499	ELSE
3500	diff_k(ip,il,iv) = diff_k_const*(un-MAX(MIN((zi_soil(il)-cryoturbation_depth(ip,iv))/ &
3501	(3.-cryoturbation_depth(ip,iv)),un),zero))
3502	ENDIF
3503	END DO
3504	!
3505	IF (printlev>=3) WRITE(,) 'cryoturb method 5: ip, iv, diff_k(ip,:,iv): ', ip, iv, diff_k(ip,:,iv)
3506	END SELECT
3507
3508	ELSE ! old cryoturbation scheme
3509	!
3510	diff_k(ip,1:altmax_ind(ip,iv),iv) = diff_k_const
3511	diff_k(ip, altmax_ind(ip,iv)+1,iv) = diff_k_const/10.
3512	diff_k(ip, altmax_ind(ip,iv)+2,iv) = diff_k_const/100.
3513	diff_k(ip,(altmax_ind(ip,iv)+3):ngrnd,iv) = zero
3514	ENDIF
3515	ELSE IF ( bioturb_location(ip,iv) ) THEN
3516	DO il = 1, ngrnd
3517	IF ( zi_soil(il) .LE. bioturbation_depth ) THEN
3518	diff_k(ip,il,iv) = bio_diff_k_const
3519	ELSE
3520	diff_k(ip,il,iv) = zero
3521	ENDIF
3522	END DO
3523	ELSE
3524	diff_k(ip,:,iv) = zero
3525	END IF
3526	END DO
3527	END DO
3528
3529	mu_soil_rev=diff_k(:,1,:)*time_step/(zf_soil(1)-zf_soil(0))/(zi_soil(2)-zi_soil(1))
3530
3531	DO il = 1,ngrnd-1
3532	WHERE ( cryoturb_location(:,:) .OR. bioturb_location(:,:) )
3533	xc_cryoturb(:,il,:) = (zf_soil(il)-zf_soil(il-1)) / time_step
3534	xd_cryoturb(:,il,:) = diff_k(:,il,:) / (zi_soil(il+1)-zi_soil(il))
3535	endwhere
3536	ENDDO
3537
3538
3539	DO iele = 1,nelements
3540	WHERE ( cryoturb_location(:,:) .OR. bioturb_location(:,:) )
3541	xc_cryoturb(:,ngrnd,:) = (zf_soil(ngrnd)-zf_soil(ngrnd-1)) / time_step
3542
3543	!bottom
3544	xe_a(:,:) = xc_cryoturb(:,ngrnd,:)+xd_cryoturb(:,ngrnd-1,:)
3545	xe_s(:,:) = xc_cryoturb(:,ngrnd,:)+xd_cryoturb(:,ngrnd-1,:)
3546	xe_p(:,:) = xc_cryoturb(:,ngrnd,:)+xd_cryoturb(:,ngrnd-1,:)
3547	alpha_a(:,ngrnd-1,:,iele) = xd_cryoturb(:,ngrnd-1,:) / xe_a(:,:)
3548	alpha_s(:,ngrnd-1,:,iele) = xd_cryoturb(:,ngrnd-1,:) / xe_s(:,:)
3549	alpha_p(:,ngrnd-1,:,iele) = xd_cryoturb(:,ngrnd-1,:) / xe_p(:,:)
3550	beta_a(:,ngrnd-1,:,iele) = xc_cryoturb(:,ngrnd,:)*deepSOM_a(:,ngrnd,:,iele) / xe_a(:,:)
3551	beta_s(:,ngrnd-1,:,iele) = xc_cryoturb(:,ngrnd,:)*deepSOM_s(:,ngrnd,:,iele) / xe_s(:,:)
3552	beta_p(:,ngrnd-1,:,iele) = xc_cryoturb(:,ngrnd,:)*deepSOM_p(:,ngrnd,:,iele) / xe_p(:,:)
3553	END WHERE
3554
3555	!other levels
3556	DO il = ngrnd-2,1,-1
3557	WHERE ( cryoturb_location(:,:) .OR. bioturb_location(:,:) )
3558	xe_a(:,:) = xc_cryoturb(:,il+1,:) + (1.-alpha_a(:,il+1,:,iele))*xd_cryoturb(:,il+1,:) + xd_cryoturb(:,il,:)
3559	xe_s(:,:) = xc_cryoturb(:,il+1,:) + (1.-alpha_s(:,il+1,:,iele))*xd_cryoturb(:,il+1,:) + xd_cryoturb(:,il,:)
3560	xe_p(:,:) = xc_cryoturb(:,il+1,:) + (1.-alpha_p(:,il+1,:,iele))*xd_cryoturb(:,il+1,:) + xd_cryoturb(:,il,:)
3561	alpha_a(:,il,:,iele) = xd_cryoturb(:,il,:) / xe_a(:,:)
3562	alpha_s(:,il,:,iele) = xd_cryoturb(:,il,:) / xe_s(:,:)
3563	alpha_p(:,il,:,iele) = xd_cryoturb(:,il,:) / xe_p(:,:)
3564	beta_a(:,il,:,iele) = (xc_cryoturb(:,il+1,:)*deepSOM_a(:,il+1,:,iele) + &
3565	xd_cryoturb(:,il+1,:)*beta_a(:,il+1,:,iele)) / xe_a(:,:)
3566	beta_s(:,il,:,iele) = (xc_cryoturb(:,il+1,:)*deepSOM_s(:,il+1,:,iele) + &
3567	xd_cryoturb(:,il+1,:)*beta_s(:,il+1,:,iele)) / xe_s(:,:)
3568	beta_p(:,il,:,iele) = (xc_cryoturb(:,il+1,:)*deepSOM_p(:,il+1,:,iele) + &
3569	xd_cryoturb(:,il+1,:)*beta_p(:,il+1,:,iele)) / xe_p(:,:)
3570	END WHERE
3571	ENDDO
3572	ENDDO !End lop over nelements
3573	ELSE
3574	!
3575	! do not know this action
3576	!
3577	CALL ipslerr_p(3, 'cryoturbate', 'DO NOT KNOW WHAT TO DO:', TRIM(action), '')
3578	!
3579	ENDIF
3580
3581	! keep last action in mind
3582	!
3583	last_action = action
3584
3585	END SUBROUTINE cryoturbate
3586
3587	!!
3588	!================================================================================================================================
3589	!! SUBROUTINE : permafrost_decomp
3590	!!
3591	!>\BRIEF This routine calculates soil organic matter decomposition
3592	!! DESCRIPTION :
3593	!!
3594	!! RECENT CHANGE(S) : None
3595	!!
3596	!! MAIN OUTPUT VARIABLE(S) :
3597	!!
3598	!! REFERENCE(S) : None
3599	!!
3600	!! FLOWCHART11 : None
3601	!! \n
3602	!_
3603	!================================================================================================================================
3604
3605	SUBROUTINE permafrost_decomp (kjpindex, time_step, tprof, airvol_soil, &
3606	oxlim, tau_CH4troph, ok_methane, fbactratio, O2m, &
3607	totporO2_soil, totporCH4_soil, hslong, clay, &
3608	deepSOM_a, deepSOM_s, deepSOM_p, &
3609	deltaCH4g, deltaCH4, deltaSOM1_a, deltaSOM1_s, deltaSOM1_p, deltaSOM2, &
3610	deltaSOM3, O2_soil, CH4_soil, fbact_out, CN_target, MG_useallCpools)
3611
3612	!! 0. Variable and parameter declaration
3613
3614	!! 0.1 Input variables
3615
3616	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
3617	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
3618	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: tprof !! deep temperature profile
3619	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: airvol_soil
3620	LOGICAL, INTENT(in) :: oxlim !! O2 limitation taken into account
3621	REAL(r_std), INTENT(in) :: tau_CH4troph !! time constant of methanetrophy (s)
3622	LOGICAL, INTENT(in) :: ok_methane !! Is Methanogenesis and -trophy taken into account?
3623	REAL(r_std), INTENT(in) :: fbactratio !! time constant of methanogenesis (ratio to that of oxic)
3624	REAL(r_std), INTENT(in) :: O2m !! oxygen concentration [g/m3] below which there is anoxy
3625	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: totporO2_soil !! total O2 porosity (Tans, 1998)
3626	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: totporCH4_soil !! total CH4 porosity (Tans, 1998)
3627	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: hslong !! deep soil humidity
3628	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: clay !! clay content
3629	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(in) :: fbact_out
3630	REAL(r_std), DIMENSION(kjpindex,nvm,ncarb), INTENT(in) :: CN_target !! C to N ratio of SOM flux from one pool to another (gN m-2 dt-1)
3631	LOGICAL, INTENT(in) :: MG_useallCpools !! Do we allow all three C pools to feed methanogenesis?
3632	!! 0.2 Output variables
3633
3634	!! 0.3 Modified variables
3635
3636	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deepSOM_a !! soil organic matter (g/m**3) active
3637	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deepSOM_s !! soil organic matter (g/m**3) slow
3638	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deepSOM_p !! soil organic matter (g/m**3) passive
3639	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: deltaCH4
3640	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: deltaCH4g
3641	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deltaSOM1_a
3642	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deltaSOM1_s
3643	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deltaSOM1_p
3644	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deltaSOM2
3645	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(inout) :: deltaSOM3
3646	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: O2_soil !! oxygen (g O2/m**3 air)
3647	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm), INTENT(inout) :: CH4_soil !! methane (g CH4/m**3 air)
3648
3649	!! 0.4 Local variables
3650
3651	INTEGER(i_std) :: ier
3652	REAL(r_std), DIMENSION(ncarb,ncarb,nelements) :: somflux !! fluxes between soil orgnanic matter reservoirs
3653	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm) :: nadd_soil !! number of moles created / m**3 of air
3654	REAL(r_std) :: fbact_a,fbact_s, fbact_p,temp
3655	REAL(r_std) :: fbactCH4_a, fbactCH4_s, fbactCH4_p
3656	REAL(r_std), DIMENSION(nelements) :: dC
3657	REAL(r_std) :: dCm
3658	REAL(r_std) :: dCH4,dCH4m,dO2
3659	INTEGER(i_std) :: il, ip, iv, iele
3660	LOGICAL, SAVE :: firstcall = .TRUE. !! first call?
3661	!$OMP THREADPRIVATE(firstcall)
3662
3663
3664	IF (firstcall) THEN
3665
3666	ALLOCATE (fc(kjpindex,3,3,nvm),stat=ier)
3667	IF (ier /= 0) CALL ipslerr_p(3,'permafrost_decomp', 'Pb in alloc for fc','','')
3668
3669	ALLOCATE (fr(kjpindex,3,nvm),stat=ier)
3670	IF (ier /= 0) CALL ipslerr_p(3,'permafrost_decomp', 'Pb in alloc for fr','','')
3671
3672	!
3673	! calculate soil organic matter flux fractions
3674	!
3675	DO iv =1,nvm
3676	fc(:,iactive,iactive,iv) = 0.0_r_std
3677	fc(:,iactive,ipassive,iv) = 0.004_r_std
3678	fc(:,iactive,islow,iv) = 1._r_std - (.85-.68*clay(:)) - fc(:,iactive,ipassive,iv)
3679	!
3680	fc(:,islow,islow,iv) = .0_r_std
3681	fc(:,islow,iactive,iv) = .42_r_std
3682	fc(:,islow,ipassive,iv) = .03_r_std
3683	!
3684	fc(:,ipassive,ipassive,iv) = .0_r_std
3685	fc(:,ipassive,iactive,iv) = .45_r_std
3686	fc(:,ipassive,islow,iv) = .0_r_std
3687	!
3688	fr(:,:,iv) = 1._r_std-fc(:,:,iactive,iv)-fc(:,:,islow,iv)-fc(:,:,ipassive,iv)
3689	firstcall = .FALSE.
3690	END DO
3691	IF (printlev>=3) THEN
3692	DO ip = 1,kjpindex
3693	WRITE(,) 'cdk: permafrost_decomp: i, fraction respired gridcell(i) :', ip, fr(ip,:,1)
3694	END DO
3695	ENDIF
3696	ENDIF
3697
3698	!
3699	! calculate som consumption
3700	!
3701	nadd_soil(:,:,:) = zero
3702	somflux(:,:,:) = zero
3703
3704	deltaSOM1_a(:,:,:,:) = zero
3705	deltaSOM1_s(:,:,:,:) = zero
3706	deltaSOM1_p(:,:,:,:) = zero
3707	deltaCH4(:,:,:) = zero
3708	deltaCH4g(:,:,:) = zero
3709	deltaSOM2(:,:,:,:) = zero
3710	deltaSOM3(:,:,:,:) = zero
3711
3712	DO ip = 1, kjpindex
3713	!
3714	DO iv = 1, nvm
3715	!
3716	IF ( veget_mask_2d(ip,iv) ) THEN
3717	!
3718	DO il = 1, ngrnd
3719	!
3720	! 1 function that gives soil organic matter residence time as a function of
3721	! soil temperature (in seconds)
3722	!
3723	temp = tprof(ip,il,iv) - ZeroCelsius
3724	fbact_a = fbact_out(ip,il,iv)
3725	fbact_a = MAX(fbact_a,time_step)
3726	!
3727	IF ( fbact_a/HUGE(1.) .GT. .1 ) THEN
3728	fbact_s = fbact_a
3729	fbact_p = fbact_a
3730	ELSE
3731	fbact_s = fbact_a * fslow
3732	fbact_p = fbact_a * fpassive
3733	ENDIF
3734	!
3735	! methanogenesis: first guess, 10 times (fbactratio) slower than oxic
3736	! decomposition
3737	IF ( fbact_a/HUGE(1.) .GT. .1 ) THEN
3738	fbactCH4_a = fbact_a
3739	fbactCH4_s = fbact_s
3740	fbactCH4_p = fbact_p
3741	ELSE
3742	fbactCH4_a = fbact_a * fbactratio
3743	IF ( MG_useallCpools ) THEN
3744	fbactCH4_s = fbact_s * fbactratio
3745	fbactCH4_p = fbact_p * fbactratio
3746	ELSE
3747	fbactCH4_s = HUGE(1.0)
3748	fbactCH4_p = HUGE(1.0)
3749	ENDIF
3750	ENDIF
3751	!
3752	! 2 oxic decomposition: carbon and oxygen consumption
3753	!
3754	! 2.1 active
3755	!
3756	DO iele = 1, nelements
3757
3758	IF (oxlim) THEN
3759	dCm = O2_soil(ip,il,iv)airvol_soil(ip,il,iv)wC/wO2
3760	dC(iele) = MIN(deepSOM_a(ip,il,iv,iele) * time_step/fbact_a,dCm)
3761	ELSE
3762	dC(iele) = deepSOM_a(ip,il,iv,iele) * time_step/fbact_a
3763	ENDIF
3764
3765	! pour actif
3766	dC(iele) = dC(iele) * ( un - som_turn_iactive_clay_frac * clay(ip) )
3767
3768	! flux vers les autres reservoirs
3769	IF (iele .EQ. icarbon) THEN
3770	somflux(iactive,ipassive,iele) = fc(ip,iactive,ipassive,iv) * dC(iele)
3771	somflux(iactive,islow,iele) = fc(ip,iactive,islow,iv) * dC(iele)
3772	ELSEIF (iele .EQ. initrogen) THEN
3773	somflux(iactive,ipassive,iele) = fc(ip,iactive,ipassive,iv) * dC(icarbon) / CN_target(ip,iv, ipassive)
3774	somflux(iactive,islow,iele) = fc(ip,iactive,islow,iv) * dC(icarbon) / CN_target(ip,iv, islow)
3775	ELSE
3776	WRITE (numout,*) ' deltaSOM. We we have an elements which is neither carbon nor nitrogen. We stop.'
3777	STOP 'stomate_soil_carbon_discretization_deep_somcycle'
3778	ENDIF
3779	!
3780	deepSOM_a(ip,il,iv,iele) = deepSOM_a(ip,il,iv,iele) - dC(iele)
3781	dO2 = wO2/wC * dC(icarbon)*fr(ip,iactive,iv) / totporO2_soil(ip,il,iv)
3782	O2_soil(ip,il,iv) = MAX( O2_soil(ip,il,iv) - dO2, zero)
3783	! keep delta C * fr in memory (generates energy)
3784	IF (iele .EQ. icarbon) THEN
3785	deltaSOM1_a(ip,il,iv,iele) = dC(iele)*fr(ip,iactive,iv) !!this line!!!
3786	ELSEIF (iele .EQ. initrogen) THEN
3787	deltaSOM1_a(ip,il,iv,iele) = dC(iele) - (somflux(iactive,ipassive,iele)+somflux(iactive,islow,iele))
3788	ELSE
3789	WRITE (numout,*) ' deltaSOM. We we have an elements which is neither carbon nor nitrogen. We stop.'
3790	STOP 'stomate_soil_carbon_discretization_deep_somcycle'
3791	ENDIF
3792	!
3793	! 2.2 slow
3794	!
3795	IF (oxlim) THEN
3796	dCm = O2_soil(ip,il,iv)airvol_soil(ip,il,iv)wC/wO2
3797	dC(iele) = MIN(deepSOM_s(ip,il,iv,iele) * time_step/fbact_s,dCm)
3798	ELSE
3799	dC(iele) = deepSOM_s(ip,il,iv,iele) * time_step/fbact_s
3800	ENDIF
3801	! flux vers les autres reservoirs
3802	IF (iele .EQ. icarbon) THEN
3803	somflux(islow,iactive,iele) = fc(ip,islow,iactive,iv) * dC(iele)
3804	somflux(islow,ipassive,iele) = fc(ip,islow,ipassive,iv) * dC(iele)
3805	ELSEIF (iele .EQ. initrogen) THEN
3806	somflux(islow,iactive,iele) = fc(ip,islow,iactive,iv) * dC(icarbon) / CN_target(ip,iv, iactive)
3807	somflux(islow,ipassive,iele) = fc(ip,islow,ipassive,iv) * dC(icarbon) / CN_target(ip,iv, ipassive)
3808	ELSE
3809	WRITE (numout,*) ' deltaSOM. We we have an elements which is neither carbon nor nitrogen. We stop.'
3810	STOP 'stomate_soil_carbon_discretization_deep_somcycle'
3811	ENDIF
3812	!
3813	deepSOM_s(ip,il,iv,iele) = deepSOM_s(ip,il,iv,iele) - dC(iele)
3814	dO2 = wO2/wC * dC(iele)*fr(ip,islow,iv) / totporO2_soil(ip,il,iv)
3815	O2_soil(ip,il,iv) = MAX( O2_soil(ip,il,iv) - dO2, zero)
3816	! keep delta C * fr in memory (generates energy)
3817	IF (iele .EQ. icarbon) THEN
3818	deltaSOM1_s(ip,il,iv,iele) = dC(iele)*fr(ip,islow,iv) !!this line!!!
3819	ELSEIF (iele .EQ. initrogen) THEN
3820	deltaSOM1_s(ip,il,iv,iele) = dC(iele) - (somflux(islow,iactive,iele)+somflux(islow,ipassive,iele))
3821	ELSE
3822	WRITE (numout,*) ' deltaSOM. We we have an elements which is neither carbon nor nitrogen. We stop.'
3823	STOP 'stomate_soil_carbon_discretization_deep_somcycle'
3824	ENDIF
3825	!
3826	! 2.3 passive
3827	!
3828	IF (oxlim) THEN
3829	dCm = O2_soil(ip,il,iv)airvol_soil(ip,il,iv)wC/wO2
3830	dC(iele) = MIN(deepSOM_p(ip,il,iv,iele) * time_step/fbact_p,dCm)
3831	ELSE
3832	dC(iele) = deepSOM_p(ip,il,iv,iele) * time_step/fbact_p
3833	ENDIF
3834	! flux vers les autres reservoirs
3835	IF (iele .EQ. icarbon) THEN
3836	somflux(ipassive,iactive,iele) = fc(ip,ipassive,iactive,iv) * dC(iele)
3837	somflux(ipassive,islow,iele) = fc(ip,ipassive,islow,iv) * dC(iele)
3838	ELSEIF (iele .EQ. initrogen) THEN
3839	somflux(ipassive,iactive,iele) = fc(ip,ipassive,iactive,iv) * dC(icarbon) / CN_target(ip,iv, iactive)
3840	somflux(ipassive,islow,iele) = fc(ip,ipassive,islow,iv) * dC(icarbon) / CN_target(ip,iv, islow)
3841	ELSE
3842	WRITE (numout,*) ' deltaSOM. We we have an elements which is neither carbon nor nitrogen. We stop.'
3843	STOP 'stomate_soil_carbon_discretization_deep_somcycle'
3844	ENDIF
3845	!
3846	deepSOM_p(ip,il,iv,iele) = deepSOM_p(ip,il,iv,iele) - dC(iele)
3847	dO2 = wO2/wC * dC(iele)*fr(ip,ipassive,iv) / totporO2_soil(ip,il,iv)
3848	O2_soil(ip,il,iv) = MAX( O2_soil(ip,il,iv) - dO2, zero)
3849	! keep delta C * fr in memory (generates energy)
3850	IF (iele .EQ. icarbon) THEN
3851	deltaSOM1_p(ip,il,iv,iele) = dC(iele)*fr(ip,ipassive,iv) !!this line!!!
3852	ELSEIF (iele .EQ. initrogen) THEN
3853	deltaSOM1_p(ip,il,iv,iele) = dC(iele)- (somflux(ipassive,iactive,iele)+somflux(ipassive,islow,iele))
3854	ELSE
3855	WRITE (numout,*) ' deltaSOM. We we have an elements which is neither carbon nor nitrogen. We stop.'
3856	STOP 'stomate_soil_carbon_discretization_deep_somcycle'
3857	ENDIF
3858	!
3859	!
3860	! 3 methanogenesis or methanotrophy
3861	!
3862	!
3863	IF ((ok_methane) .AND. (iele .EQ. icarbon)) THEN
3864	!
3865	!
3866	! 3.1 active pool methanogenesis
3867	dC(iele) = deepSOM_a(ip,il,iv,iele) * time_step / fbactCH4_a * EXP(-O2_soil(ip,il,iv)(1+hslong(ip,il,iv) &
3868	(BunsenO2-1.)) / O2m ) !DKtest: when commented, no ox lim for MG
3869	! pour actif
3870	dC(iele) = dC(iele) * ( 1. - .75 * clay(ip) )
3871	dCH4 = dC(iele)fr(ip,iactive,iv) wCH4/wC / totporCH4_soil(ip,il,iv)
3872	!
3873	!
3874	! flux vers les autres reservoirs
3875	somflux(iactive,ipassive,iele)=somflux(iactive,ipassive,iele)+fc(ip,iactive,ipassive,iv)*dC(iele)
3876	somflux(iactive,islow,iele)=somflux(iactive,islow,iele)+fc(ip,iactive,islow,iv)*dC(iele)
3877	!
3878	deepSOM_a(ip,il,iv,iele) = deepSOM_a(ip,il,iv,iele) - dC(iele)
3879	!
3880	deltaCH4g(ip,il,iv) = dCH4
3881	!
3882	CH4_soil(ip,il,iv) = CH4_soil(ip,il,iv) + dCH4
3883	! keep delta C*fr in memory (generates energy)
3884	deltaSOM2(ip,il,iv,iele) = dC(iele)*fr(ip,iactive,iv)
3885	!
3886	! how many moles of gas / m**3 of air did we generate?
3887	! (methanogenesis generates 1 molecule net if we take
3888	! B -> B' + CH4 )
3889	nadd_soil(ip,il,iv) = nadd_soil(ip,il,iv) + dCH4/wCH4
3890	!
3891	!
3892	IF ( MG_useallCpools ) THEN
3893	!
3894	! 3.2 slow pool methanogenesis cdk: adding this to allow other carbon pools to participate in MG
3895	dC(iele) = deepSOM_s(ip,il,iv,iele) * time_step / fbactCH4_s * EXP(-O2_soil(ip,il,iv)(1+hslong(ip,il,iv) &
3896	(BunsenO2-1.)) / O2m ) !DKtest: when commented, no ox lim for MG
3897	dCH4 = dC(iele)fr(ip,islow,iv) wCH4/wC / totporCH4_soil(ip,il,iv)
3898	!
3899	! flux vers les autres reservoirs
3900	somflux(islow,ipassive,iele)=somflux(islow,ipassive,iele)+fc(ip,islow,ipassive,iv)*dC(iele)
3901	somflux(islow,iactive,iele)=somflux(islow,iactive,iele)+fc(ip,islow,iactive,iv)*dC(iele)
3902	!
3903	deepSOM_s(ip,il,iv,iele) = deepSOM_s(ip,il,iv,iele) - dC(iele)
3904	!
3905	deltaCH4g(ip,il,iv) = deltaCH4g(ip,il,iv) + dCH4
3906	CH4_soil(ip,il,iv) = CH4_soil(ip,il,iv) + dCH4
3907	! keep delta C*fr in memory (generates energy)
3908	deltaSOM2(ip,il,iv,iele) = deltaSOM2(ip,il,iv,iele) + dC(iele)*fr(ip,islow,iv)
3909	!
3910	! how many moles of gas / m**3 of air did we generate?
3911	! (methanogenesis generates 1 molecule net if we take
3912	! B -> B' + CH4 )
3913	nadd_soil(ip,il,iv) = nadd_soil(ip,il,iv) + dCH4/wCH4
3914	!
3915	!
3916	!
3917	! 3.3 passive pool methanogenesis cdk: adding this to allow other carbon pools to participate in MG
3918	dC(iele) = deepSOM_p(ip,il,iv,iele) * time_step / fbactCH4_p * EXP(-O2_soil(ip,il,iv)(1+hslong(ip,il,iv) &
3919	(BunsenO2-1.)) / O2m ) !DKtest: when commented, no ox lim for MG
3920	dCH4 = dC(iele)fr(ip,ipassive,iv) wCH4/wC / totporCH4_soil(ip,il,iv)
3921	!
3922	! flux vers les autres reservoirs
3923	somflux(ipassive,islow,iele)=somflux(ipassive,islow,iele)+fc(ip,ipassive,islow,iv)*dC(iele)
3924	somflux(ipassive,iactive,iele)=somflux(ipassive,iactive,iele)+fc(ip,ipassive,iactive,iv)*dC(iele)
3925	!
3926	deepSOM_p(ip,il,iv,iele) = deepSOM_p(ip,il,iv,iele) - dC(iele)
3927	!
3928	deltaCH4g(ip,il,iv) = deltaCH4g(ip,il,iv) + dCH4
3929	CH4_soil(ip,il,iv) = CH4_soil(ip,il,iv) + dCH4
3930	! keep delta C*fr in memory (generates energy)
3931	deltaSOM2(ip,il,iv,iele) = deltaSOM2(ip,il,iv,iele) + dC(iele)*fr(ip,ipassive,iv)
3932	!
3933	! how many moles of gas / m**3 of air did we generate?
3934	! (methanogenesis generates 1 molecule net if we take
3935	! B -> B' + CH4 )
3936	nadd_soil(ip,il,iv) = nadd_soil(ip,il,iv) + dCH4/wCH4
3937	!
3938	!
3939	ENDIF
3940	!
3941	! trophy:
3942	! no temperature dependence except that T>0ï¿œï¿œC (Price et
3943	! al, GCB 2003; Koschorrek and Conrad, GBC 1993).
3944	! tau_CH4troph is such that we fall between values of
3945	! soil methane oxidation flux given by these authors.
3946	!
3947	IF ( temp .GE. zero ) THEN
3948	!
3949	dCH4m = O2_soil(ip,il,iv)/2. * wCH4/wO2 * totporO2_soil(ip,il,iv)/totporCH4_soil(ip,il,iv)
3950	! dCH4m = CH4_soil(ip,il,iv) !DKtest - no ox lim to trophy
3951	dCH4 = MIN( CH4_soil(ip,il,iv) * time_step/MAX(tau_CH4troph,time_step), dCH4m )
3952	CH4_soil(ip,il,iv) = CH4_soil(ip,il,iv) - dCH4
3953	dO2 = 2.dCH4 wO2/wCH4 * totporCH4_soil(ip,il,iv)/totporO2_soil(ip,il,iv)
3954	O2_soil(ip,il,iv) = MAX( O2_soil(ip,il,iv) - dO2, zero)
3955	! keep delta CH4 in memory (generates energy)
3956	deltaCH4(ip,il,iv) = dCH4
3957	! carbon (g/m3 soil) transformed to CO2
3958	deltaSOM3(ip,il,iv,icarbon)=dCH4/wCH4wCtotporCH4_soil(ip,il,iv)
3959	! how many moles of gas / m**3 of air did we generate?
3960	! (methanotrophy consumes 2 molecules net if we take
3961	! CH4 + 2 O2 -> CO2 + 2 H2O )
3962	nadd_soil(ip,il,iv) = nadd_soil(ip,il,iv)-2.*dCH4/wCH4
3963	!
3964	ENDIF
3965
3966	ENDIF
3967
3968	! 4 add fluxes between reservoirs
3969
3970	deepSOM_a(ip,il,iv,iele)=deepSOM_a(ip,il,iv,iele)+somflux(islow,iactive,iele)+somflux(ipassive,iactive,iele)
3971	deepSOM_s(ip,il,iv,iele)=deepSOM_s(ip,il,iv,iele)+somflux(iactive,islow,iele)+somflux(ipassive,islow,iele)
3972	deepSOM_p(ip,il,iv,iele)=deepSOM_p(ip,il,iv,iele)+somflux(iactive,ipassive,iele)+somflux(islow,ipassive,iele)
3973
3974	ENDDO
3975
3976	ENDDO ! End loop over nelements
3977
3978	ENDIF
3979
3980	ENDDO
3981
3982	ENDDO
3983	END SUBROUTINE permafrost_decomp
3984
3985
3986	!!
3987	!================================================================================================================================
3988	!! SUBROUTINE : calc_vert_int_som
3989	!!
3990	!>\BRIEF This routine calculates carbon decomposition
3991	!!
3992	!! DESCRIPTION :
3993	!!
3994	!! RECENT CHANGE(S) : None
3995	!!
3996	!! MAIN OUTPUT VARIABLE(S) :
3997	!!
3998	!! REFERENCE(S) : None
3999	!!
4000	!! FLOWCHART11 : None
4001	!! \n
4002	!_
4003	!================================================================================================================================
4004
4005	SUBROUTINE calc_vert_int_som(kjpindex, deepSOM_a, deepSOM_s, deepSOM_p, som, som_surf, zf_soil)
4006
4007	!! 0. Variable and parameter declaration
4008
4009	!! 0.1 Input variables
4010
4011	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
4012	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(in) :: deepSOM_a !! active pool deepsom
4013	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(in) :: deepSOM_s !! slow pool deepsom
4014	REAL(r_std), DIMENSION(kjpindex,ngrnd,nvm,nelements), INTENT(in) :: deepSOM_p !! passive pool deepsom
4015	REAL(r_std), DIMENSION(0:ngrnd), INTENT(in) :: zf_soil !! depths at full levels
4016
4017	!! 0.2 Output variables
4018
4019	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements), INTENT (out) :: som !! vertically-integrated som pool: active, slow, or passive, (gC/(m**2 of ground))
4020	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm,nelements), INTENT (out):: som_surf!! vertically-integrated som pool to 1 meter: active, slow, or passive,(gC/(m**2 of ground))
4021
4022	!! 0.3 Modified variables
4023
4024	!! 0.4 Local variables
4025	INTEGER(i_std) :: il,iele, ivm
4026	real(r_std), parameter :: maxdepth=2.!! depth to which we intergrate the carbon for som_surf calculation
4027
4028	som(:,:,:,:) = zero
4029	DO il = 1, ngrnd
4030	DO iele = 1, nelements
4031	WHERE ( veget_mask_2d(:,:) )
4032	som(:,iactive,:,iele) = som(:,iactive,:,iele) + deepSOM_a(:,il,:,iele)*(zf_soil(il)-zf_soil(il-1))
4033	som(:,islow,:,iele) = som(:,islow,:,iele) + deepSOM_s(:,il,:,iele)*(zf_soil(il)-zf_soil(il-1))
4034	som(:,ipassive,:,iele) = som(:,ipassive,:,iele) + deepSOM_p(:,il,:,iele)*(zf_soil(il)-zf_soil(il-1))
4035	END WHERE
4036	ENDDO
4037	ENDDO
4038
4039	som_surf(:,:,:,:) = zero
4040	DO il = 1, ngrnd
4041	DO iele = 1, nelements
4042	IF (zf_soil(il-1) .lt. maxdepth ) THEN
4043	WHERE ( veget_mask_2d(:,:) )
4044	som_surf(:,iactive,:,iele) = som_surf(:,iactive,:,iele) + &
4045	deepSOM_a(:,il,:,iele)*(min(maxdepth,zf_soil(il))-zf_soil(il-1))
4046	som_surf(:,islow,:,iele) = som_surf(:,islow,:,iele) + &
4047	deepSOM_s(:,il,:,iele)*(min(maxdepth,zf_soil(il))-zf_soil(il-1))
4048	som_surf(:,ipassive,:,iele) = som_surf(:,ipassive,:,iele) + &
4049	deepSOM_p(:,il,:,iele)*(min(maxdepth,zf_soil(il))-zf_soil(il-1))
4050	END WHERE
4051	ENDIF
4052	ENDDO
4053	ENDDO
4054
4055	END SUBROUTINE calc_vert_int_som
4056
4057	END MODULE stomate_soil_carbon_discretization

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

Download in other formats: