Context Navigation

stomate_permafrost_soilcarbon.f90 @ 6737

Last change on this file since 6737 was 5098, checked in by albert.jornet, 6 years ago

Fix: CHARACTER requires to be TRIMMED before comparing. It failed even when there was the same word. Extra spaces made the comparision failed.
Fix: include openmp threadprivate statements for stomate_permafrost_soilcarbon.
Fix: missing orch_reduce_sum_omp_rgen_scal and orch_reduce_sum_omp_igen_scal subroutine for scalar values.
Clean: change action variable name to cur_action. It's a Fortran reserved word
Clean: use ipslerr_p instead of STOP to properly signal an issue

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

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

Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/stomate_permafrost_soilcarbon.f90 @ 6737

Download in other formats: