Context Navigation

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

Last change on this file since 8367 was 7069, checked in by nicolas.vuichard, 3 years ago
corrections for avoiding warning in debug mode
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 92.5 KB

Line
1
2	! =================================================================================================================================
3	! MODULE : stomate_somdynamics
4	!
5	! CONTACT : orchidee-help _at_ listes.ipsl.fr
6	!
7	! LICENCE : IPSL (2006)
8	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
9	!
10	!>\BRIEF Calculate soil dynamics largely following the Century model
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! SVN :
17	!! $HeadURL$
18	!! $Date$
19	!! $Revision$
20	!! \n
21	!_ ================================================================================================================================
22
23	MODULE stomate_som_dynamics
24
25	! modules used:
26
27	USE ioipsl_para
28	USE stomate_data
29	USE constantes
30	USE constantes_soil
31	USE xios_orchidee
32
33	IMPLICIT NONE
34
35	! private & public routines
36
37	PRIVATE
38	PUBLIC som_dynamics,som_dynamics_clear,nitrogen_dynamics,nitrogen_dynamics_clear
39
40	! Variables shared by all subroutines in this module
41
42	LOGICAL, SAVE :: firstcall_som = .TRUE. !! Is this the first call? (true/false)
43	!$OMP THREADPRIVATE(firstcall_som)
44	LOGICAL, SAVE :: firstcall_nitrogen = .TRUE. !! Is this the first call? (true/false)
45	!$OMP THREADPRIVATE(firstcall_nitrogen)
46	! flux fractions within carbon pools
47	REAL(r_std),ALLOCATABLE,SAVE, DIMENSION(:,:,:) :: frac_carb !! Flux fractions between carbon pools
48	!! (second index=origin, third index=destination)
49	!! (unitless, 0-1)
50	!$OMP THREADPRIVATE(frac_carb)
51	REAL(r_std), ALLOCATABLE,SAVE, DIMENSION(:,:) :: frac_resp !! Flux fractions from carbon pools to the atmosphere (respiration) (unitless, 0-1)
52	!$OMP THREADPRIVATE(frac_resp)
53
54
55
56	CONTAINS
57
58
59	!! ================================================================================================================================
60	!! SUBROUTINE : som_dynamics_clear
61	!!
62	!>\BRIEF Set the flag ::firstcall_som to .TRUE. and as such activate sections 1.1.2 and 1.2 of the subroutine som_dynamics
63	!! (see below).
64	!!
65	!_ ================================================================================================================================
66
67	SUBROUTINE som_dynamics_clear
68	firstcall_som=.TRUE.
69	END SUBROUTINE som_dynamics_clear
70
71
72	!! ================================================================================================================================
73	!! SUBROUTINE : som_dynamics
74	!!
75	!>\BRIEF Computes the soil respiration and nutrient stocks, essentially
76	!! following Parton et al. (1987). Additional dynamics for the nitrogen pools are
77	!! described in the nitrogen_dynamics subroutine.
78	!!
79	!! DESCRIPTION : The soil is divided into 3 pools, with different
80	!! characteristic turnover times : active (1-5 years), slow (20-40 years)
81	!! and passive (200-1500 years).\n
82	!! There are three types of nutrient (carbon, nitrogen) transferred into the soil:\n
83	!! - input to active and slow pools from litter decomposition,\n
84	!! - nutrient fluxes between the three pools,\n
85	!! - nurtient losses from the pools to the atmosphere, i.e., soil respiration.\n
86	!!
87	!! The subroutine performs the following tasks:\n
88	!!
89	!! Section 1.\n
90	!! The flux fractions (f) between carbon pools are defined based on Parton et
91	!! al. (1987). The fractions are constants, except for the flux fraction from
92	!! the active pool to the slow pool, which depends on the clay content,\n
93	!! \latexonly
94	!! \input{soilcarbon_eq1.tex}
95	!! \endlatexonly\n
96	!! In addition, to each pool is assigned a constant turnover time. No nitrogen
97	!! is considered in this section.\n
98	!!
99	!! Section 2.\n
100	!! The carbon and nitrogen inputs, calculated in the stomate_litter module, are
101	!! added to the carbon and nitrogen stocks of the different pools.\n
102	!!
103	!! Section 3.\n
104	!! First, the outgoing flux of each pool is calculated. It is
105	!! proportional to the product of the carbon stock and the ratio between the
106	!! iteration time step and the residence time:\n
107	!! \latexonly
108	!! \input{soilcarbon_eq2.tex}
109	!! \endlatexonly
110	!! ,\n
111	!! Note that in the case of crops, the additional multiplicative factor
112	!! integrates the faster decomposition due to tillage (following Gervois et
113	!! al. (2008)).
114	!! In addition, the flux from the active pool depends on the clay content:\n
115	!! \latexonly
116	!! \input{soilcarbon_eq3.tex}
117	!! \endlatexonly
118	!! ,\n
119	!! Each pool is then cut from the carbon amount corresponding to each outgoing
120	!! flux:\n
121	!! \latexonly
122	!! \input{soilcarbon_eq4.tex}
123	!! \endlatexonly\n
124	!! Note that the total flux for both carbon and nitrogen out of the pools is
125	!! calculated first and grouped together. This total material flux is then
126	!! partitioned between the elements and the pools based on a target CN ratio.\n
127	!! Second, the flux fractions lost to the atmosphere is calculated in each pool
128	!! by subtracting from 1 the pool-to-pool flux fractions. The soil respiration
129	!! is then the summed contribution of all the pools,\n
130	!! \latexonly
131	!! \input{soilcarbon_eq5.tex}
132	!! \endlatexonly\n
133	!! Note that soil respiration only happens for the carbon pools.\n
134	!! Finally, each soil nutrient pool accumulates the contribution of the other pools:
135	!! \latexonly
136	!! \input{soilcarbon_eq6.tex}
137	!! \endlatexonly
138	!! The lefotver nitrogen flux that doesn't go into one of the four pools is
139	!! mineralized.\n
140	!! Section 4.\n
141	!! If the flag SPINUP_ANALYTIC is set to true, the matrix A is updated following
142	!! Lardy (2011).
143	!!
144	!! RECENT CHANGE(S): Merge with the Nitrogen cycle, 2018.
145	!!
146	!! MAIN OUTPUTS VARIABLE(S): carbon, resp_hetero_soil
147	!!
148	!! REFERENCE(S) :
149	!! - Parton, W.J., D.S. Schimel, C.V. Cole, and D.S. Ojima. 1987. Analysis of
150	!! factors controlling soil organic matter levels in Great Plains grasslands.
151	!! Soil Sci. Soc. Am. J., 51, 1173-1179.
152	!! - Gervois, S., P. Ciais, N. de Noblet-Ducoudre, N. Brisson, N. Vuichard,
153	!! and N. Viovy (2008), Carbon and water balance of European croplands
154	!! throughout the 20th century, Global Biogeochem. Cycles, 22, GB2022,
155	!! doi:10.1029/2007GB003018.
156	!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
157	!! Environmental Modelling & Software (2011), doi:10.1016\|j.envsoft.2011.05.016
158	!! - S. Zaehle and A. D. Friend (2010), Carbon and nitrogen cycle dynamics in the
159	!! O-CN land surface model: 1. Model description, site-scale evaluation, and
160	!! sensitivity to parameter estimates. Global Biogeochem. Cycles, 24, GB1005,
161	!! doi:10.1029/2009GB003521.
162	!!
163	!! FLOWCHART :
164	!! \latexonly
165	!! \includegraphics[scale=0.5]{soilcarbon_flowchart.jpg}
166	!! \endlatexonly
167	!! \n
168	!_ ================================================================================================================================
169
170	SUBROUTINE som_dynamics (npts, clay, silt, &
171	som_input, control_temp, control_moist, veget_cov_max, drainage_pft,&
172	CN_target,som, &
173	resp_hetero_soil, &
174	MatrixA,n_mineralisation, CN_som_litter_longterm, tau_CN_longterm)
175
176	!! 0. Variable and parameter declaration
177
178	!! 0.1 Input variables
179
180	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
181	REAL(r_std), DIMENSION(npts), INTENT(in) :: clay !! Clay fraction (unitless, 0-1)
182	REAL(r_std), DIMENSION(npts), INTENT(in) :: silt !! Silt fraction (unitless, 0-1)
183	REAL(r_std), DIMENSION(npts,ncarb,nvm,nelements), INTENT(in) :: som_input !! Amount of Organic Matter going into the SOM pools from litter decomposition \f$(gC m^{-2} day^{-1})$\f
184	REAL(r_std), DIMENSION(npts,nlevs), INTENT(in) :: control_temp !! Temperature control of heterotrophic respiration (unitless: 0->1)
185	REAL(r_std), DIMENSION(npts,nlevs), INTENT(in) :: control_moist !! Moisture control of heterotrophic respiration (unitless: 0.25->1)
186	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_cov_max !! Fractional coverage: maximum share of the pixel taken by a pft
187	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: drainage_pft !! Drainage per PFT (mm/m2 /dt_sechiba) (fraction of water content)
188	REAL(r_std), DIMENSION(npts,nvm,ncarb), INTENT(in) :: CN_target !! C to N ratio of SOM flux from one pool to another (gN m-2 dt-1)
189
190	!! 0.2 Output variables
191
192	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_hetero_soil !! Soil heterotrophic respiration \f$(gC m^{-2} (dt_sechiba one_day^{-1})^{-1})$\f
193
194	!! 0.3 Modified variables
195
196	REAL(r_std), DIMENSION(npts,ncarb,nvm,nelements), INTENT(inout) :: som !! SOM pools: active, slow, or passive, \f$(gC m^{2})$\f
197	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA !! Matrix containing the fluxes between the carbon pools
198	!! per sechiba time step
199	!! @tex $(gC.m^2.day^{-1})$ @endtex
200	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: n_mineralisation !! net nitrogen mineralisation of decomposing SOM,
201	!! (gN/m**2/day), assumed to be NH4
202	REAL(r_std), DIMENSION(npts,nvm,nbpools), INTENT(inout) :: CN_som_litter_longterm !! Longterm CN ratio of litter and som pools (gC/gN)
203	REAL(r_std), INTENT(inout) :: tau_CN_longterm !! Counter used for calculating the longterm CN ratio of SOM and litter pools (seconds)
204
205	!! 0.4 Local variables
206	REAL(r_std) :: dt !! Time step \f$(dt_sechiba one_day^{-1})$\f
207	LOGICAL :: l_error !! Diagnostic boolean for error allocation (true/false)
208	INTEGER(i_std) :: ier !! Check errors in netcdf call
209	REAL(r_std), SAVE, DIMENSION(ncarb) :: som_turn !! Residence time in SOM pools (days)
210	!$OMP THREADPRIVATE(som_turn)
211	REAL(r_std), DIMENSION(npts,ncarb,nelements) :: fluxtot !! Total flux out of carbon pools \f$(gC m^{2})$\f
212	REAL(r_std), DIMENSION(npts,ncarb,ncarb,nelements) :: flux !! Fluxes between carbon pools \f$(gC m^{2})$\f
213	CHARACTER(LEN=7), DIMENSION(ncarb) :: soilpools_str !! Name of the soil pools for informative outputs (unitless)
214	!! mineral nitrogen in the soil (gN/m**2)
215	INTEGER(i_std) :: k,kk,m,j,l, ij !! Indices (unitless)
216	REAL(r_std), DIMENSION(npts,nvm,ncarb) :: decomp_rate_soilcarbon !! Decomposition rate of the soil carbon pools (s)
217	REAL(r_std), DIMENSION(npts,ncarb) :: tsoilpools !! Diagnostic for soil carbon turnover rate by pool (1/s)
218
219	!_ ================================================================================================================================
220
221	!! printlev is the level of diagnostic information, 0 (none) to 4 (full)
222	IF (printlev>=3) WRITE(numout,*) 'Entering som_dynamics'
223
224	dt = dt_sechiba/one_day
225	IF ( firstcall_som ) THEN
226	!! 1. Initializations
227	l_error = .FALSE.
228	ALLOCATE (frac_carb(npts,ncarb,ncarb), stat=ier)
229	l_error = l_error .OR. (ier.NE.0)
230	ALLOCATE (frac_resp(npts,ncarb), stat=ier)
231	l_error = l_error .OR. (ier.NE.0)
232	IF (l_error) THEN
233	STOP 'stomate_som_dynamics: error in memory allocation'
234	ENDIF
235
236	frac_carb(:,:,:) = 0.0
237	!! 1.1 Get soil "constants"
238	!! 1.1.1 Flux fractions between carbon pools: depend on clay content, recalculated each time
239	! From active pool: depends on clay content
240	frac_carb(:,iactive,ipassive) = active_to_pass_ref_frac + active_to_pass_clay_frac*clay(:)
241	frac_carb(:,iactive,islow) = un - frac_carb(:,iactive,ipassive) - (active_to_co2_ref_frac - &
242	active_to_co2_clay_silt_frac*(clay(:)+silt(:)))
243
244	! From slow pool
245	frac_carb(:,islow,ipassive) = slow_to_pass_ref_frac + slow_to_pass_clay_frac*clay(:)
246	! OCN doesn't use Parton 1993 formulation for frac_carb(:,islow,ipassive)
247	! but the one of 1987 : ie = 0.03 ....
248	frac_carb(:,islow,iactive) = un - frac_carb(:,islow,ipassive) - slow_to_co2_ref_frac
249
250	! From passive pool
251	frac_carb(:,ipassive,iactive) = pass_to_active_ref_frac
252	frac_carb(:,ipassive,islow) = pass_to_slow_ref_frac
253
254	! From surface pool
255	frac_carb(:,isurface,islow) = surf_to_slow_ref_frac
256
257	!! 1.1.2 Determine the respiration fraction : what's left after
258	! subtracting all the 'pool-to-pool' flux fractions
259	! Diagonal elements of frac_carb are zero
260	frac_resp(:,:) = un - frac_carb(:,:,isurface) - frac_carb(:,:,iactive) - frac_carb(:,:,islow) - &
261	frac_carb(:,:,ipassive)
262
263	!! 1.1.3 Turnover in SOM pools (in days)
264	!! som_turn_ipool are the turnover (in year)
265	!! It is weighted by Temp and Humidity function later
266	som_turn(iactive) = som_turn_iactive / one_year
267	som_turn(islow) = som_turn_islow / one_year
268	som_turn(ipassive) = som_turn_ipassive / one_year
269	som_turn(isurface) = som_turn_isurface / one_year
270
271	!! 1.2 Messages : display the residence times
272	soilpools_str(iactive) = 'active'
273	soilpools_str(islow) = 'slow'
274	soilpools_str(ipassive) = 'passive'
275	soilpools_str(isurface) = 'surface'
276
277	WRITE(numout,*) 'som_dynamics:'
278
279	IF (printlev >= 2) THEN
280	WRITE(numout,*) ' > minimal SOM residence time in soil pools (d):'
281	DO k = 1, ncarb ! Loop over soil pools
282	WRITE(numout,*) '(1, ::soilpools_str(k)):',soilpools_str(k),' : (1, ::som_turn(k)):',som_turn(k)
283	WRITE(numout,*) 'FRACCARB k=',k,' ',frac_carb(1,k,isurface),frac_carb(1,k,iactive), &
284	frac_carb(1,k,islow),frac_carb(1,k,ipassive)
285	ENDDO
286
287	WRITE(numout,*) ' > flux fractions between soil pools: depend on clay content'
288	END IF
289
290	firstcall_som = .FALSE.
291
292	ENDIF
293
294	!! 1.3 Set soil respiration and decomposition rate to zero
295	resp_hetero_soil(:,:) = zero
296	decomp_rate_soilcarbon(:,:,:) = zero
297
298	!! 2. Update the SOM stocks with the different soil carbon and nitrogen input
299
300	som(:,:,:,:) = som(:,:,:,:) + som_input(:,:,:,:) * dt
301
302	!! 3. Fluxes between nutrient reservoirs, and to the atmosphere (respiration) \n
303
304
305	!! 3.2. Calculate fluxes
306
307	DO m = 1, nvm ! Loop over # PFTs
308
309	!! 3.2.1. Flux out of pools
310
311	DO k = 1, ncarb ! Loop over SOM pools from which the flux comes (active, slow, passive)
312
313	DO l = 1, nelements ! Loop over elements (Carbon, Nitrogen)
314	! Determine total flux out of pool
315	fluxtot(:,k,l) = dtsom_turn(k) som(:,k,m,l) * &
316	control_moist(:,ibelow) * control_temp(:,ibelow) * decomp_factor(m)
317
318	! Flux from active pools depends on clay content
319	IF ( k .EQ. iactive ) THEN
320	fluxtot(:,k,l) = fluxtot(:,k,l) * ( un - som_turn_iactive_clay_frac * clay(:) )
321	ENDIF
322
323	! Update the loss in each carbon pool
324	som(:,k,m,l) = som(:,k,m,l) - fluxtot(:,k,l)
325
326	! Calculate diagnostic for decomposition rate of the soil carbon pools
327	IF ( l == icarbon ) THEN
328	decomp_rate_soilcarbon(:,m,k) = dtsom_turn(k) &
329	control_moist(:,ibelow) * control_temp(:,ibelow) * decomp_factor(m)
330	IF ( k .EQ. iactive ) THEN
331	decomp_rate_soilcarbon(:,m,k) = decomp_rate_soilcarbon(:,m,k)( un - som_turn_iactive_clay_frac clay(:) )
332	ENDIF
333	END IF
334
335	ENDDO
336
337	! Fluxes towards the other pools (k -> kk)
338	DO kk = 1, ncarb ! Loop over the SOM pools where the flux goes
339	! Carbon flux
340	flux(:,k,kk,icarbon) = frac_carb(:,k,kk) * fluxtot(:,k,icarbon)
341	! Nitrogen flux - Function of the C stock of the 'departure' pool
342	! and of the C to N target ratio of the 'arrival' pool
343	flux(:,k,kk,initrogen) = frac_carb(:,k,kk) * fluxtot(:,k,icarbon) / &
344	CN_target(:,m,kk)
345	ENDDO
346	ENDDO ! End of loop over SOM pools
347
348	!! 3.2.2 respiration
349	!BE CAREFUL: Here resp_hetero_soil is divided by dt to have a value which corresponds to
350	! the sechiba time step but then in stomate.f90 resp_hetero_soil is multiplied by dt.
351	! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dtradia/one_day
352	! time step and avoid some constructions that could create bug during future developments.
353	!
354	resp_hetero_soil(:,m) = &
355	( frac_resp(:,iactive) * fluxtot(:,iactive,icarbon) + &
356	frac_resp(:,islow) * fluxtot(:,islow,icarbon) + &
357	frac_resp(:,ipassive) * fluxtot(:,ipassive,icarbon) + &
358	frac_resp(:,isurface) * fluxtot(:,isurface,icarbon) ) / dt
359
360	!! 3.2.3 add fluxes to active, slow, and passive pools
361
362	DO k = 1, ncarb ! Loop over SOM pools
363	som(:,k,m,icarbon) = som(:,k,m,icarbon) + &
364	flux(:,iactive,k,icarbon) + flux(:,ipassive,k,icarbon) + flux(:,islow,k,icarbon) + flux(:,isurface,k,icarbon)
365	som(:,k,m,initrogen) = som(:,k,m,initrogen) + &
366	flux(:,iactive,k,initrogen) + flux(:,ipassive,k,initrogen) + flux(:,islow,k,initrogen) + flux(:,isurface,k,initrogen)
367	n_mineralisation(:,m) = n_mineralisation(:,m) + fluxtot(:,k,initrogen) - &
368	(flux(:,k,iactive,initrogen)+flux(:,k,ipassive,initrogen)+&
369	flux(:,k,islow,initrogen)+flux(:,k,isurface,initrogen))
370	ENDDO ! Loop over SOM pools
371
372	ENDDO ! End loop over PFTs
373
374	!! 4. (Quasi-)Analytical Spin-up
375
376	!! 4.1.1 Finish to fill MatrixA with fluxes between soil pools
377
378	IF (spinup_analytic) THEN
379
380	DO m = 2,nvm
381
382	! flux leaving the active pool
383	MatrixA(:,m,iactive_pool,iactive_pool) = moins_un * &
384	dtsom_turn(iactive) &
385	control_moist(:,ibelow) * control_temp(:,ibelow) * &
386	( 1. - som_turn_iactive_clay_frac * clay(:)) * decomp_factor(m)
387
388	! flux received by the active pool from the slow pool
389	MatrixA(:,m,iactive_pool,islow_pool) = frac_carb(:,islow,iactive)dtsom_turn(islow) * &
390	control_moist(:,ibelow) * control_temp(:,ibelow) * decomp_factor(m)
391
392	! flux received by the active pool from the passive pool
393	MatrixA(:,m,iactive_pool,ipassive_pool) = frac_carb(:,ipassive,iactive)dtsom_turn(ipassive) * &
394	control_moist(:,ibelow) * control_temp(:,ibelow) * decomp_factor(m)
395
396	! flux leaving the slow pool
397	MatrixA(:,m,islow_pool,islow_pool) = moins_un * &
398	dtsom_turn(islow) &
399	control_moist(:,ibelow) * control_temp(:,ibelow) * decomp_factor(m)
400
401	! flux received by the slow pool from the active pool
402	MatrixA(:,m,islow_pool,iactive_pool) = frac_carb(:,iactive,islow) *&
403	dtsom_turn(iactive) &
404	control_moist(:,ibelow) * control_temp(:,ibelow) * &
405	( 1. - som_turn_iactive_clay_frac * clay(:) ) * decomp_factor(m)
406
407	! flux received by the slow pool from the surface pool
408	MatrixA(:,m,islow_pool,isurface_pool) = frac_carb(:,isurface,islow) *&
409	dtsom_turn(isurface) &
410	control_moist(:,ibelow) * control_temp(:,ibelow) * decomp_factor(m)
411
412	! flux leaving the passive pool
413	MatrixA(:,m,ipassive_pool,ipassive_pool) = moins_un * &
414	dtsom_turn(ipassive) &
415	control_moist(:,ibelow) * control_temp(:,ibelow) * decomp_factor(m)
416
417	! flux received by the passive pool from the active pool
418	MatrixA(:,m,ipassive_pool,iactive_pool) = frac_carb(:,iactive,ipassive)* &
419	dtsom_turn(iactive) &
420	control_moist(:,ibelow) * control_temp(:,ibelow) *&
421	( 1. - som_turn_iactive_clay_frac * clay(:) ) * decomp_factor(m)
422
423	! flux received by the passive pool from the slow pool
424	MatrixA(:,m,ipassive_pool,islow_pool) = frac_carb(:,islow,ipassive) * &
425	dtsom_turn(islow) &
426	control_moist(:,ibelow) * control_temp(:,ibelow) * decomp_factor(m)
427
428	! flux leaving the surface pool
429	MatrixA(:,m,isurface_pool,isurface_pool) = moins_un * &
430	dtsom_turn(isurface) &
431	control_moist(:,ibelow) * control_temp(:,ibelow) * decomp_factor(m)
432
433	WHERE (som(:,isurface,m,initrogen) .GT. min_stomate)
434	CN_som_litter_longterm(:,m,isurface_pool) = ( CN_som_litter_longterm(:,m,isurface_pool) * (tau_CN_longterm-dt) &
435	+ som(:,isurface,m,icarbon)/som(:,isurface,m,initrogen) * dt)/ (tau_CN_longterm)
436	ENDWHERE
437
438	WHERE (som(:,iactive,m,initrogen) .GT. min_stomate)
439	CN_som_litter_longterm(:,m,iactive_pool) = ( CN_som_litter_longterm(:,m,iactive_pool) * (tau_CN_longterm-dt) &
440	+ som(:,iactive,m,icarbon)/som(:,iactive,m,initrogen) * dt)/ (tau_CN_longterm)
441	ENDWHERE
442
443	WHERE(som(:,islow,m,initrogen) .GT. min_stomate)
444	CN_som_litter_longterm(:,m,islow_pool) = ( CN_som_litter_longterm(:,m,islow_pool) * (tau_CN_longterm-dt) &
445	+ som(:,islow,m,icarbon)/som(:,islow,m,initrogen) * dt)/ (tau_CN_longterm)
446	ENDWHERE
447
448	WHERE(som(:,ipassive,m,initrogen) .GT. min_stomate)
449	CN_som_litter_longterm(:,m,ipassive_pool) = ( CN_som_litter_longterm(:,m,ipassive_pool) * (tau_CN_longterm-dt) &
450	+ som(:,ipassive,m,icarbon)/som(:,ipassive,m,initrogen) * dt)/ (tau_CN_longterm)
451	ENDWHERE
452
453	IF (printlev>=4) WRITE(numout,*)'Finish to fill MatrixA'
454
455	ENDDO ! Loop over # PFTS
456
457
458	! 4.2 Add Identity for each submatrix(7,7)
459
460	DO j = 1,nbpools
461	MatrixA(:,:,j,j) = MatrixA(:,:,j,j) + un
462	ENDDO
463
464	ENDIF ! (spinup_analytic)
465
466	! Output diagnostics
467	DO k = 1, ncarb ! Loop over carbon pools
468	DO ij = 1, npts
469	IF (SUM(decomp_rate_soilcarbon(ij,:,k)*veget_cov_max(ij,:)) > min_sechiba) THEN
470	tsoilpools(ij,k) = 1./(SUM(decomp_rate_soilcarbon(ij,:,k)*veget_cov_max(ij,:))/dt_sechiba)
471	ELSE
472	tsoilpools(ij,k) = xios_default_val
473	END IF
474	END DO
475	END DO
476	CALL xios_orchidee_send_field("tSoilPools",tsoilpools)
477
478	IF (printlev>=4) WRITE(numout,*) 'Leaving som_dynamics'
479
480	END SUBROUTINE som_dynamics
481
482
483	!! ================================================================================================================================
484	!! SUBROUTINE : nitrogen_dynamics_clear
485	!!
486	!>\BRIEF Set the flag ::firstcall to .TRUE.
487	!!
488	!!
489	!_ ================================================================================================================================
490
491	SUBROUTINE nitrogen_dynamics_clear
492	firstcall_nitrogen=.TRUE.
493	END SUBROUTINE nitrogen_dynamics_clear
494
495
496
497
498	!! ================================================================================================================================
499	!! SUBROUTINE : nitrogen_dynamics
500	!!
501	!>\BRIEF : Describes mineralization dynamics of nitrogen species in the
502	!! soil. Inspired by the DNDC model of Li et al (1992,2000), but very simplified
503	!! to avoid having to calculate microbe growth for the moment. Builds on the
504	!! physico-chemical reactions with some fixed assumptions.
505	!!
506	!! DESCRIPTION : Plants can only uptake nitrogen in the form of NH4 and NO3.
507	!! In the soil, nitrogen coverts between NH3, NH4, NO3, NOX, N2O, and N2,
508	!! depending on things like the concentrations of each species, the pH of
509	!! the soil, the temperature, the amount of oxygen present, and the soil
510	!! moisture content. This subroutine describes the dynamics between
511	!! all these species, influencing the amount of nitrogen available for
512	!! both plant uptake and litter decomposition.
513	!!
514	!! Ammonium (NH4) may be the most important species, as all litter
515	!! N is assumed to decompose to NH4. NH4 can then be lost via
516	!! adsorpotion to soil clays, transformed to NH3 or NO3, or uptaken
517	!! by plants.
518	!!
519	!! The subroutine performs the following tasks:\n
520	!!
521	!! Section 2.\n
522	!! Immobilizes nitrogen, reducing the amount of nitrogen available in first
523	!! the ammonium (NH4) pool and then the nitrate (NO3) pool. This amount is
524	!! later added back into the mineralized soil NH4 pool. In essence, it
525	!! saves this quantity to be used in the decomposition routines. Without
526	!! this step, all of the nitrogen could be used up in this subroutine, thus
527	!! preventing litter decomposition from happening. \n
528	!!
529	!! If more nitrogen is needed for decomposition than is available from both
530	!! the NH4 and NO3 pools, the deficit is recorded in N_support throughout
531	!! the day. This deficit is printed to the history files, but apparently
532	!! never used.\n
533	!!
534	!! Section 3.\n
535	!! The goal of this section is to calculate the volumetric fraction of
536	!! aneorobic microsites in the soil, which gives an idea of how much
537	!! aneorobic bacteria can be transforming soil nitrogen. It depends
538	!! on the amount of oxygen in the soil, which depends on the concentration
539	!! difference of oxygen between the atmosphere and the soil, the diffusion
540	!! rate of oxygen in the soil, and how much oxygen is consumed by
541	!! respiration in the soil. For ideal gases, "concentration" is given by
542	!! the partial pressure.\n
543	!!
544	!! Section 4.\n
545	!! This section takes into account the losses of ammonium in the soil
546	!! due to adsorption onto clays, as well as the leaching of ammonium
547	!! and nitrate out of the system.\n
548	!!
549	!! Section 5.\n
550	!! Nitrification, that is, the change from ammonium (NH4+) to nitrate,
551	!! including effects of soil temperature, moisture, and pH, and
552	!! losses from NH4+ conversion to N2O and NO. \n
553	!!
554	!! Section 6.\n
555	!! Denitrification, that is the process of nitrate (NO3-) losing
556	!! oxygen atoms to eventually form N2, via the intermediate steps
557	!! of NO2-, NO, and N2O. All of this happens in the soil, and
558	!! depends on bacterial biomass, soil temperature, pH, and
559	!! soil water content.\n
560	!!
561	!! Section 7.\n
562	!! Calculates how much ammonium and nitrate is taken up by the plants.\n
563	!!
564	!! Section 8.\n
565	!! Calculates how much nitrogen is lost through emission into the
566	!! atmosphere through the surface of the soil. Depends on the rate of
567	!! diffusion through the soil. Nitrate emission is set to zero, but
568	!! the rest can pass from the soil to the atmosphere. Emission is given
569	!! the lowest priority of any loss pathway, in the sense that if the
570	!! amount of the species after changes due to leeching, nitrification, and
571	!! denitrification is less than what is calculated for the emissions, the
572	!! amount of emission is reduced to match.\n
573	!!
574	!! Section 9.\n
575	!! Update the values of the soil nitrogen pools, taking into account
576	!! all the losses and gains computed above. After this, add any
577	!! external nitrogen inputs, like fertilizer and biological nitrogen
578	!! fixation (BNF).\n
579	!!
580	!! Section 10.\n
581	!! Write output variables.\n
582	!!
583	!! RECENT CHANGE(S):
584	!!
585	!! MAIN OUTPUTS VARIABLE(S): soil_n_min, mineralisation
586	!!
587	!! REFERENCE(S) :
588	!! - S. Zaehle and A. D. Friend (2010), Carbon and nitrogen cycle dynamics in the
589	!! O-CN land surface model: 1. Model description, site-scale evaluation, and
590	!! sensitivity to parameter estimates. Global Biogeochem. Cycles, 24, GB1005,
591	!! doi:10.1029/2009GB003521.
592	!! - Li C. S., J. Aber, F. Stange, K. Butterbach-Bahl, and H. Papen (2000),
593	!! A process-oriented model of N2O and NO emissions from forest soils:
594	!! 1. Model development, Journal of Geophysical Research-Atmospheres,
595	!! 105, 4369-4384.
596	!! - Li, C., S. Frolking, and T. A. Frolking (1992), A model of nitrous
597	!! oxide evolution from soil driven by rainfall events: 1. Model
598	!! structure and sensitivity, J. Geophys. Res., 97(D9), 9759â9776,
599	!! doi:10.1029/92JD00509.
600	!! - Saxton, K.E., Rawls, W.J., Romberger, J.S., Papendick, R.I., 1986
601	!! Estimating generalized soil-water characteristics from texture.
602	!! Soil Sci. Soc. Am. J. 50, 1031-1036
603	!! - Kesik, M., Ambus, P., Baritz, R., BrÃŒggemann, N., et al.
604	!! Inventories of N2O and NO emissions from European forest soils,
605	!! Biogeosciences, 2, 353-375, https://doi.org/10.5194/bg-2-353-2005, 2005.
606	!! - Schmid, M., Neftel, A., Riedo, M. et al.
607	!! Process-based modelling of nitrous oxide emissions from different nitrogen sources in mown grassland
608	!! Nutrient Cycling in Agroecosystems 60: 177. https://doi.org/10.1023/A:1012694218748, 2001.
609	!!
610	!! FLOWCHART :
611	!!
612	!_ ================================================================================================================================
613	SUBROUTINE nitrogen_dynamics(npts, njsc, veget_cov_max,clay, sand, &
614	tsoil_decomp, tmc_pft, drainage_pft, runoff_pft, swc_pft, veget_max, resp_sol, &
615	som, biomass, input, pH, &
616	mineralisation, pb, plant_uptake, Bd,soil_n_min, &
617	p_O2,bact, N_support, &
618	cn_leaf_min_2D, cn_leaf_max_2D, cn_leaf_init_2D, &
619	mcs_hydrol, mcfc_hydrol,croot_longterm)
620	!
621	! 0 declarations
622	!
623
624	! 0.1 input
625
626	INTEGER(i_std), INTENT(in) :: npts ! Domain size
627	INTEGER(i_std),DIMENSION (npts), INTENT (in) :: njsc ! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
628	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_cov_max !! Fractional coverage: maximum share of the pixel taken by a pft
629	REAL(r_std), DIMENSION(npts), INTENT(in) :: clay ! clay fraction (between 0 and 1)
630	REAL(r_std), DIMENSION(npts), INTENT(in) :: sand ! sand fraction (between 0 and 1)
631	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsoil_decomp !! Temperature used for decompostition in soil (K)
632	REAL(r_std), DIMENSION (npts,nvm), INTENT(in) :: tmc_pft !! Total soil water per PFT (mm/m2)
633	REAL(r_std), DIMENSION (npts,nvm), INTENT(in) :: drainage_pft !! Drainage per PFT (mm/m2)
634	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: runoff_pft ! ! Runoff per PFT (mm/m2 /dt_sechiba)
635	! (fraction of water content)
636	REAL(r_std), DIMENSION (npts,nvm), INTENT(in) :: swc_pft !! Relative Soil water content [tmcr:tmcs] per pft (-)
637
638	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max! fraction of a vegetation (0-1)
639	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: resp_sol ! carbon respired from below ground
640	! (hetero+autotrophic) (gC/m**2/day)
641	REAL(r_std), DIMENSION(npts,ncarb,nvm,nelements), INTENT(in) :: som ! SOM (gC(or N)/m**2)
642	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: biomass ! Biomass pools (gC(or N)/m**2)
643	! ninput=4
644	REAL(r_std), DIMENSION(npts,nvm,ninput), INTENT(in) :: input ! nitrogen inputs into the soil (gN/m**2/day)
645	! NH4 and NOX from the atmosphere, NH4 from BNF,
646	! agricultural fertiliser as NH4/NO3
647	REAL(r_std), DIMENSION(npts), INTENT(in) :: pH ! soil pH
648	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: mineralisation ! net nitrogen mineralisation of decomposing SOM
649	! (gN/m**2/day), assumed to be NH4
650	REAL(r_std), DIMENSION(npts), INTENT(in) :: pb !! Air pressure (hPa)
651	!!!!!!!! NEVER USED?
652	REAL(r_std), DIMENSION(npts), INTENT(in) :: Bd !! Bulk density (kg/m**3)
653
654	!!!!!!!!
655	REAL(r_std),DIMENSION(npts,nvm), INTENT(in) :: cn_leaf_min_2D !! minimal leaf C/N ratio
656	REAL(r_std),DIMENSION(npts,nvm), INTENT(in) :: cn_leaf_max_2D !! maximal leaf C/N ratio
657	REAL(r_std),DIMENSION(npts,nvm), INTENT(in) :: cn_leaf_init_2D !! initial leaf C/N ratio
658	REAL(r_std),DIMENSION(npts,nvm), INTENT(in) :: croot_longterm !! "Long term" (default 3 years) root carbon mass
659	!! per ground area
660	!! @tex $(gC m^{-2} year^{-1})$ @endtex
661	REAL(r_std),DIMENSION (nscm), INTENT(in) :: mcs_hydrol !! Saturated volumetric water content output to be used in stomate_soilcarbon
662	REAL(r_std),DIMENSION (nscm), INTENT(in) :: mcfc_hydrol !! Volumetric water content at field capacity output to be used in stomate_soilcarbon
663
664	REAL(r_std), DIMENSION(npts,nvm,nnspec),INTENT(inout) :: soil_n_min !! mineral nitrogen in the soil (gN/m**2)
665	REAL(r_std), DIMENSION(npts,nvm),INTENT(inout) :: p_O2 !! partial pressure of oxygen in the soil (hPa)
666
667	REAL(r_std), DIMENSION(npts,nvm),INTENT(inout) :: bact !! denitrifier biomass (gC/m**2)
668
669
670	! 0.2 output
671
672	REAL(r_std), DIMENSION(npts,nvm,nionspec), INTENT(out) :: plant_uptake !! Uptake of soil N by plants
673	!! (gN/m**2/timestep)
674	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: N_support !! Nitrogen which is added to the ecosystem to support vegetation growth
675	! 0.3 local
676	REAL(r_std) :: dt !! Time step \f$(dt_sechiba one_day^{-1})$\f
677	LOGICAL :: l_error !! Diagnostic boolean for error allocation (true/false)
678	INTEGER(i_std) :: ier !! Check errors in netcdf call
679	REAL(r_std), DIMENSION(npts,nvm,nionspec) :: leaching !! mineral nitrogen leached from the soil
680	!! (gN/m**2/timestep)
681	REAL(r_std), DIMENSION(npts,nvm) :: immob !! N immobilized (gN/m**2/day)
682	REAL(r_std), DIMENSION(npts) :: afps_max !! maximum pore-volume of the soil
683	!! Table 2, Li et al., 2000
684	!! (fraction)
685	REAL(r_std), DIMENSION(npts,nvm) :: afps !! pore-volume of the soil
686	!! Table 2, Li et al., 2000
687	!! (fraction)
688	REAL(r_std), DIMENSION(npts,nvm) :: D_s !! Oxygen diffusion in soil (m**2/day)
689	REAL(r_std), DIMENSION(npts,nvm) :: mol_O2_resp !! Density of Moles of O2 related to respiration term
690	!! ((molesO2 m-3) (gC m-2)-1)
691	REAL(r_std), DIMENSION(npts,nvm) :: p_O2_resp !! O2 partial pressure related to the respiration term
692	!! ((hPa) (gC m-2)-1)
693	REAL(r_std), DIMENSION(npts) :: p_O2air !! Oxygen partial pressure in air (hPa)
694	REAL(r_std), DIMENSION(npts) :: g_O2 !! soil gradient of O2 partial pressure (hPa m-1)
695	REAL(r_std), DIMENSION(npts) :: d_O2 !! change in O2 partial pressure (hPa day-1)
696	REAL(r_std), DIMENSION(npts,nvm) :: anvf !! Volumetric fraction of anaerobic microsites
697	!! (fraction of pore-volume)
698	REAL(r_std), DIMENSION(npts) :: FixNH4 !! Fraction of adsorbed NH4+ (-)
699	REAL(r_std), DIMENSION(npts,nvm) :: n_adsorbed !! Ammonium adsorpted (gN/m**2)
700	!! based on Li et al. 1992, JGR, Table 4
701	REAL(r_std), DIMENSION(npts,nvm) :: fwnit,fwdenit !! Effect of soil moisture on nitrification (-)
702	!! Zhang et al. 2002, Ecological Modelling, appendix A, page 101
703	REAL(r_std), DIMENSION(npts) :: var_temp_sol!! Temperature function used for calc. ft_nit (-)
704	!! Zhang et al. 2002, Ecological Modelling, appendix A, page 101
705	REAL(r_std), DIMENSION(npts) :: ft_nit !! Effect of temperature on nitrification (-)
706	!! Zhang et al. 2002, Ecological Modelling, appendix A, page 101
707	REAL(r_std), DIMENSION(npts) :: fph !! Effect of pH on nitrification (-)
708	!! Zhang et al. 2002, Ecological Modelling, appendix A, page 101
709	REAL(r_std), DIMENSION(npts) :: ftv !! Effect of temperature on NO2 or NO production
710	!! during nitrification (-)
711	!! Zhang et al. 2002, Ecological Modelling, appendix A, page 102
712	REAL(r_std), DIMENSION(npts,nvm,3) :: nitrification !! N-compounds production (NO3, N2O, NO) related
713	!! to nitrification process (gN/m**2/tstep)
714	REAL(r_std), DIMENSION(npts) :: ft_denit !! Temperature response of relative growth rate of
715	!! total denitrifiers (-)
716	!! Eq. 2 Table 4 of Li et al., 2000
717	REAL(r_std), DIMENSION(npts) :: fph_no3 !! Soil pH response of relative growth rate of
718	!! NO3 denitrifiers (-)
719	!! Eq. 2 Table 4 of Li et al., 2000
720	REAL(r_std), DIMENSION(npts) :: fph_no !! Soil pH response of relative growth rate of
721	!! NO denitrifiers (-)
722	!! Eq. 2 Table 4 of Li et al., 2000
723	REAL(r_std), DIMENSION(npts) :: fph_n2o !! Soil pH response of relative growth rate of
724	!! N2O denitrifiers (-)
725	!! Eq. 2 Table 4 of Li et al., 2000
726	REAL(r_std), DIMENSION(npts,nvm) :: Kn_conv !! Conversion from kgN/m3 to gN/m2
727	REAL(r_std), DIMENSION(npts) :: mu_NO3 !! Relative growth rate of NO3 denitrifiers (hour**-1)
728	!! Eq.1 Table 4 Li et al., 2000
729	REAL(r_std), DIMENSION(npts) :: mu_N2O !! Relative growth rate of N2O denitrifiers (hour**-1)
730	!! Eq.1 Table 4 Li et al., 2000
731	REAL(r_std), DIMENSION(npts) :: mu_NO !! Relative growth rate of NO denitrifiers (hour**-1)
732	!! Eq.1 Table 4 Li et al., 2000
733	REAL(r_std), DIMENSION(npts) :: sum_n !! sum of all N species in the soil (gN/m**2)
734	REAL(r_std), DIMENSION(npts,nvm,3) :: denitrification !! N-compounds consumption (NO3, N2O, NO) related
735	!! to denitrificaiton process (gN/m**2/tstep)
736	REAL(r_std), DIMENSION(npts) :: dn_bact !! denitrifier biomass change (kgC/m**3/timestep)
737	REAL(r_std), DIMENSION(npts) :: ft_uptake !! Temperature response of N uptake by plants (-)
738	REAL(r_std) :: conv_fac_vmax!! Conversion from (umol (gDW)-1 h-1) to (gN (gC)-1 timestep-1)
739	REAL(r_std), DIMENSION(npts,nvm) :: nc_leaf_min !! Minimal NC ratio of leaf (gN / gC)
740	REAL(r_std), DIMENSION(npts,nvm) :: nc_leaf_max !! Maximal NC ratio of leaf (gN / gC)
741	REAL(r_std), DIMENSION(npts) :: lab_n !! Labile nitrogen in plants (gN/m**2)
742	REAL(r_std), DIMENSION(npts) :: lab_c !! Labile carbon in plants (gC/m**2)
743	REAL(r_std), DIMENSION(npts) :: NCplant !! NC ratio of the plant (gN / gC)
744	!! Eq. (9) p. 3 of SM of Zaehle & Friend, 2010
745	REAL(r_std), DIMENSION(npts) :: f_NCplant !! Response of Nitrogen uptake by plants
746	!! to N/C ratio of the labile pool
747	REAL(r_std) :: conv_fac_concent!! Conversion factor from (umol per litter) to (gN m-2)
748	REAL(r_std), DIMENSION(npts) :: frac_nh3 !! dissociation of [NH3] to [NH4+] (-)
749	REAL(r_std), DIMENSION(npts) :: F_clay !! Response of N-emissions to clay fraction (-)
750	REAL(r_std), DIMENSION(npts,nvm,nnspec) :: emission !! volatile losses of nitrogen
751	!! (gN/m**2/timestep)
752	INTEGER(i_std) :: m !! Index to loop over the nvm PFT's
753	REAL(r_std), DIMENSION(npts,nvm) :: f_drain !! fraction of tmc which has been drained in the last time step
754	REAL(r_std), DIMENSION(npts) :: temp_sol,temp_sol_K !! soil temperature (C) and soil temperature (K)
755
756	!_ ================================================================================================================================
757
758	!! printlev is the level of diagnostic information, 0 (none) to 4 (full)
759	IF (printlev>=3) WRITE(numout,*) 'Entering nitrogen_dynamics'
760	IF(printlev>=4)THEN
761	WRITE(numout,*) 'CHECK values in nitrogen dynamics'
762	WRITE(numout,*) 'soil_n_min ',soil_n_min(test_grid,test_pft,:)
763	ENDIF
764
765	!! 1. Initializations
766	dt = dt_sechiba/one_day
767
768	IF ( firstcall_nitrogen ) THEN
769	firstcall_nitrogen = .FALSE.
770	ENDIF
771
772	! Transform tsoil_decomp into degree C
773	temp_sol(:) = tsoil_decomp - tp_00
774
775	IF(printlev>=4)THEN
776	WRITE(numout,*) 'CHECK values in after init'
777	WRITE(numout,*) 'soil_n_min ',soil_n_min(test_grid,test_pft,:)
778	ENDIF
779
780	IF(printlev>=4)THEN
781	WRITE(numout,*) 'CHECK values in after leaching'
782	WRITE(numout,*) 'leaching(test_grid,test_pft,iammonium) ',leaching(test_grid,test_pft,iammonium)
783	WRITE(numout,*) 'leaching(test_grid,test_pft,initrate) ',leaching(test_grid,test_pft,initrate)
784	WRITE(numout,*) 'drainage_pft(test_grid) ',drainage_pft(test_grid,test_pft)
785	ENDIF
786
787
788	IF(printlev>=4)THEN
789	WRITE(numout,*) 'PFT=',test_pft
790	WRITE(numout,*) 'leaching(test_grid,test_pft,iammonium) ',leaching(test_grid,test_pft,iammonium)
791	WRITE(numout,*) 'leaching(test_grid,test_pft,initrate) ',leaching(test_grid,test_pft,initrate)
792	ENDIF
793
794
795	IF(printlev>=4)THEN
796	WRITE(numout,*) 'drainage_pft(test_grid,test_pft) ',drainage_pft(test_grid,test_pft)
797	ENDIF
798
799
800	!! 2. Preparation for decomposition
801	!
802	! 2.1 conservation of mass from decomposition
803	! immobilisation has absolute priority to avoid mass conservation problems
804	! the code in litter and soilcarbon has to make sure that soil ammonium can never be
805	! more than exhausted completely by immobilisation!!!
806	! As a reminder, mineralisation is the amount of litter destined to become soil NH4
807	!
808	!
809	immob(:,:) = zero
810	WHERE(mineralisation(:,:).LT.0.)
811	immob(:,:) = - mineralisation(:,:)
812	! In case, the N related to immobilisation is higher that the N
813	! available in the [NH4+] pool, we take the remaining N from the
814	! [NO3-] pool
815	soil_n_min(:,:,initrate) = soil_n_min(:,:,initrate) - &
816	MAX(0.0,immob(:,:)-(soil_n_min(:,:,iammonium)-min_stomate))
817
818	soil_n_min(:,:,iammonium) = soil_n_min(:,:,iammonium) - &
819	MIN(immob(:,:),soil_n_min(:,:,iammonium)-min_stomate)
820	ENDWHERE
821
822
823	! 2.2 Deficit
824	! In case of soil_n_min(:,:,initrate) negative, we add nitrogen and we take memory of the
825	! added nitrogen in N_support
826	! This happens when immob(:,:) is bigger than the nitrogen that can be given by the two
827	! soil_n_min pools (nitrate and ammonium) and so soil_n_min(:,:,initrate) becomes negative.
828	N_support(:,:) = 0.
829	WHERE(soil_n_min(:,:,initrate) .LT. 0.)
830	N_support(:,:) = - soil_n_min(:,:,initrate)
831	soil_n_min(:,:,initrate) = 0.
832	ENDWHERE
833
834
835	IF(printlev>=4)THEN
836	WRITE(numout,*) 'CHECK values after mineralisation'
837	WRITE(numout,*) 'soil_n_min ',soil_n_min(test_grid,test_pft,:)
838	ENDIF
839
840	!! 3. ANVF
841	!
842	! The goal of this section is to calculate the volumetric fraction of
843	! aneorobic microsites (ANVF) in the soil, which gives an idea of how much
844	! aneorobic bacteria can be transforming soil nitrogen. It depends
845	! on the oxygen diffusion and partial pressure in the soil. These euqations
846	! are taken from Table 2 of Li et al. (2000) Table 2, though some things
847	! remain unclear or undefinied. afps, for example, is never explicitly given.
848	! S. Zaehle defined it with the following equation - but did not use it.
849	! In OCN, diffusion is not of afps/afps_max ratio but accounts for soilhum and
850	! soil temperature according to Monteith & Unsworth, 1990
851	! d_ox(:) = 1.73664 * ( 0.15 * (exp(-(soilhum_av(:)*3.)/0.44)-exp(-1./0.44))) &
852	! (1.+0.007*tsoil_av(:))
853
854	! The equation for afpsmax, the maximum pore volume of the soil,
855	! is taken from Table 2 of Saxton (1986)
856	afps_max(:) = h_saxton + j_saxton * (sand(:)100.) + k_saxton log10(clay(:)*100.)
857
858	DO m=1,nvm
859	! pore volume of the soil. Unclear where this comes from.
860	afps(:,m) = MAX(0.1,( un - swc_pft(:,m) ) * afps_max(:))
861
862	! diffusion through the soil. Table 2 of Li et al. (2000)
863	D_s(:,m) = D_air * ( afps(:,m)diffusionO2_power_1 ) / ( afps_max(:)diffusionO2_power_2 )
864
865
866	! Account for the impact of frost on diffusion. Table 2 of Li et al. (2000)
867	WHERE ( temp_sol(:) .GT. zero )
868	D_s(:,m) = D_s(:,m) * F_nofrost
869	ELSEWHERE
870	D_s(:,m) = D_s(:,m) * F_frost
871	ENDWHERE
872
873	! Equation (3) - Oxygen partial pressure
874	! Written in an odd differential form in Table 2 of Li et al (2000), with the
875	! time derivative of the partial pressure being related to the depth derivative
876	! of both the partial pressure and the diffusion rate.
877
878	! So we do it this way, calculating the depth gradient afterwards.
879	! Oxygen loss from respiration has to be expressed as a partial pressure
880	! using the ideal gas PV=nRRT relationship with P in Pa, V in m-3 and T in K
881	! RR is the ideal gas constat (J mol-1 K-1) and n the number of moles
882	! Respiration, one mole of C requires two moles of oxygen (CO2)
883	! Resp_below expressed in gC m-2 day-1
884	! mol_O2_resp : Density of Moles of O2 related to respiration term (molesO2 m-3) (gC m-2)-1
885	! In OCN, afps_max is used instead of afps. We keep here the original formulation
886	mol_O2_resp(:,m) = (un / C_molar_mass * 2. ) / (zmaxh * afps(:,m))
887	! O2 partial pressure related to the respiration term (hPa) (gC m-2)-1
888	p_O2_resp(:,m) = mol_O2_resp(:,m) * RR * ( temp_sol(:) + tp_00 ) * Pa_to_hPa
889	ENDDO
890
891
892	! Change in oxygen partial pressure d_O2 - Ds x gradient
893	! Not sure that we should use z_decomp (could be a fraction of zmaxh, too)
894	! oxygen partial pressure in air - p_O2air (hPa)
895	p_O2air(:)= V_O2 * pb(:)
896
897	DO m = 1, nvm
898
899	! assume the partial pressure of oxygen in the soil is uniform over the
900	! whole depth that soil decomposers are active, in order to calculate
901	! the change in partial pressure of oxygen in the soil due to the
902	! pressure different between the air and soil and the diffusion rate
903	! of oxygen through the soil
904	g_O2(:) = ( p_O2air(:) - p_O2(:,m) ) / z_decomp
905	d_O2(:) = D_s(:,m) / z_decomp * g_O2(:)
906	! D_s / z_decomp has the unit of a conductivity (m day-1)
907
908	! The new partial pressure of oxygen in the soil depends on how
909	! much oxgen has entered or left the soil due to the pressure difference
910	! between the atmosphere and the soil, in addition to any losses from
911	! soil respiration consuming oxygen to create CO2
912	p_O2(:,m) = p_O2(:,m) + d_O2(:)dt - p_O2_resp(:,m)resp_sol(:,m)*dt
913	! Equation (4) Volumetric fraction of anaerobic microsites (ANVF)
914	! a and b constants are not specified in Li et al., 2000
915	! S. Zaehle used a=0.85 and b=1 without mention to any publication
916	! a_anvf=0.85
917	! b_anvf=1.
918	anvf(:,m) = a_anvf * ( 1 - b_anvf * p_O2(:,m) / p_O2air(:) )
919	anvf(:,m) = MAX(zero, anvf(:,m))
920	ENDDO
921
922	IF(printlev>=4 .AND. m==test_pft)THEN
923	WRITE(numout,*) 'CHECK values after nitrification nh4 to no3'
924	WRITE(numout,*) 'PFT=',test_pft
925	WRITE(numout,*) 'anvf(:,m)=',anvf(test_grid,test_pft)
926	WRITE(numout,*) 'p_O2(:,m)=',p_O2(test_grid,test_pft)
927	WRITE(numout,*) 'p_O2air(:)=',p_O2air(test_grid)
928	WRITE(numout,*) 'pb(:)=',pb(test_grid)
929	WRITE(numout,*) 'd_O2(:)=',d_O2(test_grid)
930	WRITE(numout,*) 'p_O2_resp(:)=',p_O2_resp(test_grid,:)
931	WRITE(numout,*) 'mol_O2_resp(:)=',mol_O2_resp(test_grid,:)
932	WRITE(numout,*) 'resp_soil(:,m)=',resp_sol(test_grid,test_pft)
933	WRITE(numout,*) 'afps(:)=',afps(test_grid,:)
934
935	ENDIF
936
937	!! 4. Physical removal of NH4 and NO3
938	!
939	! 4.1 Adsorption of NH4 into clay
940	!
941	! Reduce soil ammonium concentrations according to the equations in
942	! Table 4 of Li et al. (1992). This represents adsorption onto soil clays.
943	! FixNH4 : Fraction of adsorbed NH4+
944	! FixNH4=[0.41-0.47 log(NH4) clay/clay_max]
945	! NH4+ concentration in the soil liquid, gN kg-1 soil (p. 9774 of Li et al., 1992)
946	! but Zhang et al. (2002) Appendix A/B define NH4 as
947	! NH4+ in a soil layer in kgN ha-1
948	! In OCN, NH4 seems defined as kgN kg-1 or m-3 of water
949	! Comment of N. Vuichard: I don't know which definition is the good one ...
950	DO m = 1, nvm
951	WHERE(soil_n_min(:,m,iammonium).GT.min_stomate)
952	FixNH4(:) = MAX(0.,(a_FixNH4 + b_FixNH4 * &
953	MAX(0.,log10(soil_n_min(:,m,iammonium)*10 )))) &
954	* MIN(clay(:),clay_max) / clay_max
955	ELSEWHERE
956	FixNH4(:) = 0.0
957	ENDWHERE
958	! In OCN, we do not multiply by FixNH4 but by FixNH4/(1+FixNH4)
959	! Comment of N. Vuichard: It is not clear if we should keep the formulation of OCN or not
960	! So far, we keep the original formulation
961	n_adsorbed(:,m) = MIN(soil_n_min(:,m,iammonium),soil_n_min(:,m,iammonium) * FixNH4(:))
962	ENDDO
963	soil_n_min(:,:,iammonium) = soil_n_min(:,:,iammonium) - n_adsorbed(:,:)
964
965	IF(printlev>=4)THEN
966	WRITE(numout,*) 'CHECK values after adsorption'
967	WRITE(numout,*) 'soil_n_min ',soil_n_min(test_grid,test_pft,:)
968	ENDIF
969
970
971	! 4.2 Leaching of NH4 and NO3
972	! Initializations
973	leaching(:,:,:) = zero
974
975	f_drain(:,:) = zero
976
977
978	WHERE((tmc_pft(:,:)+runoff_pft(:,:)) .NE. 0)
979	f_drain(:,:) = (fracn_drainagedrainage_pft(:,:)+fracn_runoffrunoff_pft)/(tmc_pft(:,:)+runoff_pft(:,:))
980	ENDWHERE
981	f_drain(:,:) = MAX(zero, MIN(f_drain(:,:),un))
982
983	DO m = 1, nvm
984	!leaching
985	leaching(:,m,iammonium) = MIN(soil_n_min(:,m,iammonium) * f_drain(:,m), &
986	soil_n_min(:,m,iammonium) )
987	leaching(:,m,initrate) = MIN(soil_n_min(:,m,initrate) * f_drain(:,m), &
988	soil_n_min(:,m,initrate) )
989	ENDDO
990
991
992	!! 5. Nitrification of NH4 in the oxygenated part of the soil
993	!
994	! Nitrification refers to ammonium (NH4+) being oxidized to
995	! nitrate (NO3-) under aeroboic (i.e., in the prescence of
996	! oxygen) conditions. This sections thus reduces the
997	! concentration of ammonium in the soil. These equations
998	! are taken from Zhang et al. (2002), Appendix A
999	!
1000	! 5.1 Effect of environmental factors (Temp, Humidity, pH)
1001	! 5.1.1 Effect of soil moisture on nitrification
1002	fwnit(:,:) = fwnit_0 + fwnit_1 * swc_pft(:,:) + fwnit_2 * swc_pft(:,:)*2 + fwnit_3 swc_pft(:,:)**3 &
1003	+ fwnit_4 * swc_pft(:,:)**4
1004	fwnit(:,:) = MAX (0., MIN( un, fwnit(:,:) ) )
1005
1006	! 5.1.2 Effect of temperature on nitrification
1007	! Note that Zhang et al 2002 fold in the factor 0.1 with the parameter
1008	! for the term linear in temperature (ft_nit_1), while we group it
1009	! with the soil temperature as per the rest of the terms. Checking
1010	! the parameter values, this leads to a first impression that they
1011	! differ by a factor of 10, when in reality the same overall result
1012	! is calculated.
1013	var_temp_sol(:) = temp_sol(:) * 0.1
1014	ft_nit(:) = ft_nit_0 + ft_nit_1 * var_temp_sol(:) + ft_nit_2 * var_temp_sol(:)*2 + ft_nit_3 var_temp_sol(:)**3 &
1015	+ ft_nit_4 * var_temp_sol(:)**4
1016	ft_nit(:) = MAX (0.001, MIN( un, ft_nit(:) ) )
1017
1018	! 5.1.3 Effect of pH on nitrification
1019	fph(:) = fph_0 + fph_1 * pH(:) + fph_2 * ph(:)**2
1020	fph(:) = MAX(0.0, fph(:) )
1021
1022	! 5.1.4 Effect of temperature on NO2 or NO production during nitrification
1023	! I am not sure why the factor of 0.0025*NO3,N seems to be missing from
1024	! the equation reported in the literature.
1025	ftv(:) = ftv_0 **(ftv_1 - ftv_2 /(temp_sol(:) + tp_00) )
1026	ftv(:) = MAX(0.0, ftv(:) )
1027
1028	DO m = 1, nvm
1029
1030	! 5.2 Actual nitrification rate per PFT NH4 soil pool - NH4 to NO3
1031	!
1032	! Equation for the nitrification rate probably from Schmid et al., 2001, Nutr. Cycl. Agro (eq.1)
1033	! but the environmental factors used are from Zhang (2002) who used an other equation
1034	! I don't know how this can be mixed together ?
1035	! In addition, the formulation mixed the one from Schmid with the use of the anaerobic balloon defined by Li et al., 2000. I don't know how this can be mixed together ?
1036	! Last, in OCN, the default fraction of N-NH4 which is converted to N-NO3 appears equal to 2 day-1 (a factor "2" in the equation) - In Schmid et al., the nitrification rate at 20 âŠC and field capacity (knitrif,20) was set to 0.2 dâ1 (Ten time less...)
1037	! k_nitrif = 0.2
1038	nitrification(:,m,i_nh4_to_no3) = MIN(fwnit(:,m) * fph(:) * ft_nit(:) * k_nitrif * dt * &
1039	soil_n_min(:,m,iammonium) ,soil_n_min(:,m,iammonium)- leaching(:,m,iammonium))
1040
1041	IF(printlev>=4 .AND. m == test_pft)THEN
1042	WRITE(numout,*) 'CHECK values after nitrification nh4 to no3'
1043	WRITE(numout,*) 'PFT=',m
1044	WRITE(numout,*) 'nitrification(:,m,i_nh4_to_no3) ',nitrification(test_grid,m,i_nh4_to_no3)
1045	WRITE(numout,*) 'fwnit=',fwnit(test_grid,m)
1046	WRITE(numout,*) 'fph=',fph(test_grid)
1047	WRITE(numout,*) 'ft_nit=',ft_nit(test_grid)
1048	WRITE(numout,*) 'anvf=',anvf(test_grid,m)
1049	ENDIF
1050
1051	!
1052	! 5.3 Emission of N2O during nitrification - NH4 to N2O
1053	!
1054	nitrification(:,m,i_nh4_to_n2o) = ftv(:) * swc_pft(:,m) * n2o_nitrif_p * &
1055	nitrification(:,m,i_nh4_to_no3)
1056
1057	IF(printlev>=4 .AND. m==test_pft)THEN
1058	WRITE(numout,*) 'CHECK values after nitrification nh4 to n2o'
1059	WRITE(numout,*) 'nitrification(:,m,i_nh4_to_n2o) ',nitrification(test_grid,m,i_nh4_to_n2o)
1060	WRITE(numout,*) 'ftv=',ftv(test_grid)
1061	WRITE(numout,*) 'swc_pft=',swc_pft(test_grid,m)
1062	ENDIF
1063
1064	nitrification(:,m,i_nh4_to_n2o) = MIN(nitrification(:,m,i_nh4_to_no3),nitrification(:,m,i_nh4_to_n2o))
1065	!
1066	! 5.4 Production of NO during nitrification - NH4 to NO
1067	nitrification(:,m,i_nh4_to_no) = ftv(:) * no_nitrif_p * nitrification(:,m,i_nh4_to_no3)
1068
1069
1070	IF(printlev>=4 .AND. m == test_pft)THEN
1071	WRITE(numout,*) 'CHECK values after nitrification nh4 to no'
1072	WRITE(numout,*) 'nitrification(:,m,i_nh4_to_no) ',nitrification(test_grid,m,i_nh4_to_no)
1073	ENDIF
1074	! 5.5 Production of NO from chemodenitrification
1075	! Chemodenitrification is the conversion of NO2 to NO through a purely
1076	! chemical process, dependent on temperature, soil pH, and concentration of
1077	! NO2. This is the final step of the process NH4+ -> NO3- -> NO2- -> NO
1078	! Based on Chem_NO in Kesik et al.(2005)
1079	! BUT Kesik et al. used NO2 concentration and not NO3 production as it is done in OCN
1080	! and modification of a multiplicative constant from 300 (Kesik) to 30 (OCN)
1081	! without clear motivation - I don't how this is reliable
1082	! [Matt also questions this equation...as it is written below, it is the amount
1083	! of NH4 being converted to NO3 multiplied by the decomposition of NO2 to NO, but
1084	! there is no NO3 to NO2 conversion included, which implies that all NO3 is
1085	! instantaneously converted to NO2...which is nonsense based on the fact that an NO3
1086	! pool exists]
1087	nitrification(:,m,i_nh4_to_no) = nitrification(:,m,i_nh4_to_no) + &
1088	( chemo_0 * chemo_1 * exp(chemo_ph0 * pH(:) ) * &
1089	exp(chemo_t0/((temp_sol(:)+tp_00)RR))) nitrification(:,m,i_nh4_to_no3)
1090
1091	nitrification(:,m,i_nh4_to_no) = MIN(nitrification(:,m,i_nh4_to_no3)-nitrification(:,m,i_nh4_to_n2o), &
1092	nitrification(:,m,i_nh4_to_no))
1093
1094
1095	IF(printlev>=4 .AND. m == test_pft)THEN
1096	WRITE(numout,*) 'CHECK values after nitrification nh4 to no chemodenitrification'
1097	WRITE(numout,*) 'nitrification(:,m,i_nh4_to_no) ',nitrification(test_grid,m,i_nh4_to_no)
1098	WRITE(numout,*) 'pH(:)=',pH(test_grid)
1099	WRITE(numout,*) 'temp_sol(:)=',temp_sol(test_grid)
1100	ENDIF
1101
1102	! In OCN, NO production and N2O production is deducted from NO3 production due to NH4
1103	nitrification(:,m,i_nh4_to_no3) = nitrification(:,m,i_nh4_to_no3) - &
1104	(nitrification(:,m,i_nh4_to_no) + nitrification(:,m,i_nh4_to_n2o))
1105	ENDDO
1106
1107
1108	!! 6. Denitrification processes
1109	!
1110	! Denitrification is the conversion of any oxidized N comound (NO2-, NO, N2O)
1111	! to the final product of gaseous nitrogen (N2).
1112	!
1113	! Denitrifiers are organisisms responsible for denitrification (bacteria?).
1114	!
1115	! Li et al, 2000, JGR Table 4 eq 1, 2 and 4, ignoring NO2 (similar to NO)
1116	! Comment of S. Zaehle: "includes treatment of bacteria dynamics, but I do not have any confidence in this;
1117	! at least it provides rates that appear resonable."
1118
1119	! 6.1 Temperature response of relative growth rate of total denitrifiers
1120	ft_denit(:) = 2.**((temp_sol(:)-22.5)/10.)
1121
1122	! 6.2 Soil pH response of relative growth rate of total denitrifiers
1123	! See also comment about parenthesis' position in Thesis of Vincent Prieur (page 50)
1124	fph_no3(:) = 1.0 - 1.0 / ( 1.0 + EXP((pH(:)-fph_no3_0)/fph_no3_1))
1125	fph_no(:) = 1.0 - 1.0 / ( 1.0 + EXP((pH(:)-fph_no_0)/fph_no_1))
1126	fph_n2o(:) = 1.0 - 1.0 / ( 1.0 + EXP((pH(:)-fph_n2o_0)/fph_n2o_1))
1127
1128	! Half saturation value of N oxides (kgN/m3), Kn
1129	! Unit conversion factor from kgN/m3 to gN/m2
1130	! OCN used the value 0.087, but 0.083 is listed in Li et al (2000).
1131	! In OCN, we account for max_eau_eau. Not clear if this is correct or not.
1132	Kn_conv(:,:) = Kn * tmc_pft(:,:)
1133
1134	! Relative growth rate of Nox denitrifiers
1135	! Eq.1 Table 4 Li et al., 2000 - but there is an error because Eq. 1 it does not
1136	! account for Temp and ph responses.
1137	! In OCN it does not account for [DOC], though Li et al. (2000) does. We keep
1138	! the formulation from OCN. In addition, in OCN the term mu_nox(max) is missing
1139	! in the equation that defines the relative growth rate of total denitrifiers
1140	! (dn_bact) - we add the term back in here.
1141	!
1142	! 6.3 Maximum Relative growth rate of Nox denitrifiers (hour-1): mu_no3, mu_no2, mu_no
1143
1144	denitrification(:,:,:)=zero
1145
1146	DO m = 1, nvm
1147
1148	WHERE((tmc_pft(:,m)/(zmaxh*1000.)) .LT. mcfc_hydrol(njsc(:)))
1149	! verifier zah + unite m ou mm
1150	fwdenit(:,m) = fwdenitfc * &
1151	exp(kfwdenit * (mcfc_hydrol(njsc(:))-tmc_pft(:,m)/(zmaxh*1000.))/mcfc_hydrol(njsc(:)) )
1152	ELSEWHERE
1153	fwdenit(:,m) = fwdenitfc + (1 - fwdenitfc) * &
1154	(tmc_pft(:,m)/(zmaxh*1000.) - mcfc_hydrol(njsc(:))) / (mcs_hydrol(njsc(:))-mcfc_hydrol(njsc(:)))
1155	ENDWHERE
1156
1157	WHERE((soil_n_min(:,m,initrate) + Kn_conv(:,m)) .GT. min_stomate)
1158	mu_no3(:) = mu_no3_max * soil_n_min(:,m,initrate) / (soil_n_min(:,m,initrate) + Kn_conv(:,m))
1159	ELSEWHERE
1160	mu_no3(:) = zero
1161	ENDWHERE
1162	WHERE((soil_n_min(:,m,inox) + Kn_conv(:,m)*fact_kn_no) .GT. min_stomate)
1163	mu_no(:) = mu_no_max * soil_n_min(:,m,inox) / (soil_n_min(:,m,inox) + Kn_conv(:,m)*fact_kn_no)
1164	ELSEWHERE
1165	mu_no(:) = zero
1166	ENDWHERE
1167
1168	WHERE((soil_n_min(:,m,initrous) + Kn_conv(:,m)*fact_kn_n2o) .GT. min_stomate)
1169	mu_n2o(:) = mu_n2o_max * soil_n_min(:,m,initrous) / (soil_n_min(:,m,initrous) + Kn_conv(:,m)*fact_kn_n2o)
1170	ELSEWHERE
1171	mu_n2o(:) = zero
1172	ENDWHERE
1173
1174	! stocks of all N species (NO3, NO, N2O) (gN/m**2)
1175	sum_n(:) = soil_n_min(:,m,inox) + soil_n_min(:,m,initrous) + &
1176	soil_n_min(:,m,initrate)
1177
1178	WHERE(sum_n(:).GT.min_stomate)
1179
1180	! This appears to be an assumption that the denitrifier biomass is
1181	! 0.05% of the active carbon soil organic matter. Unclear where
1182	! this comes from.
1183	! bact(:,m) = 0.0005*som(:,iactive,m,icarbon)
1184	! bact(:,m) = 0.00005*som(:,iactive,m,icarbon)
1185	bact(:,m) = cte_bact*som(:,iactive,m,icarbon)
1186	! 6.4 Consumption rate of N oxides
1187	! Table 4 Li et al.(2000)
1188	! Based on the maximum growth rate on N oxides (Y_nox), maintainance
1189	! coefficient on N oxides (M_nox), denitrifier biomass (bact),
1190	! relative growth rate of NOX denitrifiers (mu_nox), concentration
1191	! of NOX in the soil (soil_n_min), total concentration of all nitrogen
1192	! species in the soil (sum_n), the amount of NOX leached out of the
1193	! soil (leaching). The fwdenit, ft_denit, and fph_nox do not come
1194	! from Table 4. Conversion from (per hour) to (per timestep) is
1195	! handled by 24*dt at the end.
1196	!
1197	! The reactions are NO3 -> (NO + NO2) -> N2O -> N2
1198	!
1199	! NO3 consumption
1200	! In OCN, multiplication by 0.1 of the NO3 consumption - no justification
1201	! This is not kept in the current version
1202	denitrification(:,m,i_no3_to_nox) = MIN(soil_n_min(:,m,initrate)-leaching(:,m,initrate), &
1203	fwdenit(:,m) * ft_denit(:) * fph_no3(:) * ( mu_no3(:) / Y_no3 + M_no3 * soil_n_min(:,m,initrate) / sum_n(:)) * &
1204	bact(:,m) * 24. * dt )
1205	!
1206	! NO consumption
1207	denitrification(:,m,i_nox_to_n2o) = MIN(soil_n_min(:,m,inox), &
1208	fwdenit(:,m) * ft_denit(:) * fph_no(:) * ( mu_no(:) / Y_no + M_no * soil_n_min(:,m,inox) / sum_n(:)) * &
1209	bact(:,m) * 24. * dt )
1210	!
1211	! N2O consumption
1212	denitrification(:,m,i_n2o_to_n2) = MIN(soil_n_min(:,m,initrous), &
1213	fwdenit(:,m) * ft_denit(:) * fph_n2o(:) * ( mu_n2o(:) / Y_n2o + M_n2o * soil_n_min(:,m,initrous) / sum_n(:)) &
1214	* bact(:,m) * 24. * dt )
1215
1216	! Dynamcis on denitrifier bacterial population change is not used, due to
1217	! lack of documention in OCN and lack of clarity in Li et al (2000).
1218	! In addition, OCN used dn_bact, but bact is used here. Not clear what
1219	! the consequences of this are.
1220
1221	ENDWHERE
1222
1223	ENDDO
1224
1225	!
1226	! 7. Plant uptake of ionic nitrogen
1227	!
1228	! This section comes from Section 3 of the supporting material of Zaehle & Friend (2010).
1229	! Only ammonium and nitrate are uptaken by plants.
1230	!
1231	! Comment from OCN : Temperature function of uptake is similar to SOM decomposition
1232	! to avoid N accumulation at low temperatures
1233	! In addition in the SM of Zaehle & Friend, 2010 P. 3 it is mentioned
1234	! "The uptake rate is observed to be sensitive to root temperature which is included as f(T),
1235	! thereby following the temperature sensitivity of net N mineralization"
1236	temp_sol_K=temp_sol(:)+tp_00
1237	ft_uptake(:) = control_temp_func (npts, temp_sol_K)
1238
1239
1240	! Comment of N. Vuichard: Plenty of parameter values used in the formulation of the plant uptake are
1241	! not traceable to anything. I tend to rely much of the param values to the values reported
1242	! in the reference publications
1243
1244	! Vmax of nitrogen uptake (in umol (g DryWeight_root)-1 h-1)
1245	!
1246	! In OCN, the same values are used both for uptake of NH4+ and NO3-
1247	! See p. 3 of the SM of Zaehle & Friend, 2010:
1248	! "As a first approximation average values for vmax, kNmin and KNmin are assumed for all PFTs
1249	! (Table S1), and for both ammonium and nitrate"
1250	! However based on the two papers of Kronzucker et al. (1995, 1996), it seems that vmax should be
1251	! much higher for NH4+ than for NO3- . We still use the same here for both NH4+ and NO3- as in OCN
1252	! The value of Vmax reported in Table S1 of Zaehle & Friend, 2010 is 5.14 ugN (g-1C) d-1
1253	! Comment of N. Vuichard : I can not relate the value of vmax in OCN expressed in
1254	! (umol (g DryWeight_root)-1 h-1) (vmax=3) to the one reported in Table S1 : 5.14 ugN (g-1C) d-1
1255	! The use of the conv_fac conversion factor used in OCN should help to convert (umol (g DryWeight_root)-1 h-1)
1256	! in (gN (g-1C) tstep-1). conv_fac is defined as 24. * dt / 2. / 1000000. * 14.
1257	! 24 conversion of hour to day
1258	! dt conversion of day to dt
1259	! 1/2 conversion of g(DryWeight) to gC
1260	! 1000000. conversion of ug to g
1261	! 14 conversion of umol to ugN
1262	! So using conv_fac, vmax expressed in ugN (g-1C) d-1 should be equal to
1263	! 324./2.14. = 504 ugN (gC)-1 d-1 not 5.14
1264	! In addition, I think there is an error in the conversion from (gDW)-1 to (gC)-1
1265	! When expressed per gC of root, the N uptake should be twice more than the one expressed per gDW of root,
1266	! not twice less. So one should multiply by 2., not divide by 2. (so vmax = 2016 ug (gC)-1 d-1 )
1267	!
1268	! Vmax of nitrogen uptake (in umol (g DryWeight_root)-1 h-1)
1269	! vmax_uptake(iammonium) = 3.
1270	! vmax_uptake(initrate) = 3.
1271
1272	! Conversion from (umol (gDW)-1 h-1) to (gN (gC)-1 timestep-1)
1273	conv_fac_vmax= 24. * dt * 2. / 1000000. * 14.
1274
1275	! minimal and maximal NC ratio of leaf (gN / gC)
1276	! used in the response of N uptake to NC ratio
1277	! see eq. 10 , p. 3 of SM of Zaehle & Friend, 2010
1278	nc_leaf_max(:,1) = 0.
1279	nc_leaf_max(:,2:nvm) = 1. / cn_leaf_min_2D(:,2:nvm)
1280
1281	nc_leaf_min(:,1) = -1./min_stomate
1282	nc_leaf_min(:,2:nvm) = 1. / cn_leaf_max_2D(:,2:nvm)
1283
1284
1285	DO m = 2, nvm
1286	! NCplant (gN / gC)
1287	! See equation (9) p. 3 of SM of Zaehle & Friend, 2010
1288	lab_n(:) = biomass(:,m,ilabile,initrogen) + &
1289	biomass(:,m,iroot,initrogen) + biomass(:,m,ileaf,initrogen)
1290	lab_c(:) = biomass(:,m,ilabile,icarbon) + &
1291	biomass(:,m,iroot,icarbon) + biomass(:,m,ileaf,icarbon)
1292	WHERE(lab_c(:).GT.min_stomate)
1293	NCplant(:) = lab_n(:) / lab_c(:)
1294	ELSEWHERE
1295	NCplant(:) = fcn_root(m)/cn_leaf_init_2D(:,m)
1296	ENDWHERE
1297	! Nitrogen demand responds to N/C ratio of the labile pool for each PFT
1298	! zero uptake at a NC in the labile pool corresponding to maximal leaf NC
1299	! and 1. for a NC corresponding to a minimal leaf NC
1300	! See equation (10) p. 3 of SM of Zaehle & Friend, 2010
1301	f_NCplant(:)=min(max( ( nc_leaf_max(:,m) - NCplant(:) ) / ( nc_leaf_max(:,m) - nc_leaf_min(:,m) ), 0. ), 1.)
1302
1303
1304	! Plant N uptake (gN m-2 tstep-1)
1305	!
1306	! See eq. (8) p. 3 of SM of Zaehle & Friend, 2010
1307	! In OCN, there was a multiplicative factor "2" that I can not explain
1308	! It may partly compensate the error mentioned above about conv_fac_vmax
1309	!
1310	! Coefficient K_N_min (umol per litter)
1311	!
1312	! It corresponds to the [NH4+] (resp. [NO3-]) for which the Nuptake equals vmax/2.
1313	! Kronzucker, 1995 reports values that range between 20 and 40 umol for NH4+ uptake
1314	! OCN seems to use 30 umol for both NH4+ and NO3-. The value used is 0.84 expressed in (gN m-2)
1315	! In Kronzucker, 1995 and 1996, the concentrations are always expressed in umol. No clear
1316	! reference about the standard volume it is related to (litter ?)
1317	! Coefficient K_N_min (umol per litter)
1318	! K_N_min(:) = 30
1319	! Conversion factor from (umol per litter) to (gN m-2)
1320	conv_fac_concent = 14.0 * 1e3 * 1e-6 * 2.0
1321	! 14 : molar mass for N (gN mol-1)
1322	! 10^3 : conversion factor (dm3 to m3)
1323	! 10-6 : conversion factor (ug to g)
1324	! 2 : m3 per m2 of soil
1325	!
1326	! Comment of N. Vuichard : I wonder why it assumes 2 m3 per m2 of soil. The depth of the soil
1327	! is 2 meter. But it doesn't mean that all the column contains only water, there is also soil, no ?
1328	! The question is : what is the volume that is considered for the concentration of NH4+ and NO3-. Is it a
1329	! volume of soil or of solution ?
1330	! If it is a volume of solution, I would suggest to multiply by a factor corresponding to the relative volume
1331	! of water within the soil. - Not done yet.
1332	!
1333	! K_N_min * conv_fac_concent = 0.84
1334	!
1335	! Coefficient low_K_N_min ((gN m-2)-1)
1336	!
1337	! See eq. 8 of SM of Zaehle et al. (2010) and Table S1
1338	! In table S1, it is defined as
1339	! "Rate of N uptake not associated with Michaelis- Menten Kinetics" with a value of 0.05
1340	! It is mentioned also (unitless) but to my opinion (N. Vuichard) it should have
1341	! the same unit that 1/K_N_min or 1/N_min, so ((gN m-2)-1)
1342	! Comment of N. Vuichard --------------------------------------------------------------
1343	! So far, I cannot relate the value of low_K_N_min (0.05) to any reference
1344	! especially Kronzucker (1996)
1345	! If I refer to Figure 4 of Kronzucker showing the NH4+ influx as a function of NH4+
1346	! concentration, the slope of the relationship could be used to define Vmax*low_K_N_min
1347	! For a concentration of NH4+ of 50 mmol, the influx equals 35 umol g-1 h-1
1348	! For a concentration of NH4+ of 20 mmol, the influx equals 17 umol g-1 h-1
1349	! slope = 10-3 * (35 - 17) / (50 - 20) = 10-3 * 18 / 30 = 0.0006 g-1 h-1
1350	! low_K_N_min = slope / Vmax = 0.0006 / 3 = 0.0002 (umol)-1
1351	! using the Conversion factor from (umol per litter) to (gN m-2)
1352	! conv_fac_concent = 14 * 1e3 * 1e-6 * 2, one should get
1353	! low_K_N_min = 0.0002 / ( 14 * 1e3 * 1e-6 * 2 ) = 0.007 ((gN m-2)-1)
1354	! The value of 0.007 does not match with the one in OCN (0.05) - This needs to be clarified
1355	! End of comment -----------------------------------------------------------------------
1356	! low_K_N_min ((umol)-1)
1357
1358	! In eq. 8, I (N. Vuichard) think there is an error. It should not be
1359	! (Nmin X KNmin) but (Nmin + KNmin). The source code of OCN is correct
1360	plant_uptake(:,m,iammonium) = vmax_uptake(iammonium) * conv_fac_vmax * croot_longterm(:,m) &
1361	* ft_uptake(:) * soil_n_min(:,m,iammonium) * ( low_K_N_min(iammonium) / conv_fac_concent &
1362	+ 1. / ( soil_n_min(:,m,iammonium) + K_N_min(iammonium) * conv_fac_concent ) ) * f_NCplant(:)
1363
1364	plant_uptake(:,m,initrate) = vmax_uptake(initrate) * conv_fac_vmax * croot_longterm(:,m) &
1365	* ft_uptake(:) * soil_n_min(:,m,initrate) * ( low_K_N_min(initrate) / conv_fac_concent &
1366	+ 1. / ( soil_n_min(:,m,initrate) + K_N_min(initrate) * conv_fac_concent ) ) * f_NCplant(:)
1367
1368
1369	IF ((printlev>=4).AND.(m==test_pft)) THEN
1370	WRITE(numout,*) 'plant_uptake ',plant_uptake(test_grid,test_pft,iammonium)
1371	WRITE(numout,*) 'vmax_uptake(iammonium) ',vmax_uptake(iammonium)
1372	WRITE(numout,*) ' biomass(:,m,iroot,icarbon)', biomass(test_grid,m,iroot,icarbon)
1373	WRITE(numout,*) ' ft_uptake(test_grid)',ft_uptake(test_grid)
1374	WRITE(numout,*) 'soil_n_min ',soil_n_min(test_grid,test_pft,iammonium)
1375	WRITE(numout,*) 'f_NCplant(test_grid)', f_NCplant(test_grid)
1376	WRITE(numout,*) 'NCplant(test_grid)', NCplant(test_grid)
1377	WRITE(numout,*) 'nc_leaf_max(m)',nc_leaf_max(test_grid,m)
1378	WRITE(numout,*) 'nc_leaf_min(m)',nc_leaf_min(test_grid,m)
1379	WRITE(numout,*) 'lab_c(test_grid)',lab_c(test_grid)
1380	WRITE(numout,*) 'lab_n(test_grid)',lab_n(test_grid)
1381	WRITE(numout,*) 'leaching ',leaching(test_grid,test_pft,iammonium)
1382	WRITE(numout,*) 'nitrification ',nitrification(test_grid,test_pft,:)
1383	WRITE(numout,*) 'tmc_pft',tmc_pft(test_grid,test_pft)
1384	ENDIF
1385
1386
1387	plant_uptake(:,m,iammonium) = MIN(soil_n_min(:,m,iammonium) - &
1388	(leaching(:,m,iammonium) + nitrification(:,m,i_nh4_to_no3) + nitrification(:,m,i_nh4_to_no) + &
1389	nitrification(:,m,i_nh4_to_n2o)), plant_uptake(:,m,iammonium))
1390
1391	plant_uptake(:,m,initrate) = MIN(soil_n_min(:,m,initrate) - &
1392	leaching(:,m,initrate) - denitrification(:,m,i_no3_to_nox) + nitrification(:,m,i_nh4_to_no3), &
1393	plant_uptake(:,m,initrate))
1394
1395	ENDDO
1396
1397	!
1398	! 8. Loss of N through drainage and gaseous emission
1399	! (the emission part should eventually be calculated in Sechiba's diffuco routines)
1400	!
1401
1402	DO m = 1,nvm
1403
1404	! 8.1 Loss of NH4 due to evaporation of NH3
1405	! NH4+ is first converted to NH3, and then it evaporates from the soil.
1406	! These equations come from Li et al. (1992), Table 4, in
1407	! addition to Appendix A of Zhang et al. (2002)
1408
1409	! Current dissociation of [NH3] to [NH4+]
1410	! See Table 4 of Li et al. 1992 and Appendix A of Zhang et al. 2002
1411	! log(K_NH4) - log(K_H20) = log(NH4/NH3) + pH
1412	! See also formula in "DISSOCIATION CONSTANTS OF INORGANIC ACIDS AND BASES" pdf file
1413	! pK_H2O = -log(K_H2O) = 14
1414	! pk_NH4 = -log(K_NH4) = 9.25
1415	! pK_H2O - pK_NH4 = log([NH4+]/[NH3]) + pH
1416	! [NH4+]/[NH3] = 1O^(pK_H2O - pK_NH4 - pH) = 10^(4.75-pH)
1417
1418	! In OCN, one makes use of frac_nh3. That should be the NH3/NH4 ratio. It is defined as:
1419	! frac_nh3(:) = 10.0(4.25-pH(:)) / (1. + 10.0(4.25-pH(:)))
1420
1421	! CommentS of N. Vuichard : I have several interrogations about the equation and value used
1422	! in OCN. But I have also concerns about the formulation of Li et al. ...
1423	! 1/
1424	! The formulation of Li et al. doesn't match with the formulas in
1425	! "DISSOCIATION CONSTANTS OF INORGANIC ACIDS AND BASES" or with the formulas
1426	! of http://www.onlinebiochemistry.com/obj-512/Chap4-StudNotes.html
1427	! To my opinion, one should replace log(K_NH4+) by log(K_NH3) in the formulation of Li et al.
1428	! with the relationship pKa + pKb = pKwater (where a and b are acid and base)
1429	! this leads to [NH4+]/[NH3] = 10^(pK_NH4 - pH) = 10^(9.25 - pH)
1430	! 2/
1431	! Wether I'm right or not about 1/, I don't understand the value used in OCN (4.25).
1432	! It should be either 4.75 or 9.25 but 4.25 looks strange
1433	! 3/
1434	! OCN used a formulation for [NH3]/[NH4+] of the type: X/(1+X) with X=[NH3]/[NH4+]
1435	! This leads to X/(1+X)=([NH3]/[NH4+])/(([NH4+]/[NH4+])+([NH3]/[NH4+]))
1436	! or X/(1+X) = [NH3] / ( [NH3] + [NH4+] )
1437	! This means that the value stored in soil_n_min(:,:,iammonium) corresponds to the total
1438	! N of both [NH4+] and [NH3]. This makes sense to my opinion. But I wonder if one should not
1439	! account for this partitioning between [NH4+] and [NH3] in other processes. To check.
1440	! 4/
1441	! The X value should relate to [NH3]/[NH4+] but to my opinion the X value used
1442	! in the equation in OCN corresponds to [NH4+]/[NH3]. Is this a bug ?
1443
1444	! In conclusion, I propose the formulation
1445	frac_nh3(:) = 10.0*(0.15(pH(:)-pk_NH4)) / (1. + 10.0*(0.15(pH(:)-pk_NH4)))
1446
1447	! This seems a patch added to OCN in order to (as mentioned in OCN)
1448	! reduced emissions at low concentration (problem of only one soil layer)
1449	! high conentrations are usually associated with fertiliser events -> top layer
1450	! and thus increased emission
1451	! NOT ACTIVATE HERE
1452	! emm_fac(:) = 0.01
1453	! WHERE(soil_n_min(:,m,iammonium).GT.0.01.AND.soil_n_min(:,m,iammonium).LE.4.)
1454	! emm_fac(:) = MAX(0.01,1-exp(-(soil_n_min(:,m,iammonium)/1.75)**8))
1455	! ENDWHERE
1456	! WHERE(soil_n_min(:,m,iammonium).GT.4.)
1457	! emm_fac(:)=1.0
1458	! ENDWHERE
1459
1460
1461	! 8.2 Volatilisation of gasous species, Table 4, Li et al. 2000 I,
1462	! using diffusivity of oxigen in air as a surrogate (from OCN)
1463	! assumes no effect of air concentration on diffusion
1464	! takes standard depth as reference
1465
1466	! Clay limitation
1467	! F_clay_0 = 0.13
1468	! F_clay_1 = -0.079
1469	F_clay(:) = F_clay_0 + F_clay_1 * clay(:)
1470
1471	! In OCN, one used formulation of the type
1472	! emission(:,m,inox-1) = &
1473	! MIN( d_ox(:) * soil_n_min(:,m,inox) * (0.13-0.079clay(:)) dt / z_decomp
1474	! It should have the unit d_ox * soil_n_min * dt / z_decomp
1475	! m2 day-1 * gN m-2 * day / m
1476	! gN / m which is not homogeneous with the unit expected (gn m-2)
1477	! To my opinion, one should not use soil_n_min (gN m-2) but a volumetric concentration (gN m-3)
1478	! But I'm not clear what is the appropriate volume to consider (volume of soil ?)
1479
1480	! NH4 emission (gN m-2 per time step)
1481	emission(:,m,iammonium) = D_s(:,m) * emm_fac * frac_nh3(:) * soil_n_min(:,m,iammonium) / zmaxh &
1482	* F_clay(:) / z_decomp * dt
1483
1484	! NO NO3 emission
1485	emission(:,m,initrate) = 0.
1486
1487	! NO2 emission (gN m-2 per time step)
1488	emission(:,m,inox) = D_s(:,m) * soil_n_min(:,m,inox) / zmaxh * F_clay(:) / z_decomp * dt
1489
1490	! N2O emission (gN m-2 per time step)
1491	emission(:,m,initrous) = D_s(:,m) * soil_n_min(:,m,initrous) / zmaxh * F_clay(:) / z_decomp * dt
1492
1493	! N2 emission (gN m-2 per time step)
1494	emission(:,m,idinitro) = D_s(:,m) * soil_n_min(:,m,idinitro) / zmaxh * F_clay(:) / z_decomp * dt
1495
1496	ENDDO
1497	emission(:,:,iammonium) = MIN(emission(:,:,iammonium), &
1498	soil_n_min(:,:,iammonium) - nitrification(:,:,i_nh4_to_no3) &
1499	- nitrification(:,:,i_nh4_to_no) - nitrification(:,:,i_nh4_to_n2o) &
1500	- leaching(:,:,iammonium))
1501
1502	emission(:,:,inox) = MIN(emission(:,:,inox), &
1503	soil_n_min(:,:,inox) + nitrification(:,:,i_nh4_to_no) &
1504	+ denitrification(:,:,i_no3_to_nox) - denitrification(:,:,i_nox_to_n2o))
1505
1506	emission(:,:,initrous) = MIN(emission(:,:,initrous), &
1507	soil_n_min(:,:,initrous) + nitrification(:,:,i_nh4_to_n2o) &
1508	+ denitrification(:,:,i_nox_to_n2o) - denitrification(:,:,i_n2o_to_n2))
1509
1510	emission(:,:,idinitro) = MIN(emission(:,:,idinitro), &
1511	soil_n_min(:,:,idinitro) + denitrification(:,:,i_n2o_to_n2))
1512
1513	!
1514	! 9. Update pools
1515	!
1516
1517	! 9.1 update pools of nitrogen in the soil from nitrification and
1518	! denitrification, plant uptake, leaching and volatile emissions,
1519	! desorption from clay, and net mineralisation
1520	!
1521	! In my opinion, for a better consistency, I would recommand to consider the leaching separately
1522	! when calculating the ammonium and nitrate budget. Leaching is calculated from sechiba
1523	! Might be better to remove leaching at the top of the routine especially due to the later calculation of
1524	! NH4+ and NO3- concentration that will vary with the soil water content
1525	! Let's imagine that from one time step to another, the change in soil water content is only due to leaching
1526	! We don't want that the NH4+ and NO3- concentration vary from one time step to the other. The best way to avoid
1527	! this is to remove first the leaching from the NH4+ and NO3- pools
1528	! THIS IS NOT DONE YET
1529
1530	IF(printlev>=4)THEN
1531	WRITE(numout,*) 'CHECK values before update'
1532	WRITE(numout,*) 'nitrification ',nitrification(test_grid,test_pft,:)
1533	WRITE(numout,*) 'denitrification ',denitrification(test_grid,test_pft,:)
1534	WRITE(numout,*) 'leaching ',leaching(test_grid,test_pft,:)
1535	WRITE(numout,*) 'emission ',emission(test_grid,test_pft,:)
1536	WRITE(numout,*) 'plant_uptake ',plant_uptake(test_grid,test_pft,:)
1537	WRITE(numout,*) 'mineralisation ',mineralisation(test_grid,test_pft)
1538	WRITE(numout,*) 'immob ',immob(test_grid,test_pft)
1539	ENDIF
1540
1541	soil_n_min(:,:,iammonium) = soil_n_min(:,:,iammonium) + n_adsorbed(:,:) &
1542	- nitrification(:,:,i_nh4_to_no3) - nitrification(:,:,i_nh4_to_no) - nitrification(:,:,i_nh4_to_n2o) &
1543	- leaching(:,:,iammonium) - emission(:,:,iammonium) &
1544	+ mineralisation(:,:) + immob(:,:) - plant_uptake(:,:,iammonium)
1545
1546	soil_n_min(:,:,initrate) = soil_n_min(:,:,initrate) + nitrification(:,:,i_nh4_to_no3) &
1547	- denitrification(:,:,i_no3_to_nox) - leaching(:,:,initrate) &
1548	- plant_uptake(:,:,initrate)
1549
1550	soil_n_min(:,:,inox) = soil_n_min(:,:,inox) + nitrification(:,:,i_nh4_to_no) &
1551	+ denitrification(:,:,i_no3_to_nox) - denitrification(:,:,i_nox_to_n2o) - emission(:,:,inox)
1552
1553	soil_n_min(:,:,initrous) = soil_n_min(:,:,initrous) + nitrification(:,:,i_nh4_to_n2o) &
1554	+ denitrification(:,:,i_nox_to_n2o) - denitrification(:,:,i_n2o_to_n2) - emission(:,:,initrous)
1555
1556	soil_n_min(:,:,idinitro) = soil_n_min(:,:,idinitro) + denitrification(:,:,i_n2o_to_n2) &
1557	- emission(:,:,idinitro)
1558
1559
1560	IF(printlev>=4)THEN
1561	WRITE(numout,*) 'CHECK values after update'
1562	WRITE(numout,*) 'soil_n_min ',soil_n_min(test_grid,test_pft,:)
1563	ENDIF
1564	! 9.2 add nitrogen either from deposition (read from fields or prescribed in run.def) as
1565	! well as BNF, needs to be revised...
1566	! iatm_ammo=1
1567	! iatm_nitr=2
1568	! ibnf=3
1569	! ifert=4
1570	DO m=1,nvm
1571	! Deposition of NHx and NOy
1572	WHERE(veget_max(:,m).GT.min_stomate.AND.som(:,iactive,m,icarbon).GT.min_stomate)
1573	soil_n_min(:,m,iammonium) = soil_n_min(:,m,iammonium) &
1574	+ input(:,m,iammonium)*dt
1575	soil_n_min(:,m,initrate) = soil_n_min(:,m,initrate) &
1576	+ input(:,m,initrate)*dt
1577	ENDWHERE
1578	WHERE(veget_max(:,m).GT.min_stomate.AND.som(:,iactive,m,icarbon).LE.min_stomate)
1579	leaching(:,m,iammonium)=leaching(:,m,iammonium) + &
1580	input(:,m,iammonium)*dt
1581	leaching(:,m,initrate)=leaching(:,m,initrate) + &
1582	input(:,m,initrate)*dt
1583	ENDWHERE
1584	! BNF
1585	IF ( natural(m) ) THEN
1586	! in the presence of a organic soil component assume that there is also BNF
1587	! as long as plant available nitrogen is not too high
1588	WHERE(som(:,iactive,m,icarbon).GT.min_stomate.AND. &
1589	soil_n_min(:,m,iammonium)+soil_n_min(:,m,initrate).LT.max_soil_n_bnf(m))
1590	soil_n_min(:,m,iammonium) = soil_n_min(:,m,iammonium) + input(:,m,ibnf)*dt
1591	ENDWHERE
1592	ENDIF
1593	! Fertiliser use for agriculture, no BNF, that's already accounted for in fertil
1594	! using the average global ratio of ammonium to nitrate to
1595	! separate the two species
1596	! ratio_nh4_fert = 7./8.
1597	WHERE(veget_max(:,m).GT.min_stomate)
1598	soil_n_min(:,m,iammonium) = soil_n_min(:,m,iammonium) &
1599	+ input(:,m,ifert)ratio_nh4_fertdt
1600	soil_n_min(:,m,initrate) = soil_n_min(:,m,initrate) &
1601	+ input(:,m,ifert)(1.-ratio_nh4_fert)dt
1602	ENDWHERE
1603	ENDDO
1604
1605	IF(printlev>=4)THEN
1606	WRITE(numout,*) 'CHECK values after N input'
1607	WRITE(numout,*) 'soil_n_min ',soil_n_min(test_grid,test_pft,:)
1608	ENDIF
1609
1610	! 10. Write output values
1611	CALL histwrite_p (hist_id_stomate, 'N_UPTAKE_NH4', itime, &
1612	plant_uptake(:,:,iammonium)/dt, npts*nvm, horipft_index)
1613	CALL histwrite_p (hist_id_stomate, 'N_UPTAKE_NO3', itime, &
1614	plant_uptake(:,:,initrate)/dt, npts*nvm, horipft_index)
1615	CALL histwrite_p (hist_id_stomate, 'N_MINERALISATION', itime, &
1616	mineralisation(:,:)/dt, npts*nvm, horipft_index)
1617	CALL histwrite_p (hist_id_stomate, 'SOIL_NH4', itime, &
1618	soil_n_min(:,:,iammonium), npts*nvm, horipft_index)
1619	CALL histwrite_p (hist_id_stomate, 'SOIL_NO3', itime, &
1620	soil_n_min(:,:,initrate), npts*nvm, horipft_index)
1621	CALL histwrite_p (hist_id_stomate, 'SOIL_NOX', itime, &
1622	soil_n_min(:,:,inox), npts*nvm, horipft_index)
1623	CALL histwrite_p (hist_id_stomate, 'SOIL_N2O', itime, &
1624	soil_n_min(:,:,initrous), npts*nvm, horipft_index)
1625	CALL histwrite_p (hist_id_stomate, 'SOIL_N2', itime, &
1626	soil_n_min(:,:,idinitro), npts*nvm, horipft_index)
1627	CALL histwrite_p (hist_id_stomate, 'SOIL_P_OX', itime, &
1628	p_O2(:,:), npts*nvm, horipft_index)
1629	CALL histwrite_p (hist_id_stomate, 'BACT', itime, &
1630	bact(:,:), npts*nvm, horipft_index)
1631	CALL histwrite_p (hist_id_stomate, 'NH3_EMISSION', itime, &
1632	emission(:,:,iammonium)/dt, npts*nvm, horipft_index)
1633	CALL histwrite_p (hist_id_stomate, 'NOX_EMISSION', itime, &
1634	emission(:,:,inox)/dt, npts*nvm, horipft_index)
1635	CALL histwrite_p (hist_id_stomate, 'N2O_EMISSION', itime, &
1636	emission(:,:,initrous)/dt, npts*nvm, horipft_index)
1637	CALL histwrite_p (hist_id_stomate, 'N2_EMISSION', itime, &
1638	emission(:,:,idinitro)/dt, npts*nvm, horipft_index)
1639	CALL histwrite_p (hist_id_stomate, 'NH4_LEACHING', itime, &
1640	leaching(:,:,iammonium)/dt, npts*nvm, horipft_index)
1641	CALL histwrite_p (hist_id_stomate, 'NO3_LEACHING', itime, &
1642	leaching(:,:,initrate)/dt, npts*nvm, horipft_index)
1643	CALL histwrite_p (hist_id_stomate, 'NITRIFICATION', itime, &
1644	nitrification(:,:,i_nh4_to_no3), npts*nvm, horipft_index)
1645	CALL histwrite_p (hist_id_stomate, 'DENITRIFICATION', itime, &
1646	denitrification(:,:,i_n2o_to_n2), npts*nvm, horipft_index)
1647	CALL histwrite_p (hist_id_stomate, 'NHX_DEPOSITION', itime, &
1648	input(:,:,iammonium), npts*nvm, horipft_index)
1649	CALL histwrite_p (hist_id_stomate, 'NOX_DEPOSITION', itime, &
1650	input(:,:,initrate), npts*nvm, horipft_index)
1651	CALL histwrite_p (hist_id_stomate, 'BNF', itime, &
1652	input(:,:,ibnf), npts*nvm, horipft_index)
1653	CALL histwrite_p (hist_id_stomate, 'N_FERTILISER', itime, &
1654	input(:,:,ifert), npts*nvm, horipft_index)
1655	CALL histwrite_p (hist_id_stomate, 'N_MANURE', itime, &
1656	input(:,:,imanure), npts*nvm, horipft_index)
1657
1658	CALL xios_orchidee_send_field("N_UPTAKE_NH4",plant_uptake(:,:,iammonium)/dt)
1659	CALL xios_orchidee_send_field("N_UPTAKE_NO3",plant_uptake(:,:,initrate)/dt)
1660	CALL xios_orchidee_send_field("N_MINERALISATION",mineralisation(:,:)/dt)
1661	CALL xios_orchidee_send_field("SOIL_NH4",soil_n_min(:,:,iammonium))
1662	CALL xios_orchidee_send_field("SOIL_NO3",soil_n_min(:,:,initrate))
1663	CALL xios_orchidee_send_field("SOIL_NOX",soil_n_min(:,:,inox))
1664	CALL xios_orchidee_send_field("SOIL_N2O",soil_n_min(:,:,initrous))
1665	CALL xios_orchidee_send_field("SOIL_N2",soil_n_min(:,:,idinitro))
1666	CALL xios_orchidee_send_field("SOIL_P_OX",p_O2(:,:))
1667	CALL xios_orchidee_send_field("BACT",bact(:,:))
1668	CALL xios_orchidee_send_field("NH3_EMISSION",emission(:,:,iammonium)/dt)
1669	CALL xios_orchidee_send_field("NOX_EMISSION",emission(:,:,inox)/dt)
1670	CALL xios_orchidee_send_field("N2O_EMISSION",emission(:,:,initrous)/dt)
1671	CALL xios_orchidee_send_field("N2_EMISSION",emission(:,:,idinitro)/dt)
1672	CALL xios_orchidee_send_field("NH4_LEACHING",leaching(:,:,iammonium)/dt)
1673	CALL xios_orchidee_send_field("NO3_LEACHING",leaching(:,:,initrate)/dt)
1674	CALL xios_orchidee_send_field("NITRIFICATION",nitrification(:,:,i_nh4_to_no3))
1675	CALL xios_orchidee_send_field("DENITRIFICATION",denitrification(:,:,i_n2o_to_n2))
1676	CALL xios_orchidee_send_field("NHX_DEPOSITION",input(:,:,iammonium))
1677	CALL xios_orchidee_send_field("NOX_DEPOSITION",input(:,:,initrate))
1678	CALL xios_orchidee_send_field("BNF",input(:,:,ibnf))
1679	CALL xios_orchidee_send_field("N_FERTILISER",input(:,:,ifert))
1680	CALL xios_orchidee_send_field("N_MANURE",input(:,:,imanure))
1681
1682	CALL xios_orchidee_send_field("fBNF",SUM(input(:,:,ibnf)*veget_cov_max,dim=2)/1e3/one_day)
1683	CALL xios_orchidee_send_field("fNdep",SUM((input(:,:,iammonium)+input(:,:,initrate))*veget_cov_max,dim=2)/1e3/one_day)
1684	CALL xios_orchidee_send_field("fNfert",SUM((input(:,:,ifert)+input(:,:,imanure))*veget_cov_max,dim=2)/1e3/one_day)
1685	CALL xios_orchidee_send_field("fNgas",SUM((emission(:,:,iammonium)+emission(:,:,inox)+emission(:,:,initrous)+emission(:,:,idinitro))*veget_cov_max,dim=2)/dt/1e3/one_day)
1686	CALL xios_orchidee_send_field("fN2O",SUM(emission(:,:,initrous)*veget_cov_max,dim=2)/dt/1e3/one_day)
1687	CALL xios_orchidee_send_field("fNOx",SUM(emission(:,:,inox)*veget_cov_max,dim=2)/dt/1e3/one_day)
1688	CALL xios_orchidee_send_field("fNleach",SUM((leaching(:,:,iammonium)+leaching(:,:,initrate))*veget_cov_max,dim=2)/dt/1e3/one_day)
1689	CALL xios_orchidee_send_field("fNnetmin",SUM((mineralisation(:,:)+immob(:,:))*veget_cov_max,dim=2)/dt/1e3/one_day)
1690	CALL xios_orchidee_send_field("fNup",SUM((plant_uptake(:,:,iammonium)+plant_uptake(:,:,initrate))*veget_cov_max,dim=2)/dt/1e3/one_day)
1691	CALL xios_orchidee_send_field("nMineral",SUM((soil_n_min(:,:,iammonium)+soil_n_min(:,:,initrate)+soil_n_min(:,:,inox)+&
1692	soil_n_min(:,:,initrous)+soil_n_min(:,:,idinitro))*veget_cov_max,dim=2)/1e3)
1693	CALL xios_orchidee_send_field("nMineralNH4",SUM(soil_n_min(:,:,iammonium)*veget_cov_max,dim=2)/1e3)
1694	CALL xios_orchidee_send_field("nMineralNO3",SUM(soil_n_min(:,:,initrate)*veget_cov_max,dim=2)/1e3)
1695
1696
1697	END SUBROUTINE nitrogen_dynamics
1698
1699
1700
1701	!! ================================================================================================================================
1702	!! FUNCTION : control_temp_func
1703	!!
1704	!>\BRIEF : Unclear.
1705	!!
1706	!! DESCRIPTION : Unable to find where this comes from.
1707	!! Referenced by Zaehle and Friend (2010) in the Appendix,
1708	!! which points to Krinner et al (2005), but the closest I found
1709	!! there was Eq. A32, which is simillar but does not
1710	!! mention Q10 at all.
1711	!!
1712	!! RECENT CHANGE(S):
1713	!!
1714	!! MAIN OUTPUTS VARIABLE(S):
1715	!!
1716	!! REFERENCE(S) :
1717	!! - S. Zaehle and A. D. Friend (2010), Carbon and nitrogen cycle dynamics in the
1718	!! O-CN land surface model: 1. Model description, site-scale evaluation, and
1719	!! sensitivity to parameter estimates. Global Biogeochem. Cycles, 24, GB1005,
1720	!! doi:10.1029/2009GB003521.
1721	!! - Krinner G., N. Viovy, N. de Noblet-Ducoudre, J. Ogee, J. Polcher, P. Friedlingstein,
1722	!! P. Ciais, S. Sitch, and I. C. Prentice (2005), A dynamic global vegetation model
1723	!! for studies of the coupled atmosphere-biosphere system, Global Biogeochemical
1724	!! Cycles, 19, doi.:10.1029/2003/GB002199.
1725	!!
1726	!! FLOWCHART :
1727	!!
1728	!_ ================================================================================================================================
1729	FUNCTION control_temp_func (npts, temp_in) RESULT (tempfunc_result)
1730
1731	!! 0. Variable and parameter declaration
1732
1733	!! 0.1 Input variables
1734	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of land pixels (unitless)
1735	REAL(r_std), DIMENSION(npts), INTENT(in) :: temp_in !! Temperature (K)
1736
1737	!! 0.2 Output variables
1738	REAL(r_std), DIMENSION(npts) :: tempfunc_result !! Temperature control factor (0-1, unitless)
1739
1740	!! 0.3 Modified variables
1741
1742	!! 0.4 Local variables
1743
1744	!_ ================================================================================================================================
1745
1746	tempfunc_result(:) = exp( soil_Q10_uptake * ( temp_in(:) - (ZeroCelsius+tsoil_ref)) / Q10 )
1747	tempfunc_result(:) = MIN( un, tempfunc_result(:) )
1748
1749	END FUNCTION control_temp_func
1750
1751
1752	END MODULE stomate_som_dynamics

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

Download in other formats: