Context Navigation

stomate_stics.f90 @ 7108

Last change on this file since 7108 was 5850, checked in by albert.jornet, 6 years ago
Fix: replace hardcoded days of the year with the global time variable year_length_in_days New: variable year_length_in_days. It is calculated at every time step to allow Orchidee run multiple years in the same run.
File size: 83.5 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_stics
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Groups the subroutines which interface with or are related with sticslai
10	!
11	!!
12	!!\n DESCRIPTION : None
13	!!
14	!! RECENT CHANGE(S) : None
15	!!
16	!! REFERENCE(S) : None
17	!!
18	!! SVN :
19	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-MICT/ORCHIDEE/src_stomate/stomate.f90 $
20	!! $Date: 2017-07-17 16:22:36 +0200 (Mon, 17 Jul 2017) $
21	!! $Revision: 4509 $
22	!! \n
23	!_ ================================================================================================================================
24
25	MODULE stomate_stics
26
27	! Modules used:
28	USE netcdf
29	USE defprec
30	USE grid
31	USE constantes
32	USE constantes_soil
33	USE pft_parameters
34	USE ioipsl
35	USE ioipsl_para
36	USE mod_orchidee_para
37	USE interpol_help ! necessary for management map input
38	USE time, ONLY: year_length_in_days
39
40	IMPLICIT NONE
41
42	! Private & public routines
43	PRIVATE
44	PUBLIC stomate_stics_ManageInput, stomate_stics_read_cycle, stomate_stics_read_rotation, stomate_stics_read_plantdate, stomate_stics_get_cmd
45	PUBLIC stomate_stics_NfertInput, stomat_stics_calc_N_limfert, stomate_stics_rotation
46
47	! Public variables (not recommended)
48	!PUBLIC
49
50	CONTAINS
51	!
52	!! ================================================================================================================================
53	!! SUBROUTINE : stomate_forcing_CropVar
54	!!
55	!>\BRIEF Interpolate (extract) Crop Variety information for each crop type
56	!!
57	!! DESCRIPTION : None
58	!!
59	!! RECENT CHANGE(S) : None
60	!!
61	!! MAIN OUTPUT VARIABLE(S): None
62	!!
63	!! REFERENCES : None
64	!!
65	!! FLOWCHART : None
66	!! \n
67	!_ ================================================================================================================================
68	SUBROUTINE slowproc_CropVar(nbpt, lalo, neighbours, resolution, contfrac,& ! In
69	plantcyc, cropvar ) ! Out
70	!
71	!
72	! INTEGER(i_std) :: iv
73	! INTEGER(i_std),DIMENSION(nvm) :: plantdate_default = (/(0,iv=1,nvm)/)
74	! !It should be defined as an array for different crops
75	! !The value should be given outside this program and subject to crop type
76	! !if being defined outside, please comment the above sentence
77
78	! 0.1 INPUT
79	!
80	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated (extracted)
81	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longtitude
82	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vectors of neighbours for each grid point
83	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
84	!REAL(r_std),INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid box in lat and lon
85	REAL(r_std) :: resolution(nbpt,2) ! The size in km of each grid box in lat and lon
86	REAL(r_std),INTENT(in) :: contfrac(nbpt) ! The fraction of land in each grid box
87	! INTEGER(I_std),INTENT(in) :: plantcyc ! The time of the rotation
88	! in this year (0 means this info unknown)
89	! INTEGER(i_std), INTENT(in) :: croptype ! The type of crop being
90	! simulated
91	! not necessary as we read planting date data for all crop types
92	!
93	! 0.2 OUTPUT
94	!
95	REAL(r_std),ALLOCATABLE, INTENT(out) :: cropvar(:,:,:) ! The crop variety matrix
96	INTEGER(i_std), INTENT(out) :: plantcyc ! times of rotation that year
97	! nvm is the number of PFTs, there may not be planting date for all the PFTs
98	!
99	! 0.3 LOCAL
100	!
101	INTEGER(i_std) :: nbvmax ! a parameter for interpolation
102	CHARACTER(LEN=80) :: filename
103	INTEGER(i_std) :: iml, jml, lml, tml, fid, fid1
104	INTEGER(i_std) :: ip, jp, ib, ilf, fopt ! for-loop variable
105	INTEGER(i_std) :: nbexp
106	REAL(r_std) :: lev(1), date, dt
107	REAL(r_std) :: missing_val
108	INTEGER(i_std) :: itau(1)
109
110	INTEGER(i_std) :: nb_dim
111	INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w
112	LOGICAL :: l_ex
113
114	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel
115	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:,:) :: cropvar_mat ! LON LAT VEGET PLANTCYC TIME_COUNTER
116	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: cropvar_mat_1 ! LON LAT VEGET, PLANTCYC it is used when input file does not have time counter
117	! INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: plntdt_mat
118	! Because flinget does not support reading integer matrix, plntdt_mat has to
119	! be real
120	! or we will have to rewrite IOIPSL
121	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
122	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: temp_var
123	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
124	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
125
126	REAL(r_std) :: sgn, sum_float
127	INTEGER(i_std) :: ivgt, icyc, pltcyc, it
128	CHARACTER(LEN=30) :: callsign
129	LOGICAL :: ok_interpol
130	INTEGER :: ALLOC_ERR
131	INTEGER(i_std) :: cps,maxcp,jlf,maxps
132	INTEGER(i_std),ALLOCATABLE,DIMENSION(:) :: hashlst1
133	REAL(r_std),ALLOCATABLE,DIMENSION(:) :: hashlst2
134	REAL(r_std) :: areamax
135	LOGICAL ok_var
136
137
138	! croptype = TRIM(croptype) !if croptype is a string
139	! else a switch expression is needed
140	filename = "/work/cont003/p529tan/WXH/CropVar_1990.nc" ! default input file
141	! String operation needed
142	CALL getin_p('CROPVAR_FILE',filename)
143
144	IF (is_root_prc) THEN
145	! ? what does is_root_prc mean?
146	CALL flininfo(filename, iml, jml, lml, tml, fid)
147	! CALL flinclo(fid)
148	ENDIF
149	CALL bcast(iml)
150	CALL bcast(jml)
151	! CALL bcast(lml)
152	! CALL bcast(tml)
153
154	! Printing information for debugging
155	WRITE(numout, *) "Xuhui's debug info for slowproc_CropVar #1:"
156	WRITE(numout, *) "filename is: ", filename
157	WRITE(numout, *) "Dimension 1, lon, iml:", iml
158	WRITE(numout, *) "Dimension 2, lat, jml:", jml
159	WRITE(numout, *) "Dimension 3, veget, lml:", lml
160	WRITE(numout, *) "Dimension 4, time, tml:", tml
161	! apparently, flinget function is not designed to take veget but levels to
162	! be the
163	! 3rd dimension, modification to lml is needed
164	CALL flioopfd(filename,fid1)
165	CALL flioinqv(fid1,v_n="CROPVAR", l_ex = l_ex, nb_dims = nb_dim, len_dims =l_d_w)
166	IF (lml == 0) THEN
167	! CALL
168	! flioinqv(fid1,v_n="PLNTDT",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
169	lml=l_d_w(3)
170	WRITE(numout, *) "len_dims: ", l_d_w
171	WRITE(numout, *) "lml AFTER revision"
172	WRITE(numout, *) "lml: ", lml
173	ENDIF
174	IF (nb_dim.EQ.5) THEN
175	plantcyc = l_d_w(4)
176	tml = l_d_w(5)
177	WRITE(numout,*) 'No. of plantcycle in the input: ', plantcyc
178	WRITE(numout,*) 'tml changed to: ', tml
179	ELSE
180	IF (nb_dim.EQ.4) THEN
181	plantcyc = 0
182	ELSE
183	WRITE(numout,*) "Error: More axis in PlantDate than expected"
184	STOP
185	ENDIF
186	ENDIF
187	IF (plantcyc>3) THEN
188	WRITE(numout,*) "Rotation cycle more than accepted: ",plantcyc
189	ENDIF
190	WRITE(numout,*) "plantcyc: ",plantcyc
191	CALL flioclo(fid1)
192	CALL bcast(lml)
193	CALL bcast(tml)
194	! CALL bcast(plantcyc)
195
196	ALLOC_ERR=-1
197	ALLOCATE(cropvar(nbpt,lml,tml),STAT=ALLOC_ERR)
198	IF (ALLOC_ERR/=0) THEN
199	WRITE(numout,*) "ERROR IN ALLOCATION OF cropvar: ", ALLOC_ERR
200	ENDIF
201	WRITE(numout,*) "cropvar ALLOCATED"
202	!
203	!! assumption: the crop variety file is a 0.5 degree planting date (julian
204	!date) file
205	!
206	ALLOC_ERR=-1
207	ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
208	IF (ALLOC_ERR/=0) THEN
209	WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
210	STOP
211	ENDIF
212	ALLOC_ERR=-1
213	ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
214	IF (ALLOC_ERR/=0) THEN
215	WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
216	STOP
217	ENDIF
218	ALLOC_ERR=-1
219	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
220	IF (ALLOC_ERR/=0) THEN
221	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
222	STOP
223	ENDIF
224
225	IF (plantcyc == 0) THEN
226	ALLOC_ERR=-1
227	ALLOCATE(cropvar_mat_1(iml,jml,lml,tml), STAT=ALLOC_ERR)
228	! !lml is supposed to be nvm (number of PFTs), if not ,change it
229	IF (ALLOC_ERR/=0) THEN
230	WRITE(numout,*) "ERROR IN ALLOCATION of cropvar_mat_1 : ",ALLOC_ERR
231	STOP
232	ENDIF
233	ALLOC_ERR=-1
234	ALLOCATE(cropvar_mat(iml,jml,lml,1,tml), STAT=ALLOC_ERR)
235	IF (ALLOC_ERR/=0) THEN
236	WRITE(numout,*) "ERROR IN ALLOCATION of cropvar_mat: ", ALLOC_ERR
237	STOP
238	ENDIF
239	ELSE
240	ALLOC_ERR=-1
241	ALLOCATE(cropvar_mat(iml,jml,lml,plantcyc,tml), STAT=ALLOC_ERR)
242	IF (ALLOC_ERR/=0) THEN
243	WRITE(numout,*) "ERROR IN ALLOCATION of cropvar_mat: ", ALLOC_ERR
244	STOP
245	ENDIF
246	ENDIF
247
248	! input of some attributes
249	! IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel,
250	! lat_rel, lev, tml, itau, date, dt, fid)
251	! CALL bcast(lon_rel)
252	! CALL bcast(lat_rel)
253	! CALL bcast(itau)
254	! CALL bcast(date)
255	! CALL bcast(dt)
256	IF (is_root_prc) THEN
257	CALL flinget(fid, 'LON', iml, jml, lml, tml, 1, 1, lon_rel)
258	CALL flinget(fid, 'LAT', iml, jml, lml, tml, 1, 1, lat_rel)
259	CALL bcast(lon_rel)
260	CALL bcast(lat_rel)
261	! WRITE (numout,*) 'lon_rel size: ', SIZE(lon_rel)
262	! WRITE (numout,*) 'lat_rel size: ', SIZE(lat_rel)
263	ENDIF
264
265	! input of the matrix
266	IF (is_root_prc) THEN
267	! CALL flinget(fid, 'CROPVAR', iml, jml, lml, tml, 1, 1, crop_mat)
268	! time_counter has to be 1, or it will not match the size of plntdt_mat
269	CALL flioopfd(filename,fid1)
270	IF (plantcyc>3) THEN
271	WRITE(numout,*) "slowproc_CropVar : Planting cycle more than model accepted"
272	STOP
273	ELSE
274	IF (plantcyc > 0) THEN
275	CALL fliogetv(fid1,'CROPVAR',cropvar_mat,start=(/1,1,1,1,1/),count=(/iml,jml,lml,plantcyc,tml/))
276	ELSEIF (plantcyc == 0) THEN
277	CALL fliogetv(fid1,'CROPVAR',cropvar_mat_1,start=(/1,1,1,1/),count=(/iml,jml,lml,tml/))
278	ELSE
279	WRITE(numout,*) "invalid plantcyc axis : ", plantcyc
280	ENDIF
281	ENDIF
282	! get missing_val
283	CALL fliogeta(fid1,'CROPVAR','missing_value',missing_val)
284	CALL flioclo(fid1)
285	ENDIF
286	IF (plantcyc == 0) THEN
287	cropvar_mat(:,:,:,1,:) = cropvar_mat_1(:,:,:,:)
288	! CALL bcast(cropvar_mat)
289	DEALLOCATE(cropvar_mat_1)
290	ELSE
291	! CALL bcast(cropvar_mat)
292	ENDIF
293	! WRITE(numout,*) 'cropvar_mat size: ',SIZE(cropvar_mat)
294	! WRITE(numout,*) 'missing value: ', missing_val
295	! WRITE(numout,*) 'lat(361,284): ',lat_rel(361,284)
296	! WRITE(numout,*) 'lon(361,284): ',lon_rel(361,284)
297	! WRITE(numout,*) 'cropvar(361,284,1,1): ',cropvar_mat(361,284,1,1)
298
299	IF (is_root_prc) CALL flinclo(fid)
300
301	cropvar(:,:,:) = zero ! nbpt veget time
302	IF (plantcyc == 0) THEN
303	pltcyc = 1
304	ELSE
305	pltcyc = plantcyc
306	ENDIF
307
308	DO it = 1,tml
309	DO icyc = 1,1 ! pltcyc ! not dealing with plant cycles for now
310	DO ivgt = 1,lml ! ? We can suppose PFTs less than 10 are natural veg without planting date, but not now
311	IF (natural(ivgt)) THEN
312	WRITE(numout,*) "veget, plantcyc, time: ", ivgt,icyc,it
313	nbexp = 0
314	! the number of exceptions
315
316	! mask of available value
317	!
318	! mask is commented because it is of no use right now.
319	mask(:,:) = zero; ! Defined in constante.f90
320	DO ip = 1,iml
321	DO jp = 1,jml
322	IF ((cropvar_mat(ip,jp,ivgt,icyc,it) .GT. min_sechiba) .AND. &
323	(cropvar_mat(ip,jp,ivgt,icyc,it) /= missing_val)) THEN
324	mask(ip,jp) = un; ! Defined in constante.f90
325	! here we assumed that for each plant cycle at each
326	! there might be missing data at different grid
327	! in this case, mask has to be generated each plant
328	! cycle each time step
329	ENDIF
330	ENDDO
331	ENDDO
332
333	! Interpolation started
334	nbvmax = 200
335	! the maximum amount of fine grids that one coarse grid may have
336
337	callsign = "Crop Variety"
338
339	ok_interpol = .FALSE.
340
341	DO WHILE ( .NOT. ok_interpol )
342	WRITE(numout,*) "Pojection arrays for ", callsign, ":"
343	WRITE(numout,*) "nbvmax = ", nbvmax
344
345	ALLOC_ERR = -1
346	ALLOCATE(temp_var(nbvmax,lml), STAT=ALLOC_ERR)
347	IF (ALLOC_ERR /=0) THEN
348	WRITE(numout,*) "ERROR IN ALLOCATION OF temp_var :",ALLOC_ERR
349	STOP
350	ENDIF
351	ALLOC_ERR = -1
352	ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
353	IF (ALLOC_ERR /=0) THEN
354	WRITE(numout,*) "ERROR IN ALLOCATION OF sub_index :", ALLOC_ERR
355	STOP
356	ENDIF
357	sub_index(:,:,:) = zero
358	ALLOC_ERR = -1
359	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
360	IF (ALLOC_ERR /=0) THEN
361	WRITE(numout,*) "ERROR IN ALLOCATION OF sub_area :", ALLOC_ERR
362	STOP
363	ENDIF
364	sub_area(:,:) = zero
365	resolution(:,:) = resolution(:,:)/1000 ! m -> km
366	write(,) "resolution updated: ", resolution(1,:), "km"
367
368	CALL aggregate_p(nbpt, lalo, neighbours, resolution,contfrac, &
369	& iml, jml, lon_rel, lat_rel, mask, callsign, &
370	& nbvmax, sub_index, sub_area, ok_interpol)
371
372	IF ( .NOT. ok_interpol ) THEN
373	DEALLOCATE(temp_var)
374	DEALLOCATE(sub_index)
375	DEALLOCATE(sub_area)
376	nbvmax = nbvmax * 2
377	ENDIF
378	ENDDO
379
380	! assign the values to plantdate
381	! values should be given to all PFTs
382	DO ib = 1, nbpt
383	! cropvar(ib,:,icyc) = zero
384
385	! examing all sub_point we found
386	fopt = COUNT(sub_area(ib,:)>zero)
387
388	! confirm that we found some points
389	IF ( fopt .EQ. 0) THEN
390	nbexp = nbexp + 1
391	cropvar(ib,ivgt,it) = val_exp ! cropvar_default(ivgt)
392	ELSE
393	DO ilf = 1,fopt
394	! !Not to get lat and lon in wrong order
395	temp_var(ilf,ivgt) = cropvar_mat(sub_index(ib,ilf,1),sub_index(ib,ilf,2),ivgt,icyc,it)
396	ENDDO
397
398
399	! sum_float = zero
400	ok_var = .FALSE.
401	maxcp = nbvmax/2;
402	DO WHILE (.NOT. ok_var)
403
404	ALLOCATE(hashlst1(maxcp), STAT=ALLOC_ERR) ! no.
405	IF (ALLOC_ERR /=0) THEN
406	WRITE(numout,*) "ERROR IN ALLOCATION OF hashlst1:", ALLOC_ERR
407	STOP
408	ENDIF
409	ALLOCATE(hashlst2(maxcp), STAT=ALLOC_ERR) ! area
410	IF (ALLOC_ERR /=0) THEN
411	WRITE(numout,*) "ERROR IN ALLOCATION OF hashlst2:", ALLOC_ERR
412	STOP
413	ENDIF
414	hashlst1(:) = 0
415	hashlst2(:) = zero
416	sgn = zero
417	cps = 1;
418	DO ilf = 1,fopt
419	sgn = sgn + sub_area(ib,ilf)
420	IF (cps==1) THEN
421	hashlst1(cps) = temp_var(ilf,ivgt)
422	hashlst2(cps) = hashlst2(cps) + sub_area(ib,ilf)
423	cps = cps + 1
424	ELSE
425	IF (cps .EQ. maxcp) THEN
426	EXIT
427	ELSE
428	DO jlf = 1,cps
429	IF (temp_var(ilf,ivgt) .EQ. hashlst1(jlf)) THEN
430	hashlst2(jlf) = hashlst2(jlf) +sub_area(ib,ilf)
431	EXIT
432	ELSEIF (jlf .EQ. cps) THEN
433	! .AND. (temp_var(ilf,ivgt) .NE.
434	! hashlst(jlf,1))
435	hashlst1(cps) = temp_var(ilf,ivgt)
436	hashlst2(cps) = hashlst2(cps) + sub_area(ib,ilf)
437	! better to multiply PFT here
438	cps = cps + 1;
439	ENDIF
440	ENDDO
441	ENDIF
442	ENDIF
443	ENDDO
444
445	WRITE(numout,*) "total variety number: ",cps
446	IF (cps .LT. maxcp) THEN
447	ok_var = .TRUE.
448	! find the variety with maximum area
449
450	ELSE
451	! ok_var = .FALSE.
452	maxcp = maxcp * 2
453	DEALLOCATE(hashlst1)
454	DEALLOCATE(hashlst2)
455	ENDIF
456	ENDDO
457
458	! Normalize the surface
459	IF ( sgn .LT. min_sechiba) THEN
460	nbexp = nbexp + 1
461	cropvar(ib,ivgt,it) = val_exp !plantdate_default(ivgt)
462	ELSE
463	areamax = zero
464	maxps = 0
465	IF (cps .LE. 1) THEN
466	WRITE(numout,*) "no sub point found, probably an error occured"
467	STOP
468	ENDIF
469	DO jlf = 1,cps-1
470	IF (hashlst2(jlf) .GT. areamax) THEN
471	areamax = hashlst2(jlf)
472	maxps = jlf
473	ENDIF
474	ENDDO
475	cropvar(ib,ivgt,it) = hashlst1(maxps)
476	DEALLOCATE(hashlst1)
477	DEALLOCATE(hashlst2)
478	! ! INMAX
479	! cropvar(ib,ivgt,icyc) = ANINT(sum_float/sgn)
480	ENDIF
481
482	ENDIF
483
484	ENDDO
485
486	IF ( nbexp .GT. 0) THEN
487	WRITE(numout,*) 'slowproc_PlntDt : default plant date was applied in ', nbexp, 'grid(s)'
488	WRITE(numout,*) 'slowproc_PlntDt : These are either coastal points or having missing data'
489	ENDIF
490	DEALLOCATE (sub_area)
491	DEALLOCATE (sub_index)
492	DEALLOCATE (temp_var)
493	! WRITE(numout,*) 'Planting Date of Site 1 veget ',ivgt,' :
494	! ',plantdate(1,ivgt,icyc)
495	ENDIF
496	ENDDO
497	! End of Veget cycle
498	ENDDO
499	! End of Rotation cycle
500	ENDDO
501	! END of Time Axis cycle
502
503	DEALLOCATE (lat_rel)
504	DEALLOCATE (lon_rel)
505	DEALLOCATE (mask)
506	! DEALLOCATE (sub_area)
507	! DEALLOCATE (sub_index)
508	! DEALLOCATE (temp_date)
509	DEALLOCATE (cropvar_mat)
510
511	WRITE (numout,*) 'Output Crop Variety:'
512	WRITE (numout,*) 'timestep 1:'
513	WRITE (numout,*) cropvar(1,:,1)
514	IF (plantcyc>1) THEN
515	WRITE (numout,*) 'timestep 2:'
516	WRITE (numout,*) cropvar(1,:,2)
517	ENDIF
518	WRITE (numout,) 'END of DEBUG INFO slowproc_CropVar*'
519	! WRITE (numout,*) 'Manual STOP after slowproc_CropVar'
520	! STOP
521	RETURN
522
523	END SUBROUTINE slowproc_CropVar
524	! End of Edition by Xuhui, Nov. 27th 010
525
526
527	!! ================================================================================================================================
528	!! SUBROUTINE : stomate_stics_ManageInput
529	!!
530	!>\BRIEF Interpolate (extract) Planting Date information for a specific crop typ
531	!!
532	!! DESCRIPTION : None
533	!!
534	!! RECENT CHANGE(S) : Xuhui Wang, Oct. 18th, 2010
535	!!
536	!! MAIN OUTPUT VARIABLE(S): None
537	!!
538	!! REFERENCES : None
539	!!
540	!! FLOWCHART : None
541	!! \n
542	!_ ================================================================================================================================
543	SUBROUTINE stomate_stics_ManageInput(nbpt, lalo, neighbours, resolution, contfrac,strIn, varname, manage, yrlen)
544
545	! INTEGER, parameter :: i_std = 4
546	! REAL, parameter :: r_std = 8
547	!
548	! 0.1 INPUT
549	!
550	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated (extracted)
551	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longtitude
552	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vectors of neighbours for each grid point
553	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
554	REAL(r_std),INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid box in lat and lon
555	REAL(r_std),INTENT(in) :: contfrac(nbpt) ! The fraction of land in each grid box
556	CHARACTER(LEN=30),INTENT(in) :: strIn ! getin parameter and Call Sign of the management data
557	CHARACTER(LEN=30),INTENT(in) :: varname ! variable name in the nc file
558	!
559	! 0.2 OUTPUT
560	!
561	REAL(r_std),ALLOCATABLE, INTENT(out) :: manage(:,:,:) ! The planting date of the crop: nbpt, veg, year
562	! nvm is the number of PFTs, there may not be planting date for all the PFTs
563	INTEGER(i_std), INTENT(out) :: yrlen ! year length of the output matrix
564	!
565	! 0.3 LOCAL
566	!
567	INTEGER(i_std) :: nbvmax ! a parameter for interpolation
568	REAL(r_std) :: myres(nbpt,2)
569	CHARACTER(LEN=80) :: filename
570	INTEGER(i_std) :: iml, jml, lml, tml, fid, fid1
571	INTEGER(i_std) :: ip, jp, ib, ilf, fopt, it ! for-loop variable
572	INTEGER(i_std) :: nbexp
573	REAL(r_std) :: lev(1), date, dt
574	REAL(r_std) :: missing_val
575	INTEGER(i_std) :: itau(1)
576
577	INTEGER(i_std) :: nb_dim
578	INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w
579	LOGICAL :: l_ex
580
581	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel
582	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: manage_mat ! LON LAT VEGET, Time
583	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
584	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: temp_data
585	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
586	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
587
588	REAL(r_std) :: sgn, sum_float
589	INTEGER(i_std) :: ivgt ! , icyc, pltcyc
590	CHARACTER(LEN=30) :: callsign
591	LOGICAL :: ok_interpol
592	INTEGER :: ALLOC_ERR
593	LOGICAL :: mydebug = .false.
594
595	! ! croptype = TRIM(croptype) !if croptype is a string
596	! ! else a switch expression is needed
597	! filename = "/work/cont003/p529tan/WXH/plt_date_modif.nc" ! default input
598	! file
599	! ! String operation needed
600	filename = "PlantingDate.nc"
601	CALL getin_p(strIn,filename)
602
603	IF (is_root_prc) THEN
604	! ? what does is_root_prc mean?
605	CALL flininfo(filename, iml, jml, lml, tml, fid)
606	ENDIF
607	CALL bcast(iml)
608	CALL bcast(jml)
609	! CALL bcast(lml)
610	! CALL bcast(tml)
611
612	! Printing information for debugging
613	IF (mydebug) THEN
614	WRITE(numout, *) "Xuhui's debug info for stomate_stics_ManageInput #1:"
615	WRITE(numout, *) "string in: ", strIn
616	WRITE(numout, *) "variable name: ", varname
617	WRITE(numout, *) "filename is: ", filename
618	WRITE(numout, *) "Dimension 1, lon, iml:", iml
619	WRITE(numout, *) "Dimension 2, lat, jml:", jml
620	WRITE(numout, *) "Dimension 3, veget, lml:", lml
621	WRITE(numout, *) "Dimension 4, time, tml:", tml
622	ENDIF
623	! apparently, flinget function is not designed to take veget but levels to
624	! be the
625	! 3rd dimension, modification to lml is needed
626
627	!JG all flio calls must be done by is_root_prc
628	IF (is_root_prc) THEN
629	CALL flioopfd(filename,fid1)
630	CALL flioinqv(fid1,v_n=varname, l_ex = l_ex, nb_dims = nb_dim, len_dims =l_d_w)
631	IF (lml == 0) THEN
632	! CALL
633	! flioinqv(fid1,v_n="PLNTDT",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
634	lml=l_d_w(3)
635	IF (mydebug) THEN
636	WRITE(numout, *) "len_dims: ", l_d_w
637	WRITE(numout, *) "lml AFTER revision"
638	WRITE(numout, *) "lml: ", lml
639	ENDIF
640	ENDIF
641	IF (mydebug) THEN
642	WRITE(numout,*) "nb_dim: ", nb_dim
643	WRITE(numout,*) "resolution: ", resolution(1,:)
644	ENDIF
645
646	IF (nb_dim .NE. 4) THEN
647	WRITE(numout,*) "dimension not supported for ", nb_dim
648	ENDIF
649	tml = l_d_w(4)
650	!yrlen = tml
651	END IF
652	IF (mydebug) THEN
653	WRITE(numout, *) "Now the tml is, :", tml
654	WRITE(numout, *) "Now the lml is:", lml
655	ENDIF
656
657	!JG REMVOVE CALL flioclo(fid1)
658	CALL bcast(lml)
659	CALL bcast(tml)
660	CALL bcast(nb_dim)
661	! CALL bcast(plantcyc)
662
663	! JG yrlen must not be done after bcast(tml)
664	yrlen = tml
665
666	ALLOC_ERR=-1
667	ALLOCATE(manage(nbpt,lml,tml),STAT=ALLOC_ERR)
668	IF (ALLOC_ERR/=0) THEN
669	WRITE(numout,*) "ERROR IN ALLOCATION OF manage: ", ALLOC_ERR
670	ENDIF
671	WRITE(numout,*) "manage ALLOCATED"
672	!CALL bcast(manage)
673
674	!
675	ALLOC_ERR=-1
676	ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
677	IF (ALLOC_ERR/=0) THEN
678	WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
679	STOP
680	ENDIF
681	! CALL bcast(lat_rel)
682
683	ALLOC_ERR=-1
684	ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
685	IF (ALLOC_ERR/=0) THEN
686	WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
687	STOP
688	ENDIF
689	! CALL bcast(lon_rel)
690
691	ALLOC_ERR=-1
692	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
693	IF (ALLOC_ERR/=0) THEN
694	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
695	STOP
696	ENDIF
697	! CALL bcast(mask)
698
699
700	ALLOC_ERR=-1
701	ALLOCATE(manage_mat(iml,jml,lml,tml), STAT=ALLOC_ERR)
702	! !lml is supposed to be nvm (number of PFTs), if not ,change it
703	IF (ALLOC_ERR/=0) THEN
704	WRITE(numout,*) "ERROR IN ALLOCATION of manage_mat : ",ALLOC_ERR
705	STOP
706	ENDIF
707	! CALL bcast(manage_mat)
708	! WRITE (numout,*) 'bcast manage_mat'
709
710	! input of some attributes
711	IF (is_root_prc) THEN
712	! JG with the flioclo, done before this was not ok. Now ok
713	CALL flinget(fid, 'LON', iml, jml, lml, tml, 1, 1, lon_rel)
714	CALL flinget(fid, 'LAT', iml, jml, lml, tml, 1, 1, lat_rel)
715	ENDIF
716	CALL bcast(lon_rel)
717	CALL bcast(lat_rel)
718	WRITE (numout,*) 'lon_rel size: ', SIZE(lon_rel)
719	WRITE (numout,*) 'lat_rel size: ', SIZE(lat_rel)
720
721
722	! input of the matrix
723	IF (is_root_prc) THEN
724	! CALL flinget(fid, 'PLNTDT', iml, jml, lml, tml, 1, 1, plntdt_mat)
725	! JG remove CALL flioopfd: already done
726	! CALL flioopfd(filename,fid1)
727	CALL fliogetv(fid1,trim(varname),manage_mat,start=(/1,1,1,1/),count=(/iml,jml,lml,tml/))
728	! get missing_val
729	CALL fliogeta(fid1,varname,'missing_value',missing_val)
730	CALL flioclo(fid1)
731	ENDIF
732	CALL bcast(manage_mat)
733	CALL bcast(missing_val)
734	WRITE (numout,*) 'bcast manage_mat'
735
736
737	! WRITE(numout,*) 'manage_mat size: ',SIZE(manage_mat)
738	! WRITE(numout,*) 'missing value: ', missing_val
739	! WRITE(numout,*) 'lat(361,284): ',lat_rel(361,284)
740	! WRITE(numout,*) 'lon(361,284): ',lon_rel(361,284)
741	! WRITE(numout,*) 'plntdt(361,284,1,1): ',plntdt_mat(361,284,1,1)
742
743	IF (is_root_prc) CALL flinclo(fid)
744
745	manage(:,:,:) = zero ! nbpt veget year
746
747	DO it = 1,tml
748	DO ivgt = 1,lml ! ? We can suppose PFTs less than 10 are natural veg without planting date, but not now
749	IF (.NOT. natural(ivgt)) THEN
750	WRITE(numout,*) "veget, time: ", ivgt,it
751	nbexp = 0
752	! the number of exceptions
753
754	! mask of available value
755	mask(:,:) = zero; ! Defined in constante.f90
756	DO ip = 1,iml
757	DO jp = 1,jml
758	IF ((manage_mat(ip,jp,ivgt,it) .GT. min_sechiba) .AND. &
759	(manage_mat(ip,jp,ivgt,it) /= missing_val)) THEN
760	mask(ip,jp) = un; ! Defined in constante.f90
761	! here we assumed that for each plant cycle at each
762	! there might be missing data at different grid
763	! in this case, mask has to be generated each plant
764	! cycle each time step
765	ENDIF
766	ENDDO
767	ENDDO
768
769	! Interpolation started
770	nbvmax = 200
771	! the maximum amount of fine grids that one coarse grid may have
772
773	callsign = strIn
774
775	ok_interpol = .FALSE.
776
777	DO WHILE ( .NOT. ok_interpol )
778	WRITE(numout,*) "Pojection arrays for ", callsign, ":"
779	WRITE(numout,*) "nbvmax = ", nbvmax
780
781	ALLOC_ERR = -1
782	ALLOCATE(temp_data(nbvmax,lml), STAT=ALLOC_ERR)
783	IF (ALLOC_ERR /=0) THEN
784	WRITE(numout,*) "ERROR IN ALLOCATION OF temp_data :", ALLOC_ERR
785	STOP
786	ENDIF
787	temp_data = zero
788	ALLOC_ERR = -1
789	ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
790	IF (ALLOC_ERR /=0) THEN
791	WRITE(numout,*) "ERROR IN ALLOCATION OF sub_index :", ALLOC_ERR
792	STOP
793	ENDIF
794	sub_index(:,:,:) = zero
795	ALLOC_ERR = -1
796	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
797	IF (ALLOC_ERR /=0) THEN
798	WRITE(numout,*) "ERROR IN ALLOCATION OF sub_area :",ALLOC_ERR
799	STOP
800	ENDIF
801	sub_area(:,:) = zero
802	myres(:,:) = resolution(:,:)/1000 !m -> km
803	write(numout,*) "resolution updated: ", myres(1,:), " km"
804	!CALL bcast(myres)
805	! CALL bcast(myres)
806
807	! write(,) "calling aggregate_p? "
808	CALL aggregate_p(nbpt, lalo, neighbours, myres, contfrac, &
809	& iml, jml, lon_rel, lat_rel, mask, callsign, &
810	& nbvmax, sub_index, sub_area, ok_interpol)
811	! write(numout,*) "wu: we finished aggregate_p:) "
812
813	IF ( .NOT. ok_interpol ) THEN
814	DEALLOCATE(temp_data)
815	DEALLOCATE(sub_index)
816	DEALLOCATE(sub_area)
817	nbvmax = nbvmax * 2
818	ENDIF
819	ENDDO
820
821	! WRITE(numout,*) "called aggregate_p"
822	! assign the values to plantdate
823	! values should be given to all PFTs
824	DO ib = 1, nbpt
825	! examing all sub_point we found
826	fopt = COUNT(sub_area(ib,:)>zero)
827
828	! confirm that we found some points
829	IF ( fopt .EQ. 0) THEN
830	nbexp = nbexp + 1
831	manage(ib,ivgt,it) = val_exp
832	ELSE
833	DO ilf = 1,fopt
834	! !Not to get lat and lon in wrong order
835	temp_data(ilf,ivgt) = manage_mat(sub_index(ib,ilf,1),sub_index(ib,ilf,2),ivgt,it)
836	ENDDO
837
838	sgn = zero
839	sum_float = zero
840	DO ilf = 1,fopt
841	! average the data weighted by area ! better to multiply
842	! PFT HERE
843	! need to add management specific judgem
844	sum_float = sum_float + temp_data(ilf,ivgt)*sub_area(ib,ilf)
845	sgn = sgn + sub_area(ib,ilf)
846	ENDDO
847
848	! Normalize the surface
849	! sgn can be a scaler, however, to prepare it for future
850	! incorporation of fraction
851	! I make it a vector with nvm values which are equal to each
852	! other
853	IF ( sgn .LT. min_sechiba) THEN
854	nbexp = nbexp + 1
855	manage(ib,ivgt,it) = val_exp ! plantdate_default(ivgt)
856	ELSE
857	manage(ib,ivgt,it) = ANINT(sum_float/sgn)
858	ENDIF
859
860	ENDIF
861
862	ENDDO ! ib
863
864	IF ( nbexp .GT. 0) THEN
865	WRITE(numout,*) 'stomate_stics_ManageInput : exp_val was applied in', nbexp, 'grid(s)'
866	WRITE(numout,*) 'stomate_stics_ManageInput : These are either coastal points or having missing data'
867	ENDIF
868	DEALLOCATE (sub_area)
869	DEALLOCATE (sub_index)
870	DEALLOCATE (temp_data)
871	! WRITE(numout,*) 'Planting Date of Site 1 veget ',ivgt,' :
872	! ',plantdate(1,ivgt,icyc)
873	ENDIF
874	ENDDO
875	! End of Veget cycle
876	ENDDO
877	! End of Time Axis cycle
878
879	DEALLOCATE (lat_rel)
880	DEALLOCATE (lon_rel)
881	DEALLOCATE (mask)
882	DEALLOCATE (manage_mat)
883
884	WRITE (numout,*) 'Output Management Date:'
885	WRITE (numout,*) 'time_step 1:'
886	WRITE (numout,*) manage(1,:,1)
887	IF (tml>1) THEN
888	WRITE (numout,*) 'time_step 2:'
889	WRITE (numout,*) manage(1,:,2)
890	ENDIF
891	WRITE (numout,) 'END of DEBUG INFO stomate_stics_ManageInput*'
892	RETURN
893
894	END SUBROUTINE stomate_stics_ManageInput
895	! End of Edition by Xuhui, Mar. 16th 2011
896
897
898
899	! Date: 22.06.2014, Xuhui Wang
900	! Interpolate (extract) Nitrogen fertilization information
901	! for a specific crop type
902	!! ================================================================================================================================
903	!! SUBROUTINE : stomate_stics_NfertInput
904	!!
905	!>\BRIEF
906	!!
907	!! DESCRIPTION : None
908	!!
909	!! RECENT CHANGE(S) : None
910	!!
911	!! MAIN OUTPUT VARIABLE(S): None
912	!!
913	!! REFERENCES : None
914	!!
915	!! FLOWCHART : None
916	!! \n
917	!_ ================================================================================================================================
918	SUBROUTINE stomate_stics_NfertInput(nbpt, lalo, neighbours, resolution, contfrac,strIn, varname, Nfert, yrlen)
919
920	! INTEGER, parameter :: i_std = 4
921	! REAL, parameter :: r_std = 8
922	!
923	! 0.1 INPUT
924	!
925	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated (extracted)
926	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longtitude
927	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vectors of neighbours for each grid point
928	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
929	REAL(r_std),INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid box in lat and lon
930	!REAL(r_std) :: resolution(nbpt,2) ! The size in km of each grid box in lat and lon
931	REAL(r_std),INTENT(in) :: contfrac(nbpt) ! The fraction of land in each grid box
932	CHARACTER(LEN=30),INTENT(in) :: strIn ! getin parameter and Call Sign of the Nfertment data
933	CHARACTER(LEN=30),INTENT(in) :: varname ! variable name in the nc file
934	!
935	! 0.2 OUTPUT
936	!
937	REAL(r_std),ALLOCATABLE, INTENT(out) :: Nfert(:,:,:) ! The planting date of the crop: nbpt, veg, year
938	! nvm is the number of PFTs, there may not be planting date for all the PFTs
939	INTEGER(i_std), INTENT(out) :: yrlen ! year length of the output matrix
940	!
941	! 0.3 LOCAL
942	!
943	INTEGER(i_std) :: nbvmax ! a parameter for interpolation
944	REAL(r_std) :: myres(nbpt,2)
945	CHARACTER(LEN=80) :: filename
946	INTEGER(i_std) :: iml, jml, lml, tml, fid, fid1
947	INTEGER(i_std) :: ip, jp, ib, ilf, fopt, it ! for-loop variable
948	INTEGER(i_std) :: nbexp
949	REAL(r_std) :: lev(1), date, dt
950	REAL(r_std) :: missing_val
951	INTEGER(i_std) :: itau(1)
952
953	INTEGER(i_std) :: nb_dim
954	INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w
955	LOGICAL :: l_ex
956
957	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel
958	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: Nfert_mat ! LON LAT VEGET, Time
959	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
960	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: temp_data
961	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
962	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
963
964	REAL(r_std) :: sgn, sum_float
965	INTEGER(i_std) :: ivgt ! , icyc, pltcyc
966	CHARACTER(LEN=30) :: callsign
967	LOGICAL :: ok_interpol
968	INTEGER :: ALLOC_ERR
969
970	! ! croptype = TRIM(croptype) !if croptype is a string
971	! ! else a switch expression is needed
972	! filename = "/work/cont003/p529tan/WXH/plt_date_modif.nc" ! default input
973	! file
974	! ! String operation needed
975	filename = "Nfert.nc"
976	CALL getin_p(strIn,filename)
977
978	IF (is_root_prc) THEN
979	! ? what does is_root_prc mean?
980	CALL flininfo(filename, iml, jml, lml, tml, fid)
981	ENDIF
982	CALL bcast(iml)
983	CALL bcast(jml)
984	! CALL bcast(lml)
985	! CALL bcast(tml)
986
987	! Printing information for debugging
988	WRITE(numout, *) "Xuhui's debug info for stomate_stics_NfertInput #1:"
989	WRITE(numout, *) "string in: ", strIn
990	WRITE(numout, *) "variable name: ", varname
991	WRITE(numout, *) "filename is: ", filename
992	WRITE(numout, *) "Dimension 1, lon, iml:", iml
993	WRITE(numout, *) "Dimension 2, lat, jml:", jml
994	WRITE(numout, *) "Dimension 3, veget, lml:", lml
995	WRITE(numout, *) "Dimension 4, time, tml:", tml
996	! apparently, flinget function is not designed to take veget but levels to
997	! be the
998	! 3rd dimension, modification to lml is needed
999
1000	!JG all flio calls must be done by is_root_prc
1001	IF (is_root_prc) THEN
1002	CALL flioopfd(filename,fid1)
1003	CALL flioinqv(fid1,v_n=varname, l_ex = l_ex, nb_dims = nb_dim, len_dims =l_d_w)
1004	IF (lml == 0) THEN
1005	! CALL
1006	! flioinqv(fid1,v_n="PLNTDT",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
1007	lml=l_d_w(3)
1008	WRITE(numout, *) "len_dims: ", l_d_w
1009	WRITE(numout, *) "lml AFTER revision"
1010	WRITE(numout, *) "lml: ", lml
1011	ENDIF
1012	WRITE(numout,*) "nb_dim: ", nb_dim
1013	WRITE(numout,*) "resolution: ", resolution(1,:)
1014
1015	IF (nb_dim .NE. 4) THEN
1016	WRITE(numout,*) "dimension not supported for ", nb_dim
1017	ENDIF
1018	tml = l_d_w(4)
1019	!yrlen = tml
1020	END IF
1021
1022	WRITE(numout, *) "Now the tml is, :", tml
1023	WRITE(numout, *) "Now the lml is:", lml
1024
1025	!JG REMVOVE CALL flioclo(fid1)
1026	CALL bcast(lml)
1027	CALL bcast(tml)
1028	CALL bcast(nb_dim)
1029	! CALL bcast(plantcyc)
1030
1031	! JG yrlen must not be done after bcast(tml)
1032	yrlen = tml
1033
1034	ALLOC_ERR=-1
1035	ALLOCATE(Nfert(nbpt,lml,tml),STAT=ALLOC_ERR)
1036	IF (ALLOC_ERR/=0) THEN
1037	WRITE(numout,*) "ERROR IN ALLOCATION OF Nfert: ", ALLOC_ERR
1038	ENDIF
1039	WRITE(numout,*) "Nfert ALLOCATED"
1040	!CALL bcast(Nfert)
1041
1042	!
1043	ALLOC_ERR=-1
1044	ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
1045	IF (ALLOC_ERR/=0) THEN
1046	WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
1047	STOP
1048	ENDIF
1049	! CALL bcast(lat_rel)
1050
1051	ALLOC_ERR=-1
1052	ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
1053	IF (ALLOC_ERR/=0) THEN
1054	WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
1055	STOP
1056	ENDIF
1057	! CALL bcast(lon_rel)
1058
1059	ALLOC_ERR=-1
1060	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
1061	IF (ALLOC_ERR/=0) THEN
1062	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
1063	STOP
1064	ENDIF
1065	! CALL bcast(mask)
1066
1067
1068	ALLOC_ERR=-1
1069	ALLOCATE(Nfert_mat(iml,jml,lml,tml), STAT=ALLOC_ERR)
1070	! !lml is supposed to be nvm (number of PFTs), if not ,change it
1071	IF (ALLOC_ERR/=0) THEN
1072	WRITE(numout,*) "ERROR IN ALLOCATION of Nfert_mat : ",ALLOC_ERR
1073	STOP
1074	ENDIF
1075	! CALL bcast(Nfert_mat)
1076	! WRITE (numout,*) 'bcast Nfert_mat'
1077
1078	! input of some attributes
1079	IF (is_root_prc) THEN
1080	! JG with the flioclo, done before this was not ok. Now ok
1081	CALL flinget(fid, 'LON', iml, jml, lml, tml, 1, 1, lon_rel)
1082	CALL flinget(fid, 'LAT', iml, jml, lml, tml, 1, 1, lat_rel)
1083	ENDIF
1084	CALL bcast(lon_rel)
1085	CALL bcast(lat_rel)
1086	WRITE (numout,*) 'lon_rel size: ', SIZE(lon_rel)
1087	WRITE (numout,*) 'lat_rel size: ', SIZE(lat_rel)
1088
1089
1090	! input of the matrix
1091	IF (is_root_prc) THEN
1092	! CALL flinget(fid, 'PLNTDT', iml, jml, lml, tml, 1, 1, plntdt_mat)
1093	! JG remove CALL flioopfd: already done
1094	! CALL flioopfd(filename,fid1)
1095	CALL fliogetv(fid1,trim(varname),Nfert_mat,start=(/1,1,1,1/),count=(/iml,jml,lml,tml/))
1096	! get missing_val
1097	CALL fliogeta(fid1,varname,'missing_value',missing_val)
1098	CALL flioclo(fid1)
1099	ENDIF
1100	CALL bcast(Nfert_mat)
1101
1102	IF (is_root_prc) CALL flinclo(fid)
1103
1104	Nfert(:,:,:) = zero ! nbpt veget year
1105
1106	DO it = 1,tml
1107	DO ivgt = 1,lml ! ? We can suppose PFTs less than 10 are natural veg without planting date, but not now
1108	IF (.NOT. natural(ivgt)) THEN
1109	WRITE(numout,*) "veget, time: ", ivgt,it
1110	nbexp = 0
1111	! the number of exceptions
1112
1113	! mask of available value
1114	mask(:,:) = zero; ! Defined in constante.f90
1115	DO ip = 1,iml
1116	DO jp = 1,jml
1117	IF ((Nfert_mat(ip,jp,ivgt,it) .GE. zero) .AND. &
1118	(Nfert_mat(ip,jp,ivgt,it) /= missing_val)) THEN
1119	mask(ip,jp) = un; ! Defined in constante.f90
1120	! here we assumed that for each plant cycle at each
1121	! there might be missing data at different grid
1122	! in this case, mask has to be generated each plant
1123	! cycle each time step
1124	ENDIF
1125	ENDDO
1126	ENDDO
1127
1128	! Interpolation started
1129	nbvmax = 200
1130	! the maximum amount of fine grids that one coarse grid may have
1131
1132	callsign = strIn
1133
1134	ok_interpol = .FALSE.
1135
1136	DO WHILE ( .NOT. ok_interpol )
1137	WRITE(numout,*) "Pojection arrays for ", callsign, ":"
1138	WRITE(numout,*) "nbvmax = ", nbvmax
1139
1140	ALLOC_ERR = -1
1141	ALLOCATE(temp_data(nbvmax,lml), STAT=ALLOC_ERR)
1142	IF (ALLOC_ERR /=0) THEN
1143	WRITE(numout,*) "ERROR IN ALLOCATION OF temp_data :", ALLOC_ERR
1144	STOP
1145	ENDIF
1146	ALLOC_ERR = -1
1147	ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
1148	IF (ALLOC_ERR /=0) THEN
1149	WRITE(numout,*) "ERROR IN ALLOCATION OF sub_index :", ALLOC_ERR
1150	STOP
1151	ENDIF
1152	sub_index(:,:,:) = zero
1153	ALLOC_ERR = -1
1154	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
1155	IF (ALLOC_ERR /=0) THEN
1156	WRITE(numout,*) "ERROR IN ALLOCATION OF sub_area :",ALLOC_ERR
1157	STOP
1158	ENDIF
1159	sub_area(:,:) = zero
1160	myres(:,:) = resolution(:,:)/1000 !m -> km
1161	write(numout,*) "resolution updated: ", myres(1,:), " km"
1162	!CALL bcast(myres)
1163
1164	! write(numout,*) "wu: stepping in aggregate_p? "
1165	CALL aggregate_p(nbpt, lalo, neighbours, myres, contfrac, &
1166	& iml, jml, lon_rel, lat_rel, mask, callsign, &
1167	& nbvmax, sub_index, sub_area, ok_interpol)
1168	! write(numout,*) "wu: we finished aggregate_p:) "
1169
1170	IF ( .NOT. ok_interpol ) THEN
1171	DEALLOCATE(temp_data)
1172	DEALLOCATE(sub_index)
1173	DEALLOCATE(sub_area)
1174	nbvmax = nbvmax * 2
1175	ENDIF
1176	ENDDO
1177
1178	! WRITE(numout,*) "called aggregate_p"
1179	! assign the values to plantdate
1180	! values should be given to all PFTs
1181	DO ib = 1, nbpt
1182	! examing all sub_point we found
1183	fopt = COUNT(sub_area(ib,:)>zero)
1184
1185	! confirm that we found some points
1186	IF ( fopt .EQ. 0) THEN
1187	nbexp = nbexp + 1
1188	! Nfert(ib,ivgt,it) = val_exp
1189	Nfert(ib,ivgt,it) = 0
1190	ELSE
1191	DO ilf = 1,fopt
1192	! !Not to get lat and lon in wrong order
1193	temp_data(ilf,ivgt) = Nfert_mat(sub_index(ib,ilf,1),sub_index(ib,ilf,2),ivgt,it)
1194	ENDDO
1195
1196	sgn = zero
1197	sum_float = zero
1198	DO ilf = 1,fopt
1199	! average the data weighted by area ! better to multiply
1200	! PFT HERE
1201	! need to add Nfertment specific judgem
1202	sum_float = sum_float + temp_data(ilf,ivgt)*sub_area(ib,ilf)
1203	sgn = sgn + sub_area(ib,ilf)
1204	ENDDO
1205
1206	! Normalize the surface
1207	! sgn can be a scaler, however, to prepare it for future
1208	! incorporation of fraction
1209	! I make it a vector with nvm values which are equal to each
1210	! other
1211	IF ( sgn .LT. min_sechiba) THEN
1212	nbexp = nbexp + 1
1213	! Nfert(ib,ivgt,it) = val_exp ! plantdate_default(ivgt)
1214	Nfert(ib,ivgt,it) = 0 ! plantdate_default(ivgt)
1215	ELSE
1216	!Nfert(ib,ivgt,it) = ANINT(sum_float/sgn)
1217	Nfert(ib,ivgt,it) = sum_float/sgn
1218	ENDIF
1219
1220	ENDIF
1221
1222	ENDDO ! ib
1223
1224	IF ( nbexp .GT. 0) THEN
1225	WRITE(numout,*) 'stomate_stics_NfertInput : 0 was applied in', nbexp, 'grid(s)'
1226	WRITE(numout,*) 'stomate_stics_NfertInput : These are either coastal points or having missing data'
1227	ENDIF
1228	DEALLOCATE (sub_area)
1229	DEALLOCATE (sub_index)
1230	DEALLOCATE (temp_data)
1231	! WRITE(numout,*) 'Planting Date of Site 1 veget ',ivgt,' :
1232	! ',plantdate(1,ivgt,icyc)
1233	ENDIF
1234	ENDDO
1235	! End of Veget cycle
1236	ENDDO
1237	! End of Time Axis cycle
1238
1239	DEALLOCATE (lat_rel)
1240	DEALLOCATE (lon_rel)
1241	DEALLOCATE (mask)
1242	DEALLOCATE (Nfert_mat)
1243
1244	END SUBROUTINE stomate_stics_NfertInput
1245	! End of Edition by Xuhui, Mar. 16th 2011
1246
1247	!! ================================================================================================================================
1248	!! SUBROUTINE : stomate_stics_calc_N_limfert
1249	!!
1250	!>\BRIEF
1251	!!
1252	!! DESCRIPTION : None
1253	!!
1254	!! RECENT CHANGE(S) : None
1255	!!
1256	!! MAIN OUTPUT VARIABLE(S): None
1257	!!
1258	!! REFERENCES : None
1259	!!
1260	!! FLOWCHART : None
1261	!! \n
1262	!_ ================================================================================================================================
1263	SUBROUTINE stomat_stics_calc_N_limfert(npts,N_fert_total,N_limfert)
1264
1265	INTEGER (i_std) , INTENT(in) :: npts
1266	!REAL(r_std), DIMENSION(npts,nvm,nstocking), INTENT(in) :: nfertamm
1267	!REAL(r_std), DIMENSION(npts,nvm,nstocking), INTENT(in) :: nfertnit
1268	!REAL(r_std), DIMENSION(npts,nvm,nstocking), INTENT(in) :: nliquidmanure
1269	!REAL(r_std), DIMENSION(npts,nvm,nstocking), INTENT(in) :: nslurry
1270	!REAL(r_std), DIMENSION(npts,nvm,nstocking), INTENT(in) :: nsolidmanure
1271	!REAL(r_std), DIMENSION(npts,nvm) , INTENT(in) :: legume_fraction
1272	!REAL(r_std), DIMENSION(npts,nvm) , INTENT(in) :: soil_fertility
1273	REAL(r_std), DIMENSION(npts,nvm) , INTENT(in) :: N_fert_total
1274	REAL(r_std), DIMENSION(:,:) , INTENT(out) :: N_limfert
1275
1276
1277	INTEGER(i_std) :: k,j,ip
1278
1279	!N_fert_total(:,:) = 0.0
1280	!DO k=1,nstocking
1281	! N_fert_total(:,:) = (N_fert_total(:,:) + nfertamm(:,:,k) + &
1282	! nfertnit(:,:,k) + nliquidmanure(:,:,k) + &
1283	! nslurry(:,:,k) + nsolidmanure(:,:,k))
1284	!ENDDO
1285	!N_fert_total(:,:) = N_fert_total(:,:) * 10000
1286	!DO j=2,nvm
1287	! IF ((management_intensity(j) .EQ. 2).AND. is_c3(j)) THEN
1288	! N_fert_total(:,mcut_C3)=N_fert_total(:,j)
1289	! N_fert_total(:,mgraze_C3)=N_fert_total(:,j)
1290	! ENDIF
1291	! IF ((management_intensity(j) .EQ. 2).AND. (.NOT.is_c3(j))) THEN
1292	! N_fert_total(:,mcut_C4)=N_fert_total(:,j)
1293	! N_fert_total(:,mgraze_C4)=N_fert_total(:,j)
1294	! ENDIF
1295	!
1296	!ENDDO
1297
1298
1299	! N_limfert(:,:) = 1. + N_effect - N_effect * (0.75 ** (N_fert_total(:,:)/30))
1300	N_limfert(:,:) = 0
1301	DO j=2, nvm
1302	IF ( (SP_neffmax(j) .lt. zero) .AND. ( ok_LAIdev(j) ) ) THEN
1303	CALL ipslerr_p(3, 'calc_limfert', 'negative N effect', '', '')
1304	ENDIF
1305	IF ( ok_LAIdev(j) ) THEN
1306	N_limfert(:,j) = 1. + SP_neffmax(j) - SP_neffmax(j)(SP_nsatrat(j) * (N_fert_total(:,j)/10))
1307	ELSE
1308	N_limfert(:,j) = 1.
1309	ENDIF
1310	ENDDO
1311
1312	! WHERE (N_limfert(:,:) .LT. 1.0) ! never reached if N_effect>0
1313	! N_limfert(:,:) = 1.0
1314	! ELSEWHERE (N_limfert(:,:) .GT. 2.0) ! will yield mistake if N_effect>1
1315	! N_limfert(:,:) = 1.+N_effect
1316	! ENDWHERE
1317
1318	END SUBROUTINE stomat_stics_calc_N_limfert
1319	! end of the nitrogen fertilization module
1320
1321
1322	!! ================================================================================================================================
1323	!! SUBROUTINE : stomate_stics_read_cycle
1324	!!
1325	!>\BRIEF
1326	!!
1327	!! DESCRIPTION : None
1328	!!
1329	!! RECENT CHANGE(S) : None
1330	!!
1331	!! MAIN OUTPUT VARIABLE(S): None
1332	!!
1333	!! REFERENCES : None
1334	!!
1335	!! FLOWCHART : None
1336	!! \n
1337	!_ ================================================================================================================================
1338	SUBROUTINE stomate_stics_read_cycle(nbpt,lalo,neighbours,resolution,contfrac, & ! in
1339	cyc_num, cyc_num_tot ) ! out
1340	! 0.1 INPUT
1341	!
1342	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated (extracted)
1343	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longtitude
1344	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vectors of neighbours for each grid point
1345	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
1346	REAL(r_std),INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid box in lat and lon
1347	!REAL(r_std) :: resolution(nbpt,2) ! The size in km of each !grid box in lat and lon
1348	REAL(r_std),INTENT(in) :: contfrac(nbpt) ! The fraction of land in each grid box
1349
1350	INTEGER(i_std),INTENT(out), DIMENSION(nbpt,nvm) :: cyc_num ! flag for starting the rotation for kjpindex, nvm
1351	INTEGER(i_std),INTENT(out), DIMENSION(nbpt) :: cyc_num_tot ! number of current rotation cycle
1352
1353	! 0.3 local variables
1354	INTEGER(i_std) :: yrlen
1355	CHARACTER(LEN=30) :: strManage
1356	CHARACTER(LEN=30) :: strVar
1357	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: manage
1358	INTEGER(i_std) :: ip
1359
1360	!!!!-------------------------------------------------------------------------------------------
1361
1362	cyc_num(:,:) = 1
1363
1364	CALL getin_p('IMPOSE_ROT',rot_1d)
1365	! by default, rot_1d is true
1366	IF (nbpt==1) THEN
1367	cyc_num_tot(1) = cyc_rot_max
1368	ELSE ! multiple point simulations
1369	IF (rot_1d) THEN
1370	cyc_num_tot(:) = cyc_rot_max
1371	ELSE
1372	strManage = 'NUMROTATE_FILE'
1373	strVar = 'CycleNum'
1374	CALL stomate_stics_ManageInput(nbpt,lalo,neighbours,resolution,contfrac,strManage,strVar,manage,yrlen)
1375	DO ip = 1,nbpt
1376	cyc_num_tot(ip) = INT(MAXVAL(manage(ip,:,1)),i_std)
1377	IF (cyc_num_tot(ip) .LT. 1) THEN
1378	cyc_num_tot(ip) = 1 ! at least one cycle in rotation
1379	WRITE(numout,*) 'WARNING xuhui: ip, cyc_num_tot(ip) ', ip, MAXVAL(manage(ip,:,1))
1380	ENDIF
1381	ENDDO
1382	IF (MAXVAL(cyc_num_tot(:))>10) THEN
1383	WRITE(numout,*) 'xuhui: rotation cycles longer than 10, likely a bug'
1384	ENDIF
1385	IF (MAXVAL(cyc_num_tot(:))>cyc_rot_max) THEN
1386	WRITE(numout,*) 'xuhui: No. of rotation updated'
1387	WRITE(numout,*) 'MAXVAL(cyc_num_tot(:)), cyc_rot_max', MAXVAL(cyc_num_tot(:)), cyc_rot_max
1388	cyc_rot_max = MAXVAL(cyc_num_tot(:))
1389	ELSEIF (MAXVAL(cyc_num_tot(:))<cyc_rot_max) THEN
1390	WRITE(numout,*) 'xuhui: MAXVAL(cyc_num_tot(:))<cyc_rot_max',MAXVAL(cyc_num_tot(:)), cyc_rot_max
1391	WRITE(numout,*) 'bad setting in cyc_rot_max'
1392	WRITE(numout,*) 'xuhui: No. of rotation updated'
1393	cyc_rot_max = MAXVAL(cyc_num_tot(:))
1394	ENDIF
1395	ENDIF
1396	ENDIF
1397	IF (ALLOCATED(manage)) DEALLOCATE(manage) ! clear manage for other input purpose
1398
1399	END SUBROUTINE stomate_stics_read_cycle
1400
1401	!! ================================================================================================================================
1402	!! SUBROUTINE : stomate_stics_read_rotation
1403	!!
1404	!>\BRIEF
1405	!!
1406	!! DESCRIPTION : None
1407	!!
1408	!! RECENT CHANGE(S) : None
1409	!!
1410	!! MAIN OUTPUT VARIABLE(S): None
1411	!!
1412	!! REFERENCES : None
1413	!!
1414	!! FLOWCHART : None
1415	!! \n
1416	!_ ================================================================================================================================
1417	SUBROUTINE stomate_stics_read_rotation(nbpt,lalo,neighbours,resolution,contfrac, & ! in
1418	rot_cmd_store ) ! out
1419	! 0.1 INPUT
1420	!
1421	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated (extracted)
1422	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longtitude
1423	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vectors of neighbours for each grid point
1424	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
1425	REAL(r_std),INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid box in lat and lon
1426	!REAL(r_std) :: resolution(nbpt,2) ! The size in km of each !grid box in lat and lon
1427	REAL(r_std),INTENT(in) :: contfrac(nbpt) ! The fraction of land in each grid box
1428
1429	INTEGER(i_std), INTENT(out), ALLOCATABLE, DIMENSION(:,:,:) :: rot_cmd_store
1430
1431	! 0.3 local variables
1432	INTEGER(i_std) :: yrlen
1433	CHARACTER(LEN=30) :: strManage
1434	CHARACTER(LEN=30) :: strVar
1435	CHARACTER(LEN=30) :: temp_varname !! temporary variable string for rot_1d only
1436	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: manage
1437	INTEGER(i_std) :: ip, kc, ier
1438	INTEGER(i_std),DIMENSION(rot_cmd_max) :: temp_cmd !! temporary command variable for rot_1d only
1439
1440	!!!!-------------------------------------------------------------------------------------------
1441
1442	!!!! rotation command
1443	IF (rot_1d) THEN
1444	DO kc = 1,cyc_rot_max
1445	WRITE(temp_varname,"('CMDROTATE_',i1)"),kc
1446	CALL getin_p(temp_varname,temp_cmd)
1447	!!!! confirm how getin_p perform when the input command is less than rot_cmd_max
1448	WRITE(numout,*) 'xuhui: kc, ', kc, 'temp_cmd, ', temp_cmd
1449	DO ip = 1,nbpt
1450	rot_cmd_store(ip,:,kc) = temp_cmd
1451	ENDDO
1452	ENDDO
1453	ELSE ! reading rotation map
1454	strManage = 'CMDROTATE_FILE'
1455	strVar = 'Command'
1456	CALL stomate_stics_ManageInput(nbpt,lalo,neighbours,resolution,contfrac,strManage,strVar,manage,yrlen)
1457	!!! dimension of manage: kjpindex, cm_max, cyc_rot_max
1458	IF ( SIZE(manage,3,i_std) .NE. cyc_rot_max ) THEN
1459	WRITE(numout,*) 'SIZE(manage,3), cyc_rot_max', SIZE(manage,3,i_std), cyc_rot_max
1460	STOP 'input rotation cycle did not match cyc_rot_max'
1461	ENDIF
1462	IF ( SIZE(manage,2,i_std) .GT. rot_cmd_max ) THEN
1463	IF (SIZE(manage,2,i_std) .GT. 10 ) THEN
1464	WRITE(numout,*) 'WARNING: more than 10 commands for one rotation'
1465	WRITE(numout,*) 'xuhui: memory issue may occur. Are you sure about so many rotation commands?'
1466	ENDIF
1467	WRITE(numout,*) 'xuhui: update rot_cmd_max'
1468	WRITE(numout,*) ' SIZE(manage,2), rot_cmd_max', SIZE(manage,2,i_std), rot_cmd_max
1469	rot_cmd_max = SIZE(manage,2,i_std)
1470	DEALLOCATE(rot_cmd_store)
1471	ALLOCATE(rot_cmd_store(nbpt,rot_cmd_max,cyc_rot_max),stat=ier)
1472	IF (ier/=0) THEN
1473	WRITE(numout,*) "ERROR IN RE-ALLOCATION OF rot_cmd_store: ",ier
1474	STOP 'stomate_init rot_cmd_store reallocate'
1475	ENDIF
1476	ENDIF
1477	IF (SIZE(manage,2,i_std) .LT. rot_cmd_max) THEN
1478	rot_cmd_store(:,1:SIZE(manage,2,i_std),:) = INT(manage,i_std)
1479	rot_cmd_store(:,SIZE(manage,2,i_std)+1:rot_cmd_max,:) = 0
1480	ELSE !!! size equivalent
1481	rot_cmd_store = INT(manage,i_std)
1482	ENDIF
1483	ENDIF
1484	IF (ALLOCATED(manage)) DEALLOCATE(manage) ! clear manage for other input purpose
1485
1486	END SUBROUTINE stomate_stics_read_rotation
1487
1488	!! ================================================================================================================================
1489	!! SUBROUTINE : stomate_stics_read_plantdate
1490	!!
1491	!>\BRIEF
1492	!!
1493	!! DESCRIPTION : None
1494	!!
1495	!! RECENT CHANGE(S) : None
1496	!!
1497	!! MAIN OUTPUT VARIABLE(S): None
1498	!!
1499	!! REFERENCES : None
1500	!!
1501	!! FLOWCHART : None
1502	!! \n
1503	!_ ================================================================================================================================
1504	SUBROUTINE stomate_stics_read_plantdate(nbpt,lalo,neighbours,resolution,contfrac, & ! in
1505	cyc_num, &
1506	plantdate, plantdate_now) ! out
1507	! 0.1 INPUT
1508	!
1509	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated (extracted)
1510	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longtitude
1511	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vectors of neighbours for each grid point
1512	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
1513	REAL(r_std),INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid box in lat and lon
1514	!REAL(r_std) :: resolution(nbpt,2) ! The size in km of each !grid box in lat and lon
1515	REAL(r_std),INTENT(in) :: contfrac(nbpt) ! The fraction of land in each grid box
1516	INTEGER(i_std),INTENT(in), DIMENSION(nbpt,nvm) :: cyc_num ! flag for starting the rotation for kjpindex, nvm
1517
1518	INTEGER(i_std), INTENT(out), DIMENSION (nbpt,nvm,cyc_rot_max) :: plantdate !! kjpindex, nvm, cyc_rot_max
1519	INTEGER(i_std), INTENT(out), DIMENSION (nbpt,nvm) :: plantdate_now !! kjpindex, nvm, plantdate of current cycle
1520
1521	! 0.3 local variables
1522	INTEGER(i_std) :: yrlen
1523	CHARACTER(LEN=30) :: strManage
1524	CHARACTER(LEN=30) :: strVar
1525	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: manage
1526	INTEGER(i_std) :: ip, kc, ier, j
1527
1528	!!!!-------------------------------------------------------------------------------------------
1529
1530	!iplt_1d = .TRUE.
1531	CALL getin_p('IMPOSE_IPLT',iplt_1d)
1532	IF (.not. iplt_1d) THEN
1533	strManage = "IPLT_FILE"
1534	strVar = "PLNTDT"
1535	CALL stomate_stics_ManageInput(nbpt,lalo,neighbours,resolution,contfrac,strManage,strVar,manage,yrlen)
1536	!manage(manage=val_exp) = 1 !when there is exceptional value, we applied 1
1537	IF ( SIZE(manage,3,i_std) .NE. cyc_rot_max ) THEN
1538	WRITE(numout,*) 'cyc_rot_max, SIZE(manage,3)', cyc_rot_max, SIZE(manage,3,i_std)
1539	STOP 'input and specified maximum rotation cycle not matched'
1540	ENDIF
1541
1542	!!!! deal with exceptional values
1543	DO ip = 1,nbpt
1544	DO j = 2,nvm
1545	IF ( ok_LAIdev(j) ) THEN
1546	DO kc = 1,cyc_rot_max
1547	IF ( (manage(ip,j,kc)==val_exp) .OR. (manage(ip,j,kc) .LT. zero) ) THEN
1548	manage(ip,j,kc) = 1 + year_length_in_days*(kc-1)
1549	ENDIF
1550	ENDDO
1551	ENDIF
1552	ENDDO
1553	ENDDO
1554
1555	DO j = 2,nvm
1556	IF (nvm_plnt) THEN
1557	IF (.NOT. (nvm .EQ. SIZE(manage,2,i_std)) ) THEN
1558	WRITE(numout,*) 'nvm_plnt: ', nvm_plnt
1559	WRITE(numout,*) 'nvm, SIZE(manage,2)', nvm, SIZE(manage,2,i_std)
1560	STOP 'PFT more than input planting date'
1561	ENDIF
1562	WHERE(manage(:,j,:) >= val_exp)
1563	manage(:,j,:) = val_exp
1564	ENDWHERE
1565	plantdate(:,j,:) = INT(manage(:,j,:),i_std)
1566	ELSE
1567	IF (MAXVAL(pft_to_mtc(:)) > SIZE(manage,2,i_std)) THEN
1568	WRITE(numout,*) 'nvm_plnt: ', nvm_plnt
1569	WRITE(numout,*) 'MAXVAL(pft_to_mtc(:)), SIZE(manage,2)', MAXVAL(pft_to_mtc(:)), SIZE(manage,2,i_std)
1570	STOP 'PFT more than input planting date'
1571	ENDIF
1572	WHERE(manage(:,j,:) >= val_exp)
1573	manage(:,j,:) = val_exp
1574	ENDWHERE
1575	plantdate(:,j,:) = INT(manage(:,pft_to_mtc(j),:),i_std)
1576	ENDIF
1577	ENDDO
1578	ELSE ! iplt_1d
1579	IF (nvm_plnt) THEN
1580	IF (nvm .GT. SIZE(SP_iplt0,1,i_std)) THEN
1581	WRITE(numout,*) 'nvm_plnt: ', nvm_plnt
1582	WRITE(numout,*) 'nvm, SIZE(SP_iplt0,1)', nvm, SIZE(SP_iplt0,1,i_std)
1583	STOP 'PFT more than input planting date'
1584	ENDIF
1585	DO j = 2,nvm
1586	plantdate(:,j,1) = SP_iplt0(j)
1587	IF ( cyc_rot_max .GT. 1) THEN
1588	plantdate(:,j,2) = SP_iplt1(j)
1589	ENDIF
1590	IF ( cyc_rot_max .GT. 2) THEN
1591	plantdate(:,j,3) = SP_iplt2(j)
1592	ENDIF
1593	ENDDO
1594	ELSE !!!! using pft_to_mtc
1595	WRITE(numout,*) 'SP_iplt0, ', SP_iplt0
1596	IF (MAXVAL(pft_to_mtc(:)) .GT. SIZE(SP_iplt0,1,i_std)) THEN
1597	WRITE(numout,*) 'nvm_plnt: ', nvm_plnt
1598	WRITE(numout,*) 'MAXVAL(pft_to_mtc(:)), SIZE(SP_iplt0,1)', MAXVAL(pft_to_mtc(:)), SIZE(SP_iplt0,1,i_std)
1599	STOP 'PFT more than input planting date'
1600	ENDIF
1601	DO j = 2,nvm
1602	plantdate(:,j,1) = SP_iplt0(pft_to_mtc(j))
1603	IF ( cyc_rot_max .GT. 1) THEN
1604	plantdate(:,j,2) = SP_iplt1(pft_to_mtc(j))
1605	ENDIF
1606	IF ( cyc_rot_max .GT. 2) THEN
1607	plantdate(:,j,3) = SP_iplt2(pft_to_mtc(j))
1608	ENDIF
1609	ENDDO
1610	ENDIF
1611	ENDIF ! iplt_1d
1612
1613	IF (ok_rotate) THEN
1614	DO ip = 1,nbpt
1615	DO j = 2,nvm
1616	IF (ok_LAIdev(j)) THEN
1617	plantdate_now(ip,j) = plantdate(ip,j,cyc_num(ip,j))
1618	ENDIF
1619	ENDDO
1620	ENDDO
1621	ELSE
1622	plantdate_now = plantdate(:,:,1)
1623	ENDIF
1624	END SUBROUTINE stomate_stics_read_plantdate
1625
1626	!! ================================================================================================================================
1627	!! SUBROUTINE : stomate_stics_get_cmd
1628	!!
1629	!>\BRIEF
1630	!!
1631	!! DESCRIPTION : None
1632	!!
1633	!! RECENT CHANGE(S) : None
1634	!!
1635	!! MAIN OUTPUT VARIABLE(S): None
1636	!!
1637	!! REFERENCES : None
1638	!!
1639	!! FLOWCHART : None
1640	!! \n
1641	!_ ================================================================================================================================
1642	SUBROUTINE stomate_stics_get_cmd(cmdin, src_rot, tgt_rot, prc_rot)
1643	! 0.1 Input
1644	INTEGER(i_std),INTENT(in) :: cmdin
1645	! 0.2 Output
1646	INTEGER(i_std),INTENT(out) :: src_rot, tgt_rot
1647	REAL(r_std),INTENT(out) :: prc_rot
1648
1649	! 0.3 Local variables
1650
1651	!!!! --------------------------------------------------------------
1652	IF (cmdin > 1010000 .OR. cmdin < 0 ) THEN
1653	WRITE(numout,*) 'cmdin, ',cmdin
1654	STOP 'cmd error in stomate_stics_get_cmd'
1655	ENDIF
1656	IF (cmdin == 0) THEN
1657	tgt_rot = 0
1658	src_rot = 0
1659
1660	prc_rot = 0.0
1661	ELSE
1662	tgt_rot = MOD(cmdin, 100)
1663	src_rot = MOD(FLOOR(FLOAT(cmdin)/100), 100)
1664
1665	prc_rot = FLOAT(FLOOR(FLOAT(cmdin)/10000))/100.0
1666	IF (printlev >=4) THEN
1667	WRITE(numout,*) 'xuhui: cmdin, tgt_rot, src_rot, prc_rot', cmdin, src_rot, tgt_rot, prc_rot
1668	ENDIF
1669	IF (prc_rot .GT. 1.0 .AND. prc_rot .LT. 1.0+0.01) THEN ! resolve potential precision issues
1670	prc_rot = 1.0
1671	ENDIF
1672	!!! consistency check
1673	IF (prc_rot .GT. 1.0) THEN
1674	WRITE(numout,*) 'percent change larger than 1..., prc_rot',prc_rot
1675	STOP 'incorrect percent rotation, stomate_stics_get_cmd'
1676	ENDIF
1677	IF ( (tgt_rot .GT. nvm) .OR. ( .NOT. (ok_LAIdev(tgt_rot) .OR. (tgt_rot .EQ. 1)) ) ) THEN
1678	WRITE(numout,*) 'rotation target error: tgt_rot ', tgt_rot
1679	WRITE(numout,*) 'nvm, ok_LAIdev', nvm, ok_LAIdev
1680	STOP 'incorrect rotation target, stomate_stics_get_cmd'
1681	ENDIF
1682	IF ( (src_rot .GT. nvm) .OR. ( .NOT. (ok_LAIdev(src_rot) .OR. (src_rot .EQ. 1)) ) ) THEN
1683	WRITE(numout,*) 'rotation target error: src_rot ', src_rot
1684	WRITE(numout,*) 'nvm, ok_LAIdev', nvm, ok_LAIdev
1685	STOP 'incorrect rotation source, stomate_stics_get_cmd'
1686	ENDIF
1687	ENDIF
1688	END SUBROUTINE stomate_stics_get_cmd
1689
1690
1691	!! ================================================================================================================================
1692	!! SUBROUTINE : stomate_stics_rotation
1693	!!
1694	!>\BRIEF
1695	!!
1696	!! DESCRIPTION : None
1697	!!
1698	!! RECENT CHANGE(S) : None
1699	!!
1700	!! MAIN OUTPUT VARIABLE(S): None
1701	!!
1702	!! REFERENCES : None
1703	!!
1704	!! FLOWCHART : None
1705	!! \n
1706	!_ ================================================================================================================================
1707	SUBROUTINE stomate_stics_rotation(kjpindex, ip, matrix_prc, maxfrac, & ! in
1708	in_cycle, cyc_num_tot, plantdate, nrec, nlev, &
1709	plantdate_now, & ! inout
1710	soilc_mict, cyc_num, litter, turnover_daily, &
1711	deepC_a, deepC_s, deepC_p, carbon, &
1712	age, ind, PFTpresent, senescence, &
1713	when_growthinit, everywhere, leaf_frac )
1714	!!! This subroutine transfer litter and soil carbon when the cropland is rotated
1715	!!! update plantdate_now
1716	!!! veget_max and soil water and thermo properties will be transferred in sechiba_main
1717	! 0.1 Input
1718	INTEGER(i_std),INTENT(in) :: kjpindex ! number of land points
1719	INTEGER(i_std),INTENT(in) :: ip ! target point
1720	REAL(r_std),DIMENSION(nvm,nvm),INTENT(in) :: matrix_prc ! rotation matrix in % of veget_max
1721	REAL(r_std),DIMENSION(nvm),INTENT(in) :: maxfrac ! veget_max before rotation
1722	LOGICAL, DIMENSION(kjpindex, nvm), INTENT(in) :: in_cycle
1723	INTEGER(i_std), DIMENSION(kjpindex), INTENT(in) :: cyc_num_tot ! number of current rotation cycle
1724	INTEGER(i_std), DIMENSION(kjpindex,nvm,cyc_rot_max), INTENT(in) :: plantdate !! kjpindex, nvm, cyc_rot_max
1725	INTEGER(i_std), DIMENSION(kjpindex, nvm), INTENT(in) :: nrec
1726	INTEGER(i_std), DIMENSION(kjpindex, nvm), INTENT(in) :: nlev
1727
1728	! 0.2 Output
1729	INTEGER(i_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: plantdate_now !! kjpindex, nvm, plantdate of current cycle
1730	REAL(r_std), DIMENSION(ndeep,nvm), INTENT(out) :: soilc_mict
1731	INTEGER(i_std),INTENT(inout), DIMENSION(kjpindex,nvm) :: cyc_num ! flag for starting the rotation for kjpindex, nvm
1732	REAL(r_std), INTENT(inout), DIMENSION(:,:,:,:,:) :: litter ! Above and below ground metabolic and structural litter per ground area
1733	REAL(r_std), INTENT(inout), DIMENSION(:,:,:,:) :: turnover_daily ! Senescence-driven turnover (better: mortality) of leaves and roots
1734	REAL(r_std), INTENT(inout), DIMENSION(:,:,:) :: deepC_a ! deep active carbon profile
1735	REAL(r_std), INTENT(inout), DIMENSION(:,:,:) :: deepC_s ! deep slow carbon profile
1736	REAL(r_std), INTENT(inout), DIMENSION(:,:,:) :: deepC_p ! deep passive carbon profile
1737	REAL(r_std), INTENT(inout), DIMENSION(:,:,:) :: carbon ! Soil carbon pools per ground area: active, slow, or passive
1738	REAL(r_std), INTENT(inout), DIMENSION(:,:) :: age !! Age of PFT it normalized by biomass - can increase and decrease - (years)
1739	REAL(r_std), INTENT(inout), DIMENSION(:,:) :: ind !! Vegetation density, number of individuals per unit
1740	LOGICAL, INTENT(inout), DIMENSION(:,:) :: PFTpresent !! PFT exists (equivalent to veget > 0 for natural PFTs)
1741	LOGICAL, INTENT(inout), DIMENSION(:,:) :: senescence !! The PFT is senescent
1742	REAL(r_std), INTENT(inout), DIMENSION(:,:) :: when_growthinit !! Days since beginning of growing season (days)
1743	REAL(r_std), INTENT(inout), DIMENSION(:,:) :: everywhere !! Is the PFT everywhere in the grid box or very localized
1744	REAL(r_std), INTENT(inout), DIMENSION(:,:,:) :: leaf_frac !! PFT fraction of leaf mass in leaf age class (0-1,
1745
1746
1747	! 0.3 Local Variables
1748	REAL(r_std),DIMENSION(nvm,nlitt,nlevs,nelements) :: dilu_lit !! Litter dilution
1749	REAL(r_std),DIMENSION(nvm,nparts,nelements) :: dilu_turn !! Turnover dilution
1750	REAL(r_std),DIMENSION(nvm,ncarb) :: dilu_soil_carbon
1751	REAL(r_std),DIMENSION(nvm,ndeep,ncarb) :: dilu_soil_carbon_vertres !!vertically-resolved Soil C arbondilution (gC/mÂ²)
1752
1753	REAL(r_std),DIMENSION(ncarb,nvm) :: carbon_old
1754	REAL(r_std),DIMENSION(nlitt,nvm,nlevs,nelements) :: litter_old
1755	REAL(r_std),DIMENSION(nvm,nparts,nelements) :: turnover_old
1756	REAL(r_std), DIMENSION(ndeep,nvm) :: deepC_a_old
1757	REAL(r_std), DIMENSION(ndeep,nvm) :: deepC_s_old
1758	REAL(r_std), DIMENSION(ndeep,nvm) :: deepC_p_old
1759
1760	REAL(r_std),DIMENSION(nvm) :: maxfrac_new ! veget_max after rotation
1761	INTEGER(i_std) :: i,jtar,jsrc,tempcyc
1762	INTEGER(i_std),DIMENSION(nvm) :: cyc_num_old
1763	INTEGER(i_std) :: rngTol = 14
1764	LOGICAL :: f_convert
1765
1766	!!!! -------------------------------------------------------
1767
1768	!1.0 update cyc_num plantdate_now
1769	! temp_cyc = cyc_num
1770	cyc_num_old = cyc_num(ip,:)
1771	!!! source crops
1772	DO jsrc = 2,nvm
1773	IF (SUM(matrix_prc(jsrc,:)) .GT. min_stomate) THEN
1774	!!!! rotation from the jsrc, clean the rotation cycle of this PFT
1775	cyc_num(ip,jsrc) = 1
1776	ENDIF
1777	ENDDO
1778
1779	maxfrac_new = maxfrac
1780	litter_old = litter(ip,:,:,:,:)
1781	turnover_old = turnover_daily(ip,:,:,:)
1782	IF (ok_pc) THEN
1783	deepC_a_old = deepC_a(ip,:,:)
1784	deepC_s_old = deepC_s(ip,:,:)
1785	deepC_p_old = deepC_p(ip,:,:)
1786	ELSE
1787	carbon_old = carbon(ip,:,:)
1788	ENDIF
1789	!!! target crops
1790	DO jtar = 1, nvm
1791	dilu_lit(:,:,:,:) = zero
1792	dilu_turn(:,:,:) = zero
1793	dilu_soil_carbon(:,:) = zero
1794	dilu_soil_carbon_vertres(:,:,:) = zero
1795	tempcyc = cyc_num_old(jtar)
1796	IF (SUM(matrix_prc(:,jtar)) .GT. min_stomate) THEN
1797	DO jsrc = 2, nvm
1798	IF (matrix_prc(jsrc,jtar) .GT. min_stomate) THEN
1799	f_convert = .FALSE.
1800	IF (in_cycle(ip,jtar)) THEN
1801	! a complex case when the target PFT is still growing
1802	! this is a bad case we do not allow
1803	WRITE(numout,*) 'rotate to a growing PFT: jsrc, jtar',jsrc, jtar
1804	WRITE(numout,*) 'in_cycle(ip,jtar) ', in_cycle(ip,jtar)
1805	STOP 'non-permitted rotation case occurred'
1806	ELSE
1807	IF ( cyc_num(ip,jtar) == 1 ) THEN
1808	!! simple case when target PFT is awating start
1809	IF (cyc_num_old(jsrc) .LT. cyc_num_tot(ip)) THEN
1810	cyc_num(ip,jtar) = cyc_num_old(jsrc) + 1
1811	ELSE ! complete of one rotation cycle, back to 1
1812	cyc_num(ip,jtar) = 1
1813	ENDIF
1814	plantdate_now(ip,jtar) = plantdate(ip,jtar,cyc_num(ip,jtar))
1815	IF ( (nrec(ip,jsrc) < nlev(ip,jsrc)) .AND. (nrec(ip,jsrc) .GT. 0) ) THEN
1816	!harvest occured in the next-year of planting, adjusting planting date next cycle
1817	IF (plantdate_now(ip,jtar) > year_length_in_days) plantdate_now(ip,jtar) = plantdate_now(ip,jtar) - year_length_in_days
1818	ENDIF
1819	!!! when next planting date passed but within tolerance range
1820	IF ( (nrec(ip,jsrc) .GE. plantdate_now(ip,jtar)) .AND. (nrec(ip,jsrc) .GT. 0) .AND. &
1821	(nrec(ip,jsrc) .LT. rngTol+plantdate_now(ip,jtar)) ) THEN
1822	plantdate_now(ip,jtar) = nrec(ip,jsrc)+2 ! we still grow it the day after tomorrow
1823	ENDIF
1824	maxfrac_new(jtar) = maxfrac_new(jtar) + matrix_prc(jsrc,jtar) * maxfrac(jsrc)
1825	maxfrac_new(jsrc) = maxfrac_new(jsrc) - matrix_prc(jsrc,jtar) * maxfrac(jsrc)
1826	f_convert = .TRUE.
1827	ELSEIF ( (tempcyc-1 .EQ. cyc_num_old(jsrc)) .OR. (tempcyc==1 .AND. cyc_num_old(jsrc)==cyc_num_tot(ip)) ) THEN
1828	! one cycle ahead waiting to start
1829	! no need to change cyc_num & planting date
1830	maxfrac_new(jtar) = maxfrac_new(jtar) + matrix_prc(jsrc,jtar) * maxfrac(jsrc)
1831	maxfrac_new(jsrc) = maxfrac_new(jsrc) - matrix_prc(jsrc,jtar) * maxfrac(jsrc)
1832	f_convert = .TRUE.
1833	ELSE
1834	WRITE(numout,*) 'rotation situation unexpected:'
1835	WRITE(numout,*) 'cyc_num_old(ip,jsrc), cyc_num_old(ip,jtar)',cyc_num_old(jsrc), cyc_num_old(jtar)
1836	WRITE(numout,*) 'cyc_num(ip,jtar)',cyc_num(ip,jtar)
1837	WRITE(numout,*) 'maxfrac',maxfrac
1838	STOP 'stomate_stics_rotation error with unexpected command'
1839	ENDIF ! tempcyc == 1
1840	ENDIF ! in_cycle(jtar)
1841
1842	IF (f_convert) THEN
1843	!!!! transfering soil and litter/turnover carbon pools
1844	!!! this is actually the better place to modify f_rot_sech & rot_cmd
1845
1846	!!! first is source crops
1847	dilu_lit(jsrc,:,:,:) = litter_old(:,jsrc,:,:)
1848	dilu_turn(jsrc,:,:) = turnover_old(jsrc,:,:)
1849	IF (ok_pc) THEN
1850	dilu_soil_carbon_vertres(jsrc,:,iactive) = deepC_a_old(:,jsrc)
1851	dilu_soil_carbon_vertres(jsrc,:,islow) = deepC_s_old(:,jsrc)
1852	dilu_soil_carbon_vertres(jsrc,:,ipassive) = deepC_p_old(:,jsrc)
1853	ELSE
1854	dilu_soil_carbon(jsrc,:) = carbon_old(:,jsrc)
1855	ENDIF
1856	ENDIF ! we indeed converting jsrc to jtar
1857	ENDIF !jsrc -> jtar
1858	ENDDO !jsrc
1859	!!! then is the conversion of litter and carbon to target crop
1860	! by design, SUM(matrix_prc(jtar,:)), should be either 0 or 1, this
1861	! is simply in case for more complex situations...
1862	litter(ip,:,jtar,:,:) = litter_old(:,jtar,:,:) * maxfrac(jtar) * (1.0 - SUM(matrix_prc(jtar,:)))
1863	turnover_daily(ip,jtar,:,:) = turnover_old(jtar,:,:) * maxfrac(jtar) * (1.0 - SUM(matrix_prc(jtar,:)))
1864	IF (ok_pc) THEN
1865	deepC_a(ip,:,jtar) = deepC_a_old(:,jtar) * maxfrac(jtar) * (1.0 - SUM(matrix_prc(jtar,:)))
1866	deepC_s(ip,:,jtar) = deepC_s_old(:,jtar) * maxfrac(jtar) * (1.0 - SUM(matrix_prc(jtar,:)))
1867	deepC_p(ip,:,jtar) = deepC_p_old(:,jtar) * maxfrac(jtar) * (1.0 - SUM(matrix_prc(jtar,:)))
1868	ELSE
1869	carbon(ip,:,jtar) = carbon_old(:,jtar) * maxfrac(jtar) * (1.0 - SUM(matrix_prc(jtar,:)))
1870	ENDIF
1871	DO jsrc = 1,nvm
1872	litter(ip,:,jtar,:,:) = litter(ip,:,jtar,:,:) + (maxfrac(jsrc) * matrix_prc(jsrc,jtar) * dilu_lit(jsrc,:,:,:))
1873	turnover_daily(ip,jtar,:,:) = turnover_daily(ip,jtar,:,:) + (maxfrac(jsrc) * matrix_prc(jsrc,jtar) * dilu_turn(jsrc,:,:))
1874	IF (ok_pc) THEN
1875	deepC_a(ip,:,jtar) = deepC_a(ip,:,jtar) + (maxfrac(jsrc) * matrix_prc(jsrc,jtar) * dilu_soil_carbon_vertres(jsrc,:,iactive))
1876	deepC_s(ip,:,jtar) = deepC_s(ip,:,jtar) + (maxfrac(jsrc) * matrix_prc(jsrc,jtar) * dilu_soil_carbon_vertres(jsrc,:,islow))
1877	deepC_p(ip,:,jtar) = deepC_p(ip,:,jtar) + (maxfrac(jsrc) * matrix_prc(jsrc,jtar) * dilu_soil_carbon_vertres(jsrc,:,ipassive))
1878	ELSE
1879	carbon(ip,:,jtar) = carbon(ip,:,jtar) + (maxfrac(jsrc) * matrix_prc(jsrc,jtar) * dilu_soil_carbon(jsrc,:))
1880	ENDIF
1881	ENDDO
1882	litter(ip,:,jtar,:,:) = litter(ip,:,jtar,:,:) / maxfrac_new(jtar)
1883	turnover_daily(ip,jtar,:,:) = turnover_daily(ip,jtar,:,:) / maxfrac_new(jtar)
1884	IF (ok_pc) THEN
1885	deepC_a(ip,:,jtar) = deepC_a(ip,:,jtar) / maxfrac_new(jtar)
1886	deepC_s(ip,:,jtar) = deepC_s(ip,:,jtar) / maxfrac_new(jtar)
1887	deepC_p(ip,:,jtar) = deepC_p(ip,:,jtar) / maxfrac_new(jtar)
1888	ELSE
1889	carbon(ip,:,jtar) = carbon(ip,:,jtar) / maxfrac_new(jtar)
1890	ENDIF
1891	age(ip,jtar) = zero
1892
1893	IF ( (maxfrac(jtar) .LT. min_stomate) .AND. (maxfrac_new(jtar) .GT. min_stomate) ) THEN !initialize a previously non-existed PFT
1894	ind(ip,jtar) = maxfrac_new(jtar)
1895	PFTpresent(ip,jtar) = .TRUE.
1896	!everywhere(ip,jtar) = 1.
1897	senescence(ip,jtar) = .TRUE.
1898	age(ip,jtar) = zero
1899	when_growthinit(ip,jtar) = large_value
1900	leaf_frac(ip,jtar,1) = 1.0
1901	ENDIF
1902
1903	ENDIF ! to jtar > 0
1904	ENDDO !jtar
1905	IF (printlev>=4) THEN
1906	WRITE(numout,*) 'maxfrac', maxfrac
1907	WRITE(numout,*) 'maxfrac_new', maxfrac_new
1908	DO i = 1,3
1909	WRITE(numout,*) 'i, carbon_old(i,:)', i, carbon_old(i,:)
1910	WRITE(numout,*) 'i, carbon(ip,i,:)', i, carbon(ip,i,:)
1911	ENDDO
1912	WRITE(numout,*) 'SUM(litter_old(:,:,:,icarbon))',SUM(litter_old(:,:,:,icarbon))
1913	WRITE(numout,*) 'SUM(litter(ip,:,:,:,icarbon))',SUM(litter(ip,:,:,:,icarbon))
1914	ENDIF
1915
1916	DO jsrc = 1,nvm
1917	IF ( maxfrac_new(jsrc) .LT. min_stomate ) THEN
1918	! clean the litter and soil carbon pool of the PFT no longer existed
1919	cyc_num(ip,jsrc) = 1
1920	litter(ip,:,jsrc,:,:) = zero
1921	turnover_daily(ip,jsrc,:,:) = zero
1922	IF (ok_pc) THEN
1923	deepC_a(ip,:,jsrc) = zero
1924	deepC_s(ip,:,jsrc) = zero
1925	deepC_p(ip,:,jsrc) = zero
1926
1927	ELSE
1928	carbon(ip,:,jsrc) = zero
1929	ENDIF
1930	ind(ip,jsrc) = zero
1931	PFTpresent(ip,jsrc) = .FALSE.
1932	senescence(ip,jsrc) = .FALSE.
1933	age(ip,jsrc) = zero
1934	when_growthinit(ip,jsrc) = large_value
1935	everywhere(ip,jsrc) = zero
1936	ENDIF
1937	ENDDO
1938
1939	soilc_mict(:,:) = deepC_a(ip,:,:) + deepC_s(ip,:,:) + deepC_p(ip,:,:)
1940
1941	!!! now check the maxfrac_new for consistency
1942	IF ( SUM(maxfrac_new(:)) .GT. 1.0 ) THEN
1943	WRITE(numout,*) 'vegetation fraction greater than 1.0 after rotation'
1944	WRITE(numout,*) 'ip, maxfrac_new,', ip, maxfrac_new
1945	STOP 'stomate_stics_rotation: wrong vegetation fraction'
1946	ENDIF
1947
1948	END SUBROUTINE stomate_stics_rotation
1949
1950	END MODULE stomate_stics

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

Context Navigation

source: branches/publications/ORCHIDEE-MICT-OP-r6850/src_stomate/stomate_stics.f90 @ 7108

Download in other formats: