Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-261/src_stomate/stomate_soilcarbon.f90 @ 8692

Last change on this file since 8692 was 4998, checked in by nicolas.vuichard, 7 years ago
rev29012018
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 80.1 KB

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

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

Download in other formats: