Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_stomate/stomate_permafrost_soilcarbon.f90

Last change on this file was 5488, checked in by chunjing.qiu, 6 years ago
C balance checked
File size: 196.5 KB

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

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

Download in other formats: