Context Navigation

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

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

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

Download in other formats: