Context Navigation

stomate_permafrost_soilcarbon.f90 @ 6940

Last change on this file since 6940 was 6940, checked in by jinfeng.chang, 4 years ago
add missing files for ORCHIDEE-GMv3.2
File size: 188.9 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, febul, kjpindex*nvm, horipft_index)
1098	CALL histwrite_p (hist_id_stomate, 'flupmt', itime, flupmt, kjpindex*nvm, horipft_index)
1099	CALL histwrite_p (hist_id_stomate, 'alt', itime, alt, kjpindex*nvm, horipft_index)
1100	CALL histwrite_p (hist_id_stomate, 'altmax', itime, altmax, kjpindex*nvm, horipft_index)
1101	CALL histwrite_p (hist_id_stomate, 'sfluxCH4_deep', itime, sfluxCH4_deep, kjpindex, index)
1102	CALL histwrite_p (hist_id_stomate, 'sfluxCO2_deep', itime, sfluxCO2_deep, kjpindex, index)
1103	CALL histwrite_p (hist_id_stomate, 'pb', itime, pb, kjpindex, index)
1104	call histwrite_p (hist_id_stomate, 'deepC_a_pftmean', itime, deepC_pftmean(:,:,iactive), kjpindex*ndeep, horideep_index)
1105	call histwrite_p (hist_id_stomate, 'deepC_s_pftmean', itime, deepC_pftmean(:,:,islow), kjpindex*ndeep, horideep_index)
1106	call histwrite_p (hist_id_stomate, 'deepC_p_pftmean', itime, deepC_pftmean(:,:,ipassive), kjpindex*ndeep, horideep_index)
1107
1108	DO jv = 1, nvm
1109	IF (permafrost_veg_exists(jv)) THEN !don't bother to write if there are pfts that don't exist in our domain
1110	WRITE(part_str,'(I2)') jv
1111	IF (jv < 10) part_str(1:1) = '0'
1112	IF (writehist_deepC) THEN
1113	CALL histwrite_p (hist_id_stomate, 'deepC_a_'//part_str(1:LEN_TRIM(part_str)), &
1114	itime, deepC_a(:,:,jv), kjpindex*ndeep, horideep_index)
1115	CALL histwrite_p (hist_id_stomate, 'deepC_s_'//part_str(1:LEN_TRIM(part_str)), &
1116	itime, deepC_s(:,:,jv), kjpindex*ndeep, horideep_index)
1117	CALL histwrite_p (hist_id_stomate, 'deepC_p_'//part_str(1:LEN_TRIM(part_str)), &
1118	itime, deepC_p(:,:,jv), kjpindex*ndeep, horideep_index)
1119	ENDIF
1120	IF (writehist_soilgases) THEN
1121	CALL histwrite_p (hist_id_stomate, 'O2_soil_'//part_str(1:LEN_TRIM(part_str)), &
1122	itime, O2_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1123	CALL histwrite_p (hist_id_stomate, 'CH4_soil_'//part_str(1:LEN_TRIM(part_str)), &
1124	itime, CH4_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1125	CALL histwrite_p (hist_id_stomate, 'O2_snow_'//part_str(1:LEN_TRIM(part_str)), &
1126	itime, O2_snow(:,:,jv), kjpindex*nsnow, horisnow_index)
1127	CALL histwrite_p (hist_id_stomate, 'CH4_snow_'//part_str(1:LEN_TRIM(part_str)), &
1128	itime, CH4_snow(:,:,jv), kjpindex*nsnow, horisnow_index)
1129	ENDIF
1130	IF (writehist_deltaC) THEN
1131	CALL histwrite_p (hist_id_stomate, 'deltaCH4g_'//part_str(1:LEN_TRIM(part_str)), &
1132	itime, deltaCH4g(:,:,jv), kjpindex*ndeep, horideep_index)
1133	CALL histwrite_p (hist_id_stomate, 'deltaCH4_'//part_str(1:LEN_TRIM(part_str)), &
1134	itime, deltaCH4(:,:,jv), kjpindex*ndeep, horideep_index)
1135	CALL histwrite_p (hist_id_stomate, 'deltaC1_'//part_str(1:LEN_TRIM(part_str)), &
1136	itime, deltaC1_a(:,:,jv)+deltaC1_s(:,:,jv)+deltaC1_p(:,:,jv), kjpindex*ndeep, horideep_index)
1137	CALL histwrite_p (hist_id_stomate, 'deltaC2_'//part_str(1:LEN_TRIM(part_str)), &
1138	itime, deltaC2(:,:,jv), kjpindex*ndeep, horideep_index)
1139	CALL histwrite_p (hist_id_stomate, 'deltaC3_'//part_str(1:LEN_TRIM(part_str)), &
1140	itime, deltaC3(:,:,jv), kjpindex*ndeep, horideep_index)
1141	ENDIF
1142
1143	IF (writehist_zimovheat) THEN
1144	CALL histwrite_p (hist_id_stomate, 'heat_Zimov_'//part_str(1:LEN_TRIM(part_str)), &
1145	itime, heat_Zimov(:,:,jv), kjpindex*ndeep, horideep_index)
1146	ENDIF
1147
1148	IF (writehist_deltaC_litter) THEN
1149	CALL histwrite_p (hist_id_stomate, 'deltaC_litter_act_'//part_str(1:LEN_TRIM(part_str)), &
1150	itime, dc_litter_z(:,iactive,:,jv)/ time_step, kjpindex*ndeep, horideep_index)
1151	CALL histwrite_p (hist_id_stomate, 'deltaC_litter_slo_'//part_str(1:LEN_TRIM(part_str)), &
1152	itime, dc_litter_z(:,islow,:,jv)/ time_step, kjpindex*ndeep, horideep_index)
1153	ENDIF
1154	!------------------------------ further output for debugging/diagnosing
1155
1156	IF (writehist_gascoeff) THEN
1157	CALL histwrite_p (hist_id_stomate, 'totporO2_soil_'//part_str(1:LEN_TRIM(part_str)), &
1158	itime, totporO2_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1159	CALL histwrite_p (hist_id_stomate, 'diffO2_soil_'//part_str(1:LEN_TRIM(part_str)), &
1160	itime, diffO2_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1161	CALL histwrite_p (hist_id_stomate, 'alphaO2_soil_'//part_str(1:LEN_TRIM(part_str)), &
1162	itime, alphaO2_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1163	CALL histwrite_p (hist_id_stomate, 'betaO2_soil_'//part_str(1:LEN_TRIM(part_str)), &
1164	itime, betaO2_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1165
1166	CALL histwrite_p (hist_id_stomate, 'totporCH4_soil_'//part_str(1:LEN_TRIM(part_str)), &
1167	itime, totporCH4_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1168	CALL histwrite_p (hist_id_stomate, 'diffCH4_soil_'//part_str(1:LEN_TRIM(part_str)), &
1169	itime, diffCH4_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1170	CALL histwrite_p (hist_id_stomate, 'alphaCH4_soil_'//part_str(1:LEN_TRIM(part_str)), &
1171	itime, alphaCH4_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1172	CALL histwrite_p (hist_id_stomate, 'betaCH4_soil_'//part_str(1:LEN_TRIM(part_str)), &
1173	itime, betaCH4_soil(:,:,jv), kjpindex*ndeep, horideep_index)
1174	ENDIF
1175	END IF
1176	END DO
1177
1178	ENDIF
1179
1180	! XIOS history output
1181	IF ( .NOT. soilc_isspinup ) THEN
1182
1183	CALL xios_orchidee_send_field ('tsurf', tsurf)
1184	CALL xios_orchidee_send_field ('fluxCH4', sfluxCH4)
1185	CALL xios_orchidee_send_field ('febul', febul)
1186	CALL xios_orchidee_send_field ('flupmt', flupmt)
1187	CALL xios_orchidee_send_field ( 'alt', alt )
1188	CALL xios_orchidee_send_field ( 'altmax', altmax)
1189	! JC add for ISIMIP2b
1190	CALL xios_orchidee_send_field ( 'altmax_annual', altmax(:,1))
1191	! End JC add
1192	CALL xios_orchidee_send_field ( 'sfluxCH4_deep', sfluxCH4_deep)
1193	CALL xios_orchidee_send_field ( 'sfluxCO2_deep', sfluxCO2_deep)
1194	CALL xios_orchidee_send_field ( 'pb', pb)
1195	call xios_orchidee_send_field ( 'deepC_a_pftmean', deepC_pftmean(:,:,iactive))
1196	call xios_orchidee_send_field ( 'deepC_s_pftmean', deepC_pftmean(:,:,islow))
1197	call xios_orchidee_send_field ( 'deepC_p_pftmean', deepC_pftmean(:,:,ipassive))
1198
1199	IF (writehist_deepC) THEN
1200	CALL xios_orchidee_send_field ( 'deepC_a', deepC_a)
1201	CALL xios_orchidee_send_field ( 'deepC_s', deepC_s)
1202	CALL xios_orchidee_send_field ( 'deepC_p', deepC_p)
1203	ENDIF
1204
1205	IF (writehist_soilgases) THEN
1206	CALL xios_orchidee_send_field ( 'O2_soil', O2_soil)
1207	CALL xios_orchidee_send_field ( 'CH4_soil', CH4_soil)
1208	CALL xios_orchidee_send_field ('O2_snow', O2_snow)
1209	CALL xios_orchidee_send_field ( 'CH4_snow', CH4_snow)
1210	ENDIF
1211
1212	IF (writehist_deltaC) THEN
1213	CALL xios_orchidee_send_field ( 'deltaCH4g', deltaCH4g)
1214	CALL xios_orchidee_send_field ( 'deltaCH4', deltaCH4)
1215	CALL xios_orchidee_send_field ( 'deltaC1', deltaC1_a+deltaC1_s+deltaC1_p)
1216	CALL xios_orchidee_send_field ( 'deltaC2', deltaC2)
1217	CALL xios_orchidee_send_field ( 'deltaC3', deltaC3)
1218	ENDIF
1219
1220	IF (writehist_zimovheat) THEN
1221	CALL xios_orchidee_send_field ( 'heat_Zimov', heat_Zimov)
1222	ENDIF
1223
1224	IF (writehist_deltaC_litter) THEN
1225	CALL xios_orchidee_send_field ( 'deltaC_litter_act', dc_litter_z(:,iactive,:,:)/ time_step)
1226	CALL xios_orchidee_send_field ( 'deltaC_litter_slo', dc_litter_z(:,islow,:,:)/ time_step)
1227	ENDIF
1228
1229	IF (writehist_gascoeff) THEN
1230	CALL xios_orchidee_send_field ( 'totporO2_soil', totporO2_soil)
1231	CALL xios_orchidee_send_field ( 'diffO2_soil', diffO2_soil)
1232	CALL xios_orchidee_send_field ( 'alphaO2_soil', alphaO2_soil)
1233	CALL xios_orchidee_send_field ( 'betaO2_soil', betaO2_soil)
1234
1235	CALL xios_orchidee_send_field ( 'totporCH4_soil', totporCH4_soil)
1236	CALL xios_orchidee_send_field ( 'diffCH4_soil', diffCH4_soil)
1237	CALL xios_orchidee_send_field ('alphaCH4_soil', alphaCH4_soil)
1238	CALL xios_orchidee_send_field ( 'betaCH4_soil', betaCH4_soil)
1239	ENDIF
1240
1241	ENDIF
1242
1243	IF (printlev>=3) WRITE(,) 'cdk: leaving deep_carbcycle'
1244
1245	IF ( firstcall ) firstcall = .FALSE.
1246
1247
1248	END SUBROUTINE deep_carbcycle
1249
1250	!!
1251	!================================================================================================================================
1252	!! SUBROUTINE : altcalc
1253	!!
1254	!>\BRIEF This routine calculate active layer thickness
1255	!!
1256	!! DESCRIPTION :
1257	!!
1258	!! RECENT CHANGE(S) : None
1259	!!
1260	!! MAIN OUTPUT VARIABLE(S) : alt
1261	!!
1262	!! REFERENCE(S) : None
1263	!!
1264	!! FLOWCHART11 : None
1265	!! \n
1266	!_
1267	!================================================================================================================================
1268	SUBROUTINE altcalc (kjpindex,time_step,dayno,scnd, temp, zprof, alt, alt_ind, altmax, altmax_ind, &
1269	altmax_lastyear, altmax_ind_lastyear)
1270
1271	!! 0. Variable and parameter declaration
1272
1273	!! 0.1 Input variables
1274
1275	INTEGER(i_std), INTENT(in) :: kjpindex
1276	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
1277	INTEGER(i_std), INTENT(in) :: dayno !! number of the day in the current year
1278	REAL(r_std), INTENT(in) :: scnd !! model time & time step
1279	REAL(r_std), DIMENSION(kjpindex,ndeep, nvm), INTENT(in) :: temp !! soil temperature
1280	REAL(r_std), DIMENSION(ndeep), INTENT(in) :: zprof !! soil depths (m)
1281
1282	!! 0.2 Output variables
1283
1284	REAL(r_std), DIMENSION(kjpindex, nvm), INTENT(out) :: alt !! active layer thickness
1285	INTEGER, DIMENSION(kjpindex, nvm), INTENT(out) :: alt_ind !! active layer index
1286
1287	!! 0.3 Modified variables
1288
1289	REAL(r_std), DIMENSION(kjpindex, nvm),INTENT(inout) :: altmax_lastyear !! Maximum active-layer thickness
1290	REAL(r_std), DIMENSION(kjpindex, nvm),INTENT(inout) :: altmax !! Maximum active-layer thickness
1291	INTEGER(i_std), DIMENSION(kjpindex, nvm),INTENT(inout) :: altmax_ind !! Maximum over the year active-layer index
1292	INTEGER(i_std), DIMENSION(kjpindex, nvm),INTENT(inout) :: altmax_ind_lastyear !! Maximum over the year active-layer index
1293
1294	!! 0.4 Local variables
1295
1296	INTEGER :: ix,iz,il,iv !! grid indices
1297	LOGICAL, SAVE :: firstcall = .TRUE.
1298	!$OMP THREADPRIVATE(firstcall)
1299	INTEGER, save :: tcounter
1300	INTEGER(i_std), SAVE :: id, id2
1301	!$OMP THREADPRIVATE(id)
1302	!$OMP THREADPRIVATE(id2)
1303	LOGICAL, SAVE :: check = .FALSE.
1304	!$OMP THREADPRIVATE(check)
1305	LOGICAL, SAVE :: newaltcalc = .FALSE.
1306	!$OMP THREADPRIVATE(newaltcalc)
1307	LOGICAL, DIMENSION(kjpindex,nvm) :: inalt, bottomlevelthawed
1308	CHARACTER(LEN=16) :: buf
1309	INTEGER :: lev
1310
1311
1312	IF ( firstcall ) THEN
1313
1314	! calculate altmax_ind from altmax
1315	altmax_ind(:,:) = 0
1316	DO ix = 1, kjpindex
1317	DO iv = 1, nvm
1318	IF ( veget_mask_2d(ix,iv) ) THEN
1319	DO il=1,ndeep
1320	IF ( altmax(ix,iv) .GE. zprof(il) ) THEN
1321	altmax_ind(ix,iv) = altmax_ind(ix,iv) + 1
1322	END IF
1323	END DO
1324	END IF
1325	END DO
1326	END DO
1327	altmax_lastyear(:,:) = altmax(:,:)
1328	altmax_ind_lastyear(:,:) = altmax_ind(:,:)
1329	firstcall = .FALSE.
1330
1331	!Config Key = newaltcalc
1332	!Config Desc = calculate alt ?
1333	!Config Def = n
1334	!Config If = OK_PC
1335	!Config Help =
1336	!Config Unit = [flag]
1337	CALL getin_p('newaltcalc', newaltcalc)
1338
1339	ELSE
1340	! all other timesteps
1341	IF ( .NOT. newaltcalc ) THEN
1342	DO ix = 1, kjpindex
1343	DO iv = 1, nvm
1344	IF ( veget_mask_2d(ix,iv) ) THEN
1345	iz = 1
1346	DO WHILE( temp(ix,iz,iv) > ZeroCelsius .AND. iz < ndeep )
1347	iz = iz + 1
1348	END DO
1349	IF( iz == 1 ) THEN
1350	! it means that all is frozen
1351	alt(ix,iv) = zero
1352	ELSE
1353	alt(ix,iv) = zprof(iz-1)
1354	END IF
1355	alt_ind(ix,iv) = iz-1
1356	END IF
1357	END DO
1358	END DO
1359	ELSE
1360	! initialize for pfts that don't exist
1361	alt(:,:) = zprof(ndeep)
1362	bottomlevelthawed(:,:) = .FALSE.
1363	! start from bottom and work up instead
1364	WHERE (temp(:,ndeep,:) > ZeroCelsius )
1365	bottomlevelthawed(:,:) = .TRUE.
1366	alt(:,:) = zprof(ndeep)
1367	alt_ind(:,:) = ndeep
1368	END WHERE
1369	inalt(:,:) = .FALSE.
1370	DO iz = 1, ndeep - 1
1371	lev = ndeep - iz
1372	WHERE ( temp(:,lev,:) > ZeroCelsius .AND. .NOT. inalt(:,:) .AND. .NOT. bottomlevelthawed(:,:) )
1373	inalt(:,:) = .TRUE.
1374	alt(:,:) = zprof(lev)
1375	alt_ind(:,:) = lev
1376	ELSEWHERE ( temp(:,lev,:) <= ZeroCelsius .AND. inalt(:,:) .AND. .NOT. bottomlevelthawed(:,:) )
1377	inalt(:,:) = .FALSE.
1378	END WHERE
1379	END DO
1380	WHERE ( .NOT. inalt .AND. .NOT. bottomlevelthawed(:,:) )
1381	alt(:,:) = zero
1382	alt_ind(:,:) = 0
1383	END WHERE
1384	ENDIF
1385
1386	! debug
1387	IF ( check ) THEN
1388	IF (ANY(alt(:,:) .GT. zprof(ndeep))) THEN
1389	WRITE(,) 'error: alt greater than soil depth.'
1390	ENDIF
1391	ENDIF
1392
1393	! Maximum over the year active layer thickness
1394	WHERE ( ( alt(:,:) .GT. altmax(:,:) ) .AND. veget_mask_2d(:,:) )
1395	altmax(:,:) = alt(:,:)
1396	altmax_ind(:,:) = alt_ind(:,:)
1397	ENDWHERE
1398
1399	IF ( .NOT. soilc_isspinup ) THEN
1400	! do it on the second timestep, that way when we are writing restart files it is not done before that!
1401	! now we are doing daily permafrost calcs, so just run it on the second day.
1402	IF ( ( dayno .EQ. 2) ) THEN
1403	! Reinitialize ALT_max
1404	altmax_lastyear(:,:) = altmax(:,:)
1405	altmax_ind_lastyear(:,:) = altmax_ind(:,:)
1406	altmax(:,:) = alt(:,:)
1407	altmax_ind(:,:) = alt_ind(:,:)
1408	END IF
1409	ELSE
1410
1411	! for spinup, best to set altmax_lastyear to altmax, and not boter to reset since every year is the same,
1412	! and if you try to do so, it doesn't work properly -- 06 may 2010
1413	altmax_lastyear(:,:) = altmax(:,:)
1414	altmax_ind_lastyear(:,:) = altmax_ind(:,:)
1415	END IF
1416	END IF
1417
1418	IF (printlev>=3) WRITE(,) 'leaving altcalc'
1419	END SUBROUTINE altcalc
1420
1421	!!
1422	!================================================================================================================================
1423	!! SUBROUTINE : soil_gasdiff_main
1424	!!
1425	!>\BRIEF This routine calculate oxygen and methane in the snow/soil medium
1426	!!
1427	!! DESCRIPTION :
1428	!!
1429	!! RECENT CHANGE(S) : None
1430	!!
1431	!! MAIN OUTPUT VARIABLE(S) :
1432	!!
1433	!! REFERENCE(S) : None
1434	!!
1435	!! FLOWCHART11 : None
1436	!! \n
1437	!_
1438	!================================================================================================================================
1439	SUBROUTINE soil_gasdiff_main( kjpindex,time_step,index,cur_action, &
1440	psol,tsurf,tprof,diffO2_snow,diffCH4_snow, &
1441	totporO2_snow,totporCH4_snow,O2_snow,CH4_snow,diffO2_soil,diffCH4_soil, &
1442	totporO2_soil,totporCH4_soil,O2_soil,CH4_soil, zi_snow, zf_snow)
1443
1444	!! 0. Variable and parameter declaration
1445
1446	!! 0.1 Input variables
1447
1448	INTEGER(i_std), INTENT(in) :: kjpindex !! number of grid points
1449	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
1450	CHARACTER(LEN=*), INTENT(in) :: cur_action !! what to do
1451	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: psol !! surface pressure (Pa)
1452	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: tsurf !! Surface temperature (K)
1453	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: tprof !! Soil temperature (K)
1454	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: diffO2_snow !! oxygen diffusivity (m**2/s)
1455	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: diffCH4_snow !! methane diffusivity (m**2/s)
1456	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: totporO2_snow !! total O2 porosity (Tans, 1998)
1457	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: totporCH4_snow !! total CH4 porosity (Tans, 1998)
1458	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: diffO2_soil !! oxygen diffusivity (m**2/s)
1459	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: diffCH4_soil !! methane diffusivity (m**2/s)
1460	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: totporO2_soil !! total O2 porosity (Tans, 1998)
1461	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: totporCH4_soil !! total CH4 porosity (Tans, 1998)
1462	INTEGER(i_std),DIMENSION(kjpindex),INTENT(in) :: index !! Indeces of permafrost points on the map
1463
1464	!! 0.2 Output variables
1465
1466	!! 0.3 Modified variables
1467
1468	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: O2_snow !! oxygen (g O2/m**3 air)
1469	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: CH4_snow !! methane (g CH4/m**3 air)
1470	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: O2_soil !! oxygen (g O2/m**3 air)
1471	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: CH4_soil !! methane (g CH4/m**3 air)
1472	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm), intent(inout):: zf_snow !! depths of full levels (m)
1473	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), intent(inout) :: zi_snow !! depths of intermediate levels (m)
1474
1475	!! 0.4 local variables
1476
1477	CHARACTER(LEN=50), SAVE :: last_action = 'not called'
1478	!$OMP THREADPRIVATE(last_action)
1479
1480
1481	! 1. ensure that we do not repeat actions
1482	!
1483	IF ( TRIM(cur_action) .EQ. TRIM(last_action) ) THEN
1484	!
1485	CALL ipslerr_p(3, 'soil_gasdiff_main','CANNOT TAKE THE SAME ACTION TWICE: ',cur_action, last_action)
1486	!
1487	ENDIF
1488	!
1489	! 2. decide what to do
1490	!
1491	IF ( TRIM(cur_action) .EQ. 'initialize' ) THEN
1492	!
1493	! 2.1 initialize
1494	!
1495	IF ( TRIM(last_action) .NE. 'not called' ) THEN
1496	!
1497	CALL ipslerr_p(3, 'soil_gasdiff_main','SOIL MODEL CANNOT BE INITIALIZED TWICE.', '', '')
1498	!
1499	ENDIF
1500	!
1501	CALL soil_gasdiff_alloc( kjpindex )
1502	!
1503	ELSEIF ( TRIM(cur_action) .EQ. 'diffuse' ) THEN
1504	!
1505	! 2.2 calculate soil temperatures
1506	!
1507	CALL soil_gasdiff_diff( kjpindex,time_step,index,psol,tsurf, O2_snow, CH4_snow, O2_soil, CH4_soil)
1508	!
1509	ELSEIF ( TRIM(cur_action) .EQ. 'coefficients' ) THEN
1510	!
1511	! 2.3 calculate coefficients (heat flux and apparent surface heat capacity)
1512	!
1513	CALL soil_gasdiff_coeff( kjpindex,time_step,tprof,O2_snow,CH4_snow, &
1514	diffO2_snow,diffCH4_snow,totporO2_snow,totporCH4_snow,O2_soil,CH4_soil, &
1515	diffO2_soil,diffCH4_soil,totporO2_soil,totporCH4_soil, zi_snow, zf_snow)
1516	!
1517	ELSE
1518	!
1519	! 2.4 do not know this action
1520	!
1521	CALL ipslerr_p(3,'soil_gasdiff_main', 'This action does not exists and must be implemented', &
1522	'or its wrong:',cur_action)
1523	!
1524	ENDIF
1525	!
1526	! 2.5 keep last action in mind
1527	!
1528	last_action = TRIM(cur_action)
1529
1530	IF (printlev>=3) WRITE(,) 'leaving soil_gasdiff_main'
1531	END SUBROUTINE soil_gasdiff_main
1532
1533	!!
1534	!================================================================================================================================
1535	!! SUBROUTINE : soil_gasdiff_alloc
1536	!!
1537	!>\BRIEF This routine allocate arrays related to oxygen and methane in the snow/soil medium
1538	!!
1539	!! DESCRIPTION :
1540	!!
1541	!! RECENT CHANGE(S) : None
1542	!!
1543	!! MAIN OUTPUT VARIABLE(S) :
1544	!!
1545	!! REFERENCE(S) : None
1546	!!
1547	!! FLOWCHART11 : None
1548	!! \n
1549	!_
1550	!================================================================================================================================
1551	SUBROUTINE soil_gasdiff_alloc( kjpindex )
1552
1553	!! 0. Variable and parameter declaration
1554
1555	!! 0.1 Input variables
1556
1557	INTEGER(i_std), INTENT(in) :: kjpindex
1558
1559	!! 0.2 Output variables
1560
1561	!! 0.3 Modified variables
1562
1563	!! 0.4 local variables
1564
1565	INTEGER(i_std) :: ier
1566
1567	! Allocate the variables that need to be saved after soil_gasdiff_coeff
1568
1569	ALLOCATE (alphaO2_soil(kjpindex,ndeep,nvm),stat=ier)
1570	IF (ier.NE.0) THEN
1571	WRITE (numout,*) ' error in alphaO2_soil allocation. We stop. We need', kjpindex, ' fois ',ndeep, ' fois ',nvm,' words = '&
1572	& , kjpindexndeepnvm
1573	STOP 'deep_carbcycle'
1574	END IF
1575
1576	ALLOCATE (betaO2_soil(kjpindex,ndeep,nvm),stat=ier)
1577	IF (ier.NE.0) THEN
1578	WRITE (numout,*) ' error in betaO2_soil allocation. We stop. We need', kjpindex, ' fois ',ndeep, ' fois ',nvm,' words = '&
1579	& , kjpindexndeepnvm
1580	STOP 'deep_carbcycle'
1581	END IF
1582
1583	ALLOCATE (alphaCH4_soil(kjpindex,ndeep,nvm),stat=ier)
1584	IF (ier.NE.0) THEN
1585	WRITE (numout,*) ' error in alphaCH4_soil allocation. We stop. We need', kjpindex, ' fois ',ndeep, ' fois ',nvm,' words = '&
1586	& , kjpindexndeepnvm
1587	STOP 'deep_carbcycle'
1588	END IF
1589
1590	ALLOCATE (betaCH4_soil(kjpindex,ndeep,nvm),stat=ier)
1591	IF (ier.NE.0) THEN
1592	WRITE (numout,*) ' error in betaCH4_soil allocation. We stop. We need', kjpindex, ' fois ',ndeep, ' fois ',nvm,' words = '&
1593	& , kjpindexndeepnvm
1594	STOP 'deep_carbcycle'
1595	END IF
1596
1597	ALLOCATE (alphaO2_snow(kjpindex,nsnow,nvm),stat=ier)
1598	IF (ier.NE.0) THEN
1599	WRITE (numout,*) ' error in alphaO2_snow allocation. We stop. We need', kjpindex, ' fois ',nsnow, ' fois ',nvm,' words = '&
1600	& , kjpindexnsnownvm
1601	STOP 'deep_carbcycle'
1602	END IF
1603
1604	ALLOCATE (betaO2_snow(kjpindex,nsnow,nvm),stat=ier)
1605	IF (ier.NE.0) THEN
1606	WRITE (numout,*) ' error in betaO2_snow allocation. We stop. We need', kjpindex, ' fois ',nsnow, ' fois ',nvm,' words = '&
1607	& , kjpindexnsnownvm
1608	STOP 'deep_carbcycle'
1609	END IF
1610
1611	ALLOCATE (alphaCH4_snow(kjpindex,nsnow,nvm),stat=ier)
1612	IF (ier.NE.0) THEN
1613	WRITE (numout,*) ' error in alphaCH4_snow allocation. We stop. We need', kjpindex, ' fois ',nsnow, ' fois ',nvm,' words = '&
1614	& , kjpindexnsnownvm
1615	STOP 'deep_carbcycle'
1616	END IF
1617
1618	ALLOCATE (betaCH4_snow(kjpindex,nsnow,nvm),stat=ier)
1619	IF (ier.NE.0) THEN
1620	WRITE (numout,*) ' error in betaCH4_snow allocation. We stop. We need', kjpindex, ' fois ',nsnow, ' fois ',nvm,' words = '&
1621	& , kjpindexnsnownvm
1622	STOP 'deep_carbcycle'
1623	END IF
1624
1625	ALLOCATE (zf_coeff_snow(kjpindex,0:nsnow,nvm),stat=ier)
1626	IF (ier.NE.0) THEN
1627	WRITE (numout,*) ' error in zf_coeff_snow allocation. We stop. We need', kjpindex, ' fois ',nsnow+1, ' fois ',nvm,' words = '&
1628	& , kjpindex(nsnow+1)nvm
1629	STOP 'deep_carbcycle'
1630	END IF
1631
1632	ALLOCATE (zi_coeff_snow(kjpindex,nsnow,nvm),stat=ier)
1633	IF (ier.NE.0) THEN
1634	WRITE (numout,*) ' error in zi_coeff_snow allocation. We stop. We need', kjpindex, ' fois ',nsnow, ' fois ',nvm,' words = '&
1635	& , kjpindexnsnownvm
1636	STOP 'deep_carbcycle'
1637	END IF
1638
1639	ALLOCATE (mu_snow(kjpindex,nvm),stat=ier)
1640	IF (ier.NE.0) THEN
1641	WRITE (numout,*) ' error in mu_snow allocation. We stop. We need', kjpindex, ' fois ',nvm,' words = '&
1642	& , kjpindex*nvm
1643	STOP 'deep_carbcycle'
1644	END IF
1645
1646	alphaO2_soil(:,:,:) = zero
1647	betaO2_soil(:,:,:) = zero
1648	alphaCH4_soil(:,:,:) = zero
1649	betaCH4_soil(:,:,:) = zero
1650	alphaO2_snow(:,:,:) = zero
1651	betaO2_snow(:,:,:) = zero
1652	alphaCH4_snow(:,:,:) = zero
1653	betaCH4_snow(:,:,:) = zero
1654	zf_coeff_snow(:,:,:) = zero
1655	zi_coeff_snow(:,:,:) = zero
1656	mu_snow(:,:) = zero
1657
1658	END SUBROUTINE soil_gasdiff_alloc
1659
1660	!!
1661	!================================================================================================================================
1662	!! SUBROUTINE : soil_gasdiff_coeff
1663	!!
1664	!>\BRIEF This routine calculate coeff related to gas diffuvisity
1665	!!
1666	!! DESCRIPTION :
1667	!!
1668	!! RECENT CHANGE(S) : None
1669	!!
1670	!! MAIN OUTPUT VARIABLE(S) :
1671	!!
1672	!! REFERENCE(S) : None
1673	!!
1674	!! FLOWCHART11 : None
1675	!! \n
1676	!_
1677	!================================================================================================================================
1678
1679	SUBROUTINE soil_gasdiff_coeff( kjpindex,time_step,tprof,O2_snow,CH4_snow, &
1680	diffO2_snow,diffCH4_snow,totporO2_snow,totporCH4_snow,O2_soil,CH4_soil, &
1681	diffO2_soil,diffCH4_soil,totporO2_soil,totporCH4_soil, zi_snow, zf_snow)
1682
1683
1684	!! 0. Variable and parameter declaration
1685
1686	!! 0.1 Input variables
1687
1688	INTEGER(i_std), INTENT(in) :: kjpindex !! number of grid points
1689	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
1690	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: tprof !! Soil temperature (K)
1691	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: diffO2_snow !! oxygen diffusivity (m**2/s)
1692	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: diffCH4_snow !! methane diffusivity (m**2/s)
1693	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: totporO2_snow !! total O2 porosity (Tans, 1998)
1694	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: totporCH4_snow !! total CH4 porosity (Tans, 1998)
1695	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: diffO2_soil !! oxygen diffusivity (m**2/s)
1696	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: diffCH4_soil !! methane diffusivity (m**2/s)
1697	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: totporO2_soil !! total O2 porosity (Tans, 1998)
1698	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: totporCH4_soil !! total CH4 porosity (Tans, 1998)
1699	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: O2_snow !! oxygen (g O2/m**3 air)
1700	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: CH4_snow !! methane (g CH4/m**3 air)
1701	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: O2_soil !! oxygen (g O2/m**3 air)
1702	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: CH4_soil !! methane (g CH4/m**3 air)
1703	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm), INTENT(in) :: zf_snow !! depths of full levels (m)
1704	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(in) :: zi_snow !! depths of intermediate levels (m)
1705
1706	!! 0.2 Output variables
1707
1708	!! 0.3 Modified variables
1709
1710	!! 0.4 local variables
1711
1712	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: xcO2_snow,xdO2_snow
1713	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: xcCH4_snow,xdCH4_snow
1714	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm) :: xcO2_soil,xdO2_soil
1715	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm) :: xcCH4_soil,xdCH4_soil
1716	INTEGER(i_std) :: il
1717	REAL(r_std), DIMENSION(kjpindex,nvm) :: xeO2,xeCH4
1718	LOGICAL, DIMENSION(kjpindex,nvm) :: snow_height_mask_2d
1719	LOGICAL, SAVE :: firstcall = .true.
1720	!$OMP THREADPRIVATE(firstcall)
1721
1722	! loop over materials (soil, snow), beginning at the bottom
1723	!
1724	! 1. define useful variables linked to geometry and physical properties
1725	!
1726	! 1.1 normal levels
1727	!
1728	! default value if inexistent
1729	xcO2_snow(:,1,:) = xcO2_soil(:,1,:)
1730	xdO2_snow(:,1,:) = xdO2_soil(:,1,:)
1731	xcCH4_snow(:,1,:) = xcCH4_soil(:,1,:)
1732	xdCH4_snow(:,1,:) = xdCH4_soil(:,1,:)
1733	!
1734	snow_height_mask_2d(:,:) = ( heights_snow(:,:) .GT. hmin_tcalc )
1735	!
1736	DO il = 1,nsnow-1
1737	!
1738	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1739	!
1740	xcO2_snow(:,il,:) = ( zf_snow(:,il,:) - zf_snow(:,il-1,:) ) * &
1741	totporO2_snow(:,il,:) / time_step
1742	xcCH4_snow(:,il,:) = ( zf_snow(:,il,:) - zf_snow(:,il-1,:) ) * &
1743	totporCH4_snow(:,il,:) / time_step
1744	!
1745	xdO2_snow(:,il,:) = diffO2_snow(:,il,:) / &
1746	(zi_snow(:,il+1,:)-zi_snow(:,il,:))
1747	xdCH4_snow(:,il,:) = diffCH4_snow(:,il,:) / &
1748	(zi_snow(:,il+1,:)-zi_snow(:,il,:))
1749	!
1750	ENDWHERE
1751	END DO
1752	!
1753	DO il = 1,ndeep-1
1754	!
1755	WHERE ( veget_mask_2d(:,:) )
1756	!
1757	xcO2_soil(:,il,:) = ( zf_soil(il) - zf_soil(il-1) ) * &
1758	totporO2_soil(:,il,:) / time_step
1759	xcCH4_soil(:,il,:) = ( zf_soil(il) - zf_soil(il-1) ) * &
1760	totporCH4_soil(:,il,:) / time_step
1761	!
1762	xdO2_soil(:,il,:) = diffO2_soil(:,il,:) / &
1763	(zi_soil(il+1)-zi_soil(il))
1764	xdCH4_soil(:,il,:) = diffCH4_soil(:,il,:) / &
1765	(zi_soil(il+1)-zi_soil(il))
1766	!
1767	ENDWHERE
1768	!
1769	ENDDO
1770	!
1771	! 1.2 for the lower boundary, define a similar geometric variable.
1772	!
1773	!snow
1774	!
1775	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1776	xcO2_snow(:,nsnow,:) = ( zf_snow(:,nsnow,:) - &
1777	zf_snow(:,nsnow-1,:) ) * &
1778	totporO2_snow(:,nsnow,:) / time_step
1779	xdO2_snow(:,nsnow,:) = diffO2_snow(:,nsnow,:) / &
1780	( zi_soil(1) + &
1781	zf_snow(:,nsnow,:) - zi_snow(:,nsnow,:) )
1782	xcCH4_snow(:,nsnow,:) = ( zf_snow(:,nsnow,:) - &
1783	zf_snow(:,nsnow-1,:) ) * &
1784	totporCH4_snow(:,nsnow,:) / time_step
1785	xdCH4_snow(:,nsnow,:) = diffCH4_snow(:,nsnow,:) / &
1786	( zi_soil(1) + &
1787	zf_snow(:,nsnow,:) - zi_snow(:,nsnow,:) )
1788	ENDWHERE
1789	!
1790	! soil
1791	!
1792	WHERE ( veget_mask_2d(:,:) ) ! removed heights_soil logic
1793	xcO2_soil(:,ndeep,:) = &
1794	( zf_soil(ndeep) - zf_soil(ndeep-1) ) * &
1795	totporO2_soil(:,ndeep,:) / time_step
1796	xdO2_soil(:,ndeep,:) = diffO2_soil(:,ndeep,:) / &
1797	( zf_soil(ndeep) - zi_soil(ndeep) )
1798	xcCH4_soil(:,ndeep,:) = &
1799	( zf_soil(ndeep) - zf_soil(ndeep-1) ) * &
1800	totporCH4_soil(:,ndeep,:) / time_step
1801	xdCH4_soil(:,ndeep,:) = diffCH4_soil(:,ndeep,:) / &
1802	( zf_soil(ndeep) - zi_soil(ndeep) )
1803	ENDWHERE
1804	!
1805	! 1.3 extrapolation factor from first levels to surface
1806	!
1807	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1808	mu_snow(:,:) = zi_snow(:,1,:) / ( zi_snow(:,2,:) - zi_snow(:,1,:) )
1809	ELSEWHERE ( veget_mask_2d(:,:) )
1810	mu_snow(:,:) = .5 ! any value
1811	ENDWHERE
1812	!
1813	mu_soil = zi_soil(1) / ( zi_soil(2) - zi_soil(1) )
1814	!
1815	! 2. bottom level: treatment depends on lower boundary condition
1816	!
1817	! soil
1818	!
1819	WHERE ( veget_mask_2d(:,:) ) ! removed heights_soil logic
1820	!
1821	xeO2(:,:) = xcO2_soil(:,ndeep,:) + xdO2_soil(:,ndeep-1,:)
1822	xeCH4(:,:) = xcCH4_soil(:,ndeep,:) + xdCH4_soil(:,ndeep-1,:)
1823	!
1824	alphaO2_soil(:,ndeep-1,:) = xdO2_soil(:,ndeep-1,:) / xeO2(:,:)
1825	alphaCH4_soil(:,ndeep-1,:) = xdCH4_soil(:,ndeep-1,:) &
1826	/ xeCH4(:,:)
1827	!
1828	betaO2_soil(:,ndeep-1,:) = &
1829	(xcO2_soil(:,ndeep,:)*O2_soil(:,ndeep,:))/xeO2(:,:)
1830	betaCH4_soil(:,ndeep-1,:) = &
1831	(xcCH4_soil(:,ndeep,:)*CH4_soil(:,ndeep,:))/xeCH4(:,:)
1832	!
1833	ENDWHERE
1834	!
1835	!snow
1836	!
1837	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1838	!
1839	! dernier niveau
1840	!
1841	xeO2(:,:) = xcO2_soil(:,1,:) + &
1842	(1.-alphaO2_soil(:,1,:))*xdO2_soil(:,1,:) + &
1843	xdO2_snow(:,nsnow,:)
1844	xeCH4(:,:) = xcCH4_soil(:,1,:) + &
1845	(1.-alphaCH4_soil(:,1,:))*xdCH4_soil(:,1,:) + &
1846	xdCH4_snow(:,nsnow,:)
1847	!
1848	alphaO2_snow(:,nsnow,:) = xdO2_snow(:,nsnow,:)/xeO2(:,:)
1849	alphaCH4_snow(:,nsnow,:) = xdCH4_snow(:,nsnow,:) &
1850	/xeCH4(:,:)
1851	!
1852	betaO2_snow(:,nsnow,:) = &
1853	( xcO2_soil(:,1,:)*O2_soil(:,1,:) + &
1854	xdO2_soil(:,1,:)*betaO2_soil(:,1,:) ) &
1855	/ xeO2(:,:)
1856	betaCH4_snow(:,nsnow,:) = &
1857	( xcCH4_soil(:,1,:)*CH4_soil(:,1,:) + &
1858	xdCH4_soil(:,1,:)*betaCH4_soil(:,1,:) ) &
1859	/ xeCH4(:,:)
1860	!
1861	! avant-dernier niveau
1862	!
1863	xeO2(:,:) = xcO2_snow(:,nsnow,:) + &
1864	(1.-alphaO2_snow(:,nsnow,:))*xdO2_snow(:,nsnow,:) + &
1865	xdO2_snow(:,nsnow-1,:)
1866	xeCH4(:,:) = xcCH4_snow(:,nsnow,:) + &
1867	(1.-alphaCH4_snow(:,nsnow,:))*xdCH4_snow(:,nsnow,:) &
1868	+ xdCH4_snow(:,nsnow-1,:)
1869	!
1870	alphaO2_snow(:,nsnow-1,:) = &
1871	xdO2_snow(:,nsnow-1,:) / xeO2(:,:)
1872	alphaCH4_snow(:,nsnow-1,:) = &
1873	xdCH4_snow(:,nsnow-1,:) / xeCH4(:,:)
1874	!
1875	betaO2_snow(:,nsnow-1,:) = &
1876	( xcO2_snow(:,nsnow,:)*O2_snow(:,nsnow,:) + &
1877	xdO2_snow(:,nsnow,:)*betaO2_snow(:,nsnow,:) ) &
1878	/ xeO2(:,:)
1879	betaCH4_snow(:,nsnow-1,:) = &
1880	( xcCH4_snow(:,nsnow,:)*CH4_snow(:,nsnow,:) + &
1881	xdCH4_snow(:,nsnow,:)*betaCH4_snow(:,nsnow,:) ) &
1882	/ xeCH4(:,:)
1883	!
1884	ELSEWHERE ( veget_mask_2d(:,:) )
1885	!
1886	alphaO2_snow(:,nsnow,:) = 1.
1887	alphaCH4_snow(:,nsnow,:) = 1.
1888	betaO2_snow(:,nsnow,:) = zero
1889	betaCH4_snow(:,nsnow,:) = zero
1890	!
1891	alphaO2_snow(:,nsnow-1,:) = 1.
1892	alphaCH4_snow(:,nsnow-1,:) = 1.
1893	betaO2_snow(:,nsnow-1,:) = zero
1894	betaCH4_snow(:,nsnow-1,:) = zero
1895	!
1896	ENDWHERE
1897	!
1898
1899	!
1900	! 3. the other levels
1901	!
1902	DO il = nsnow-2,1,-1 !snow
1903	!
1904	WHERE ( snow_height_mask_2d(:,:) .AND. veget_mask_2d(:,:) )
1905	!
1906	xeO2(:,:) = xcO2_snow(:,il+1,:) + &
1907	(1.-alphaO2_snow(:,il+1,:))*xdO2_snow(:,il+1,:) + xdO2_snow(:,il,:)
1908	xeCH4(:,:) = xcCH4_snow(:,il+1,:) + &
1909	(1.-alphaCH4_snow(:,il+1,:))*xdCH4_snow(:,il+1,:) + &
1910	xdCH4_snow(:,il,:)
1911	!
1912	alphaO2_snow(:,il,:) = xdO2_snow(:,il,:) / xeO2(:,:)
1913	alphaCH4_snow(:,il,:) = xdCH4_snow(:,il,:) / xeCH4(:,:)
1914	!
1915	betaO2_snow(:,il,:) = &
1916	( xcO2_snow(:,il+1,:)*O2_snow(:,il+1,:) + &
1917	xdO2_snow(:,il+1,:)*betaO2_snow(:,il+1,:) ) / xeO2(:,:)
1918	betaCH4_snow(:,il,:) = &
1919	( xcCH4_snow(:,il+1,:)*CH4_snow(:,il+1,:) + &
1920	xdCH4_snow(:,il+1,:)*betaCH4_snow(:,il+1,:) ) / xeCH4(:,:)
1921	!
1922	ELSEWHERE ( veget_mask_2d(:,:) )
1923	!
1924	alphaO2_snow(:,il,:) = 1.
1925	alphaCH4_snow(:,il,:) = 1.
1926	!
1927	betaO2_snow(:,il,:) = zero
1928	betaCH4_snow(:,il,:) = zero
1929	!
1930	ENDWHERE
1931	!
1932	ENDDO
1933	!
1934	DO il = ndeep-2,1,-1 !soil
1935	!
1936	WHERE ( veget_mask_2d(:,:) ) !removed heights_soil logic
1937	!
1938	xeO2(:,:) = xcO2_soil(:,il+1,:) + &
1939	(1.-alphaO2_soil(:,il+1,:))*xdO2_soil(:,il+1,:) + xdO2_soil(:,il,:)
1940	xeCH4(:,:) = xcCH4_soil(:,il+1,:) + &
1941	(1.-alphaCH4_soil(:,il+1,:))*xdCH4_soil(:,il+1,:) + &
1942	xdCH4_soil(:,il,:)
1943	!
1944	alphaO2_soil(:,il,:) = xdO2_soil(:,il,:) / xeO2(:,:)
1945	alphaCH4_soil(:,il,:) = xdCH4_soil(:,il,:) / xeCH4(:,:)
1946	!
1947	betaO2_soil(:,il,:) = &
1948	( xcO2_soil(:,il+1,:)*O2_soil(:,il+1,:) + &
1949	xdO2_soil(:,il+1,:)*betaO2_soil(:,il+1,:) ) / xeO2(:,:)
1950	betaCH4_soil(:,il,:) = &
1951	( xcCH4_soil(:,il+1,:)*CH4_soil(:,il+1,:) + &
1952	xdCH4_soil(:,il+1,:)*betaCH4_soil(:,il+1,:) ) / xeCH4(:,:)
1953	!
1954	ENDWHERE
1955	!
1956	ENDDO
1957	!
1958	! 4. store thickness of the different levels for all soil types (for security)
1959	!
1960	zf_coeff_snow(:,:,:) = zf_snow(:,:,:)
1961	zi_coeff_snow(:,:,:) = zi_snow(:,:,:)
1962
1963	!--hist out for keeping track of these
1964	IF (firstcall) THEN
1965	firstcall = .false.
1966	ELSE
1967	ENDIF
1968
1969	END SUBROUTINE soil_gasdiff_coeff
1970
1971	!!
1972	!================================================================================================================================
1973	!! SUBROUTINE : soil_gasdiff_diff
1974	!!
1975	!>\BRIEF This routine update oxygen and methane in the snow and soil
1976	!!
1977	!! DESCRIPTION :
1978	!!
1979	!! RECENT CHANGE(S) : None
1980	!!
1981	!! MAIN OUTPUT VARIABLE(S) :
1982	!!
1983	!! REFERENCE(S) : None
1984	!!
1985	!! FLOWCHART11 : None
1986	!! \n
1987	!_
1988	!================================================================================================================================
1989
1990	SUBROUTINE soil_gasdiff_diff( kjpindex,time_step,index,pb,tsurf, O2_snow, CH4_snow, O2_soil, CH4_soil)
1991
1992	!! 0. Variable and parameter declaration
1993
1994	!! 0.1 Input variables
1995
1996	INTEGER(i_std), INTENT(in) :: kjpindex !! number of grid points
1997	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
1998	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure
1999	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: tsurf !! Surface temperature
2000	INTEGER(i_std),DIMENSION(kjpindex),INTENT(in) :: index !! Indeces of the points on the map
2001	!! 0.2 Output variables
2002
2003	!! 0.3 Modified variables
2004
2005	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: O2_snow !! oxygen (g O2/m**3 air)
2006	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: CH4_snow !! methane (g CH4/m**3 air)
2007	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: O2_soil !! oxygen (g O2/m**3 air)
2008	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: CH4_soil !! methane (g CH4/m**3 air)
2009
2010	!! 0.4 local variables
2011
2012	INTEGER(i_std) :: it, ip, il, iv
2013	LOGICAL, DIMENSION(kjpindex,nvm) :: snowtop
2014	REAL(r_std), DIMENSION(kjpindex,nvm) :: O2sa, CH4sa
2015
2016	!
2017	! 1.1 Determine which is the first existing soil type.
2018	!
2019	snowtop(:,:) = .FALSE.
2020	!
2021	!ignore snow for now...
2022	WHERE ( heights_snow(:,:) .GT. hmin_tcalc )
2023	snowtop(:,:) = .TRUE.
2024	ENDWHERE
2025	!
2026	! 2.gas diffusion
2027	!
2028	! 2.1 top level
2029	!
2030	! 2.1.1 non-existing
2031	!
2032	DO iv = 1, nvm
2033	O2sa(:,iv) = pb(:)/(RRtsurf(:)) O2_surf * wO2
2034	CH4sa(:,iv) = pb(:)/(RRtsurf(:)) CH4_surf * wCH4
2035	ENDDO
2036	!
2037	WHERE ( (.NOT. snowtop(:,:)) .AND. veget_mask_2d(:,:) ) ! it equals 1 (snow) but there is no snow...
2038	!
2039	O2_snow(:,1,:) = O2sa(:,:)
2040	CH4_snow(:,1,:) = CH4sa(:,:)
2041	!
2042	O2_soil(:,1,:) = ( O2sa(:,:) + mu_soil*betaO2_soil(:,1,:) ) / &
2043	( 1. + mu_soil*(1.-alphaO2_soil(:,1,:)) )
2044	CH4_soil(:,1,:) = ( CH4sa(:,:) + mu_soil*betaCH4_soil(:,1,:) ) / &
2045	( 1. + mu_soil*(1.-alphaCH4_soil(:,1,:)) )
2046	!
2047	ENDWHERE
2048	!
2049	! 2.1.2 first existing soil type
2050	!
2051	WHERE ( snowtop(:,:) .AND. veget_mask_2d(:,:) )
2052	!
2053	O2_snow(:,1,:) = ( O2sa(:,:) + mu_snow(:,:)*betaO2_snow(:,1,:) ) / &
2054	( 1. + mu_snow(:,:)*(1.-alphaO2_snow(:,1,:)) )
2055	CH4_snow(:,1,:) = ( CH4sa(:,:) + mu_snow(:,:)*betaCH4_snow(:,1,:) ) / &
2056	( 1. + mu_snow(:,:)*(1.-alphaCH4_snow(:,1,:)) )
2057	!
2058	O2_soil(:,1,:) = &
2059	alphaO2_snow(:,nsnow,:) * O2_snow(:,nsnow,:) + &
2060	betaO2_snow(:,nsnow,:)
2061	CH4_soil(:,1,:) = &
2062	alphaCH4_snow(:,nsnow,:) * CH4_snow(:,nsnow,:) + &
2063	betaCH4_snow(:,nsnow,:)
2064	! debug: need to check for weird numbers here!
2065	ENDWHERE
2066	!
2067	! 2.2 other levels
2068	!
2069	DO il = 2, nsnow
2070
2071	WHERE ( veget_mask_2d(:,:) )
2072	!
2073	O2_snow(:,il,:) = &
2074	alphaO2_snow(:,il-1,:) * O2_snow(:,il-1,:) + &
2075	betaO2_snow(:,il-1,:)
2076	CH4_snow(:,il,:) = &
2077	alphaCH4_snow(:,il-1,:) * CH4_snow(:,il-1,:) + &
2078	betaCH4_snow(:,il-1,:)
2079	END WHERE
2080	ENDDO
2081	DO il = 2, ndeep
2082
2083	WHERE ( veget_mask_2d(:,:) )
2084	!
2085	O2_soil(:,il,:) = &
2086	alphaO2_soil(:,il-1,:) * O2_soil(:,il-1,:) + &
2087	betaO2_soil(:,il-1,:)
2088	CH4_soil(:,il,:) = &
2089	alphaCH4_soil(:,il-1,:) * CH4_soil(:,il-1,:) + &
2090	betaCH4_soil(:,il-1,:)
2091	END WHERE
2092	ENDDO
2093
2094	END SUBROUTINE soil_gasdiff_diff
2095
2096	!!
2097	!================================================================================================================================
2098	!! SUBROUTINE : get_gasdiff
2099	!!
2100	!>\BRIEF This routine update oxygen and methane in the snow and soil
2101	!!
2102	!! DESCRIPTION :
2103	!!
2104	!! RECENT CHANGE(S) : None
2105	!!
2106	!! MAIN OUTPUT VARIABLE(S) :
2107	!!
2108	!! REFERENCE(S) : None
2109	!!
2110	!! FLOWCHART11 : None
2111	!! \n
2112	!_
2113	!================================================================================================================================
2114	SUBROUTINE get_gasdiff (kjpindex,hslong,tprof,snow,airvol_snow, &
2115	totporO2_snow,totporCH4_snow,diffO2_snow,diffCH4_snow, &
2116	airvol_soil,totporO2_soil,totporCH4_soil,diffO2_soil,diffCH4_soil, z_organic, snowrho)
2117
2118	!! 0. Variable and parameter declaration
2119
2120	!! 0.1 Input variables
2121
2122	INTEGER(i_std), INTENT(in) :: kjpindex !! number of grid points
2123	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: hslong !! deep long term soil humidity profile
2124	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: tprof !! Soil temperature (K)
2125	REAL(r_std), DIMENSION(kjpindex,nsnow), INTENT(in) :: snowrho !! snow density
2126	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: snow !! Snow mass [Kg/m^2]
2127	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: z_organic !! depth to organic soil
2128
2129	!! 0.2 Output variables
2130
2131	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(out) :: airvol_soil
2132	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(out) :: totporO2_soil !! total O2 porosity (Tans, 1998)
2133	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(out) :: totporCH4_soil !! total CH4 porosity
2134	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(out) :: diffO2_soil !! oxygen diffusivity (m**2/s)
2135	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(out) :: diffCH4_soil !! methane diffusivity (m**2/s)
2136	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: airvol_snow
2137	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: totporO2_snow !! total O2 porosity (Tans, 1998)
2138	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: totporCH4_snow !! total CH4 porosity (Tans, 1998)
2139	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: diffO2_snow !! oxygen diffusivity (m**2/s)
2140	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(out) :: diffCH4_snow !! methane diffusivity (m**2/s)
2141
2142	!! 0.3 Modified variables
2143
2144	!! 0.4 local variables
2145
2146	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: density_snow
2147	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: porosity_snow
2148	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: tortuosity_snow
2149	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm) :: density_soil
2150	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm) :: porosity_soil
2151	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm) :: tortuosity_soil
2152	INTEGER(i_std) :: it,ip, il, iv
2153	REAL(r_std) :: x, rho_iw
2154	REAL(r_std) :: csat, fng
2155	REAL(r_std), SAVE :: cond_fact
2156	!$OMP THREADPRIVATE(cond_fact)
2157	LOGICAL, SAVE :: pr_fois=.TRUE.
2158	!$OMP THREADPRIVATE(pr_fois)
2159
2160	IF (pr_fois) THEN
2161	cond_fact=1.
2162	CALL getin_p('COND_FACT',cond_fact)
2163	WRITE(,) 'COND_FACT=',cond_fact
2164	pr_fois=.FALSE.
2165	ENDIF
2166
2167	!
2168	! 1. Three-layers snow model with snow density resolved at each snow layer
2169	!
2170	DO iv = 1, nvm
2171	density_snow(:,:,iv) = snowrho(:,:)
2172	ENDDO
2173	porosity_snow(:,:,:) = (1. - density_snow(:,:,:)/rho_ice )
2174	tortuosity_snow(:,:,:) = porosity_snow(:,:,:)**(1./3.) ! based on Sommerfeld et al., GBC, 1996
2175	diffO2_snow(:,:,:) = diffO2_air * porosity_snow(:,:,:) * tortuosity_snow(:,:,:)
2176	diffCH4_snow(:,:,:) = diffCH4_air * porosity_snow(:,:,:) * tortuosity_snow(:,:,:)
2177	airvol_snow(:,:,:) = MAX(porosity_snow(:,:,:),avm)
2178	totporO2_snow(:,:,:) = airvol_snow(:,:,:)
2179	totporCH4_snow(:,:,:) = airvol_snow(:,:,:)
2180	!
2181	! 2. soil: depends on temperature and soil humidity
2182	!
2183	DO ip = 1, kjpindex
2184	!
2185	DO iv = 1, nvm
2186	!
2187	IF ( veget_mask_2d(ip,iv) ) THEN
2188	!
2189	DO il = 1, ndeep
2190	!
2191	! 2.1 soil dry density, porosity, and dry heat capacity
2192	!
2193	porosity_soil(ip,il,iv) = tetasat
2194	!
2195	!
2196	! 2.2 heat capacity and density as a function of
2197	! ice and water content
2198	! removed these as we are calculating thermal evolution in the sechiba subroutines
2199
2200	!
2201	! 2.3 oxygen diffusivity: soil can get waterlogged,
2202	! therefore take soil humidity into account
2203	!
2204	tortuosity_soil(ip,il,iv) = 2./3. ! Hillel, 1980
2205	airvol_soil(ip,il,iv) = porosity_soil(ip,il,iv)*(1.-hslong(ip,il,iv))
2206	totporO2_soil(ip,il,iv) = airvol_soil(ip,il,iv) + porosity_soil(ip,il,iv)BunsenO2hslong(ip,il,iv)
2207	totporCH4_soil(ip,il,iv) = airvol_soil(ip,il,iv) + porosity_soil(ip,il,iv)BunsenCH4hslong(ip,il,iv)
2208	diffO2_soil(ip,il,iv) = (diffO2_air*airvol_soil(ip,il,iv) + &
2209	diffO2_wBunsenO2hslong(ip,il,iv)porosity_soil(ip,il,iv))tortuosity_soil(ip,il,iv)
2210	diffCH4_soil(ip,il,iv) = (diffCH4_air*airvol_soil(ip,il,iv) + &
2211	diffCH4_wBunsenCH4hslong(ip,il,iv)porosity_soil(ip,il,iv))tortuosity_soil(ip,il,iv)
2212	!
2213	END DO
2214	ELSE
2215	tortuosity_soil(ip,:,iv) = EPSILON(0.)
2216	airvol_soil(ip,:,iv) = EPSILON(0.)
2217	totporO2_soil(ip,:,iv) = EPSILON(0.)
2218	totporCH4_soil(ip,:,iv) = EPSILON(0.)
2219	diffO2_soil(ip,:,iv) = EPSILON(0.)
2220	diffCH4_soil(ip,:,iv) = EPSILON(0.)
2221	END IF
2222	ENDDO
2223	ENDDO
2224
2225	END SUBROUTINE get_gasdiff
2226
2227	!!
2228	!================================================================================================================================
2229	!! SUBROUTINE : traMplan
2230	!!
2231	!>\BRIEF This routine calculates plant-mediated transport of methane
2232	!!
2233	!! DESCRIPTION :
2234	!!
2235	!! RECENT CHANGE(S) : None
2236	!!
2237	!! MAIN OUTPUT VARIABLE(S) :
2238	!!
2239	!! REFERENCE(S) : None
2240	!!
2241	!! FLOWCHART11 : None
2242	!! \n
2243	!_
2244	!================================================================================================================================
2245	SUBROUTINE traMplan(CH4,O2,kjpindex,time_step,totporCH4,totporO2,z_root,rootlev,Tgr,Tref,hslong,flupmt, &
2246	refdep, zi_soil, tprof)
2247
2248	!! 0. Variable and parameter declaration
2249
2250	!! 0.1 Input variables
2251
2252	INTEGER(i_std), INTENT(in) :: kjpindex
2253	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
2254	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: totporO2 !! total oxygen porosity
2255	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: totporCH4 !! total methane porosity
2256	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm),INTENT(in) :: tprof !! soil temperature (K)
2257	INTEGER(i_std),DIMENSION(kjpindex,nvm),INTENT(in) :: rootlev !! the deepest model level within the rooting depth
2258	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: z_root !! the rooting depth
2259	REAL(r_std), INTENT(in) :: Tgr !! Temperature at which plants begin to grow (C)
2260	REAL(r_std), DIMENSION(ndeep), INTENT(in) :: zi_soil !! depths at intermediate levels
2261	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: hslong !! deep soil humidity
2262
2263	!! 0.2 Output variables
2264
2265	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: flupmt !! plant-mediated methane flux (g m-2 s-1)
2266
2267	!! 0.3 Modified variables
2268
2269	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(inout) :: Tref !! Ref. temperature for growing season caluculation (C)
2270	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: O2
2271	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: CH4
2272
2273	!! 0.4 local variables
2274	REAL(r_std), DIMENSION(kjpindex,nvm) :: CH4atm !! CH4 atm concentration
2275	REAL(r_std), DIMENSION(kjpindex,nvm) :: dCH4 !! delta CH4 per m3 air
2276	REAL(r_std), DIMENSION(kjpindex,nvm) :: dO2 !! O2 change
2277	REAL(r_std), DIMENSION(kjpindex,nvm) :: fgrow !! Plant growing state (maturity index)
2278	REAL(r_std) :: froot !! vertical distribution of roots
2279	REAL(r_std) :: Tmat !! Temperature at which plants reach maturity (C)
2280	REAL(r_std), PARAMETER :: La_min = zero
2281	REAL(r_std), PARAMETER :: La = 4.
2282	REAL(r_std), PARAMETER :: La_max = La_min + La
2283	REAL(r_std), PARAMETER :: Tveg = 10 !! Vegetation type control on the plant-mediated transport, Adjustable parameter,
2284	!! but we start from 10 following Walter et al (2001) tundra value
2285	REAL(r_std), PARAMETER :: Pox = 0.5 !! fraction of methane oxydized near the roots
2286	LOGICAL, SAVE :: firstcall=.TRUE.
2287	!$OMP THREADPRIVATE(firstcall)
2288	INTEGER(i_std) :: il,ip, iv
2289	LOGICAL, SAVE :: check = .FALSE.
2290	!$OMP THREADPRIVATE(check)
2291	REAL(r_std), INTENT(in) :: refdep !! Depth to compute reference temperature for the growing season (m)
2292	INTEGER(i_std), SAVE :: reflev = 0 !! Level closest to reference depth refdep
2293	!$OMP THREADPRIVATE(reflev)
2294
2295
2296	IF (firstcall) THEN
2297	firstcall = .FALSE.
2298
2299	! Find the level closest to refdep
2300	DO il=1,ndeep
2301	IF (zi_soil(il) .GT. refdep .AND. reflev.EQ.0) reflev = il-1
2302	ENDDO
2303	IF (reflev.EQ.0) reflev = ndeep
2304
2305
2306	IF (check) THEN
2307	OPEN (28,file='pmt.dat',status='unknown')
2308	OPEN (29,file='pmtf.dat',status='unknown')
2309	ENDIF
2310	ENDIF
2311
2312	! Update seasonal reference temperature trace record
2313	WHERE ( veget_mask_2d(:,:) )
2314	Tref(:,:) = tprof(:,reflev,:) - ZeroCelsius
2315	END WHERE
2316
2317	Tmat = Tgr + 10._r_std
2318	flupmt(:,:) = zero
2319	CH4atm(:,:) = zero
2320
2321
2322	! Plant growing state (maturity index)
2323	WHERE (Tref(:,:).LE.Tgr .AND. veget_mask_2d(:,:) )
2324	fgrow(:,:) = La_min
2325	ELSEWHERE (Tref(:,:).GE.Tmat .AND. veget_mask_2d(:,:) )
2326	fgrow(:,:) = La_max
2327	ELSEWHERE ( veget_mask_2d(:,:))
2328	fgrow(:,:) = La_min + La * (1 - ((Tmat - Tref(:,:))/(Tmat - Tgr))**2)
2329	ENDWHERE
2330
2331	DO ip=1,kjpindex
2332	DO iv = 1, nvm
2333	IF ( (z_root(ip,iv) .GT. 0.) .AND. veget_mask_2d(ip,iv) ) THEN ! added this to prevent pmt calcs when soil frozen
2334	DO il=1,rootlev(ip,iv)
2335	! vertical distribution of roots
2336	froot = MAX( 2 * (z_root(ip,iv) - REAL( zi_soil(il) )) / z_root(ip,iv), zero)
2337	! Methane removal from a given depth. We assume that the methane
2338	! in air pores is always in equilibrium with that dissolved
2339	! in water-filled pores. If soil humidity is low,
2340	! with root water as well
2341	! We assume that PMT is proportional to soil humidity
2342	dCH4(ip,iv) = 0.01_r_std * Tveg * froot * fgrow(ip,iv) * hslong(ip,il,iv) * (CH4(ip,il,iv) - CH4atm(ip,iv))
2343	! No transport if soil concentration is less than atmospheric
2344	IF (dCH4(ip,iv).LT.CH4atm(ip,iv)) dCH4(ip,iv) = zero
2345	! Strange thing in WH 2001: 0.01Tvegfroot*fgrow > 1
2346	! at Tveg=15, froot&fgrow=max, i.e. more CH4 is taken than available
2347	! So need to impose a limitation:
2348	IF (dCH4(ip,iv).GT.CH4(ip,il,iv)) dCH4(ip,iv) = CH4(ip,il,iv)
2349	! Methane concentration is decreased within the root layer:
2350
2351	CH4(ip,il,iv) = CH4(ip,il,iv) - dCH4(ip,iv)
2352	! O2 concentration is decreased in reaction with
2353	! dCH4Poxtime_step
2354	dO2(ip,iv) = dCH4(ip,iv)Pox wO2/wCH4 * totporCH4(ip,il,iv)/totporO2(ip,il,iv)
2355	IF ( dO2(ip,iv).LT.O2(ip,il,iv) ) O2(ip,il,iv) = O2(ip,il,iv) - dO2(ip,iv)
2356
2357	! CO2 concentration is increased by dCH4(:)*Pox
2358
2359	! Integration
2360	flupmt(ip,iv) = flupmt(ip,iv) + dCH4(ip,iv)totporCH4(ip,il,iv)/time_step (1 - Pox) * &
2361	( zf_soil(il) - zf_soil(il-1) )
2362	ENDDO
2363	END IF
2364	ENDDO
2365	ENDDO
2366
2367	IF (check) THEN
2368	WRITE(29,*) flupmt(:,:)
2369	CALL flush(28)
2370	CALL flush(29)
2371	END IF
2372
2373	END SUBROUTINE traMplan
2374
2375	!!
2376	!================================================================================================================================
2377	!! SUBROUTINE : ebullition
2378	!!
2379	!>\BRIEF This routine calculates CH4 ebullition
2380	!!
2381	!! DESCRIPTION :
2382	!!
2383	!! RECENT CHANGE(S) : None
2384	!!
2385	!! MAIN OUTPUT VARIABLE(S) :
2386	!!
2387	!! REFERENCE(S) : None
2388	!!
2389	!! FLOWCHART11 : None
2390	!! \n
2391	!_
2392	!================================================================================================================================
2393	SUBROUTINE ebullition (kjpindex,time_step,tprof,totporCH4_soil,hslong,Ch4_soil,febul)
2394
2395	!! 0. Variable and parameter declaration
2396
2397	!! 0.1 Input variables
2398
2399	INTEGER(i_std), INTENT(in) :: kjpindex
2400	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
2401	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm),INTENT(in) :: tprof !! soil temperature (K)
2402	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: totporCH4_soil !! total methane porosity
2403	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: hslong !! deep soil humidity
2404
2405	!! 0.2 Output variables
2406
2407	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: febul !! CH4 ebullition
2408
2409	!! 0.3 Modified variables
2410
2411	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: Ch4_soil !! methane
2412
2413	!! 0.4 Local variables
2414	REAL(r_std) :: dCH4, CH4d
2415	INTEGER(i_std) :: ip, il, iv
2416	REAL(r_std) :: dz
2417	REAL(r_std), PARAMETER :: tortuosity=2./3.
2418	REAL(r_std), PARAMETER :: wsize=0.01
2419	REAL(r_std), PARAMETER :: CH4wm = 12.E-3 !! CH4 concentration threshold for ebullition (8-16 mg/m3 in Walter&Heimann 2000)
2420	REAL(r_std) :: hum
2421
2422	DO ip=1,kjpindex
2423	DO iv = 1, nvm
2424	IF ( veget_mask_2d(ip,iv) ) THEN
2425	febul(ip,iv) = zero
2426	IF (hslong(ip,1,iv).GT.ebuthr) THEN
2427	DO il = ndeep, 1, -1
2428	CH4d = Ch4_soil(ip,il,iv) - CH4wm/BunsenCH4
2429	IF (CH4d .GT. EPSILON(0.)) THEN
2430	IF (il.GT.1) THEN
2431	dz = zi_soil(il) - zi_soil(il-1)
2432	hum = ( hslong(ip,il,iv) + hslong(ip,il-1,iv) ) / 2
2433	ELSE
2434	dz = zi_soil(1)
2435	hum = hslong(ip,1,iv)
2436	ENDIF
2437
2438	dCH4 = hum*( dz/wsize/tortuosity ) CH4d
2439	dCH4 = CH4d
2440
2441	Ch4_soil(ip,il,iv) = Ch4_soil(ip,il,iv) - dCH4
2442
2443
2444	febul(ip,iv) = febul(ip,iv) + dCH4 * totporCH4_soil(ip,il,iv) * &
2445	( zf_soil(il) - zf_soil(il-1) ) / time_step
2446
2447	ENDIF
2448	ENDDO
2449	ENDIF
2450	END IF
2451	ENDDO
2452	ENDDO
2453	END SUBROUTINE ebullition
2454
2455	!!
2456	!================================================================================================================================
2457	!! SUBROUTINE : microactem
2458	!!
2459	!>\BRIEF This routine calculates parameters describing bacterial activity (time constant tau[s]) as a function of temperature
2460	!!
2461	!! DESCRIPTION :
2462	!!
2463	!! RECENT CHANGE(S) : None
2464	!!
2465	!! MAIN OUTPUT VARIABLE(S) :
2466	!!
2467	!! REFERENCE(S) : None
2468	!!
2469	!! FLOWCHART11 : None
2470	!! \n
2471	!_
2472	!================================================================================================================================
2473	FUNCTION microactem ( temp_kelvin, frozen_respiration_func, moist_in, i_ind, j_ind, k_ind, zi_soil ) RESULT ( fbact )
2474
2475	!! 0. Variable and parameter declaration
2476
2477	!! 0.1 Input variables
2478
2479	INTEGER(i_std), INTENT(in) :: i_ind
2480	INTEGER(i_std), INTENT(in) :: j_ind
2481	INTEGER(i_std), INTENT(in) :: k_ind
2482	INTEGER(i_std), INTENT(in) :: frozen_respiration_func
2483	REAL, DIMENSION(i_ind, j_ind, k_ind), INTENT(in) :: moist_in
2484	REAL, DIMENSION(i_ind, j_ind, k_ind), INTENT(in) :: temp_kelvin
2485	REAL, DIMENSION(j_ind), INTENT(in) :: zi_soil
2486
2487	!! 0.2 Output variables
2488
2489	!! 0.3 Modified variables
2490
2491	!! 0.4 Local variables
2492
2493	REAL, DIMENSION(i_ind, j_ind, k_ind) :: fbact
2494	REAL, DIMENSION(i_ind, j_ind, k_ind) :: tempfunc_result
2495	INTEGER(i_std), PARAMETER :: ntconfun = 7
2496	REAL(r_std), DIMENSION(ntconfun) :: tconfun
2497	REAL(r_std), DIMENSION(ntconfun) :: tauconfun
2498	INTEGER :: itz
2499	INTEGER :: ii, ij, ik
2500	REAL, DIMENSION(i_ind, j_ind, k_ind) :: moistfunc_result
2501	REAL(r_std), parameter :: q10 = 2.0
2502	REAL(r_std), PARAMETER :: stomate_tau = 4.699E6 !4.7304E7 !4.699E6
2503	logical, parameter :: limit_decomp_moisture = .true.
2504	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
2505	!$OMP THREADPRIVATE(depth_modifier)
2506	LOGICAL, SAVE :: is_firstcall = .TRUE.
2507	!$OMP THREADPRIVATE(is_firstcall)
2508
2509	IF (is_firstcall) THEN
2510	!
2511	!Config Key = depth_modifier
2512	!Config Desc = e-folding depth of turnover rates,following Koven et
2513	!Config DEsc al.,2013,Biogeosciences. A very large value means no depth modification
2514	!Config Def = 1.E6
2515	!Config If =
2516	!Config Help =
2517	!Config
2518	!Config Units = [-]
2519	!
2520	CALL getin_p('depth_modifier',depth_modifier)
2521	is_firstcall = .FALSE.
2522	ENDIF
2523
2524	SELECT CASE(frozen_respiration_func)
2525
2526	CASE(0) ! this is the standard ORCHIDEE state
2527
2528	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2529	tempfunc_result(:,:,:) = MIN( 1._r_std, tempfunc_result(:,:,:) )
2530
2531	CASE(1) ! cutoff respiration when T < -1C
2532	WHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius ) ! normal as above
2533	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2534	ELSEWHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius - 1. ) ! linear dropoff to zero
2535	tempfunc_result(:,:,:) = (temp_kelvin(:,:,:) - (ZeroCelsius - 1.)) * &
2536	EXP( log(q10) * ( ZeroCelsius - (ZeroCelsius+30.) ) / 10. )
2537	ELSEWHERE ! zero
2538	tempfunc_result(:,:,:) = EPSILON(0.)
2539	endwhere
2540
2541	tempfunc_result(:,:,:) = MAX(MIN( 1._r_std, tempfunc_result(:,:,:) ), EPSILON(0.))
2542
2543	CASE(2) ! cutoff respiration when T < -3C
2544	WHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius ) ! normal as above
2545	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2546	ELSEWHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius - 3. ) ! linear dropoff to zero
2547	tempfunc_result(:,:,:) = ((temp_kelvin(:,:,:) - (ZeroCelsius - 3.))/3.) * &
2548	EXP( log(q10) * ( ZeroCelsius - (ZeroCelsius+30.) ) / 10. )
2549	ELSEWHERE ! zero
2550	tempfunc_result(:,:,:) = EPSILON(0.)
2551	endwhere
2552
2553	CASE(3) ! q10 = 100 when below zero
2554	WHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius ) ! normal as above
2555	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2556	ELSEWHERE
2557	tempfunc_result(:,:,:) = EXP( log(100.) * ( temp_kelvin(:,:,:) - (ZeroCelsius) ) / 10. ) * &
2558	EXP( log(q10) * ( -30. ) / 10. )
2559	endwhere
2560
2561	CASE(4) ! q10 = 1000 when below zero
2562	WHERE (temp_kelvin(:,:,:) .GT. ZeroCelsius ) ! normal as above
2563	tempfunc_result(:,:,:) = EXP( log(q10) * ( temp_kelvin(:,:,:) - (ZeroCelsius+30.) ) / 10. )
2564	ELSEWHERE
2565	tempfunc_result(:,:,:) = EXP( log(1000.) * ( temp_kelvin(:,:,:) - (ZeroCelsius) ) / 10. ) * &
2566	EXP( log(q10) * ( -30. ) / 10. )
2567	endwhere
2568
2569	CASE DEFAULT
2570	WRITE(,) 'microactem ERROR: frozen_respiration_func not in list: ', frozen_respiration_func
2571	STOP
2572
2573	END SELECT
2574	tempfunc_result(:,:,:) = MAX(MIN( 1._r_std, tempfunc_result(:,:,:) ), EPSILON(0.))
2575
2576	!---- stomate residence times: -----!
2577	! residence times in carbon pools (days)
2578	!carbon_tau(iactive) = .149 * one_year !!!!???? 1.5 years
2579	!carbon_tau(islow) = 5.48 * one_year !!!!???? 25 years
2580	!carbon_tau(ipassive) = 241. * one_year !!!!???? 1000 years
2581	!-----------------------------------!
2582	IF ( limit_decomp_moisture ) THEN
2583	! stomate moisture control function
2584	moistfunc_result(:,:,:) = -1.1 * moist_in(:,:,:) * moist_in(:,:,:) + 2.4 * moist_in(:,:,:) - 0.29
2585	moistfunc_result(:,:,:) = max( 0.25_r_std, min( 1._r_std, moistfunc_result(:,:,:) ) )
2586	ELSE
2587	moistfunc_result(:,:,:) = 1._r_std
2588	ENDIF
2589
2590	DO ij = 1, ndeep
2591	fbact(:,ij,:) = stomate_tau/(moistfunc_result(:,ij,:) * tempfunc_result(:,ij,:)) / EXP(-zi_soil(ij)/depth_modifier)
2592	ENDDO
2593
2594	! chaoyue: We tentatively increase the turnover of soil C in croplands,
2595	! as shown here to decrease its tau -- the residence time.
2596	DO ik = 1,nvm
2597	IF ( (.NOT. natural(ik)) .AND. (.NOT. is_c4(ik)) .AND. (.NOT. pasture(ik)) ) THEN
2598	fbact(:,:,ik) = fbact(:,:,ik)/flux_tot_coeff(1)
2599	ELSEIF ( (.NOT. natural(ik)) .AND. is_c4(ik) .AND. (.NOT. pasture(ik)) ) THEN
2600	fbact(:,:,ik) = fbact(:,:,ik)/flux_tot_coeff(2)
2601	ENDIF
2602	ENDDO
2603
2604	END FUNCTION microactem
2605
2606
2607	!!
2608	!================================================================================================================================
2609	!! SUBROUTINE : snowlevels
2610	!!
2611	!>\BRIEF This routine calculates depths of full levels and intermediate
2612	!! levels related to snow pack
2613	!!
2614	!! DESCRIPTION :
2615	!!
2616	!! RECENT CHANGE(S) : None
2617	!!
2618	!! MAIN OUTPUT VARIABLE(S) :
2619	!!
2620	!! REFERENCE(S) : None
2621	!!
2622	!! FLOWCHART11 : None
2623	!! \n
2624	!_
2625	!================================================================================================================================
2626
2627	SUBROUTINE snowlevels( kjpindex, snowdz, zi_snow, zf_snow, veget_max )
2628
2629	!! 0. Variable and parameter declaration
2630
2631	!! 0.1 Input variables
2632
2633	INTEGER(i_std), INTENT(in) :: kjpindex
2634	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: veget_max !! maximum vegetation fraction
2635	REAL(r_std), DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz !! snow depth
2636
2637	!! 0.2 Output variables
2638
2639	!! 0.3 Modified variables
2640
2641	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm), INTENT(inout) :: zf_snow !! depths of full levels (m)
2642	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: zi_snow !! depths of intermediate levels (m)
2643
2644	!! 0.4 Local variables
2645
2646	REAL(r_std), DIMENSION(kjpindex,nvm) :: z_alpha !! parameter of the geometric series
2647	INTEGER(i_std) :: il,it, ix, iv
2648	INTEGER(i_std) :: it_beg,it_end
2649	INTEGER(i_std), PARAMETER :: niter = 10
2650	REAL(r_std), DIMENSION(kjpindex) :: dxmin
2651	INTEGER(i_std), DIMENSION(kjpindex) :: imin
2652	INTEGER(i_std) :: i,j
2653	REAL(r_std), DIMENSION(kjpindex,nvm) :: xi, xf
2654	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: snowdz_pft
2655
2656	snowdz_pft(:,:,:) = 0.0
2657	DO il = 1,nsnow
2658	DO iv = 1, nvm
2659	WHERE ( veget_mask_2d(:,iv) )
2660	snowdz_pft(:,il,iv) = snowdz(:,il)
2661	ENDWHERE
2662	ENDDO
2663	ENDDO
2664	!
2665	! calculate snow discretisation
2666	!
2667	WHERE ( veget_mask_2d(:,:) )
2668	zf_snow(:,0,:) = 0.
2669	END WHERE
2670	!
2671	DO il = 1, nsnow
2672	IF ( il .EQ. 1 ) THEN
2673	WHERE ( veget_mask_2d(:,:) )
2674
2675	zi_snow(:,il,:) = snowdz_pft(:,1,:) / 2.
2676
2677	zf_snow(:,il,:) = snowdz_pft(:,1,:)
2678
2679	END WHERE
2680	ENDIF
2681
2682	IF ( il .GT. 1 ) THEN
2683	WHERE ( veget_mask_2d(:,:) )
2684
2685	zi_snow(:,il,:) = zf_snow(:,il-1,:) + snowdz_pft(:,il,:) / 2
2686
2687	zf_snow(:,il,:) = SUM(snowdz_pft(:,1:il,:),2)
2688
2689	END WHERE
2690	ENDIF
2691
2692	ENDDO
2693
2694	DO ix = 1, kjpindex
2695	DO il = 1, nsnow
2696	zi_snow_nopftdim(ix,il) = SUM(zi_snow(ix,il,:)*veget_max(ix,:))
2697	zf_snow_nopftdim(ix,il) = SUM(zf_snow(ix,il,:)*veget_max(ix,:))
2698	END DO
2699	END DO
2700
2701	END SUBROUTINE snowlevels
2702
2703	!!
2704	!================================================================================================================================
2705	!! SUBROUTINE : snow_interpol
2706	!!
2707	!>\BRIEF This routine interpolates oxygen and methane into snow layers
2708	!!
2709	!! DESCRIPTION :
2710	!!
2711	!! RECENT CHANGE(S) : None
2712	!!
2713	!! MAIN OUTPUT VARIABLE(S) :
2714	!!
2715	!! REFERENCE(S) : None
2716	!!
2717	!! FLOWCHART11 : None
2718	!! \n
2719	!_
2720	!================================================================================================================================
2721
2722	SUBROUTINE snow_interpol (kjpindex,snowO2, snowCH4, zi_snow, zf_snow, veget_max, snowdz)
2723
2724	!! 0. Variable and parameter declaration
2725
2726	!! 0.1 Input variables
2727
2728	INTEGER(i_std), INTENT(in) :: kjpindex
2729	REAL(r_std), DIMENSION(kjpindex,nsnow), INTENT(in) :: snowdz !! snow depth at each layer
2730	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: veget_max !! maximum vegetation fraction
2731
2732	!! 0.2 Output variables
2733
2734	!! 0.3 Modified variables
2735
2736	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: snowO2 !! snow oxygen (g O2/m**3 air)
2737	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: snowCH4 !! snow methane (g CH4/m**3 air), needed just for num. scheme
2738	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm), INTENT(inout) :: zf_snow !! depths at full levels
2739	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm), INTENT(inout) :: zi_snow !! depths at intermediate levels
2740
2741	!! 0.4 Local variables
2742	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: isnow !! index of first old layer that is deeper
2743	INTEGER(i_std), DIMENSION(kjpindex,nsnow,nvm) :: i1,i2 !! indices of the layers used for the inter- or extrapolation
2744	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: snowO2o !! initial snow oxygen (g O2/m**3 air)
2745	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: snowCH4o !! initial snow methane (g CH4/m**3 air)
2746	REAL(r_std), DIMENSION(kjpindex,nvm) :: dzio !! initial distance between two levels
2747	INTEGER(i_std) :: il, it, ip, ill, iv !! indices
2748	REAL(r_std), DIMENSION(kjpindex,0:nsnow,nvm) :: zfo !! initial depths at full levels
2749	REAL(r_std), DIMENSION(kjpindex,nsnow,nvm) :: zio !! initial depths at intermediate levels
2750
2751
2752
2753	! 1. save old discretisation and temperatures
2754
2755	zio(:,:,:) = zi_snow(:,:,:)
2756
2757	zfo(:,:,:) = zf_snow(:,:,:)
2758
2759	snowO2o(:,:,:) = snowO2(:,:,:)
2760	snowCH4o(:,:,:) = snowCH4(:,:,:)
2761
2762	! 2. new discretisation
2763
2764	CALL snowlevels( kjpindex, snowdz, zi_snow, zf_snow, veget_max)
2765
2766	! 3. for each new intermediate layer, look for the first old intermediate
2767	! layer that is deeper
2768
2769	DO il = 1, nsnow
2770
2771	isnow(:,il,:) = -1
2772
2773	DO ill = nsnow,1,-1
2774
2775	WHERE ( zio(:,ill,:) .GT. zi_snow(:,il,:) .AND. veget_mask_2d(:,:) )
2776
2777	isnow(:,il,:) = ill
2778
2779	ENDWHERE
2780
2781	ENDDO
2782
2783	ENDDO
2784
2785	! 4. determine which levels to take for the inter- or extrapolation
2786
2787
2788	DO ip = 1, kjpindex
2789	DO iv = 1, nvm
2790	IF ( veget_mask_2d(ip,iv) ) THEN
2791	DO il = 1, nsnow
2792	!
2793	IF ( isnow(ip,il,iv) .EQ. 1 ) THEN
2794	!
2795	! 4.1 first old layer is below new layer:
2796	! extrapolation from layers 1 and 2
2797	!
2798	i1(ip,il,iv) = 1
2799	i2(ip,il,iv) = 2
2800	!
2801	ELSEIF ( isnow(ip,il,iv) .EQ. -1 ) THEN
2802	!
2803	! 4.2 new layer is below last old layer:
2804	! extrapolation from layers nsnow-1 and nsnow
2805	!
2806	i1(ip,il,iv) = nsnow-1
2807	i2(ip,il,iv) = nsnow
2808	!
2809	ELSE
2810	!
2811	! 4.3 new layer is between two old layers: interpolation
2812	!
2813	i1(ip,il,iv) = isnow(ip,il,iv)-1
2814	i2(ip,il,iv) = isnow(ip,il,iv)
2815	!
2816	ENDIF
2817
2818	ENDDO
2819	ENDIF
2820	ENDDO
2821	ENDDO
2822
2823	! 5. inter- or extrapolate
2824
2825	DO ip = 1, kjpindex
2826	DO iv = 1, nvm
2827	IF ( veget_mask_2d(ip,iv) ) THEN
2828	DO il = 1, nsnow
2829	dzio(ip,iv) = zio(ip,i2(ip,il,iv),iv) - zio(ip,i1(ip,il,iv),iv)
2830
2831	IF ( dzio(ip,iv) .GT. min_stomate ) THEN
2832
2833	snowO2(ip,il,iv) = snowO2o(ip,i1(ip,il,iv),iv) + &
2834	( zi_snow(ip,il,iv) - zio(ip,i1(ip,il,iv),iv) ) / dzio(ip,iv) * &
2835	( snowO2o(ip,i2(ip,il,iv),iv) - snowO2o(ip,i1(ip,il,iv),iv) )
2836	snowCH4(ip,il,iv) = snowCH4o(ip,i1(ip,il,iv),iv) + &
2837	( zi_snow(ip,il,iv) - zio(ip,i1(ip,il,iv),iv) ) / dzio(ip,iv) * &
2838	( snowCH4o(ip,i2(ip,il,iv),iv) - snowCH4o(ip,i1(ip,il,iv),iv) )
2839
2840	ELSE
2841
2842	snowO2(ip,il,iv) = snowO2o(ip,i1(ip,il,iv),iv)
2843	snowCH4(ip,il,iv) = snowCH4o(ip,i1(ip,il,iv),iv)
2844
2845	ENDIF
2846
2847	ENDDO
2848	ENDIF
2849	ENDDO
2850
2851	ENDDO
2852	END SUBROUTINE snow_interpol
2853
2854	!!
2855	!================================================================================================================================
2856	!! SUBROUTINE : permafrost_carbon_clear
2857	!!
2858	!>\BRIEF
2859	!!
2860	!! DESCRIPTION :
2861	!!
2862	!! RECENT CHANGE(S) : None
2863	!!
2864	!! MAIN OUTPUT VARIABLE(S) :
2865	!!
2866	!! REFERENCE(S) : None
2867	!!
2868	!! FLOWCHART11 : None
2869	!! \n
2870	!_
2871	!================================================================================================================================
2872	SUBROUTINE permafrost_carbon_clear()
2873	IF (ALLOCATED(veget_mask_2d)) DEALLOCATE(veget_mask_2d)
2874	IF (ALLOCATED(heights_snow)) DEALLOCATE(heights_snow)
2875	IF (ALLOCATED(zf_soil)) DEALLOCATE(zf_soil)
2876	IF (ALLOCATED(zi_soil)) DEALLOCATE(zi_soil)
2877	IF (ALLOCATED(zf_snow)) DEALLOCATE(zf_snow)
2878	IF (ALLOCATED(zi_snow)) DEALLOCATE(zi_snow)
2879	IF (ALLOCATED(alphaO2_soil )) DEALLOCATE(alphaO2_soil )
2880	IF (ALLOCATED(betaO2_soil )) DEALLOCATE(betaO2_soil )
2881	IF (ALLOCATED(alphaCH4_soil )) DEALLOCATE(alphaCH4_soil )
2882	IF (ALLOCATED(betaCH4_soil )) DEALLOCATE(betaCH4_soil )
2883	IF (ALLOCATED(alphaO2_snow )) DEALLOCATE(alphaO2_snow )
2884	IF (ALLOCATED(betaO2_snow )) DEALLOCATE(betaO2_snow )
2885	IF (ALLOCATED(alphaCH4_snow )) DEALLOCATE(alphaCH4_snow )
2886	IF (ALLOCATED(betaCH4_snow )) DEALLOCATE(betaCH4_snow )
2887	IF (ALLOCATED(zf_coeff_snow )) DEALLOCATE(zf_coeff_snow )
2888	IF (ALLOCATED(zi_coeff_snow )) DEALLOCATE(zi_coeff_snow )
2889	IF (ALLOCATED(mu_snow )) DEALLOCATE(mu_snow )
2890	IF (ALLOCATED(deepc_pftmean )) DEALLOCATE(deepc_pftmean )
2891
2892	END SUBROUTINE permafrost_carbon_clear
2893
2894	!!
2895	!================================================================================================================================
2896	!! SUBROUTINE : initialize_yedoma_carbonstocks
2897	!!
2898	!>\BRIEF This routine intialize soil carbon in yedoma region
2899	!!
2900	!! DESCRIPTION :
2901	!!
2902	!! RECENT CHANGE(S) : None
2903	!!
2904	!! MAIN OUTPUT VARIABLE(S) :
2905	!!
2906	!! REFERENCE(S) : None
2907	!!
2908	!! FLOWCHART11 : None
2909	!! \n
2910	!_
2911	!================================================================================================================================
2912
2913	SUBROUTINE initialize_yedoma_carbonstocks(kjpindex, lalo, soilc_a, soilc_s, soilc_p, zz_deep, &
2914	yedoma_map_filename, yedoma_depth, yedoma_cinit_act, yedoma_cinit_slo, yedoma_cinit_pas, altmax_ind)
2915
2916	!! 0. Variable and parameter declaration
2917
2918	!! 0.1 Input variables
2919
2920	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
2921	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo !! geographic lat/lon
2922	REAL(r_std), DIMENSION(ndeep), INTENT (in) :: zz_deep !! deep vertical profile
2923	CHARACTER(LEN=80), INTENT (in) :: yedoma_map_filename !! yedoma map
2924	REAL(r_std), INTENT(in) :: yedoma_depth !! depth of yedoma carbon stock
2925	REAL(r_std), INTENT(in) :: yedoma_cinit_act !! initial active soil C concentration
2926	REAL(r_std), INTENT(in) :: yedoma_cinit_slo !! initial slow soil C concentration
2927	REAL(r_std), INTENT(in) :: yedoma_cinit_pas !! initial passive soil C concentration
2928	INTEGER(i_std), DIMENSION(kjpindex,nvm),INTENT(in) :: altmax_ind !! Maximum over the year active-layer index
2929
2930	!! 0.2 Output variables
2931
2932	!! 0.3 Modified variables
2933
2934	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: soilc_a !! active soil C concentration
2935	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: soilc_s !! slow soil C concentration
2936	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: soilc_p !! passive soil C concentration
2937
2938	!! 0.4 Local variables
2939	REAL(r_std), DIMENSION(kjpindex) :: yedoma
2940	INTEGER(i_std) :: il, ils, ip, ix, iy, imin, jmin, ier, iv
2941	REAL(r_std) :: dlon, dlonmin, dlat, dlatmin
2942	INTEGER(i_std) :: iml, jml, lml, tml, fid
2943	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: xx,yy, yedoma_file
2944	REAL(r_std),ALLOCATABLE,DIMENSION(:) :: x,y
2945	REAL(r_std) :: lev(1), date, dt
2946	INTEGER(i_std) :: itau(1)
2947	INTEGER(i_std) :: yedoma_depth_index, iz
2948
2949	! plus bas, on prend la temperature lue dans un fichier climato si celui-ci existe
2950
2951	IF ( yedoma_map_filename .EQ. "NONE" ) THEN
2952	yedoma(:) = zero
2953	ELSE IF ( yedoma_map_filename .EQ. "EVERYWHERE" ) THEN
2954	yedoma(:) = 1.
2955	ELSE
2956	CALL flininfo(yedoma_map_filename,iml, jml, lml, tml, fid)
2957
2958	ALLOCATE (yy(iml,jml),stat=ier)
2959	IF (ier.NE.0) THEN
2960	WRITE (numout,*) ' error in yy allocation. We stop. We need ',iml,' fois ',jml,' words = '&
2961	& , iml*jml
2962	STOP 'deep_carbcycle'
2963	END IF
2964
2965	ALLOCATE (xx(iml,jml),stat=ier)
2966	IF (ier.NE.0) THEN
2967	WRITE (numout,*) ' error in xx allocation. We stop. We need ',iml,'fois ',jml,' words = '&
2968	& , iml*jml
2969	STOP 'deep_carbcycle'
2970	END IF
2971
2972	ALLOCATE (x(iml),stat=ier)
2973	IF (ier.NE.0) THEN
2974	WRITE (numout,*) ' error in x allocation. We stop. We need',iml,' words = '&
2975	& , iml
2976	STOP 'deep_carbcycle'
2977	END IF
2978
2979	ALLOCATE (y(jml),stat=ier)
2980	IF (ier.NE.0) THEN
2981	WRITE (numout,*) ' error in y allocation. We stop. We need',jml,'words = '&
2982	& , jml
2983	STOP 'deep_carbcycle'
2984	END IF
2985
2986	ALLOCATE (yedoma_file(iml,jml),stat=ier)
2987	IF (ier.NE.0) THEN
2988	WRITE (numout,*) ' error in yedoma_file allocation. We stop. We need ',iml,'fois ',jml,' words = '&
2989	& , iml*jml
2990	STOP 'deep_carbcycle'
2991	END IF
2992
2993	CALL flinopen (yedoma_map_filename, .FALSE., iml, jml, lml, &
2994	xx, yy, lev, tml, itau, date, dt, fid)
2995	CALL flinget (fid, 'yedoma', iml, jml, lml, tml, &
2996	1, 1, yedoma_file)
2997	CALL flinclo (fid)
2998	! On suppose que le fichier est regulier.
2999	! Si ce n'est pas le cas, tant pis. Les temperatures seront mal
3000	! initialisees et puis voila. De toute maniere, il faut avoir
3001	! l'esprit mal tourne pour avoir l'idee de faire un fichier de
3002	! climatologie avec une grille non reguliere.
3003	x(:) = xx(:,1)
3004	y(:) = yy(1,:)
3005	! prendre la valeur la plus proche
3006	DO ip = 1, kjpindex
3007	dlonmin = HUGE(1.)
3008	DO ix = 1,iml
3009	dlon = MIN( ABS(lalo(ip,2)-x(ix)), ABS(lalo(ip,2)+360.-x(ix)), ABS(lalo(ip,2)-360.-x(ix)) )
3010	IF ( dlon .LT. dlonmin ) THEN
3011	imin = ix
3012	dlonmin = dlon
3013	ENDIF
3014	ENDDO
3015	dlatmin = HUGE(1.)
3016	DO iy = 1,jml
3017	dlat = ABS(lalo(ip,1)-y(iy))
3018	IF ( dlat .LT. dlatmin ) THEN
3019	jmin = iy
3020	dlatmin = dlat
3021	ENDIF
3022	ENDDO
3023	yedoma(ip) = yedoma_file(imin,jmin)
3024	ENDDO
3025	DEALLOCATE (yy)
3026	DEALLOCATE (xx)
3027	DEALLOCATE (x)
3028	DEALLOCATE (y)
3029	DEALLOCATE (yedoma_file)
3030	ENDIF
3031
3032	yedoma_depth_index = 0
3033	DO iz = 1, ndeep
3034	IF (zz_deep(iz) .LE. yedoma_depth ) yedoma_depth_index = yedoma_depth_index + 1
3035	END DO
3036	WRITE(,) 'yedoma_depth_index ', yedoma_depth_index, ' at depth ', yedoma_depth
3037
3038	IF ( yedoma_depth_index .GT. 0) THEN
3039	DO ix = 1, kjpindex
3040	DO iv = 2, nvm !!! no yedoma carbon for PFT zero.
3041	IF ( veget_mask_2d(ix,iv) ) THEN
3042	DO iz = 1, yedoma_depth_index
3043	IF (yedoma(ix) .GT. 0.) THEN
3044	IF ( iz .GE. altmax_ind(ix,iv) ) THEN !!! only put yedoma carbon at base of and below the active layer
3045	soilc_a(ix, iz,iv) = yedoma_cinit_act
3046	soilc_s(ix, iz,iv) = yedoma_cinit_slo
3047	soilc_p(ix, iz,iv) = yedoma_cinit_pas
3048	ELSE
3049	soilc_a(ix, iz,iv) = zero
3050	soilc_s(ix, iz,iv) = zero
3051	soilc_p(ix, iz,iv) = zero
3052	ENDIF
3053	ELSE
3054	soilc_a(ix, iz,iv) = zero
3055	soilc_s(ix, iz,iv) = zero
3056	soilc_p(ix, iz,iv) = zero
3057	END IF
3058	END DO
3059	ENDIF
3060	ENDDO
3061	ENDDO
3062	ENDIF
3063
3064	END SUBROUTINE initialize_yedoma_carbonstocks
3065	!!
3066	!================================================================================================================================
3067	!! SUBROUTINE : carbinput
3068	!!
3069	!>\BRIEF This routine calculate carbon input to the soil
3070	!!
3071	!! DESCRIPTION :
3072	!!
3073	!! RECENT CHANGE(S) : None
3074	!!
3075	!! MAIN OUTPUT VARIABLE(S) :
3076	!!
3077	!! REFERENCE(S) : None
3078	!!
3079	!! FLOWCHART11 : None
3080	!! \n
3081	!_
3082	!================================================================================================================================
3083	SUBROUTINE carbinput(kjpindex,time_step,time,no_pfrost_decomp,tprof,tsurf,hslong, dayno,z_root,altmax, &
3084	soilc_a, soilc_s, soilc_p, soilc_in, dc_litter_z, z_organic, veget_max, rprof)
3085
3086	!! 0. Variable and parameter declaration
3087
3088	!! 0.1 Input variables
3089
3090	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
3091	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
3092	REAL(r_std), INTENT(in) :: time
3093	LOGICAL, INTENT(in) :: no_pfrost_decomp
3094	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: tprof !! Soil temperature (K)
3095	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: tsurf !! Surface temperature (K)
3096	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: hslong !! deep soil humidity
3097	INTEGER(i_std), INTENT(in) :: dayno !! current day of year
3098	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: z_root !! the rooting depth
3099	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: altmax !! Maximum over the year active-layer thickness
3100	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: z_organic !! depth to organic soil
3101	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in) :: veget_max !! Maximum fraction of vegetation type
3102	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)
3103
3104	!! 0.2 Output variables
3105
3106	REAL(r_std), DIMENSION(kjpindex,ncarb,ndeep,nvm), INTENT(out) :: dc_litter_z !! depth_dependent carbon input due to litter
3107
3108	!! 0.3 Modified variables
3109
3110	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: soilc_a !! active soil C
3111	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: soilc_s !! slow soil C
3112	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: soilc_p !! passive soil C
3113
3114	!! 0.4 Local variables
3115
3116	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm) :: dc_litter !! depth-integrated carbon input due to litter decomposition
3117	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm) :: soilc_in_finite
3118	REAL(r_std), DIMENSION(kjpindex,nvm) :: intdep !! integral depth of carbon deposition
3119	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm) :: carbinp_correction
3120	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm) :: soilc_in_TS
3121	LOGICAL, SAVE :: firstcall = .TRUE.
3122	!$OMP THREADPRIVATE(firstcall)
3123	REAL(r_std), DIMENSION(kjpindex,nvm) :: z_lit !! litter input e-folding depth
3124	INTEGER :: il,ic,iv, ix
3125	LOGICAL, SAVE :: check = .FALSE.
3126	!$OMP THREADPRIVATE(check)
3127	REAL(r_std), PARAMETER :: dgyrst = 96.
3128	INTEGER(i_std), SAVE :: id, id2, id3, id4
3129	!$OMP THREADPRIVATE(id)
3130	!$OMP THREADPRIVATE(id2)
3131	!$OMP THREADPRIVATE(id3)
3132	!$OMP THREADPRIVATE(id4)
3133	CHARACTER(LEN=16) :: buf
3134	INTEGER :: recn
3135	LOGICAL, SAVE :: correct_carboninput_vertprof = .TRUE.
3136	!$OMP THREADPRIVATE(correct_carboninput_vertprof)
3137	LOGICAL, SAVE :: new_carbinput_intdepzlit = .FALSE.
3138	!$OMP THREADPRIVATE(new_carbinput_intdepzlit)
3139	REAL(r_std), DIMENSION(ndeep), SAVE :: z_thickness
3140	!$OMP THREADPRIVATE(z_thickness)
3141	REAL(r_std), DIMENSION(ndeep) :: root_prof
3142	REAL(r_std), SAVE :: minaltmax = 0.1
3143	!$OMP THREADPRIVATE(minaltmax)
3144	REAL(r_std), SAVE :: maxaltmax = 2.
3145	!$OMP THREADPRIVATE(maxaltmax)
3146	REAL(r_std), SAVE :: finerootdepthratio = 0.5 !! the ratio of fine root to overall root e-folding depth (for C inputs)
3147	!$OMP THREADPRIVATE(finerootdepthratio)
3148	REAL(r_std), SAVE :: altrootratio = 0.5 !! the maximum ratio of fine root depth to active layer thickness (for C inputs)
3149	!$OMP THREADPRIVATE(altrootratio)
3150	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: rprof !! root depth (m)
3151	INTEGER, save :: tcounter
3152
3153
3154
3155	IF (no_pfrost_decomp) THEN
3156	!
3157	! no carbon input during spinup
3158	!
3159	dc_litter(:,:,:) = 0.
3160	!
3161	ELSE
3162	!
3163	IF (firstcall) THEN
3164
3165	DO il = 1, ndeep
3166	z_thickness(il) = zf_soil(il) - zf_soil(il-1)
3167	END DO
3168
3169	!
3170	!Config Key = new_carbinput_intdepzlit
3171	!Config Desc =
3172	!Config Def = n
3173	!Config If = OK_PC
3174	!Config Help =
3175	!Config Units = [flag]
3176	CALL getin_p('new_carbinput_intdepzlit', new_carbinput_intdepzlit)
3177
3178	!
3179	!Config Key = correct_carboninput_vertprof
3180	!Config Desc =
3181	!Config Def = n
3182	!Config If = OK_PC
3183	!Config Help =
3184	!Config Units = [flag]
3185	CALL getin_p('correct_carboninput_vertprof', correct_carboninput_vertprof)
3186
3187
3188	! Diagnostic output init
3189
3190	IF (check) THEN
3191	tcounter = 1
3192	WRITE(buf,'(I3)') yr_len
3193	id2 = 0
3194	CALL fliocrfd ('alt.nc', (/'geo ','veg ','time'/), (/kjpindex, nvm, -1/), id, id2, 'REPLACE')
3195	CALL fliodefv (id,'time',(/ 3 /),units='seconds since 0000-01-01 00:00:00',v_t=flio_r8)
3196	CALL flioputa (id,'time','title','time')
3197	CALL flioputa (id,'time','calendar',TRIM(buf)//'d')
3198	CALL fliodefv (id,'alt',(/ 1,2,3 /),units='m',v_t=flio_r8)
3199
3200	CALL fliocrfd ('soilc_litterinput.nc', (/'geo ','carb','veg ','time'/), (/kjpindex,ncarb,nvm,-1/), id3, id4, 'REPLACE')
3201	CALL fliodefv (id3,'time',(/ 4 /),units='seconds since 0000-01-01 00:00:00',v_t=flio_r8)
3202	CALL flioputa (id3,'time','title','time')
3203	CALL flioputa (id3,'time','calendar',TRIM(buf)//'d')
3204	CALL fliodefv (id3,'dc_litter',(/ 1,2,3,4 /),units='g C / ts',v_t=flio_r8)
3205	CALL fliodefv (id3,'soilc_in_TS',(/ 1,2,3,4 /),units='g C / ts',v_t=flio_r8)
3206
3207
3208	ENDIF ! check
3209
3210	firstcall = .FALSE.
3211	!
3212	ENDIF ! firstcall
3213
3214	!
3215	! 1. Litter input and decomposition
3216	!
3217	! 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
3218	soilc_in_TS(:,:,:) = soilc_in(:,:,:)*time_step/one_day
3219
3220
3221	! 2. Carbon input e-folding depth. We distribute with e-depth = min(z_root,intdep)
3222	! and integral depth = min(altmax,z_org)
3223	! ! e-folding depth cannot be greater than integral depth
3224
3225	! change to make intdep equal to z_root alone
3226	IF ( .NOT. new_carbinput_intdepzlit ) THEN
3227	z_lit(:,:) = z_root(:,:)
3228	intdep(:,:) = z_root(:,:)
3229	ELSE
3230	!change to separate e-folding depths for roots from total depth over which to integrate
3231	z_lit(:,:) = MIN(rprof(:,:)finerootdepthratio, altmax(:,:)altrootratio) ! z_lit is the e-folding depth
3232	intdep(:,:) = MIN(altmax(:,:), maxaltmax) ! intdep is the maximum depth of integration;
3233	ENDIF
3234
3235	! Litter is decomposed somehow (?) even when alt == 0. To avoid carbon loss,
3236	! we distribute this carbon within the first 2 soil layers when alt == 0
3237	WHERE ( intdep(:,:) .LT. zi_soil(2) ) intdep(:,:) = zi_soil(2) +EPSILON(0.)
3238	WHERE ( z_lit(:,:) .LT. zi_soil(2) ) z_lit(:,:) = zi_soil(2)
3239
3240	!
3241	! 3. Carbon input.
3242	!
3243	dc_litter_z(:,:,:,:) = zero
3244
3245	dc_litter(:,:,:)=zero
3246
3247
3248	DO il = 1, ndeep
3249	DO ic = 1, ncarb
3250
3251	! 3.1. from litter.
3252
3253	WHERE ( zi_soil(il) .LT. intdep(:,:) .AND. veget_mask_2d(:,:) )
3254	dc_litter_z(:,ic,il,:) = soilc_in_TS(:,ic,:) / z_lit(:,:) / ( 1. - EXP( -intdep(:,:) / z_lit(:,:) ) ) &
3255	* EXP( -zi_soil(il) / z_lit(:,:) )
3256	ELSEWHERE
3257	dc_litter_z(:,ic,il,:) = zero
3258	ENDWHERE
3259
3260	dc_litter(:,ic,:) = dc_litter(:,ic,:) + dc_litter_z(:,ic,il,:) * (zf_soil(il)-zf_soil(il-1))
3261
3262	ENDDO
3263
3264	ENDDO
3265
3266
3267	IF ( correct_carboninput_vertprof ) THEN
3268	!! correct for the truncated carbon adddition profile here by multiplying by a scalar
3269	DO ic = 1, ncarb
3270	WHERE ( dc_litter(:,ic,:) .GT. EPSILON(0.) .AND. veget_mask_2d(:,:) )
3271	carbinp_correction(:,ic,:) = soilc_in_TS(:,ic,:)/dc_litter(:,ic,:)
3272	ELSEWHERE
3273	carbinp_correction(:,ic,:) = 0.
3274	END WHERE
3275	END DO
3276
3277	dc_litter(:,:,:)=0.
3278	DO ic = 1, ncarb
3279	DO il = 1, ndeep
3280	WHERE ( veget_mask_2d(:,:) )
3281	dc_litter_z(:,ic,il,:) = carbinp_correction(:,ic,:)*dc_litter_z(:,ic,il,:)
3282	END WHERE
3283	dc_litter(:,ic,:) = dc_litter(:,ic,:) + dc_litter_z(:,ic,il,:) * (zf_soil(il)-zf_soil(il-1)) !! check again
3284	END DO
3285	END DO
3286
3287
3288	ENDIF
3289
3290	DO il = 1, ndeep
3291	WHERE ( veget_mask_2d(:,:) )
3292	soilc_a(:,il,:) = soilc_a(:,il,:) + dc_litter_z(:,iactive,il,:)
3293	soilc_s(:,il,:) = soilc_s(:,il,:) + dc_litter_z(:,islow,il,:)
3294	soilc_p(:,il,:) = soilc_p(:,il,:) + dc_litter_z(:,ipassive,il,:)
3295	END WHERE
3296	END DO
3297
3298	! Diagnostic output
3299
3300	IF (check) THEN
3301	recn = NINT(time/time_step)
3302	tcounter = tcounter + 1
3303	WRITE(,) 'carbinput check: output to .nc number',recn
3304	WRITE(,) 'time',time
3305	WRITE(,) 'time_step',time_step
3306
3307	CALL flioputv (id,'time', time, (/ tcounter /) )
3308	CALL flioputv (id,'alt', altmax(:,:), start = (/ 1, 1, tcounter /), count = (/ kjpindex, nvm, 1 /) )
3309	CALL fliosync(id)
3310
3311	CALL flioputv (id3,'time', time, (/ tcounter /) )
3312	CALL flioputv (id3,'soilc_in_TS', soilc_in_TS(:,:,:), start = (/ 1, 1, 1, tcounter /), &
3313	count = (/ kjpindex, ncarb, nvm, 1 /) )
3314	CALL flioputv (id3,'dc_litter', dc_litter(:,:,:), start = (/ 1, 1, 1, tcounter /), &
3315	count = (/ kjpindex, ncarb, nvm, 1 /) )
3316	CALL fliosync(id3)
3317	ENDIF
3318
3319	ENDIF
3320
3321	END SUBROUTINE carbinput
3322
3323	!!
3324	!================================================================================================================================
3325	!! SUBROUTINE : cryoturbate
3326	!!
3327	!>\BRIEF This routine calculates cryoturbation process
3328	!!
3329	!! DESCRIPTION :
3330	!!
3331	!! RECENT CHANGE(S) : None
3332	!!
3333	!! MAIN OUTPUT VARIABLE(S) :
3334	!!
3335	!! REFERENCE(S) : None
3336	!!
3337	!! FLOWCHART11 : None
3338	!! \n
3339	!_
3340	!================================================================================================================================
3341
3342	SUBROUTINE cryoturbate(kjpindex, time_step, dayno, altmax_ind, deepC_a, deepC_s, deepC_p, &
3343	action, diff_k_const, bio_diff_k_const, altmax_lastyear, fixed_cryoturbation_depth)
3344
3345	!! 0. Variable and parameter declaration
3346
3347	!! 0.1 Input variables
3348
3349	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
3350	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
3351	INTEGER(i_std), INTENT(in) :: dayno !! number of the day in the current year
3352	INTEGER(i_std), DIMENSION(kjpindex,nvm),INTENT(in) :: altmax_ind !! Maximum over the year active-layer index
3353	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(in) :: altmax_lastyear !! Maximum over the year active-layer thickness
3354	CHARACTER(LEN=*), INTENT(in) :: action !! what to do
3355	REAL(r_std), INTENT(in) :: diff_k_const
3356	REAL(r_std), INTENT(in) :: bio_diff_k_const
3357
3358	!! 0.2 Output variables
3359
3360	!! 0.3 Modified variables
3361	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deepC_a !! soil carbon (g/m**3) active
3362	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deepC_s !! soil carbon (g/m**3) slow
3363	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deepC_p !! soil carbon (g/m**3) passive
3364	REAL(r_std), DIMENSION(kjpindex,nvm),INTENT(inout) :: fixed_cryoturbation_depth !! depth to hold cryoturbation to for fixed runs
3365
3366	!! 0.4 Local variables
3367	LOGICAL, SAVE :: firstcall = .TRUE.
3368	!$OMP THREADPRIVATE(firstcall)
3369	LOGICAL, SAVE :: use_new_cryoturbation
3370	!$OMP THREADPRIVATE(use_new_cryoturbation)
3371	INTEGER, SAVE :: cryoturbation_method
3372	!$OMP THREADPRIVATE(cryoturbation_method)
3373	REAL(r_std), DIMENSION(kjpindex,nvm) :: altC_a_old !! soil carbon (g/m**2) active integrated over active layer before cryoturbation
3374	REAL(r_std), DIMENSION(kjpindex,nvm) :: altC_s_old !! soil carbon (g/m**2) slow integrated over active layer before cryoturbation
3375	REAL(r_std), DIMENSION(kjpindex,nvm) :: altC_p_old !! soil carbon (g/m**2) passive integrated over active layer before cryoturbation
3376	REAL(r_std), DIMENSION(kjpindex,nvm) :: altC_a
3377	REAL(r_std), DIMENSION(kjpindex,nvm) :: altC_s
3378	REAL(r_std), DIMENSION(kjpindex,nvm) :: altC_p
3379	INTEGER(i_std), PARAMETER :: n_totakefrom = 3 !! how many surface layers to subtract from in mass balance
3380	REAL(r_std), DIMENSION(kjpindex,nvm) :: surfC_totake_a !! active soil carbon to subtract from surface layers to maintain mass balance (g/m**3)
3381	REAL(r_std), DIMENSION(kjpindex,nvm) :: surfC_totake_s !! slow soil carbon to subtract from surface layers to maintain mass balance (g/m**3)
3382	REAL(r_std), DIMENSION(kjpindex,nvm) :: surfC_totake_p !! passive soil carbon to subtract from surface layers to maintain mass balance (g/m**3)
3383	REAL(r_std), DIMENSION(kjpindex,nvm) :: error_a
3384	REAL(r_std), DIMENSION(kjpindex,nvm) :: error_s
3385	REAL(r_std), DIMENSION(kjpindex,nvm) :: error_p
3386	INTEGER(i_std) :: ip, il, ier, iv
3387	CHARACTER(LEN=20), SAVE :: last_action = 'not called'
3388	!$OMP THREADPRIVATE(last_action)
3389	INTEGER(i_std) :: cryoturb_date
3390	REAL(r_std), SAVE :: max_cryoturb_alt
3391	!$OMP THREADPRIVATE(max_cryoturb_alt)
3392	REAL(r_std), SAVE :: min_cryoturb_alt
3393	!$OMP THREADPRIVATE(min_cryoturb_alt)
3394	REAL(r_std), SAVE :: bioturbation_depth
3395	!$OMP THREADPRIVATE(bioturbation_depth)
3396	LOGICAL, SAVE :: reset_fixed_cryoturbation_depth = .FALSE.
3397	!$OMP THREADPRIVATE(reset_fixed_cryoturbation_depth)
3398	LOGICAL, SAVE :: use_fixed_cryoturbation_depth = .FALSE.
3399	!$OMP THREADPRIVATE(use_fixed_cryoturbation_depth)
3400	REAL(r_std), DIMENSION(kjpindex,nvm) :: cryoturbation_depth
3401
3402
3403	! 1. ensure that we do not repeat actions
3404	!
3405	IF ( action .EQ. last_action ) THEN
3406	!
3407	WRITE(,) 'CANNOT TAKE THE SAME ACTION TWICE: ',TRIM(action)
3408	STOP
3409	!
3410	ENDIF
3411
3412	IF (firstcall) THEN
3413
3414	! 2. faire les trucs du debut
3415
3416	! 2.1 allocation des variables
3417	ALLOCATE (xe_a(kjpindex,nvm),stat=ier)
3418	IF (ier.NE.0) THEN
3419	WRITE (numout,*) ' error in xe_a allocation. We stop. We need ',kjpindex,' fois ',nvm,' words = '&
3420	& , kjpindex*nvm
3421	STOP 'deep_carbcycle'
3422	END IF
3423
3424	ALLOCATE (xe_s(kjpindex,nvm),stat=ier)
3425	IF (ier.NE.0) THEN
3426	WRITE (numout,*) ' error in xe_s allocation. We stop. We need ',kjpindex,' fois ',nvm,' words = '&
3427	& , kjpindex*nvm
3428	STOP 'deep_carbcycle'
3429	END IF
3430
3431	ALLOCATE (xe_p(kjpindex,nvm),stat=ier)
3432	IF (ier.NE.0) THEN
3433	WRITE (numout,*) ' error in xe_p allocation. We stop. We need ',kjpindex,' fois ',nvm,' words = '&
3434	& , kjpindex*nvm
3435	STOP 'deep_carbcycle'
3436	END IF
3437
3438	ALLOCATE (xc_cryoturb(kjpindex,ndeep,nvm),stat=ier)
3439	IF (ier.NE.0) THEN
3440	WRITE (numout,*) ' error in xc_cryoturb allocation. We stop. We need ',kjpindex,' fois ',ndeep,' fois ',nvm,' words = '&
3441	& , kjpindexndeepnvm
3442	STOP 'deep_carbcycle'
3443	END IF
3444
3445	ALLOCATE (xd_cryoturb(kjpindex,ndeep,nvm),stat=ier)
3446	IF (ier.NE.0) THEN
3447	WRITE (numout,*) ' error in xd_cryoturb allocation. We stop. We need ',kjpindex,' fois ',ndeep,' fois ',nvm,' words = '&
3448	& , kjpindexndeepnvm
3449	STOP 'deep_carbcycle'
3450	END IF
3451
3452	ALLOCATE (alpha_a(kjpindex,ndeep,nvm),stat=ier)
3453	IF (ier.NE.0) THEN
3454	WRITE (numout,*) ' error in alpha_a allocation. We stop. We need ',kjpindex,' fois ',ndeep,' fois ',nvm,' words = '&
3455	& , kjpindexndeepnvm
3456	STOP 'deep_carbcycle'
3457	END IF
3458	alpha_a(:,:,:)=0.
3459
3460	ALLOCATE (alpha_s(kjpindex,ndeep,nvm),stat=ier)
3461	IF (ier.NE.0) THEN
3462	WRITE (numout,*) ' error in alpha_s allocation. We stop. We need ',kjpindex,' fois ',ndeep,' fois ',nvm,' words = '&
3463	& , kjpindexndeepnvm
3464	STOP 'deep_carbcycle'
3465	END IF
3466	alpha_s(:,:,:)=0.
3467
3468	ALLOCATE (alpha_p(kjpindex,ndeep,nvm),stat=ier)
3469	IF (ier.NE.0) THEN
3470	WRITE (numout,*) ' error in alpha_p allocation. We stop. We need ',kjpindex,' fois ',ndeep,' fois ',nvm,' words = '&
3471	& , kjpindexndeepnvm
3472	STOP 'deep_carbcycle'
3473	END IF
3474	alpha_p(:,:,:)=0.
3475
3476	ALLOCATE (mu_soil_rev(kjpindex,nvm),stat=ier)
3477	IF (ier.NE.0) THEN
3478	WRITE (numout,*) ' error in mu_soil_rev allocation. We stop. We need ',kjpindex,' fois ',nvm,' words = '&
3479	& , kjpindex*nvm
3480	STOP 'deep_carbcycle'
3481	END IF
3482	mu_soil_rev(:,:)=0.
3483
3484	ALLOCATE (beta_a(kjpindex,ndeep,nvm),stat=ier)
3485	IF (ier.NE.0) THEN
3486	WRITE (numout,*) ' error in beta_a allocation. We stop. We need ',kjpindex,' fois ',ndeep,' fois ',nvm,' words = '&
3487	& , kjpindexndeepnvm
3488	STOP 'deep_carbcycle'
3489	END IF
3490	beta_a(:,:,:)=0.
3491
3492	ALLOCATE (beta_s(kjpindex,ndeep,nvm),stat=ier)
3493	IF (ier.NE.0) THEN
3494	WRITE (numout,*) ' error in beta_s allocation. We stop. We need ',kjpindex,' fois ',ndeep,' fois ',nvm,' words = '&
3495	& , kjpindexndeepnvm
3496	STOP 'deep_carbcycle'
3497	END IF
3498	beta_s(:,:,:)=0.
3499
3500	ALLOCATE (beta_p(kjpindex,ndeep,nvm),stat=ier)
3501	IF (ier.NE.0) THEN
3502	WRITE (numout,*) ' error in beta_p allocation. We stop. We need ',kjpindex,' fois ',ndeep,' fois ',nvm,' words = '&
3503	& , kjpindexndeepnvm
3504	STOP 'deep_carbcycle'
3505	END IF
3506	beta_p(:,:,:)=0.
3507
3508	ALLOCATE (diff_k(kjpindex,ndeep,nvm),stat=ier)
3509	IF (ier.NE.0) THEN
3510	WRITE (numout,*) ' error in diff_k allocation. We stop. We need ',kjpindex,' fois ',ndeep,' fois ',nvm,' words = '&
3511	& , kjpindexndeepnvm
3512	STOP 'deep_carbcycle'
3513	END IF
3514
3515	ALLOCATE (cryoturb_location(kjpindex,nvm),stat=ier)
3516	IF (ier.NE.0) THEN
3517	WRITE (numout,*) ' error in cryoturb_location allocation. We stop. We need ',kjpindex,' fois ',nvm,' words = '&
3518	& , kjpindex*nvm
3519	STOP 'deep_carbcycle'
3520	END IF
3521
3522	ALLOCATE (bioturb_location(kjpindex,nvm),stat=ier)
3523	IF (ier.NE.0) THEN
3524	WRITE (numout,*) ' error in bioturb_location allocation. We stop. We need ',kjpindex,' fois ',nvm,' words = '&
3525	& , kjpindex*nvm
3526	STOP 'deep_carbcycle'
3527	END IF
3528
3529
3530	cryoturb_location(:,:) = .false.
3531	use_new_cryoturbation = .false.
3532	!
3533	!Config Key = use_new_cryoturbation
3534	!Config Desc =
3535	!Config Def = n
3536	!Config If = OK_PC
3537	!Config Help =
3538	!Config Units = [flag]
3539	CALL getin_p('use_new_cryoturbation', use_new_cryoturbation)
3540	!
3541	!Config Key = cryoturbation_method
3542	!Config Desc =
3543	!Config Def = 1
3544	!Config If = OK_PC
3545	!Config Help =
3546	!Config Units = []
3547	cryoturbation_method = 4
3548	CALL getin_p('cryoturbation_method', cryoturbation_method)
3549	!
3550	!Config Key = max_cryoturb_alt
3551	!Config Desc =
3552	!Config Def = 1
3553	!Config If = OK_PC
3554	!Config Help =
3555	!Config Units = []
3556	max_cryoturb_alt = 3.
3557	CALL getin_p('max_cryoturb_alt',max_cryoturb_alt)
3558	!
3559	!Config Key = min_cryoturb_alt
3560	!Config Desc =
3561	!Config Def = 1
3562	!Config If = OK_PC
3563	!Config Help =
3564	!Config Units = []
3565	min_cryoturb_alt = 0.01
3566	CALL getin_p('min_cryoturb_alt',min_cryoturb_alt)
3567	!
3568	!Config Key = reset_fixed_cryoturbation_depth
3569	!Config Desc =
3570	!Config Def = n
3571	!Config If = OK_PC
3572	!Config Help =
3573	!Config Units = [flag]
3574	CALL getin_p('reset_fixed_cryoturbation_depth',reset_fixed_cryoturbation_depth)
3575	IF (reset_fixed_cryoturbation_depth) THEN
3576	fixed_cryoturbation_depth = altmax_lastyear
3577	ENDIF
3578	!
3579	!Config Key = use_fixed_cryoturbation_depth
3580	!Config Desc =
3581	!Config Def = n
3582	!Config If = OK_PC
3583	!Config Help =
3584	!Config Units = [flag]
3585	CALL getin_p('use_fixed_cryoturbation_depth',use_fixed_cryoturbation_depth)
3586	bioturb_location(:,:) = .false.
3587	!
3588	!Config Key = bioturbation_depth
3589	!Config Desc = maximum bioturbation depth
3590	!Config Def = 2
3591	!Config If = ok_pc
3592	!Config Help =
3593	!Config Units = m
3594	bioturbation_depth = 2.
3595	CALL getin_p('bioturbation_depth',bioturbation_depth)
3596
3597	firstcall = .FALSE.
3598	ENDIF
3599
3600	IF ( action .EQ. 'diffuse' ) THEN
3601	! 1. calculate the total soil carbon in the active layer
3602	altC_a_old(:,:) = zero
3603	altC_s_old(:,:) = zero
3604	altC_p_old(:,:) = zero
3605	altC_a(:,:) = zero
3606	altC_s(:,:) = zero
3607	altC_p(:,:) = zero
3608
3609	DO ip = 1, kjpindex
3610	DO iv = 1, nvm
3611	IF ( cryoturb_location(ip,iv) .OR. bioturb_location(ip,iv) )THEN
3612	! 1. calculate the total soil carbon
3613	DO il = 1, ndeep
3614	altC_a_old(ip,iv) = altC_a_old(ip,iv) + deepC_a(ip,il,iv)*(zf_soil(il)-zf_soil(il-1))
3615	altC_s_old(ip,iv) = altC_s_old(ip,iv) + deepC_s(ip,il,iv)*(zf_soil(il)-zf_soil(il-1))
3616	altC_p_old(ip,iv) = altC_p_old(ip,iv) + deepC_p(ip,il,iv)*(zf_soil(il)-zf_soil(il-1))
3617	ENDDO
3618
3619	! 2. diffuse the soil carbon
3620	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)))
3621	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)))
3622	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)))
3623
3624	DO il = 2, ndeep
3625	deepC_a(ip,il,iv) = alpha_a(ip,il-1,iv)*deepC_a(ip,il-1,iv) + beta_a(ip,il-1,iv)
3626	deepC_s(ip,il,iv) = alpha_s(ip,il-1,iv)*deepC_s(ip,il-1,iv) + beta_s(ip,il-1,iv)
3627	deepC_p(ip,il,iv) = alpha_p(ip,il-1,iv)*deepC_p(ip,il-1,iv) + beta_p(ip,il-1,iv)
3628	ENDDO
3629
3630	! 3. recalculate the total soil carbon
3631	DO il = 1, ndeep
3632	altC_a(ip,iv) = altC_a(ip,iv) + deepC_a(ip,il,iv)*(zf_soil(il)-zf_soil(il-1))
3633	altC_s(ip,iv) = altC_s(ip,iv) + deepC_s(ip,il,iv)*(zf_soil(il)-zf_soil(il-1))
3634	altC_p(ip,iv) = altC_p(ip,iv) + deepC_p(ip,il,iv)*(zf_soil(il)-zf_soil(il-1))
3635	ENDDO
3636
3637
3638	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
3639	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)
3640	CALL ipslerr_p (3,'cryoturbate','','','')
3641	ENDIF
3642
3643	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
3644	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)
3645	CALL ipslerr_p (3,'cryoturbate','','','')
3646	ENDIF
3647
3648	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
3649	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)
3650	CALL ipslerr_p (3,'cryoturbate','','','')
3651	ENDIF
3652
3653	! 4. subtract the soil carbon in the top layer(s) so that the total carbon content of the active layer is conserved.
3654	! for now remove this correction term...
3655	! surfC_totake_a(ip,iv) = (altC_a(ip,iv)-altC_a_old(ip,iv))/(zf_soil(altmax_ind(ip,iv))-zf_soil(0))
3656	! surfC_totake_s(ip,iv) = (altC_s(ip,iv)-altC_s_old(ip,iv))/(zf_soil(altmax_ind(ip,iv))-zf_soil(0))
3657	! surfC_totake_p(ip,iv) = (altC_p(ip,iv)-altC_p_old(ip,iv))/(zf_soil(altmax_ind(ip,iv))-zf_soil(0))
3658	! deepC_a(ip,1:altmax_ind(ip,iv),iv) = deepC_a(ip,1:altmax_ind(ip,iv),iv) - surfC_totake_a(ip,iv)
3659	! deepC_s(ip,1:altmax_ind(ip,iv),iv) = deepC_s(ip,1:altmax_ind(ip,iv),iv) - surfC_totake_s(ip,iv)
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	!
3662	! ! if negative values appear, we don't subtract the delta-C from top layers
3663	! IF (ANY(deepC_a(ip,1:altmax_ind(ip,iv),iv) .LT. zero) ) THEN
3664	! deepC_a(ip,1:altmax_ind(ip,iv),iv)=deepC_a(ip,1:altmax_ind(ip,iv),iv)+surfC_totake_a(ip,iv)
3665	! IF (altC_a(ip,iv) .GT. zero) THEN
3666	! deepC_a(ip,:,iv)=deepC_a(ip,:,iv)*altC_a_old(ip,iv)/altC_a(ip,iv)
3667	! ENDIF
3668	! ENDIF
3669	! IF (ANY(deepC_s(ip,1:altmax_ind(ip,iv),iv) .LT. zero) ) THEN
3670	! deepC_s(ip,1:altmax_ind(ip,iv),iv)=deepC_s(ip,1:altmax_ind(ip,iv),iv)+surfC_totake_s(ip,iv)
3671	! IF (altC_s(ip,iv) .GT. zero) THEN
3672	! deepC_s(ip,:,iv)=deepC_s(ip,:,iv)*altC_s_old(ip,iv)/altC_s(ip,iv)
3673	! ENDIF
3674	! ENDIF
3675	! IF (ANY(deepC_p(ip,1:altmax_ind(ip,iv),iv) .LT. zero) ) THEN
3676	! deepC_p(ip,1:altmax_ind(ip,iv),iv)=deepC_p(ip,1:altmax_ind(ip,iv),iv)+surfC_totake_p(ip,iv)
3677	! IF (altC_p(ip,iv) .GT. zero) THEN
3678	! deepC_p(ip,:,iv)=deepC_p(ip,:,iv)*altC_p_old(ip,iv)/altC_p(ip,iv)
3679	! ENDIF
3680	! ENDIF
3681
3682	! Consistency check. Potentially add to STRICT_CHECK flag
3683	IF ( ANY(deepC_a(ip,:,iv) .LT. zero) ) THEN
3684	WRITE (numout,*) 'cryoturbate: deepC_a<0','ip=',ip,'iv=',iv,'deepC_a=',deepC_a(ip,:,iv)
3685	CALL ipslerr_p (3,'cryoturbate','','','')
3686	ENDIF
3687	IF ( ANY(deepC_s(ip,:,iv) .LT. zero) ) THEN
3688	WRITE (numout,*) 'cryoturbate: deepC_s<0','ip=',ip,'iv=',iv,'deepC_s=',deepC_s(ip,:,iv)
3689	CALL ipslerr_p (3,'cryoturbate','','','')
3690	ENDIF
3691	IF ( ANY(deepC_p(ip,:,iv) .LT. zero) ) THEN
3692	WRITE (numout,*) 'cryoturbate: deepC_p<0','ip=',ip,'iv=',iv,'deepC_p=',deepC_p(ip,:,iv)
3693	CALL ipslerr_p (3,'cryoturbate','','','')
3694	ENDIF
3695
3696	ENDIF
3697	ENDDO
3698	ENDDO
3699
3700
3701	ELSEIF ( action .EQ. 'coefficients' ) THEN
3702	IF (firstcall) THEN
3703	WRITE(,) 'error: initilaizations have to happen before coefficients calculated. we stop.'
3704	STOP
3705	ENDIF
3706
3707	cryoturb_location(:,:) = ( altmax_lastyear(:,:) .LT. max_cryoturb_alt ) &
3708	!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,
3709	.AND. ( altmax_lastyear(:,:) .GE. min_cryoturb_alt ) .AND. veget_mask_2d(:,:)
3710	IF (use_fixed_cryoturbation_depth) THEN
3711	cryoturbation_depth(:,:) = fixed_cryoturbation_depth(:,:)
3712	ELSE
3713	cryoturbation_depth(:,:) = altmax_lastyear(:,:)
3714	ENDIF
3715
3716	bioturb_location(:,:) = ( ( altmax_lastyear(:,:) .GE. max_cryoturb_alt ) .AND. veget_mask_2d(:,:) )
3717
3718	DO ip = 1, kjpindex
3719	DO iv = 1,nvm
3720	IF ( cryoturb_location(ip,iv) ) THEN
3721	!
3722	IF (use_new_cryoturbation) THEN
3723	SELECT CASE(cryoturbation_method)
3724	CASE(1)
3725	!
3726	DO il = 1, ndeep ! linear dropoff to zero between alt and 2*alt
3727	IF ( zi_soil(il) .LE. cryoturbation_depth(ip,iv) ) THEN
3728	diff_k(ip,il,iv) = diff_k_const
3729	ELSE
3730	diff_k(ip,il,iv) = diff_k_const*(un-MAX(MIN((zi_soil(il)/cryoturbation_depth(ip,iv))-un,un),zero))
3731	ENDIF
3732	END DO
3733	!
3734	CASE(2)
3735	!
3736	DO il = 1, ndeep ! exponential dropoff with e-folding distace = alt, below the active layer
3737	IF ( zi_soil(il) .LE. cryoturbation_depth(ip,iv) ) THEN
3738	diff_k(ip,il,iv) = diff_k_const
3739	ELSE
3740	diff_k(ip,il,iv) = diff_k_const*(EXP(-MAX((zi_soil(il)/cryoturbation_depth(ip,iv)-un),zero)))
3741	ENDIF
3742	END DO
3743	!
3744	CASE(3)
3745	!
3746	! exponential dropoff with e-folding distace = alt, starting at surface
3747	diff_k(ip,:,iv) = diff_k_const*(EXP(-(zi_soil(:)/cryoturbation_depth(ip,iv))))
3748	!
3749	CASE(4)
3750	!
3751	DO il = 1, ndeep ! linear dropoff to zero between alt and 3*alt
3752	IF ( zi_soil(il) .LE. cryoturbation_depth(ip,iv) ) THEN
3753	diff_k(ip,il,iv) = diff_k_const
3754	ELSE
3755	diff_k(ip,il,iv) = diff_k_const*(un-MAX(MIN((zi_soil(il)-cryoturbation_depth(ip,iv))/ &
3756	(2.*cryoturbation_depth(ip,iv)),un),zero))
3757	ENDIF
3758	IF ( zf_soil(il) .GT. max_cryoturb_alt ) THEN
3759	diff_k(ip,il,iv) = zero
3760	ENDIF
3761	END DO
3762	!
3763	IF (printlev>=3) WRITE(,) 'cryoturb method 4: ip, iv, diff_k(ip,:,iv): ', ip, iv, diff_k(ip,:,iv)
3764	CASE(5)
3765	!
3766	DO il = 1, ndeep ! linear dropoff to zero between alt and 3m
3767	IF ( zi_soil(il) .LE. cryoturbation_depth(ip,iv) ) THEN
3768	diff_k(ip,il,iv) = diff_k_const
3769	ELSE
3770	diff_k(ip,il,iv) = diff_k_const*(un-MAX(MIN((zi_soil(il)-cryoturbation_depth(ip,iv))/ &
3771	(3.-cryoturbation_depth(ip,iv)),un),zero))
3772	ENDIF
3773	END DO
3774	!
3775	IF (printlev>=3) WRITE(,) 'cryoturb method 5: ip, iv, diff_k(ip,:,iv): ', ip, iv, diff_k(ip,:,iv)
3776	END SELECT
3777
3778	ELSE ! old cryoturbation scheme
3779	!
3780	diff_k(ip,1:altmax_ind(ip,iv),iv) = diff_k_const
3781	diff_k(ip, altmax_ind(ip,iv)+1,iv) = diff_k_const/10.
3782	diff_k(ip, altmax_ind(ip,iv)+2,iv) = diff_k_const/100.
3783	diff_k(ip,(altmax_ind(ip,iv)+3):ndeep,iv) = zero
3784	ENDIF
3785	ELSE IF ( bioturb_location(ip,iv) ) THEN
3786	DO il = 1, ndeep
3787	IF ( zi_soil(il) .LE. bioturbation_depth ) THEN
3788	diff_k(ip,il,iv) = bio_diff_k_const
3789	ELSE
3790	diff_k(ip,il,iv) = zero
3791	ENDIF
3792	END DO
3793	ELSE
3794	diff_k(ip,:,iv) = zero
3795	END IF
3796	END DO
3797	END DO
3798
3799	mu_soil_rev=diff_k(:,1,:)*time_step/(zf_soil(1)-zf_soil(0))/(zi_soil(2)-zi_soil(1))
3800
3801	DO il = 1,ndeep-1
3802	WHERE ( cryoturb_location(:,:) .OR. bioturb_location(:,:) )
3803	xc_cryoturb(:,il,:) = (zf_soil(il)-zf_soil(il-1)) / time_step
3804	xd_cryoturb(:,il,:) = diff_k(:,il,:) / (zi_soil(il+1)-zi_soil(il))
3805	endwhere
3806	ENDDO
3807
3808	WHERE ( cryoturb_location(:,:) .OR. bioturb_location(:,:) )
3809	xc_cryoturb(:,ndeep,:) = (zf_soil(ndeep)-zf_soil(ndeep-1)) / time_step
3810
3811	!bottom
3812	xe_a(:,:) = xc_cryoturb(:,ndeep,:)+xd_cryoturb(:,ndeep-1,:)
3813	xe_s(:,:) = xc_cryoturb(:,ndeep,:)+xd_cryoturb(:,ndeep-1,:)
3814	xe_p(:,:) = xc_cryoturb(:,ndeep,:)+xd_cryoturb(:,ndeep-1,:)
3815	alpha_a(:,ndeep-1,:) = xd_cryoturb(:,ndeep-1,:) / xe_a(:,:)
3816	alpha_s(:,ndeep-1,:) = xd_cryoturb(:,ndeep-1,:) / xe_s(:,:)
3817	alpha_p(:,ndeep-1,:) = xd_cryoturb(:,ndeep-1,:) / xe_p(:,:)
3818	beta_a(:,ndeep-1,:) = xc_cryoturb(:,ndeep,:)*deepC_a(:,ndeep,:) / xe_a(:,:)
3819	beta_s(:,ndeep-1,:) = xc_cryoturb(:,ndeep,:)*deepC_s(:,ndeep,:) / xe_s(:,:)
3820	beta_p(:,ndeep-1,:) = xc_cryoturb(:,ndeep,:)*deepC_p(:,ndeep,:) / xe_p(:,:)
3821	END WHERE
3822
3823	!other levels
3824	DO il = ndeep-2,1,-1
3825	WHERE ( cryoturb_location(:,:) .OR. bioturb_location(:,:) )
3826	xe_a(:,:) = xc_cryoturb(:,il+1,:) + (1.-alpha_a(:,il+1,:))*xd_cryoturb(:,il+1,:) + xd_cryoturb(:,il,:)
3827	xe_s(:,:) = xc_cryoturb(:,il+1,:) + (1.-alpha_s(:,il+1,:))*xd_cryoturb(:,il+1,:) + xd_cryoturb(:,il,:)
3828	xe_p(:,:) = xc_cryoturb(:,il+1,:) + (1.-alpha_p(:,il+1,:))*xd_cryoturb(:,il+1,:) + xd_cryoturb(:,il,:)
3829	alpha_a(:,il,:) = xd_cryoturb(:,il,:) / xe_a(:,:)
3830	alpha_s(:,il,:) = xd_cryoturb(:,il,:) / xe_s(:,:)
3831	alpha_p(:,il,:) = xd_cryoturb(:,il,:) / xe_p(:,:)
3832	beta_a(:,il,:) = (xc_cryoturb(:,il+1,:)deepC_a(:,il+1,:)+xd_cryoturb(:,il+1,:)beta_a(:,il+1,:)) / xe_a(:,:)
3833	beta_s(:,il,:) = (xc_cryoturb(:,il+1,:)deepC_s(:,il+1,:)+xd_cryoturb(:,il+1,:)beta_s(:,il+1,:)) / xe_s(:,:)
3834	beta_p(:,il,:) = (xc_cryoturb(:,il+1,:)deepC_p(:,il+1,:)+xd_cryoturb(:,il+1,:)beta_p(:,il+1,:)) / xe_p(:,:)
3835	END WHERE
3836	ENDDO
3837
3838	ELSE
3839	!
3840	! do not know this action
3841	!
3842	CALL ipslerr_p(3, 'cryoturbate', 'DO NOT KNOW WHAT TO DO:', TRIM(action), '')
3843	!
3844	ENDIF
3845
3846	! keep last action in mind
3847	!
3848	last_action = action
3849
3850	END SUBROUTINE cryoturbate
3851
3852	!!
3853	!================================================================================================================================
3854	!! SUBROUTINE : permafrost_decomp
3855	!!
3856	!>\BRIEF This routine calculates carbon decomposition
3857	!! DESCRIPTION :
3858	!!
3859	!! RECENT CHANGE(S) : None
3860	!!
3861	!! MAIN OUTPUT VARIABLE(S) :
3862	!!
3863	!! REFERENCE(S) : None
3864	!!
3865	!! FLOWCHART11 : None
3866	!! \n
3867	!_
3868	!================================================================================================================================
3869
3870	SUBROUTINE permafrost_decomp (kjpindex, time_step, tprof, airvol_soil, &
3871	oxlim, tau_CH4troph, ok_methane, fbactratio, O2m, &
3872	totporO2_soil, totporCH4_soil, hslong, clay, &
3873	no_pfrost_decomp, deepC_a, deepC_s, deepC_p, deltaCH4g, deltaCH4, deltaC1_a, deltaC1_s, deltaC1_p, deltaC2, &
3874	deltaC3, O2_soil, CH4_soil, fbact_out, MG_useallCpools)
3875
3876	!! 0. Variable and parameter declaration
3877
3878	!! 0.1 Input variables
3879
3880	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
3881	REAL(r_std), INTENT(in) :: time_step !! time step in seconds
3882	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: tprof !! deep temperature profile
3883	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: airvol_soil
3884	LOGICAL, INTENT(in) :: oxlim !! O2 limitation taken into account
3885	REAL(r_std), INTENT(in) :: tau_CH4troph !! time constant of methanetrophy (s)
3886	LOGICAL, INTENT(in) :: ok_methane !! Is Methanogenesis and -trophy taken into account?
3887	REAL(r_std), INTENT(in) :: fbactratio !! time constant of methanogenesis (ratio to that of oxic)
3888	REAL(r_std), INTENT(in) :: O2m !! oxygen concentration [g/m3] below which there is anoxy
3889	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: totporO2_soil !! total O2 porosity (Tans, 1998)
3890	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: totporCH4_soil !! total CH4 porosity (Tans, 1998)
3891	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: hslong !! deep soil humidity
3892	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: clay !! clay content
3893	LOGICAL, INTENT(in) :: no_pfrost_decomp!! Whether this is a spinup run
3894	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: fbact_out
3895	LOGICAL, INTENT(in) :: MG_useallCpools !! Do we allow all three C pools to feed methanogenesis?
3896	!! 0.2 Output variables
3897
3898	!! 0.3 Modified variables
3899
3900	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deepC_a !! soil carbon (g/m**3) active
3901	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deepC_s !! soil carbon (g/m**3) slow
3902	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deepC_p !! soil carbon (g/m**3) passive
3903	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deltaCH4
3904	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deltaCH4g
3905	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deltaC1_a
3906	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deltaC1_s
3907	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deltaC1_p
3908	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deltaC2
3909	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: deltaC3
3910	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: O2_soil !! oxygen (g O2/m**3 air)
3911	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(inout) :: CH4_soil !! methane (g CH4/m**3 air)
3912
3913	!! 0.4 Local variables
3914
3915	INTEGER(i_std) :: ier
3916	REAL(r_std), DIMENSION(3,3) :: cflux !! fluxes between soil carbon reservoirs
3917	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm) :: nadd_soil !! number of moles created / m**3 of air
3918	REAL(r_std) :: fbact_a,fbact_s, fbact_p,temp
3919	REAL(r_std) :: fbactCH4_a, fbactCH4_s, fbactCH4_p
3920	REAL(r_std) :: dC,dCm
3921	REAL(r_std) :: dCH4,dCH4m,dO2
3922	INTEGER(i_std) :: il, ip, iv
3923	LOGICAL, SAVE :: firstcall = .TRUE. !! first call?
3924	!$OMP THREADPRIVATE(firstcall)
3925
3926
3927	IF (firstcall) THEN
3928
3929	ALLOCATE (fc(kjpindex,3,3,nvm),stat=ier)
3930	IF (ier.NE.0) THEN
3931	WRITE (numout,*) ' error in fc allocation. We stop. We need ',kjpindex,' fois ',3,' fois ',3,' fois ',nvm,' words = '&
3932	& , kjpindex33*nvm
3933	STOP 'deep_carbcycle'
3934	END IF
3935	ALLOCATE (fr(kjpindex,3,nvm),stat=ier)
3936	IF (ier.NE.0) THEN
3937	WRITE (numout,*) ' error in fc allocation. We stop. We need ',kjpindex,' fois ',3,' fois ',nvm,' words = '&
3938	& , kjpindex3nvm
3939	STOP 'deep_carbcycle'
3940	END IF
3941
3942	!
3943	! calculate carbon flux fractions
3944	!
3945	DO iv =1,nvm
3946	fc(:,iactive,iactive,iv) = 0.0_r_std
3947	fc(:,iactive,ipassive,iv) = 0.004_r_std
3948	fc(:,iactive,islow,iv) = 1._r_std - (.85-.68*clay(:)) - fc(:,iactive,ipassive,iv)
3949	!
3950	fc(:,islow,islow,iv) = .0_r_std
3951	fc(:,islow,iactive,iv) = .42_r_std
3952	fc(:,islow,ipassive,iv) = .03_r_std
3953	!
3954	fc(:,ipassive,ipassive,iv) = .0_r_std
3955	fc(:,ipassive,iactive,iv) = .45_r_std
3956	fc(:,ipassive,islow,iv) = .0_r_std
3957	!
3958	fr(:,:,iv) = 1._r_std-fc(:,:,iactive,iv)-fc(:,:,islow,iv)-fc(:,:,ipassive,iv)
3959	firstcall = .FALSE.
3960	END DO
3961	IF (printlev>=3) THEN
3962	DO ip = 1,kjpindex
3963	WRITE(,) 'cdk: permafrost_decomp: i, fraction respired gridcell(i) :', ip, fr(ip,:,1)
3964	END DO
3965	ENDIF
3966	ENDIF
3967
3968	!
3969	! calculate carbon consumption
3970	!
3971	nadd_soil(:,:,:) = zero
3972	cflux(:,:) = zero
3973
3974	deltaC1_a(:,:,:) = zero
3975	deltaC1_s(:,:,:) = zero
3976	deltaC1_p(:,:,:) = zero
3977	deltaCH4(:,:,:) = zero
3978	deltaCH4g(:,:,:) = zero
3979	deltaC2(:,:,:) = zero
3980	deltaC3(:,:,:) = zero
3981	DO ip = 1, kjpindex
3982	!
3983	DO iv = 1, nvm
3984	!
3985	IF ( veget_mask_2d(ip,iv) ) THEN
3986	!
3987	DO il = 1, ndeep
3988	!
3989	! 1 function that gives carbon residence time as a function of
3990	! soil temperature (in seconds)
3991	!
3992	temp = tprof(ip,il,iv) - ZeroCelsius
3993	IF (no_pfrost_decomp) THEN
3994	! no decomposition during spinup
3995	fbact_a = HUGE(1.0)
3996	ELSE
3997	fbact_a = fbact_out(ip,il,iv)
3998	fbact_a = MAX(fbact_a,time_step)
3999	ENDIF
4000	!
4001	IF ( fbact_a/HUGE(1.) .GT. .1 ) THEN
4002	fbact_s = fbact_a
4003	fbact_p = fbact_a
4004	ELSE
4005	fbact_s = fbact_a * fslow
4006	fbact_p = fbact_a * fpassive
4007	ENDIF
4008	!
4009	! methanogenesis: first guess, 10 times (fbactratio) slower than oxic
4010	! decomposition
4011	IF ( fbact_a/HUGE(1.) .GT. .1 ) THEN
4012	fbactCH4_a = fbact_a
4013	fbactCH4_s = fbact_s
4014	fbactCH4_p = fbact_p
4015	ELSE
4016	fbactCH4_a = fbact_a * fbactratio
4017	IF ( MG_useallCpools ) THEN
4018	fbactCH4_s = fbact_s * fbactratio
4019	fbactCH4_p = fbact_p * fbactratio
4020	ELSE
4021	fbactCH4_s = HUGE(1.0)
4022	fbactCH4_p = HUGE(1.0)
4023	ENDIF
4024	ENDIF
4025	!
4026	! 2 oxic decomposition: carbon and oxygen consumption
4027	!
4028	! 2.1 active
4029	!
4030	IF (oxlim) THEN
4031	dCm = O2_soil(ip,il,iv)airvol_soil(ip,il,iv)wC/wO2
4032	dC = MIN(deepC_a(ip,il,iv) * time_step/fbact_a,dCm)
4033	ELSE
4034	dC = deepC_a(ip,il,iv) * time_step/fbact_a
4035	ENDIF
4036
4037	! pour actif
4038	dC = dC * ( 1. - .75 * clay(ip) )
4039
4040	! flux vers les autres reservoirs
4041	cflux(iactive,ipassive) = fc(ip,iactive,ipassive,iv) * dC
4042	cflux(iactive,islow) = fc(ip,iactive,islow,iv) * dC
4043	!
4044	deepC_a(ip,il,iv) = deepC_a(ip,il,iv) - dC
4045	dO2 = wO2/wC * dC*fr(ip,iactive,iv) / totporO2_soil(ip,il,iv)
4046	O2_soil(ip,il,iv) = MAX( O2_soil(ip,il,iv) - dO2, zero)
4047	! keep delta C * fr in memory (generates energy)
4048	deltaC1_a(ip,il,iv) = dC*fr(ip,iactive,iv) !!this line!!!
4049	!
4050	! 2.2 slow
4051	!
4052	IF (oxlim) THEN
4053	dCm = O2_soil(ip,il,iv)airvol_soil(ip,il,iv)wC/wO2
4054	dC = MIN(deepC_s(ip,il,iv) * time_step/fbact_s,dCm)
4055	ELSE
4056	dC = deepC_s(ip,il,iv) * time_step/fbact_s
4057	ENDIF
4058	! flux vers les autres reservoirs
4059	cflux(islow,iactive) = fc(ip,islow,iactive,iv) * dC
4060	cflux(islow,ipassive) = fc(ip,islow,ipassive,iv) * dC
4061	!
4062	deepC_s(ip,il,iv) = deepC_s(ip,il,iv) - dC
4063	dO2 = wO2/wC * dC*fr(ip,islow,iv) / totporO2_soil(ip,il,iv)
4064	O2_soil(ip,il,iv) = MAX( O2_soil(ip,il,iv) - dO2, zero)
4065	! keep delta C * fr in memory (generates energy)
4066	deltaC1_s(ip,il,iv) = dC*fr(ip,islow,iv)
4067	!
4068	! 2.3 passive
4069	!
4070	IF (oxlim) THEN
4071	dCm = O2_soil(ip,il,iv)airvol_soil(ip,il,iv)wC/wO2
4072	dC = MIN(deepC_p(ip,il,iv) * time_step/fbact_p,dCm)
4073	ELSE
4074	dC = deepC_p(ip,il,iv) * time_step/fbact_p
4075	ENDIF
4076	! flux vers les autres reservoirs
4077	cflux(ipassive,iactive) = fc(ip,ipassive,iactive,iv) * dC
4078	cflux(ipassive,islow) = fc(ip,ipassive,islow,iv) * dC
4079	!
4080	deepC_p(ip,il,iv) = deepC_p(ip,il,iv) - dC
4081	dO2 = wO2/wC * dC*fr(ip,ipassive,iv) / totporO2_soil(ip,il,iv)
4082	O2_soil(ip,il,iv) = MAX( O2_soil(ip,il,iv) - dO2, zero)
4083	! keep delta C * fr in memory (generates energy)
4084	deltaC1_p(ip,il,iv) = dC*fr(ip,ipassive,iv)
4085	!
4086	!
4087	! 3 methanogenesis or methanotrophy
4088	!
4089	!
4090	IF (ok_methane) THEN
4091	!
4092	!
4093	! 3.1 active pool methanogenesis
4094	dC = deepC_a(ip,il,iv) * time_step / fbactCH4_a * EXP(-O2_soil(ip,il,iv)(1+hslong(ip,il,iv) &
4095	(BunsenO2-1.)) / O2m ) !DKtest: when commented, no ox lim for MG
4096	! pour actif
4097	dC = dC * ( 1. - .75 * clay(ip) )
4098	dCH4 = dcfr(ip,iactive,iv) wCH4/wC / totporCH4_soil(ip,il,iv)
4099	!
4100	!
4101	! flux vers les autres reservoirs
4102	cflux(iactive,ipassive)=cflux(iactive,ipassive)+fc(ip,iactive,ipassive,iv)*dC
4103	cflux(iactive,islow)=cflux(iactive,islow)+fc(ip,iactive,islow,iv)*dC
4104	!
4105	deepC_a(ip,il,iv) = deepC_a(ip,il,iv) - dC
4106	!
4107	deltaCH4g(ip,il,iv) = dCH4
4108	!
4109	CH4_soil(ip,il,iv) = CH4_soil(ip,il,iv) + dCH4
4110	! keep delta C*fr in memory (generates energy)
4111	deltaC2(ip,il,iv) = dC*fr(ip,iactive,iv)
4112	!
4113	! how many moles of gas / m**3 of air did we generate?
4114	! (methanogenesis generates 1 molecule net if we take
4115	! B -> B' + CH4 )
4116	nadd_soil(ip,il,iv) = nadd_soil(ip,il,iv) + dCH4/wCH4
4117	!
4118	!
4119	IF ( MG_useallCpools ) THEN
4120	!
4121	! 3.2 slow pool methanogenesis cdk: adding this to allow other carbon pools to participate in MG
4122	dC = deepC_s(ip,il,iv) * time_step / fbactCH4_s * EXP(-O2_soil(ip,il,iv)(1+hslong(ip,il,iv) &
4123	(BunsenO2-1.)) / O2m ) !DKtest: when commented, no ox lim for MG
4124	dCH4 = dcfr(ip,islow,iv) wCH4/wC / totporCH4_soil(ip,il,iv)
4125	!
4126	! flux vers les autres reservoirs
4127	cflux(islow,ipassive)=cflux(islow,ipassive)+fc(ip,islow,ipassive,iv)*dC
4128	cflux(islow,iactive)=cflux(islow,iactive)+fc(ip,islow,iactive,iv)*dC
4129	!
4130	deepC_s(ip,il,iv) = deepC_s(ip,il,iv) - dC
4131	!
4132	deltaCH4g(ip,il,iv) = deltaCH4g(ip,il,iv) + dCH4
4133	CH4_soil(ip,il,iv) = CH4_soil(ip,il,iv) + dCH4
4134	! keep delta C*fr in memory (generates energy)
4135	deltaC2(ip,il,iv) = deltaC2(ip,il,iv) + dC*fr(ip,islow,iv)
4136	!
4137	! how many moles of gas / m**3 of air did we generate?
4138	! (methanogenesis generates 1 molecule net if we take
4139	! B -> B' + CH4 )
4140	nadd_soil(ip,il,iv) = nadd_soil(ip,il,iv) + dCH4/wCH4
4141	!
4142	!
4143	!
4144	! 3.3 passive pool methanogenesis cdk: adding this to allow other carbon pools to participate in MG
4145	dC = deepC_p(ip,il,iv) * time_step / fbactCH4_p * EXP(-O2_soil(ip,il,iv)(1+hslong(ip,il,iv) &
4146	(BunsenO2-1.)) / O2m ) !DKtest: when commented, no ox lim for MG
4147	dCH4 = dcfr(ip,ipassive,iv) wCH4/wC / totporCH4_soil(ip,il,iv)
4148	!
4149	! flux vers les autres reservoirs
4150	cflux(ipassive,islow)=cflux(ipassive,islow)+fc(ip,ipassive,islow,iv)*dC
4151	cflux(ipassive,iactive)=cflux(ipassive,iactive)+fc(ip,ipassive,iactive,iv)*dC
4152	!
4153	deepC_p(ip,il,iv) = deepC_p(ip,il,iv) - dC
4154	!
4155	deltaCH4g(ip,il,iv) = deltaCH4g(ip,il,iv) + dCH4
4156	CH4_soil(ip,il,iv) = CH4_soil(ip,il,iv) + dCH4
4157	! keep delta C*fr in memory (generates energy)
4158	deltaC2(ip,il,iv) = deltaC2(ip,il,iv) + dC*fr(ip,ipassive,iv)
4159	!
4160	! how many moles of gas / m**3 of air did we generate?
4161	! (methanogenesis generates 1 molecule net if we take
4162	! B -> B' + CH4 )
4163	nadd_soil(ip,il,iv) = nadd_soil(ip,il,iv) + dCH4/wCH4
4164	!
4165	!
4166	ENDIF
4167	!
4168	! trophy:
4169	! no temperature dependence except that T>0ï¿œï¿œC (Price et
4170	! al, GCB 2003; Koschorrek and Conrad, GBC 1993).
4171	! tau_CH4troph is such that we fall between values of
4172	! soil methane oxidation flux given by these authors.
4173	!
4174	IF ( temp .GE. zero ) THEN
4175	!
4176	dCH4m = O2_soil(ip,il,iv)/2. * wCH4/wO2 * totporO2_soil(ip,il,iv)/totporCH4_soil(ip,il,iv)
4177	! dCH4m = CH4_soil(ip,il,iv) !DKtest - no ox lim to trophy
4178	dCH4 = MIN( CH4_soil(ip,il,iv) * time_step/MAX(tau_CH4troph,time_step), dCH4m )
4179	CH4_soil(ip,il,iv) = CH4_soil(ip,il,iv) - dCH4
4180	dO2 = 2.dCH4 wO2/wCH4 * totporCH4_soil(ip,il,iv)/totporO2_soil(ip,il,iv)
4181	O2_soil(ip,il,iv) = MAX( O2_soil(ip,il,iv) - dO2, zero)
4182	! keep delta CH4 in memory (generates energy)
4183	deltaCH4(ip,il,iv) = dCH4
4184	! carbon (g/m3 soil) transformed to CO2
4185	deltaC3(ip,il,iv)=dCH4/wCH4wCtotporCH4_soil(ip,il,iv)
4186	! how many moles of gas / m**3 of air did we generate?
4187	! (methanotrophy consumes 2 molecules net if we take
4188	! CH4 + 2 O2 -> CO2 + 2 H2O )
4189	nadd_soil(ip,il,iv) = nadd_soil(ip,il,iv)-2.*dCH4/wCH4
4190	!
4191	ENDIF
4192
4193	ENDIF
4194
4195	! 4 add fluxes between reservoirs
4196
4197	deepC_a(ip,il,iv)=deepC_a(ip,il,iv)+cflux(islow,iactive)+cflux(ipassive,iactive)
4198	deepC_s(ip,il,iv)=deepC_s(ip,il,iv)+cflux(iactive,islow)+cflux(ipassive,islow)
4199	deepC_p(ip,il,iv)=deepC_p(ip,il,iv)+cflux(iactive,ipassive)+cflux(islow,ipassive)
4200
4201	ENDDO
4202
4203	ELSE
4204
4205	ENDIF
4206
4207	ENDDO
4208
4209	ENDDO
4210	END SUBROUTINE permafrost_decomp
4211
4212
4213	!!
4214	!================================================================================================================================
4215	!! SUBROUTINE : calc_vert_int_soil_carbon
4216	!!
4217	!>\BRIEF This routine calculates carbon decomposition
4218	!!
4219	!! DESCRIPTION :
4220	!!
4221	!! RECENT CHANGE(S) : None
4222	!!
4223	!! MAIN OUTPUT VARIABLE(S) :
4224	!!
4225	!! REFERENCE(S) : None
4226	!!
4227	!! FLOWCHART11 : None
4228	!! \n
4229	!_
4230	!================================================================================================================================
4231
4232	SUBROUTINE calc_vert_int_soil_carbon(kjpindex, deepC_a, deepC_s, deepC_p, carbon, carbon_surf, zf_soil)
4233
4234	!! 0. Variable and parameter declaration
4235
4236	!! 0.1 Input variables
4237
4238	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
4239	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: deepC_a !! active pool deepc
4240	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: deepC_s !! slow pool deepc
4241	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT(in) :: deepC_p !! passive pool deepc
4242	REAL(r_std), DIMENSION(0:ndeep), INTENT(in) :: zf_soil !! depths at full levels
4243
4244	!! 0.2 Output variables
4245
4246	REAL(r_std), DIMENSION(kjpindex,ncarb,nvm), INTENT (out) :: carbon !! vertically-integrated carbon pool: active, slow, or passive, (gC/(m**2 of ground))
4247	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))
4248
4249	!! 0.3 Modified variables
4250
4251	!! 0.4 Local variables
4252	INTEGER(i_std) :: il
4253	real(r_std), parameter :: maxdepth=2.!! depth to which we intergrate the carbon for carbon_surf calculation
4254
4255	carbon(:,:,:) = zero
4256	DO il = 1, ndeep
4257	WHERE ( veget_mask_2d(:,:) )
4258	carbon(:,iactive,:) = carbon(:,iactive,:) + deepC_a(:,il,:)*(zf_soil(il)-zf_soil(il-1))
4259	carbon(:,islow,:) = carbon(:,islow,:) + deepC_s(:,il,:)*(zf_soil(il)-zf_soil(il-1))
4260	carbon(:,ipassive,:) = carbon(:,ipassive,:) + deepC_p(:,il,:)*(zf_soil(il)-zf_soil(il-1))
4261	END WHERE
4262	ENDDO
4263
4264	carbon_surf(:,:,:) = zero
4265	DO il = 1, ndeep
4266	if (zf_soil(il-1) .lt. maxdepth ) then
4267	where ( veget_mask_2d(:,:) )
4268	carbon_surf(:,iactive,:) = carbon_surf(:,iactive,:) + deepC_a(:,il,:)*(min(maxdepth,zf_soil(il))-zf_soil(il-1))
4269	carbon_surf(:,islow,:) = carbon_surf(:,islow,:) + deepC_s(:,il,:)*(min(maxdepth,zf_soil(il))-zf_soil(il-1))
4270	carbon_surf(:,ipassive,:) = carbon_surf(:,ipassive,:) + deepC_p(:,il,:)*(min(maxdepth,zf_soil(il))-zf_soil(il-1))
4271	end where
4272	endif
4273	ENDDO
4274
4275	END SUBROUTINE calc_vert_int_soil_carbon
4276
4277	END MODULE stomate_permafrost_soilcarbon

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

Context Navigation

source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_stomate/stomate_permafrost_soilcarbon.f90 @ 6940

Download in other formats: