Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90 @ 6289

Last change on this file since 6289 was 6145, checked in by josefine.ghattas, 5 years ago
Correction related to VEGET_UPDATE and small fractions as done in the trunk [6034]: Separated into new subroutine slowproc_veget_max_limit the code were veget_max and frac_nobio is set to 0 if smaller than min_vegfrac. This subroutine is now called in slowproc_veget (as before) and also called after the new map has been read in slowproc_readvegetmax. See ticket #456
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 184.9 KB

Rev	Line
[5364]	1
[947]	2	! =================================================================================================================================
	3	! MODULE : slowproc
[8]	4	!
[4470]	5	! CONTACT : orchidee-help _at_ listes.ipsl.fr
[8]	6	!
[947]	7	! LICENCE : IPSL (2006)
	8	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
[8]	9	!
[947]	10	!>\BRIEF Groups the subroutines that: (1) initialize all variables used in
	11	!! slowproc_main, (2) prepare the restart file for the next simulation, (3) Update the
	12	!! vegetation cover if needed, and (4) handle all slow processes if the carbon
	13	!! cycle is activated (call STOMATE) or update the vegetation properties (LAI and
	14	!! fractional cover) in the case of a run with only SECHIBA.
	15	!!
	16	!!\n DESCRIPTION: None
	17	!!
[2419]	18	!! RECENT CHANGE(S): Allowed reading of USDA map, Nov 2014, ADucharne
[947]	19	!!
	20	!! REFERENCE(S) :
	21	!!
	22	!! SVN :
	23	!! $HeadURL$
	24	!! $Date$
	25	!! $Revision$
	26	!! \n
	27	!_ ================================================================================================================================
	28
[8]	29	MODULE slowproc
	30
	31	USE defprec
	32	USE constantes
[947]	33	USE constantes_soil
[511]	34	USE pft_parameters
[8]	35	USE ioipsl
[1788]	36	USE xios_orchidee
[1392]	37	USE ioipsl_para
[4281]	38	USE sechiba_io_p
[8]	39	USE interpol_help
	40	USE stomate
[511]	41	USE stomate_data
[8]	42	USE grid
[4646]	43	USE time, ONLY : dt_sechiba, dt_stomate, one_day, FirstTsYear, LastTsDay
	44	USE time, ONLY : year_start, month_start, day_start, sec_start
	45	USE time, ONLY : month_end, day_end
[1392]	46	USE mod_orchidee_para
[8]	47
	48	IMPLICIT NONE
	49
[947]	50	! Private & public routines
	51
[8]	52	PRIVATE
[5364]	53	PUBLIC slowproc_main, slowproc_clear, slowproc_initialize, slowproc_finalize, slowproc_change_frac, slowproc_xios_initialize
[8]	54
	55	!
	56	! variables used inside slowproc module : declaration and initialisation
	57	!
[2718]	58	REAL(r_std), SAVE :: slope_default = 0.1
[1078]	59	!$OMP THREADPRIVATE(slope_default)
[2718]	60	INTEGER, SAVE :: printlev_loc !! Local printlev in slowproc module
[2348]	61	!$OMP THREADPRIVATE(printlev_loc)
[2718]	62	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: clayfraction !! Clayfraction (0-1, unitless)
[1078]	63	!$OMP THREADPRIVATE(clayfraction)
[4808]	64	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: sandfraction !! Sandfraction (0-1, unitless)
	65	!$OMP THREADPRIVATE(sandfraction)
	66	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: siltfraction !! Siltfraction (0-1, unitless)
	67	!$OMP THREADPRIVATE(siltfraction)
[2718]	68	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: laimap !! LAI map when the LAI is prescribed and not calculated by STOMATE
[1078]	69	!$OMP THREADPRIVATE(laimap)
[2718]	70	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: soilclass_default
[1078]	71	!$OMP THREADPRIVATE(soilclass_default)
[3094]	72	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_max_new !! New year fraction of vegetation type (0-1, unitless)
	73	!$OMP THREADPRIVATE(veget_max_new)
[4657]	74	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: woodharvest !! New year wood harvest
	75	!$OMP THREADPRIVATE(woodharvest)
[3094]	76	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio_new !! New year fraction of ice+lakes+cities+... (0-1, unitless)
	77	!$OMP THREADPRIVATE(frac_nobio_new)
[2718]	78	INTEGER(i_std), SAVE :: lcanop !! canopy levels used for LAI
[2581]	79	!$OMP THREADPRIVATE(lcanop)
[2718]	80	INTEGER(i_std) , SAVE :: veget_year !! year for vegetation update
[2581]	81	!$OMP THREADPRIVATE(veget_year)
	82
[8]	83	CONTAINS
	84
[5364]	85
	86
	87
	88	!! =============================================================================================================================
	89	!! SUBROUTINE: slowproc_xios_initialize
	90	!!
	91	!>\BRIEF Initialize xios dependant defintion before closing context defintion
	92	!!
	93	!! DESCRIPTION: Initialize xios dependant defintion before closing context defintion
	94	!!
	95	!! \n
	96	!_ ==============================================================================================================================
	97
	98	SUBROUTINE slowproc_xios_initialize
	99
	100	CHARACTER(LEN=255) :: filename, name
	101	LOGICAL :: lerr
	102	REAL(r_std) :: slope_noreinf
	103	LOGICAL :: get_slope
	104	CHARACTER(LEN=30) :: veget_str !! update frequency for landuse
	105	INTEGER :: l
	106
	107	IF (printlev>=3) WRITE(numout,*) 'In slowproc_xios_initialize'
	108	!! 1. Prepare for reading of soils_param file
	109	! Get the file name from run.def file and set file attributes accordingly
	110	filename = 'soils_param.nc'
	111	CALL getin_p('SOILCLASS_FILE',filename)
	112	name = filename(1:LEN_TRIM(FILENAME)-3)
	113	CALL xios_orchidee_set_file_attr("soils_param_file",name=name)
	114
	115	! Determine if soils_param_file will be read. If not, deactivate the file.
	116	IF (xios_interpolation .AND. restname_in=='NONE' .AND. .NOT. impsoilt) THEN
	117	! Reading will be done with XIOS later
	118	IF (printlev>=2) WRITE(numout,*) 'Reading of soils_param file will be done later using XIOS. The filename is ', filename
	119	ELSE
	120	! No reading, deactivate soils_param_file
	121	IF (printlev>=2) WRITE(numout,*) 'Reading of soils_param file will not be done with XIOS.'
	122	CALL xios_orchidee_set_file_attr("soils_param_file",enabled=.FALSE.)
	123	CALL xios_orchidee_set_fieldgroup_attr("soil_text",enabled=.FALSE.)
	124	END IF
	125
	126
	127	!! 2. Prepare for reading of PFTmap file
	128	filename = 'PFTmap.nc'
	129	CALL getin_p('VEGETATION_FILE',filename)
	130	name = filename(1:LEN_TRIM(FILENAME)-3)
	131	CALL xios_orchidee_set_file_attr("PFTmap_file",name=name)
	132
	133	! Get veget_update from run.def needed to know if the file needs to be read
	134	veget_update=0
	135	WRITE(veget_str,'(a)') '0Y'
	136	CALL getin_p('VEGET_UPDATE', veget_str)
	137	l=INDEX(TRIM(veget_str),'Y')
	138	READ(veget_str(1:(l-1)),"(I2.2)") veget_update
	139
	140
	141	! Check if PFTmap file will be read by XIOS in this execution
	142	IF ( xios_interpolation .AND. .NOT. impveg .AND. &
	143	(veget_update>0 .OR. restname_in=='NONE')) THEN
	144	! PFTmap will not be read if impveg=TRUE
	145	! PFTmap file will be read each year if veget_update>0
	146	! PFTmap is read if the restart file do not exist and if impveg=F
	147
	148	! Reading will be done
	149	IF (printlev>=2) WRITE(numout,*) 'Reading of PFTmap file will be done later using XIOS. The filename is ', filename
	150	ELSE
	151	! No reading, deactivate PFTmap file
	152	IF (printlev>=2) WRITE(numout,*) 'Reading of PFTmap file will not be done with XIOS.'
	153
	154	CALL xios_orchidee_set_file_attr("PFTmap_file",enabled=.FALSE.)
	155	CALL xios_orchidee_set_field_attr("frac_veget",enabled=.FALSE.)
	156	CALL xios_orchidee_set_field_attr("frac_veget_frac",enabled=.FALSE.)
	157	ENDIF
	158
	159
	160	!! 3. Prepare for reading of topography file
	161	filename = 'cartepente2d_15min.nc'
	162	CALL getin_p('TOPOGRAPHY_FILE',filename)
	163	name = filename(1:LEN_TRIM(FILENAME)-3)
	164	CALL xios_orchidee_set_file_attr("topography_file",name=name)
	165
	166	! Set default values used by XIOS for the interpolation
	167	slope_noreinf = 0.5
	168	CALL getin_p('SLOPE_NOREINF',slope_noreinf)
	169	lerr=xios_orchidee_setvar('slope_noreinf',slope_noreinf)
	170	lerr=xios_orchidee_setvar('slope_default',slope_default)
	171
	172	get_slope = .FALSE.
	173	CALL getin_p('GET_SLOPE',get_slope)
[5454]	174	IF (xios_interpolation .AND. (restname_in=='NONE' .OR. get_slope)) THEN
[5364]	175	! The slope file will be read using XIOS
	176	IF (printlev>=2) WRITE(numout,*) 'Reading of albedo file will be done later using XIOS. The filename is ', filename
	177	ELSE
	178	! Deactivate slope reading
	179	IF (printlev>=2) WRITE(numout,*) 'The slope file will not be read by XIOS'
	180	CALL xios_orchidee_set_file_attr("topography_file",enabled=.FALSE.)
	181	CALL xios_orchidee_set_field_attr("frac_slope_interp",enabled=.FALSE.)
	182	CALL xios_orchidee_set_field_attr("reinf_slope_interp",enabled=.FALSE.)
	183	END IF
	184
	185	!! 4. Prepare for reading of lai file
	186	filename = 'lai2D.nc'
	187	CALL getin_p('LAI_FILE',filename)
	188	name = filename(1:LEN_TRIM(FILENAME)-3)
	189	CALL xios_orchidee_set_file_attr("lai_file",name=name)
	190	! Determine if lai file will be read. If not, deactivate the file.
	191	IF (xios_interpolation .AND. restname_in=='NONE' .AND. read_lai) THEN
	192	! Reading will be done
	193	IF (printlev>=2) WRITE(numout,*) 'Reading of lai file will be done later using XIOS. The filename is ', filename
	194	ELSE
	195	! No reading, deactivate lai file
	196	IF (printlev>=2) WRITE(numout,*) 'Reading of lai file will not be done with XIOS.'
	197	CALL xios_orchidee_set_file_attr("lai_file",enabled=.FALSE.)
	198	CALL xios_orchidee_set_field_attr("frac_lai_interp",enabled=.FALSE.)
	199	CALL xios_orchidee_set_field_attr("lai_interp",enabled=.FALSE.)
	200	END IF
	201
	202	!! 5. Prepare for reading of woodharvest file
	203	filename = 'woodharvest.nc'
	204	CALL getin_p('WOODHARVEST_FILE',filename)
	205	name = filename(1:LEN_TRIM(FILENAME)-3)
	206	CALL xios_orchidee_set_file_attr("woodharvest_file",name=name)
	207
	208	IF (xios_interpolation .AND. do_wood_harvest .AND. &
	209	(veget_update>0 .OR. restname_in=='NONE' )) THEN
	210	! Woodharvest file will be read each year if veget_update>0 or if no restart file exists
	211
	212	! Reading will be done
	213	IF (printlev>=2) WRITE(numout,*) 'Reading of woodharvest file will be done later using XIOS. The filename is ', filename
	214	ELSE
	215	! No reading, deactivate woodharvest file
	216	IF (printlev>=2) WRITE(numout,*) 'Reading of woodharvest file will not be done with XIOS.'
	217	CALL xios_orchidee_set_file_attr("woodharvest_file",enabled=.FALSE.)
	218	CALL xios_orchidee_set_field_attr("woodharvest_interp",enabled=.FALSE.)
	219	ENDIF
	220
	221	IF (printlev_loc>=3) WRITE(numout,*) 'End slowproc_xios_intialize'
	222
	223	END SUBROUTINE slowproc_xios_initialize
	224
	225
[947]	226	!! ================================================================================================================================
[2581]	227	!! SUBROUTINE : slowproc_initialize
	228	!!
	229	!>\BRIEF Initialize slowproc module and call initialization of stomate module
	230	!!
	231	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values
	232	!! Call initialization of stomate module.
	233	!!
	234	!! MAIN OUTPUT VARIABLE(S) :
	235	!!
	236	!! REFERENCE(S) :
	237	!!
	238	!! FLOWCHART : None
	239	!! \n
	240	!_ ================================================================================================================================
	241
[3998]	242	SUBROUTINE slowproc_initialize (kjit, kjpij, kjpindex, &
[2718]	243	rest_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
	244	IndexLand, indexveg, lalo, neighbours, &
[5573]	245	resolution, contfrac, temp_air, &
[2718]	246	soiltile, reinf_slope, deadleaf_cover, assim_param, &
	247	lai, frac_age, height, veget, &
[4723]	248	frac_nobio, njsc, veget_max, fraclut, &
[4882]	249	nwdfraclut, tot_bare_soil,totfrac_nobio, qsintmax, &
[4919]	250	co2_flux, co2_to_bm, fco2_lu, temp_growth)
[2581]	251
	252	!! 0.1 Input variables
	253	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
	254	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
	255	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
[2653]	256	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
[2581]	257	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
	258	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _history_ file identifier
	259	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier
	260	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indices of the points on the land map
	261	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in):: indexveg !! Indices of the points on the vegetation (3D map ???)
	262	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
[3447]	263	INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours !! neighbouring grid points if land.
[2581]	264	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
	265	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid (0-1, unitless)
[5573]	266	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: temp_air !! Air temperature at first atmospheric model layer (K)
[2581]	267
	268	!! 0.2 Output variables
[2653]	269	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: co2_flux !! CO2 flux per average ground area (gC m^{-2} dt_stomate^{-1})
[4919]	270	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: co2_to_bm !! Virtual gpp per average ground area (gC m^{-2} dt_stomate^{-1})
[2653]	271	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_lu !! CO2 flux from land-use (without forest management) (gC m^{-2} dt_stomate^{-1})
[5364]	272	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
[2653]	273	INTEGER(i_std), DIMENSION(kjpindex), INTENT(out) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
	274	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai !! Leaf area index (m^2 m^{-2})
	275	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! height of vegetation (m)
	276	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT(out):: frac_age !! Age efficacity from STOMATE for isoprene
[3969]	277	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget !! Fraction of vegetation type in the mesh (unitless)
	278	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh (unitless)
	279	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_max !! Maximum fraction of vegetation type in the mesh (unitless)
	280	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_bare_soil !! Total evaporating bare soil fraction in the mesh (unitless)
	281	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh (unitless)
	282	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
[4882]	283	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
	284	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: nwdFraclut !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
[2653]	285	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: reinf_slope !! slope coef for reinfiltration
	286	REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (out):: assim_param !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
	287	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: deadleaf_cover !! Fraction of soil covered by dead leaves (unitless)
	288	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintmax !! Maximum water storage on vegetation from interception (mm)
[2581]	289
[4808]	290	!! 0.3 Local variables
	291	INTEGER(i_std) :: jsl
	292	REAL(r_std),DIMENSION (kjpindex,nslm) :: land_frac !! To ouput the clay/sand/silt fractions with a vertical dim
[2581]	293
	294	!_ ================================================================================================================================
	295
[2718]	296	!! 1. Perform the allocation of all variables, define some files and some flags.
[2581]	297	! Restart file read for Sechiba.
[3998]	298	CALL slowproc_init (kjit, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, &
[4882]	299	rest_id, lai, frac_age, veget, frac_nobio, totfrac_nobio, soiltile, fraclut, nwdfraclut, reinf_slope, &
[2718]	300	veget_max, tot_bare_soil, njsc, &
[2581]	301	height, lcanop, veget_update, veget_year)
	302
[2718]	303
	304	!! 2. Define Time step in days for stomate
[2591]	305	dt_days = dt_stomate / one_day
[2581]	306
[2718]	307
	308	!! 3. check time step coherence between slow processes and fast processes
[2591]	309	IF ( dt_stomate .LT. dt_sechiba ) THEN
[4646]	310	WRITE(numout,*) 'slow_processes: time step smaller than forcing time step, dt_sechiba=',dt_sechiba,' dt_stomate=',dt_stomate
[2591]	311	CALL ipslerr_p(3,'slowproc_initialize','Coherence problem between dt_stomate and dt_sechiba',&
[2581]	312	'Time step smaller than forcing time step','')
	313	ENDIF
	314
[2718]	315	!! 4. Call stomate to initialize all variables manadged in stomate,
[2656]	316	IF ( ok_stomate ) THEN
[2581]	317
	318	CALL stomate_initialize &
[2591]	319	(kjit, kjpij, kjpindex, &
	320	rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
[2581]	321	indexLand, lalo, neighbours, resolution, &
[5573]	322	contfrac, totfrac_nobio, clayfraction, temp_air, &
[2718]	323	lai, veget, veget_max, &
[4919]	324	co2_flux, co2_to_bm, fco2_lu, &
	325	deadleaf_cover, assim_param, temp_growth )
[5941]	326	ELSE
	327	!! ok_stomate is not activated
	328	!! Define the CO2 fluxes to zero (no carbone cycle)
	329	co2_flux(:,:) = zero
	330	co2_to_bm(:,:) = zero
[2581]	331	ENDIF
	332
[2718]	333	!! 5. Specific run without the carbon cycle (STOMATE not called):
	334	!! Need to initialize some variables that will be used in SECHIBA:
[2861]	335	!! height, deadleaf_cover, assim_param, qsintmax.
[2581]	336	IF (.NOT. ok_stomate ) THEN
	337	CALL slowproc_derivvar (kjpindex, veget, lai, &
	338	qsintmax, deadleaf_cover, assim_param, height, temp_growth)
	339	ELSE
	340	qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:)
	341	qsintmax(:,1) = zero
	342	ENDIF
[4808]	343
	344
	345	!! 6. Output with XIOS for variables done only once per run
	346
	347	DO jsl=1,nslm
	348	land_frac(:,jsl) = clayfraction(:)
	349	ENDDO
	350	CALL xios_orchidee_send_field("clayfraction",land_frac) ! mean fraction of clay in grid-cell
	351	DO jsl=1,nslm
	352	land_frac(:,jsl) = sandfraction(:)
	353	ENDDO
	354	CALL xios_orchidee_send_field("sandfraction",land_frac) ! mean fraction of sand in grid-cell
	355	DO jsl=1,nslm
	356	land_frac(:,jsl) = siltfraction(:)
	357	ENDDO
	358	CALL xios_orchidee_send_field("siltfraction",land_frac) ! mean fraction of silt in grid-cell
[2581]	359
	360	END SUBROUTINE slowproc_initialize
	361
	362
	363	!! ================================================================================================================================
[947]	364	!! SUBROUTINE : slowproc_main
	365	!!
	366	!>\BRIEF Main routine that manage variable initialisation (slowproc_init),
	367	!! prepare the restart file with the slowproc variables, update the time variables
	368	!! for slow processes, and possibly update the vegetation cover, before calling
	369	!! STOMATE in the case of the carbon cycle activated or just update LAI (and possibly
	370	!! the vegetation cover) for simulation with only SECHIBA
	371	!!
	372	!!
	373	!! DESCRIPTION : (definitions, functional, design, flags): The subroutine manages
	374	!! diverses tasks:
	375	!! (1) Initializing all variables of slowproc (first call)
	376	!! (2) Preparation of the restart file for the next simulation with all prognostic variables
	377	!! (3) Compute and update time variable for slow processes
	378	!! (4) Update the vegetation cover if there is some land use change (only every years)
	379	!! (5) Call STOMATE for the runs with the carbone cycle activated (ok_stomate) and compute the respiration
	380	!! and the net primary production
	381	!! (6) Compute the LAI and possibly update the vegetation cover for run without STOMATE
	382	!!
	383	!! RECENT CHANGE(S): None
	384	!!
	385	!! MAIN OUTPUT VARIABLE(S): ::co2_flux, ::fco2_lu, ::lai, ::height, ::veget, ::frac_nobio,
[4657]	386	!! ::veget_max, ::woodharvest, ::totfrac_nobio, ::soiltype, ::assim_param, ::deadleaf_cover, ::qsintmax,
[947]	387	!! and resp_maint, resp_hetero, resp_growth, npp that are calculated and stored
	388	!! in stomate is activated.
	389	!!
	390	!! REFERENCE(S) : None
	391	!!
	392	!! FLOWCHART :
	393	! \latexonly
	394	! \includegraphics(scale=0.5){SlowprocMainFlow.eps} !PP to be finalize!!)
	395	! \endlatexonly
	396	!! \n
	397	!_ ================================================================================================================================
	398
[3998]	399	SUBROUTINE slowproc_main (kjit, kjpij, kjpindex, &
[4825]	400	IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltile, fraclut, nwdFraclut, &
[4384]	401	temp_air, temp_sol, stempdiag, &
[8]	402	humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
	403	deadleaf_cover, &
	404	assim_param, &
[3687]	405	lai, frac_age, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
[8]	406	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
[4919]	407	co2_flux, fco2_lu, co2_to_bm, temp_growth, tot_bare_soil)
[947]	408
	409	!! INTERFACE DESCRIPTION
[8]	410
[947]	411	!! 0.1 Input variables
[8]	412
[947]	413	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
	414	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
	415	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
	416	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier
	417	INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier
	418	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
	419	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _history_ file identifier
	420	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier
	421	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indices of the points on the land map
	422	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in):: indexveg !! Indices of the points on the vegetation (3D map ???)
	423	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
[3447]	424	INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in) :: neighbours !! neighbouring grid points if land
[947]	425	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
	426	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid (0-1, unitless)
	427	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Relative humidity ("moisture stress") (0-1, unitless)
[4384]	428	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: temp_air !! Temperature of first model layer (K)
[947]	429	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Surface temperature (K)
[4631]	430	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: stempdiag !! Soil temperature (K)
	431	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag !! Relative soil moisture (0-1, unitless)
[947]	432	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: litterhumdiag !! Litter humidity (0-1, unitless)
[2591]	433	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation (mm dt_stomate^{-1})
	434	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation (mm dt_stomate^{-1})
[947]	435	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: gpp !! GPP of total ground area (gC m^{-2} time step^{-1}).
	436	!! Calculated in sechiba, account for vegetation cover and
	437	!! effective time step to obtain gpp_d
	438
	439	!! 0.2 Output variables
[2591]	440	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: co2_flux !! CO2 flux per average ground area (gC m^{-2} dt_stomate^{-1})
[4919]	441	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: co2_to_bm !! virtual gpp flux per average ground area (gC m^{-2} dt_stomate^{-1})
[2591]	442	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_lu !! CO2 flux from land-use (without forest management) (gC m^{-2} dt_stomate^{-1})
[5364]	443	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
[3969]	444	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_bare_soil !! Total evaporating bare soil fraction in the mesh
[947]	445
	446	!! 0.3 Modified variables
	447	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: lai !! Leaf area index (m^2 m^{-2})
	448	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: height !! height of vegetation (m)
	449	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT(inout):: frac_age !! Age efficacity from STOMATE for isoprene
[3969]	450	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: veget !! Fraction of vegetation type including none biological fractionin the mesh (unitless)
[947]	451	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
[3969]	452	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: veget_max !! Maximum fraction of vegetation type in the mesh (unitless)
[947]	453	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh
[3969]	454	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(inout) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
[4723]	455	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(inout) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
[4825]	456	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(inout) :: nwdFraclut !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
[947]	457	REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (inout):: assim_param !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
	458	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: deadleaf_cover !! Fraction of soil covered by dead leaves (unitless)
	459	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintmax !! Maximum water storage on vegetation from interception (mm)
[8]	460
[947]	461	!! 0.4 Local variables
[2718]	462	INTEGER(i_std) :: j, jv, ji !! indices
[2591]	463	REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_maint !! Maitanance component of autotrophic respiration in (gC m^{-2} dt_stomate^{-1})
[947]	464	REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_hetero !! heterotrophic resp. (gC/(m**2 of total ground)/time step)
[2591]	465	REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_growth !! Growth component of autotrophic respiration in gC m^{-2} dt_stomate^{-1})
[947]	466	REAL(r_std), DIMENSION(kjpindex,nvm) :: npp !! Net Ecosystem Exchange (gC/(m**2 of total ground)/time step)
[3831]	467	REAL(r_std),DIMENSION (kjpindex) :: totfrac_nobio_new !! Total fraction for the next year
[4823]	468	REAL(r_std),DIMENSION (kjpindex) :: histvar !! Temporary variable for output
[4646]	469
[947]	470	!_ ================================================================================================================================
[8]	471
[2718]	472	!! 1. Compute and update all variables linked to the date and time
[4646]	473	IF (printlev_loc>=5) WRITE(numout,*) 'Entering slowproc_main, year_start, month_start, day_start, sec_start=',&
	474	year_start, month_start,day_start,sec_start
	475
[2718]	476	!! 2. Activate slow processes if it is the end of the day
	477	IF ( LastTsDay ) THEN
	478	! 3.2.2 Activate slow processes in the end of the day
[8]	479	do_slow = .TRUE.
[2718]	480
[947]	481	! 3.2.3 Count the number of days
[4646]	482	days_since_beg = days_since_beg + 1
	483	IF (printlev_loc>=4) WRITE(numout,*) "New days_since_beg : ",days_since_beg
[8]	484	ELSE
	485	do_slow = .FALSE.
	486	ENDIF
[3831]	487
[2718]	488	!! 3. Update the vegetation if it is time to do so.
	489	!! This is done at the first sechiba time step on a new year and only every "veget_update" years.
	490	!! veget_update correspond to a number of years between each vegetation updates.
	491	!! Nothing is done if veget_update=0.
[4193]	492	!! Update will never be done if impveg=true because veget_update=0.
[5518]	493	IF ( FirstTsYear ) THEN
[4657]	494	IF (veget_update > 0) THEN
	495	veget_year = veget_year + 1
	496
	497	! Update of the vegetation cover with Land Use only if
	498	! the current year match the requested condition (a multiple of "veget_update")
	499	IF ( MOD(veget_year - veget_year_orig, veget_update) == 0 ) THEN
	500	IF (printlev_loc>=1) WRITE(numout,*) 'We are updating the vegetation map for year =' , veget_year
	501
	502	! Read the new the vegetation from file. Output is veget_max_new and frac_nobio_new
	503	CALL slowproc_readvegetmax(kjpindex, lalo, neighbours, resolution, contfrac, &
	504	veget_max, veget_max_new, frac_nobio_new, veget_year, .FALSE.)
	505
[4677]	506	IF (do_wood_harvest) THEN
[4657]	507	! Read the new the wood harvest map from file. Output is wood harvest
	508	CALL slowproc_woodharvest(kjpindex, lalo, neighbours, resolution, contfrac, woodharvest)
	509	ENDIF
	510
	511	! Set the flag do_now_stomate_lcchange to activate stomate_lcchange.
	512	! This flag will be kept to true until stomate_lcchange has been done.
	513	! The variable totfrac_nobio_new will only be used in stomate when this flag is activated
	514	do_now_stomate_lcchange=.TRUE.
	515	IF ( .NOT. ok_stomate ) THEN
	516	! Special case if stomate is not activated : set the variable done_stomate_lcchange=true
	517	! so that the subroutine slowproc_change_frac will be called in the end of sechiba_main.
	518	done_stomate_lcchange=.TRUE.
	519	END IF
	520	ENDIF
[8]	521	ENDIF
[4677]	522	IF ( do_wood_harvest) THEN
[4657]	523	! Set the flag do_now_stomate_woodharvest to activate stomate_woodharvest.
	524	! This flag will be kept to true until stomate_woodharvest has been done.
	525	do_now_stomate_woodharvest=.TRUE.
	526	ENDIF
[8]	527	ENDIF
	528
[2718]	529	!! 4. Main call to STOMATE
[2656]	530	IF ( ok_stomate ) THEN
[3831]	531
	532	! Calculate totfrac_nobio_new only for the case when the land use map has been read previously
	533	IF (do_now_stomate_lcchange) THEN
	534	totfrac_nobio_new(:) = zero
	535	DO jv = 1, nnobio
[3094]	536	totfrac_nobio_new(:) = totfrac_nobio_new(:) + frac_nobio_new(:,jv)
[3831]	537	ENDDO
	538	ELSE
	539	totfrac_nobio_new(:) = zero
	540	END IF
	541
[2718]	542	!! 4.1 Call stomate main routine that will call all c-cycle routines !
[2591]	543	CALL stomate_main (kjit, kjpij, kjpindex, &
[8]	544	IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, &
[4384]	545	temp_air, temp_sol, stempdiag, &
[8]	546	humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
	547	deadleaf_cover, &
	548	assim_param, &
[947]	549	lai, frac_age, height, veget, veget_max, &
[4723]	550	veget_max_new, woodharvest, totfrac_nobio_new, fraclut, &
	551	rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
[4919]	552	co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth,co2_to_bm,temp_growth)
[947]	553
[2718]	554	!! 4.2 Output the respiration terms and the net primary
	555	!! production (NPP) that are calculated in STOMATE
[1788]	556
[2718]	557	! 4.2.1 Output the 3 respiration terms
[4823]	558	! These variables could be output from stomate.
	559	! Variables per pft
[2718]	560	CALL xios_orchidee_send_field("maint_resp",resp_maint/dt_sechiba)
	561	CALL xios_orchidee_send_field("hetero_resp",resp_hetero/dt_sechiba)
	562	CALL xios_orchidee_send_field("growth_resp",resp_growth/dt_sechiba)
[4723]	563
[4823]	564	! Variables on grid-cell
	565	CALL xios_orchidee_send_field("rh_ipcc2",SUM(resp_hetero,dim=2)/dt_sechiba)
	566	histvar(:)=zero
	567	DO jv = 2, nvm
	568	IF ( .NOT. is_tree(jv) .AND. natural(jv) ) THEN
	569	histvar(:) = histvar(:) + resp_hetero(:,jv)
	570	ENDIF
	571	ENDDO
	572	CALL xios_orchidee_send_field("rhGrass",histvar/dt_sechiba)
[4958]	573
[4823]	574	histvar(:)=zero
	575	DO jv = 2, nvm
[4958]	576	IF ( (.NOT. is_tree(jv)) .AND. (.NOT. natural(jv)) ) THEN
	577	histvar(:) = histvar(:) + resp_hetero(:,jv)
	578	ENDIF
	579	ENDDO
	580	CALL xios_orchidee_send_field("rhCrop",histvar/dt_sechiba)
	581
	582	histvar(:)=zero
	583	DO jv = 2, nvm
[4823]	584	IF ( is_tree(jv) ) THEN
	585	histvar(:) = histvar(:) + resp_hetero(:,jv)
	586	ENDIF
	587	ENDDO
	588	CALL xios_orchidee_send_field("rhTree",histvar/dt_sechiba)
	589
	590	! Output with IOIPSL
[2718]	591	CALL histwrite_p(hist_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg)
	592	CALL histwrite_p(hist_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg)
	593	CALL histwrite_p(hist_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg)
	594
	595	! 4.2.2 Compute the net primary production as the diff from
[4723]	596	! Gross primary productin and the growth and maintenance respirations
[2718]	597	npp(:,1)=zero
	598	DO j = 2,nvm
	599	npp(:,j) = gpp(:,j) - resp_growth(:,j) - resp_maint(:,j)
	600	ENDDO
	601
	602	CALL xios_orchidee_send_field("npp",npp/dt_sechiba)
	603
	604	CALL histwrite_p(hist_id, 'npp', kjit, npp, kjpindex*nvm, indexveg)
	605
[8]	606	IF ( hist2_id > 0 ) THEN
[2718]	607	CALL histwrite_p(hist2_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg)
	608	CALL histwrite_p(hist2_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg)
	609	CALL histwrite_p(hist2_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg)
	610	CALL histwrite_p(hist2_id, 'npp', kjit, npp, kjpindex*nvm, indexveg)
[8]	611	ENDIF
[2718]	612
	613	ELSE
	614	!! ok_stomate is not activated
	615	!! Define the CO2 flux from the grid point to zero (no carbone cycle)
	616	co2_flux(:,:) = zero
[4919]	617	co2_to_bm(:,:) = zero
[8]	618	ENDIF
[947]	619
[2718]	620
	621	!! 5. Do daily processes if necessary
	622	!!
	623	IF ( do_slow ) THEN
[8]	624
[2718]	625	!! 5.1 Calculate the LAI if STOMATE is not activated
	626	IF ( .NOT. ok_stomate ) THEN
[8]	627	CALL slowproc_lai (kjpindex, lcanop,stempdiag, &
[4646]	628	lalo,resolution,lai,laimap)
[2718]	629
[890]	630	frac_age(:,:,1) = un
	631	frac_age(:,:,2) = zero
	632	frac_age(:,:,3) = zero
	633	frac_age(:,:,4) = zero
[947]	634	ENDIF
[8]	635
[2718]	636	!! 5.2 Update veget
[4825]	637	CALL slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
[947]	638
[2718]	639	!! 5.3 updates qsintmax and other derived variables
[2548]	640	IF ( .NOT. ok_stomate ) THEN
[8]	641	CALL slowproc_derivvar (kjpindex, veget, lai, &
[2031]	642	qsintmax, deadleaf_cover, assim_param, height, temp_growth)
[8]	643	ELSE
[947]	644	qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:)
	645	qsintmax(:,1) = zero
	646	ENDIF
	647	END IF
[8]	648
[3969]	649	!! 6. Calculate tot_bare_soil needed in hydrol, diffuco and condveg (fraction in the mesh)
[2718]	650	tot_bare_soil(:) = veget_max(:,1)
	651	DO jv = 2, nvm
	652	DO ji =1, kjpindex
	653	tot_bare_soil(ji) = tot_bare_soil(ji) + (veget_max(ji,jv) - veget(ji,jv))
	654	ENDDO
	655	END DO
	656
[8]	657
[3049]	658	!! 7. Do some basic tests on the surface fractions updated above, only if
	659	!! slowproc_veget has been done (do_slow). No change of the variables.
	660	IF (do_slow) THEN
	661	CALL slowproc_checkveget(kjpindex, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
	662	END IF
[1788]	663
[2718]	664	!! 8. Write output fields
	665	CALL xios_orchidee_send_field("tot_bare_soil",tot_bare_soil)
	666
	667	IF ( .NOT. almaoutput) THEN
	668	CALL histwrite_p(hist_id, 'tot_bare_soil', kjit, tot_bare_soil, kjpindex, IndexLand)
[947]	669	END IF
[8]	670
[2718]	671
[2348]	672	IF (printlev_loc>=3) WRITE (numout,*) ' slowproc_main done '
[8]	673
	674	END SUBROUTINE slowproc_main
[947]	675
[2718]	676
[2581]	677	!! ================================================================================================================================
	678	!! SUBROUTINE : slowproc_finalize
	679	!!
	680	!>\BRIEF Write to restart file variables for slowproc module and call finalization of stomate module
	681	!!
	682	!! DESCRIPTION :
	683	!!
	684	!! MAIN OUTPUT VARIABLE(S) :
	685	!!
	686	!! REFERENCE(S) :
	687	!!
	688	!! FLOWCHART : None
	689	!! \n
	690	!_ ================================================================================================================================
[947]	691
[2581]	692	SUBROUTINE slowproc_finalize (kjit, kjpindex, rest_id, IndexLand, &
	693	njsc, lai, height, veget, &
[3221]	694	frac_nobio, veget_max, reinf_slope, &
[4919]	695	co2_to_bm, assim_param, frac_age )
[2581]	696
	697	!! 0.1 Input variables
[2653]	698	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
	699	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
	700	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
	701	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indices of the points on the land map
	702	INTEGER(i_std), DIMENSION(kjpindex), INTENT(in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
	703	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
	704	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: height !! height of vegetation (m)
	705	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type including none biological fraction (unitless)
	706	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
	707	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
	708	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: reinf_slope !! slope coef for reinfiltration
[4919]	709	REAL(r_std),DIMENSION (kjpindex,nvm),INTENT(in) :: co2_to_bm !! virtual gpp flux between atmosphere and biosphere
[2868]	710	REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (in):: assim_param !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
[3221]	711	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT(in):: frac_age !! Age efficacity from STOMATE for isoprene
[2581]	712	!! 0.4 Local variables
[2653]	713	REAL(r_std) :: tmp_day(1) !! temporary variable for I/O
[3221]	714	INTEGER :: jf !! Indice
	715	CHARACTER(LEN=4) :: laistring !! Temporary character string
	716	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
[2581]	717	!_ ================================================================================================================================
	718
	719	IF (printlev_loc>=3) WRITE (numout,*) 'Write restart file with SLOWPROC variables '
	720
	721	! 2.1 Write a series of variables controled by slowproc: day
	722	! counter, vegetation fraction, max vegetation fraction, LAI
	723	! variable from stomate, fraction of bare soil, soiltype
	724	! fraction, clay fraction, height of vegetation, map of LAI
	725
[2654]	726	CALL restput_p (rest_id, 'veget', nbp_glo, nvm, 1, kjit, veget, 'scatter', nbp_glo, index_g)
[2581]	727
	728	CALL restput_p (rest_id, 'veget_max', nbp_glo, nvm, 1, kjit, veget_max, 'scatter', nbp_glo, index_g)
	729
[4677]	730	IF (do_wood_harvest) THEN
[4657]	731	CALL restput_p (rest_id, 'woodharvest', nbp_glo, 1, 1, kjit, woodharvest, 'scatter', nbp_glo, index_g)
	732	END IF
	733
[2581]	734	CALL restput_p (rest_id, 'lai', nbp_glo, nvm, 1, kjit, lai, 'scatter', nbp_glo, index_g)
	735
	736	CALL restput_p (rest_id, 'frac_nobio', nbp_glo, nnobio, 1, kjit, frac_nobio, 'scatter', nbp_glo, index_g)
[3221]	737
	738
	739	DO jf = 1, nleafages
	740	! variable name is somewhat complicated as ioipsl does not allow 3d variables for the moment...
	741	WRITE(laistring,'(i4)') jf
	742	laistring=ADJUSTL(laistring)
	743	var_name='frac_age_'//laistring(1:LEN_TRIM(laistring))
	744	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, frac_age(:,:,jf), 'scatter', nbp_glo, index_g)
	745	ENDDO
	746
[4428]	747	! Add the soil_classif as suffix for the variable name of njsc when it is stored in the restart file.
	748	IF (soil_classif == 'zobler') THEN
	749	var_name= 'njsc_zobler'
	750	ELSE IF (soil_classif == 'usda') THEN
	751	var_name= 'njsc_usda'
	752	END IF
	753	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, REAL(njsc, r_std), 'scatter', nbp_glo, index_g)
[5454]	754	CALL restput_p (rest_id, 'reinf_slope', nbp_glo, 1, 1, kjit, reinf_slope, 'scatter', nbp_glo, index_g)
[2581]	755	CALL restput_p (rest_id, 'clay_frac', nbp_glo, 1, 1, kjit, clayfraction, 'scatter', nbp_glo, index_g)
[4871]	756	CALL restput_p (rest_id, 'sand_frac', nbp_glo, 1, 1, kjit, sandfraction, 'scatter', nbp_glo, index_g)
[2581]	757	!
	758	! The height of the vegetation could in principle be recalculated at the beginning of the run.
	759	! However, this is very tedious, as many special cases have to be taken into account. This variable
	760	! is therefore saved in the restart file.
	761	CALL restput_p (rest_id, 'height', nbp_glo, nvm, 1, kjit, height, 'scatter', nbp_glo, index_g)
	762	!
	763	! Specific case where the LAI is read and not calculated by STOMATE: need to be saved
	764	IF (read_lai) THEN
	765	CALL restput_p (rest_id, 'laimap', nbp_glo, nvm, 12, kjit, laimap)
	766	ENDIF
	767	!
	768	! If there is some land use change, write the year for the land use ???
[5518]	769	tmp_day(1) = REAL(veget_year,r_std)
	770	IF (is_root_prc) CALL restput (rest_id, 'veget_year', 1 , 1 , 1, kjit, tmp_day)
	771
[2656]	772	! 2.2 Write restart variables managed by STOMATE
	773	IF ( ok_stomate ) THEN
[4919]	774	CALL stomate_finalize (kjit, kjpindex, indexLand, clayfraction, co2_to_bm, assim_param)
[2581]	775	ENDIF
	776
	777	END SUBROUTINE slowproc_finalize
	778
	779
[947]	780	!! ================================================================================================================================
	781	!! SUBROUTINE : slowproc_init
[8]	782	!!
[947]	783	!>\BRIEF Initialisation of all variables linked to SLOWPROC
[8]	784	!!
[947]	785	!! DESCRIPTION : (definitions, functional, design, flags): The subroutine manages
	786	!! diverses tasks:
[8]	787	!!
[947]	788	!! RECENT CHANGE(S): None
	789	!!
	790	!! MAIN OUTPUT VARIABLE(S): ::lcanop, ::veget_update, ::veget_year,
	791	!! ::lai, ::veget, ::frac_nobio, ::totfrac_nobio, ::veget_max, ::height, ::soiltype
	792	!!
	793	!! REFERENCE(S) : None
	794	!!
	795	!! FLOWCHART : None
	796	!! \n
	797	!_ ================================================================================================================================
[8]	798
[3998]	799	SUBROUTINE slowproc_init (kjit, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, &
[4882]	800	rest_id, lai, frac_age, veget, frac_nobio, totfrac_nobio, soiltile, fraclut, nwdfraclut, reinf_slope, &
[2740]	801	veget_max, tot_bare_soil, njsc, &
	802	height, lcanop, veget_update, veget_year)
	803
	804	!! INTERFACE DESCRIPTION
[8]	805
[2740]	806	!! 0.1 Input variables
	807	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
	808	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - Terrestrial pixels only
	809	INTEGER(i_std), INTENT (in) :: rest_id !! Restart file identifier
[947]	810
[2740]	811	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indices of the land points on the map
	812	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
[3447]	813	INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours !! Vector of neighbours for each grid point
	814	!! (1=North and then clockwise)
[2740]	815	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x and y of the grid (m)
	816	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid (unitless)
	817
	818	!! 0.2 Output variables
	819	INTEGER(i_std), INTENT(out) :: lcanop !! Number of Canopy level used to compute LAI
	820	INTEGER(i_std), INTENT(out) :: veget_update !! update frequency in timesteps (years) for landuse
	821	INTEGER(i_std), INTENT(out) :: veget_year !! first year for landuse (year or index ???)
	822
	823	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai !! Leaf Area index (m^2 / m^2)
[3969]	824	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget !! Fraction of vegetation type in the mesh (unitless)
	825	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: frac_nobio !! Fraction of ice,lakes,cities, ... in the mesh (unitless)
	826	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... in the mesh (unitless)
	827	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_max !! Max fraction of vegetation type in the mesh (unitless)
	828	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_bare_soil !! Total evaporating bare soil fraction in the mesh
[2740]	829	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! Height of vegetation or surface in genral ??? (m)
	830	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT (out):: frac_age !! Age efficacity from STOMATE for isoprene
[3969]	831	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
[4723]	832	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: fraclut !! Fraction of each landuse tile
[4882]	833	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: nwdfraclut !! Fraction of non woody vegetation in each landuse tile
[2740]	834	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: reinf_slope !! slope coef for reinfiltration
	835	INTEGER(i_std), DIMENSION(kjpindex), INTENT(out) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
	836
	837	!! 0.3 Local variables
	838	REAL(r_std) :: tmp_veget_year(1) !! temporary variable
	839	REAL(r_std) :: zcanop !! ???? soil depth taken for canopy
	840	INTEGER(i_std) :: vtmp(1) !! temporary variable
[4631]	841	REAL(r_std), DIMENSION(nslm) :: zsoil !! soil depths at diagnostic levels
[3221]	842	CHARACTER(LEN=4) :: laistring !! Temporary character string
[3831]	843	INTEGER(i_std) :: l, jf !! Indices
[2740]	844	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
	845	INTEGER(i_std) :: ji, jv, ier,jst !! Indices
	846	LOGICAL :: get_slope
	847	REAL(r_std) :: frac_nobio1 !! temporary variable for frac_nobio(see above)
	848	REAL(r_std), DIMENSION(kjpindex) :: tmp_real
[4631]	849	REAL(r_std), DIMENSION(kjpindex,nslm) :: stempdiag2_bid !! matrix to store stempdiag_bid
[2740]	850	REAL(r_std), DIMENSION (kjpindex,nscm) :: soilclass !! Fractions of each soil textural class in the grid cell (0-1, unitless)
	851	CHARACTER(LEN=30), SAVE :: veget_str !! update frequency for landuse
	852	!$OMP THREADPRIVATE(veget_str)
	853	REAL(r_std), DIMENSION(kjpindex) :: frac_crop_tot !! Total fraction occupied by crops (0-1, unitless)
	854	LOGICAL :: found_restart !! found_restart=true if all 3 variables veget_max, veget and
	855	!! frac_nobio are read from restart file
	856	!_ ================================================================================================================================
	857
	858	!! 0. Initialize local printlev
[2348]	859	printlev_loc=get_printlev('slowproc')
	860	IF (printlev_loc>=3) WRITE (numout,*) "In slowproc_init"
[2740]	861
	862
	863	!! 1. Allocation
[8]	864
	865	ALLOCATE (clayfraction(kjpindex),stat=ier)
[2613]	866	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable clayfraction','','')
[8]	867	clayfraction(:)=undef_sechiba
[4808]	868
	869	ALLOCATE (sandfraction(kjpindex),stat=ier)
	870	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable sandfraction','','')
	871	sandfraction(:)=undef_sechiba
	872
	873	ALLOCATE (siltfraction(kjpindex),stat=ier)
	874	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable siltfraction','','')
	875	siltfraction(:)=undef_sechiba
[947]	876
	877	! Initialisation of the fraction of the different vegetation: Start with 100% of bare soil
	878	ALLOCATE (soilclass_default(nscm),stat=ier)
[2613]	879	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable soilclass_default','','')
[947]	880	soilclass_default(:)=undef_sechiba
[2613]	881
[2718]	882	! Allocation of last year vegetation fraction in case of land use change
[3094]	883	ALLOCATE(veget_max_new(kjpindex, nvm), STAT=ier)
	884	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable veget_max_new','','')
[947]	885
[4657]	886	! Allocation of wood harvest
	887	ALLOCATE(woodharvest(kjpindex), STAT=ier)
	888	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable woodharvest','','')
	889
[947]	890	! Allocation of the fraction of non biospheric areas
[3094]	891	ALLOCATE(frac_nobio_new(kjpindex, nnobio), STAT=ier)
	892	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable frac_nobio_new','','')
[947]	893
[2740]	894	! Allocate laimap
	895	IF (read_lai)THEN
	896	ALLOCATE (laimap(kjpindex,nvm,12),stat=ier)
	897	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable laimap','','')
	898	ELSE
	899	ALLOCATE (laimap(1,1,1), stat=ier)
	900	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable laimap(1,1,1)','','')
	901	ENDIF
	902
	903
	904	!! 2. Read variables from restart file
	905
[2654]	906	found_restart=.TRUE.
	907	var_name= 'veget'
	908	CALL ioconf_setatt_p('UNITS', '-')
	909	CALL ioconf_setatt_p('LONG_NAME','Vegetation fraction')
	910	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget, "gather", nbp_glo, index_g)
	911	IF ( ALL( veget(:,:) .EQ. val_exp ) ) found_restart=.FALSE.
	912
[8]	913	var_name= 'veget_max'
[1078]	914	CALL ioconf_setatt_p('UNITS', '-')
	915	CALL ioconf_setatt_p('LONG_NAME','Maximum vegetation fraction')
[8]	916	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget_max, "gather", nbp_glo, index_g)
[2654]	917	IF ( ALL( veget_max(:,:) .EQ. val_exp ) ) found_restart=.FALSE.
	918
[4677]	919	IF (do_wood_harvest) THEN
[4657]	920	var_name= 'woodharvest'
	921	CALL ioconf_setatt_p('UNITS', 'gC m-2 yr-1')
	922	CALL ioconf_setatt_p('LONG_NAME','Harvest wood biomass')
	923	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., woodharvest, "gather", nbp_glo, index_g)
	924	IF ( ALL( woodharvest(:) .EQ. val_exp ) ) woodharvest(:)=zero
	925	END IF
	926
[8]	927	var_name= 'frac_nobio'
[1078]	928	CALL ioconf_setatt_p('UNITS', '-')
	929	CALL ioconf_setatt_p('LONG_NAME','Special soil type fraction')
[8]	930	CALL restget_p (rest_id, var_name, nbp_glo, nnobio, 1, kjit, .TRUE., frac_nobio, "gather", nbp_glo, index_g)
[2654]	931	IF ( ALL( frac_nobio(:,:) .EQ. val_exp ) ) found_restart=.FALSE.
[947]	932
[5518]	933	IF (.NOT. impveg) THEN
[8]	934	var_name= 'veget_year'
[1078]	935	CALL ioconf_setatt_p('UNITS', '-')
	936	CALL ioconf_setatt_p('LONG_NAME','Last year get in Land Use file.')
[8]	937	IF (is_root_prc) THEN
	938	CALL restget (rest_id, var_name, 1 , 1 , 1, kjit, .TRUE., tmp_veget_year)
[3434]	939	IF (veget_reinit) THEN
	940	! Do not take the value read from restart file
[8]	941	veget_year=veget_year_orig
[3434]	942	ELSE IF (tmp_veget_year(1) == val_exp) THEN
	943	! veget_year was not found in restart file
	944	veget_year=veget_year_orig
[8]	945	ELSE
[3434]	946	! veget_year was found in restart file, transform to integer
	947	veget_year=INT(tmp_veget_year(1))
[8]	948	ENDIF
	949	ENDIF
	950	CALL bcast(veget_year)
[947]	951
[8]	952	!
[566]	953	!Config Key = VEGET_UPDATE
[5605]	954	!Config Desc = Update vegetation frequency: 0Y or 1Y
[5518]	955	!Config If =
[566]	956	!Config Def = 0Y
[5518]	957	!Config Help = The veget datas will be update each this time step. Must be 0Y if IMPOSE_VEG=y.
[843]	958	!Config Units = [years]
[8]	959	!
	960	veget_update=0
	961	WRITE(veget_str,'(a)') '0Y'
	962	CALL getin_p('VEGET_UPDATE', veget_str)
	963	l=INDEX(TRIM(veget_str),'Y')
	964	READ(veget_str(1:(l-1)),"(I2.2)") veget_update
[4693]	965	IF (printlev_loc >= 2) WRITE(numout,*) "Update frequency for land use in years :",veget_update
[3165]	966
	967	! Coherence test
	968	IF (veget_update > 0 .AND. ok_dgvm .AND. .NOT. agriculture) THEN
	969	CALL ipslerr_p(3,'slowproc_init',&
	970	'The combination DGVM=TRUE, AGRICULTURE=FALSE and VEGET_UPDATE>0 is not possible', &
	971	'Set VEGET_UPDATE=0Y in run.def','')
	972	END IF
[2718]	973	ELSE
[5518]	974	! impveg=TRUE: there can not be any land use change, veget_update must be =0
[4193]	975	! Read VEGET_UPDATE from run.def and exit if it is different from 0Y
[2718]	976	veget_update=0
[4193]	977	WRITE(veget_str,'(a)') '0Y'
	978	CALL getin_p('VEGET_UPDATE', veget_str)
	979	l=INDEX(TRIM(veget_str),'Y')
	980	READ(veget_str(1:(l-1)),"(I2.2)") veget_update
	981	IF (veget_update /= 0) THEN
[5518]	982	WRITE(numout,*) 'veget_update=',veget_update,' is not coeherent with impveg=',impveg
	983	CALL ipslerr_p(3,'slowproc_init','Incoherent values between impveg and veget_update', &
	984	'VEGET_UPDATE must be equal to 0Y if IMPOSE_VEG=y (impveg=true)','')
[4193]	985	END IF
	986
[8]	987	ENDIF
[947]	988
[5454]	989	var_name= 'reinf_slope'
	990	CALL ioconf_setatt_p('UNITS', '-')
	991	CALL ioconf_setatt_p('LONG_NAME','Slope coef for reinfiltration')
	992	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., reinf_slope, "gather", nbp_glo, index_g)
[2718]	993
[5454]	994
[2718]	995	! Below we define the soil texture of the grid-cells
[4428]	996	! Add the soil_classif as suffix for the variable name of njsc when it is stored in the restart file.
	997	IF (soil_classif == 'zobler') THEN
	998	var_name= 'njsc_zobler'
	999	ELSE IF (soil_classif == 'usda') THEN
	1000	var_name= 'njsc_usda'
	1001	ELSE
	1002	CALL ipslerr_p(3,'slowproc_init','Non supported soil type classification','','')
	1003	END IF
	1004
[2718]	1005	CALL ioconf_setatt_p('UNITS', '-')
	1006	CALL ioconf_setatt_p('LONG_NAME','Index of soil type')
	1007	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., tmp_real, "gather", nbp_glo, index_g)
[2740]	1008	IF ( ALL( tmp_real(:) .EQ. val_exp) ) THEN
	1009	njsc (:) = undef_int
	1010	ELSE
[2718]	1011	njsc = NINT(tmp_real)
[2740]	1012	END IF
[2718]	1013
[8]	1014	var_name= 'clay_frac'
[1078]	1015	CALL ioconf_setatt_p('UNITS', '-')
	1016	CALL ioconf_setatt_p('LONG_NAME','Fraction of clay in each mesh')
[8]	1017	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., clayfraction, "gather", nbp_glo, index_g)
[2740]	1018
[4871]	1019	var_name= 'sand_frac'
	1020	CALL ioconf_setatt_p('UNITS', '-')
	1021	CALL ioconf_setatt_p('LONG_NAME','Fraction of sand in each mesh')
	1022	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., sandfraction, "gather", nbp_glo, index_g)
	1023
	1024	! Calculate siltfraction not needed to be in restart file
	1025	IF ( ALL( sandfraction(:) .EQ. val_exp) ) THEN
	1026	siltfraction(:) = val_exp
	1027	ELSE
	1028	siltfraction(:) = 1. - clayfraction(:) - sandfraction(:)
	1029	END IF
	1030
[8]	1031	var_name= 'lai'
[1078]	1032	CALL ioconf_setatt_p('UNITS', '-')
	1033	CALL ioconf_setatt_p('LONG_NAME','Leaf area index')
[8]	1034	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., lai, "gather", nbp_glo, index_g)
[2740]	1035
[8]	1036	! The height of the vegetation could in principle be recalculated at the beginning of the run.
	1037	! However, this is very tedious, as many special cases have to be taken into account. This variable
	1038	! is therefore saved in the restart file.
	1039	var_name= 'height'
[1078]	1040	CALL ioconf_setatt_p('UNITS', 'm')
	1041	CALL ioconf_setatt_p('LONG_NAME','Height of vegetation')
[8]	1042	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., height, "gather", nbp_glo, index_g)
[2740]	1043
[8]	1044	IF (read_lai)THEN
	1045	var_name= 'laimap'
[1078]	1046	CALL ioconf_setatt_p('UNITS', '-')
	1047	CALL ioconf_setatt_p('LONG_NAME','Leaf area index read')
[8]	1048	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 12, kjit, .TRUE., laimap)
	1049	ENDIF
[2740]	1050
[3221]	1051	CALL ioconf_setatt_p('UNITS', '-')
	1052	CALL ioconf_setatt_p('LONG_NAME','Fraction of leaves in leaf age class ')
	1053	DO jf = 1, nleafages
	1054	! variable name is somewhat complicated as ioipsl does not allow 3d variables for the moment...
	1055	WRITE(laistring,'(i4)') jf
	1056	laistring=ADJUSTL(laistring)
	1057	var_name='frac_age_'//laistring(1:LEN_TRIM(laistring))
	1058	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE.,frac_age(:,:,jf), "gather", nbp_glo, index_g)
	1059	ENDDO
[2740]	1060
	1061	!! 3. Some other initializations
	1062
[566]	1063	!Config Key = SECHIBA_ZCANOP
[846]	1064	!Config Desc = Soil level used for canopy development (if STOMATE disactivated)
	1065	!Config If = OK_SECHIBA and .NOT. OK_STOMATE
[566]	1066	!Config Def = 0.5
	1067	!Config Help = The temperature at this soil depth is used to determine the LAI when
	1068	!Config STOMATE is not activated.
[846]	1069	!Config Units = [m]
[8]	1070	zcanop = 0.5_r_std
	1071	CALL setvar_p (zcanop, val_exp, 'SECHIBA_ZCANOP', 0.5_r_std)
[2740]	1072
[8]	1073	! depth at center of the levels
	1074	zsoil(1) = diaglev(1) / 2.
[4631]	1075	DO l = 2, nslm
[8]	1076	zsoil(l) = ( diaglev(l) + diaglev(l-1) ) / 2.
	1077	ENDDO
	1078
	1079	! index of this level
	1080	vtmp = MINLOC ( ABS ( zcanop - zsoil(:) ) )
	1081	lcanop = vtmp(1)
	1082
	1083	!
[947]	1084	! Interception reservoir coefficient
[8]	1085	!
[947]	1086	!Config Key = SECHIBA_QSINT
	1087	!Config Desc = Interception reservoir coefficient
	1088	!Config If = OK_SECHIBA
	1089	!Config Def = 0.1
	1090	!Config Help = Transforms leaf area index into size of interception reservoir
	1091	!Config for slowproc_derivvar or stomate
	1092	!Config Units = [m]
	1093	CALL getin_p('SECHIBA_QSINT', qsintcst)
[4693]	1094	IF (printlev >= 2) WRITE(numout, *)' SECHIBA_QSINT, qsintcst = ', qsintcst
[947]	1095
[8]	1096
	1097
[2740]	1098
	1099	!! 4. Initialization of variables not found in restart file
	1100
[8]	1101	IF ( impveg ) THEN
[2740]	1102
	1103	!! 4.1.a Case impveg=true: Initialization of variables by reading run.def
	1104	!! The routine setvar_p will only initialize the variable if it was not found in restart file.
	1105	!! We are on a point and thus we can read the information from the run.def
	1106
[566]	1107	!Config Key = SECHIBA_VEGMAX
	1108	!Config Desc = Maximum vegetation distribution within the mesh (0-dim mode)
	1109	!Config If = IMPOSE_VEG
	1110	!Config Def = 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0
	1111	!Config Help = The fraction of vegetation is read from the restart file. If
	1112	!Config it is not found there we will use the values provided here.
	1113	!Config Units = [-]
[8]	1114	CALL setvar_p (veget_max, val_exp, 'SECHIBA_VEGMAX', veget_ori_fixed_test_1)
	1115
[566]	1116	!Config Key = SECHIBA_FRAC_NOBIO
	1117	!Config Desc = Fraction of other surface types within the mesh (0-dim mode)
	1118	!Config If = IMPOSE_VEG
	1119	!Config Def = 0.0
	1120	!Config Help = The fraction of ice, lakes, etc. is read from the restart file. If
	1121	!Config it is not found there we will use the values provided here.
	1122	!Config For the moment, there is only ice.
	1123	!Config Units = [-]
[8]	1124	frac_nobio1 = frac_nobio(1,1)
	1125	CALL setvar_p (frac_nobio1, val_exp, 'SECHIBA_FRAC_NOBIO', frac_nobio_fixed_test_1)
	1126	frac_nobio(:,:) = frac_nobio1
	1127
[2740]	1128	IF (.NOT. found_restart) THEN
	1129	! Call slowproc_veget to correct veget_max and to calculate veget and soiltiles
[4825]	1130	CALL slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
[2740]	1131	END IF
	1132
[566]	1133	!Config Key = SECHIBA_LAI
	1134	!Config Desc = LAI for all vegetation types (0-dim mode)
	1135	!Config Def = 0., 8., 8., 4., 4.5, 4.5, 4., 4.5, 4., 2., 2., 2., 2.
	1136	!Config If = IMPOSE_VEG
	1137	!Config Help = The maximum LAI used in the 0dim mode. The values should be found
	1138	!Config in the restart file. The new values of LAI will be computed anyway
	1139	!Config at the end of the current day. The need for this variable is caused
	1140	!Config by the fact that the model may stop during a day and thus we have not
	1141	!Config yet been through the routines which compute the new surface conditions.
	1142	!Config Units = [-]
[8]	1143	CALL setvar_p (lai, val_exp, 'SECHIBA_LAI', llaimax)
	1144
[44]	1145	IF (impsoilt) THEN
[8]	1146
[1079]	1147	! If njsc is not in restart file, then initialize soilclass from values
	1148	! from run.def file and recalculate njsc
	1149	IF ( ALL(njsc(:) .EQ. undef_int )) THEN
	1150	!Config Key = SOIL_FRACTIONS
	1151	!Config Desc = Fraction of the 3 soil types (0-dim mode)
	1152	!Config Def = undef_sechiba
	1153	!Config If = IMPOSE_VEG and IMPOSE_SOILT
	1154	!Config Help = Determines the fraction for the 3 soil types
	1155	!Config in the mesh in the following order : sand loam and clay.
	1156	!Config Units = [-]
	1157
	1158	soilclass(1,:) = soilclass_default(:)
	1159	CALL getin_p('SOIL_FRACTIONS',soilclass(1,:))
[2419]	1160	! Assign for each grid-cell the % of the different textural classes (up to 12 if 'usda')
[1079]	1161	DO ji=2,kjpindex
[2589]	1162	! here we read, for the prescribed grid-cell, the % occupied by each of the soil texture classes
[1079]	1163	soilclass(ji,:) = soilclass(1,:)
	1164	ENDDO
	1165
[2419]	1166	! Simplify an heterogeneous grid-cell into an homogeneous one with the dominant texture
[1079]	1167	njsc(:) = 0
	1168	DO ji = 1, kjpindex
[2589]	1169	! here we reduce to the dominant texture class
[1079]	1170	njsc(ji) = MAXLOC(soilclass(ji,:),1)
	1171	ENDDO
	1172	END IF
	1173
[566]	1174	!Config Key = CLAY_FRACTION
	1175	!Config Desc = Fraction of the clay fraction (0-dim mode)
	1176	!Config Def = 0.2
	1177	!Config If = IMPOSE_VEG and IMPOSE_SOIL
	1178	!Config Help = Determines the fraction of clay in the grid box.
	1179	!Config Units = [-]
[2419]	1180
	1181	! If clayfraction was not in restart file it will be read fro run.def file instead of deduced
	1182	! based on fractions of each textural class
[44]	1183	CALL setvar_p (clayfraction, val_exp, 'CLAY_FRACTION', clayfraction_default)
[4871]	1184
	1185	!Config Key = SAND_FRACTION
	1186	!Config Desc = Fraction of the clay fraction (0-dim mode)
	1187	!Config Def = 0.4
	1188	!Config If = IMPOSE_VEG and IMPOSE_SOIL
	1189	!Config Help = Determines the fraction of clay in the grid box.
	1190	!Config Units = [-]
	1191
	1192	! If sand fraction was not in restart file it will be read fro run.def file
	1193	CALL setvar_p (sandfraction, val_exp, 'SAND_FRACTION', sandfraction_default)
	1194
	1195	! Calculate silt fraction
	1196	siltfraction(:) = 1. - clayfraction(:) - sandfraction(:)
	1197
	1198	ELSE ! impveg=T and impsoil=F
[4193]	1199	! Case impsoilt=false and impveg=true
[1079]	1200	IF ( MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp .OR. &
[4871]	1201	MINVAL(sandfraction) .EQ. MAXVAL(sandfraction) .AND. MAXVAL(sandfraction) .EQ. val_exp .OR. &
[1079]	1202	MINVAL(njsc) .EQ. MAXVAL(njsc) .AND. MAXVAL(njsc) .EQ. undef_int ) THEN
	1203
[4808]	1204	CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, &
	1205	clayfraction, sandfraction, siltfraction)
[1079]	1206	njsc(:) = 0
	1207	DO ji = 1, kjpindex
	1208	njsc(ji) = MAXLOC(soilclass(ji,:),1)
	1209	ENDDO
[44]	1210	ENDIF
	1211	ENDIF
[947]	1212
	1213	!Config Key = REINF_SLOPE
	1214	!Config Desc = Slope coef for reinfiltration
	1215	!Config Def = 0.1
	1216	!Config If = IMPOSE_VEG
	1217	!Config Help = Determines the reinfiltration ratio in the grid box due to flat areas
	1218	!Config Units = [-]
	1219	!
	1220	slope_default=0.1
	1221	CALL setvar_p (reinf_slope, val_exp, 'SLOPE', slope_default)
	1222
[566]	1223	!Config Key = SLOWPROC_HEIGHT
[843]	1224	!Config Desc = Height for all vegetation types
[566]	1225	!Config Def = 0., 30., 30., 20., 20., 20., 15., 15., 15., .5, .6, 1.0, 1.0
[843]	1226	!Config If = OK_SECHIBA
[566]	1227	!Config Help = The height used in the 0dim mode. The values should be found
	1228	!Config in the restart file. The new values of height will be computed anyway
	1229	!Config at the end of the current day. The need for this variable is caused
	1230	!Config by the fact that the model may stop during a day and thus we have not
	1231	!Config yet been through the routines which compute the new surface conditions.
[843]	1232	!Config Units = [m]
[8]	1233	CALL setvar_p (height, val_exp, 'SLOWPROC_HEIGHT', height_presc)
	1234
	1235
[5605]	1236	ELSE IF ( .NOT. found_restart .OR. vegetmap_reset ) THEN
[4193]	1237
[2740]	1238	!! 4.1.b Case impveg=false and no restart files: Initialization by reading vegetation map
	1239
	1240	! Initialize veget_max and frac_nobio
[5518]	1241	! Case without restart file
	1242	IF (printlev_loc>=3) WRITE(numout,*) 'Before call slowproc_readvegetmax in initialization phase without restart files'
	1243	IF (printlev_loc>=3) WRITE(numout,*) 'veget_year=', veget_year
[2740]	1244
[5518]	1245	! Call the routine to read the vegetation from file (output is veget_max_new)
	1246	CALL slowproc_readvegetmax(kjpindex, lalo, neighbours, resolution, contfrac, &
	1247	veget_max, veget_max_new, frac_nobio_new, veget_year, .TRUE.)
	1248	IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_readvegetmax in initialization phase'
[2740]	1249
[5518]	1250	! Update vegetation with values read from the file
	1251	veget_max = veget_max_new
	1252	frac_nobio = frac_nobio_new
	1253
	1254	IF (do_wood_harvest) THEN
	1255	! Read the new the wood harvest map from file. Output is wood harvest
	1256	CALL slowproc_woodharvest(kjpindex, lalo, neighbours, resolution, contfrac, woodharvest)
	1257	ENDIF
	1258
	1259
[2740]	1260	!! Reset totaly or partialy veget_max if using DGVM
	1261	IF ( ok_dgvm ) THEN
	1262	! If we are dealing with dynamic vegetation then all natural PFTs should be set to veget_max = 0
	1263	! In case no agriculture is desired, agriculture PFTS should be set to 0 as well
	1264	IF (agriculture) THEN
	1265	DO jv = 2, nvm
	1266	IF (natural(jv)) THEN
	1267	veget_max(:,jv)=zero
	1268	ENDIF
	1269	ENDDO
[8]	1270
[2740]	1271	! Calculate the fraction of crop for each point.
	1272	! Sum only on the indexes corresponding to the non_natural pfts
	1273	frac_crop_tot(:) = zero
	1274	DO jv = 2, nvm
	1275	IF(.NOT. natural(jv)) THEN
[8]	1276	DO ji = 1, kjpindex
[2740]	1277	frac_crop_tot(ji) = frac_crop_tot(ji) + veget_max(ji,jv)
[8]	1278	ENDDO
	1279	ENDIF
[2740]	1280	END DO
[4193]	1281
[2740]	1282	! Calculate the fraction of bare soil
	1283	DO ji = 1, kjpindex
	1284	veget_max(ji,1) = un - frac_crop_tot(ji) - SUM(frac_nobio(ji,:))
	1285	ENDDO
[8]	1286	ELSE
[2740]	1287	veget_max(:,:) = zero
	1288	DO ji = 1, kjpindex
	1289	veget_max(ji,1) = un - SUM(frac_nobio(ji,:))
	1290	ENDDO
	1291	END IF
	1292	END IF ! end ok_dgvm
	1293
[2419]	1294
[2740]	1295	! Call slowproc_veget to correct veget_max and to calculate veget and soiltiles
[4825]	1296	CALL slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
[2740]	1297
	1298	END IF ! end impveg
[2718]	1299
[2740]	1300	!! 4.2 Continue initializing variables not found in restart file. Case for both impveg=true and false.
[2718]	1301
[2740]	1302	! Initialize laimap for the case read_lai if not found in restart file
	1303	IF (read_lai) THEN
	1304	IF ( ALL( laimap(:,:,:) .EQ. val_exp) ) THEN
	1305	! Interpolation of LAI
	1306	CALL slowproc_interlai (kjpindex, lalo, resolution, neighbours, contfrac, laimap)
[8]	1307	ENDIF
[2740]	1308	ENDIF
	1309
	1310	! Initialize lai if not found in restart file and not already initialized using impveg
	1311	IF ( MINVAL(lai) .EQ. MAXVAL(lai) .AND. MAXVAL(lai) .EQ. val_exp) THEN
[8]	1312	IF (read_lai) THEN
[4631]	1313	stempdiag2_bid(1:kjpindex,1:nslm) = stempdiag_bid
[8]	1314	CALL slowproc_lai (kjpindex, lcanop, stempdiag2_bid, &
[4646]	1315	lalo,resolution,lai,laimap)
[2740]	1316	ELSE
[890]	1317	! If we start from scratch, we set lai to zero for consistency with stomate
[2740]	1318	lai(:,:) = zero
[8]	1319	ENDIF
[947]	1320
[2740]	1321	frac_age(:,:,1) = un
	1322	frac_age(:,:,2) = zero
	1323	frac_age(:,:,3) = zero
	1324	frac_age(:,:,4) = zero
	1325	ENDIF
	1326
	1327	! Initialize heigth if not found in restart file and not already initialized using impveg
	1328	IF ( MINVAL(height) .EQ. MAXVAL(height) .AND. MAXVAL(height) .EQ. val_exp) THEN
	1329	! Impose height
	1330	DO jv = 1, nvm
	1331	height(:,jv) = height_presc(jv)
	1332	ENDDO
	1333	ENDIF
	1334
	1335	! Initialize clayfraction and njsc if not found in restart file and not already initialized using impveg
	1336	IF ( MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp .OR. &
[4871]	1337	MINVAL(sandfraction) .EQ. MAXVAL(sandfraction) .AND. MAXVAL(sandfraction) .EQ. val_exp .OR. &
[2740]	1338	MINVAL(njsc) .EQ. MAXVAL(njsc) .AND. MAXVAL(njsc) .EQ. undef_int ) THEN
	1339
	1340	IF (printlev_loc>=4) WRITE (numout,*) 'clayfraction or njcs were not in restart file, call slowproc_soilt'
[4808]	1341	CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, &
	1342	clayfraction, sandfraction, siltfraction)
[2740]	1343	IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_soilt'
	1344	njsc(:) = 0
	1345	DO ji = 1, kjpindex
	1346	njsc(ji) = MAXLOC(soilclass(ji,:),1)
	1347	ENDDO
	1348	ENDIF
[8]	1349
[2740]	1350	!Config Key = GET_SLOPE
	1351	!Config Desc = Read slopes from file and do the interpolation
	1352	!Config Def = n
	1353	!Config If =
	1354	!Config Help = Needed for reading the slopesfile and doing the interpolation. This will be
	1355	! used by the re-infiltration parametrization
	1356	!Config Units = [FLAG]
	1357	get_slope = .FALSE.
	1358	CALL getin_p('GET_SLOPE',get_slope)
	1359
[5454]	1360	IF ( MINVAL(reinf_slope) .EQ. MAXVAL(reinf_slope) .AND. MAXVAL(reinf_slope) .EQ. val_exp .OR. get_slope) THEN
	1361	IF (printlev_loc>=4) WRITE (numout,*) 'reinf_slope was not in restart file. Now call slowproc_slope'
	1362
	1363	CALL slowproc_slope(kjpindex, lalo, neighbours, resolution, contfrac, reinf_slope)
	1364	IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_slope'
	1365
	1366	ENDIF
[2749]	1367
[5454]	1368
[2740]	1369
	1370	!! 5. Some calculations always done, with and without restart files
[947]	1371
[2740]	1372	! The variables veget, veget_max and frac_nobio were all read from restart file or initialized above.
	1373	! Calculate now totfrac_nobio and soiltiles using these variables.
	1374
	1375	! Calculate totfrac_nobio
	1376	totfrac_nobio(:) = zero
	1377	DO jv = 1, nnobio
	1378	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv)
	1379	ENDDO
	1380
	1381	! Calculate soiltile. This variable do not need to be in the restart file.
[3969]	1382	! The sum of all soiltiles makes one, and corresponds to the bio fraction
	1383	! of the grid cell (called vegtot in hydrol)
[2740]	1384	soiltile(:,:) = zero
	1385	DO jv = 1, nvm
	1386	jst = pref_soil_veg(jv)
	1387	DO ji = 1, kjpindex
	1388	soiltile(ji,jst) = soiltile(ji,jst) + veget_max(ji,jv)
	1389	ENDDO
	1390	ENDDO
[3588]	1391	DO ji = 1, kjpindex
	1392	IF (totfrac_nobio(ji) .LT. (1-min_sechiba)) THEN
	1393	soiltile(ji,:)=soiltile(ji,:)/(1-totfrac_nobio(ji))
	1394	ENDIF
	1395	ENDDO
[2740]	1396
	1397	! Always calculate tot_bare_soil
[3969]	1398	! Fraction of bare soil in the mesh (bio+nobio)
[2740]	1399	tot_bare_soil(:) = veget_max(:,1)
	1400	DO jv = 2, nvm
	1401	DO ji =1, kjpindex
	1402	tot_bare_soil(ji) = tot_bare_soil(ji) + (veget_max(ji,jv) - veget(ji,jv))
	1403	ENDDO
	1404	END DO
	1405
[4882]	1406
	1407	!! Calculate fraction of landuse tiles to be used only for diagnostic variables
[4723]	1408	fraclut(:,:)=0
[4882]	1409	nwdFraclut(:,id_psl)=0
	1410	nwdFraclut(:,id_crp)=1.
	1411	nwdFraclut(:,id_urb)=xios_default_val
	1412	nwdFraclut(:,id_pst)=xios_default_val
[4723]	1413	DO jv=1,nvm
	1414	IF (natural(jv)) THEN
	1415	fraclut(:,id_psl) = fraclut(:,id_psl) + veget_max(:,jv)
[4882]	1416	IF(.NOT. is_tree(jv)) THEN
	1417	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl) + veget_max(:,jv)
	1418	ENDIF
[4723]	1419	ELSE
	1420	fraclut(:,id_crp) = fraclut(:,id_crp) + veget_max(:,jv)
	1421	ENDIF
	1422	END DO
[4882]	1423
	1424	WHERE (fraclut(:,id_psl) > min_sechiba)
	1425	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl)/fraclut(:,id_psl)
	1426	ELSEWHERE
	1427	nwdFraclut(:,id_psl) = xios_default_val
	1428	END WHERE
[4723]	1429
[4882]	1430
[2348]	1431	IF (printlev_loc>=3) WRITE (numout,*) ' slowproc_init done '
[2740]	1432
[8]	1433	END SUBROUTINE slowproc_init
[947]	1434
	1435	!! ================================================================================================================================
	1436	!! SUBROUTINE : slowproc_clear
	1437	!!
	1438	!>\BRIEF Clear all variables related to slowproc and stomate modules
	1439	!!
	1440	!_ ================================================================================================================================
	1441
[8]	1442	SUBROUTINE slowproc_clear
	1443
[2718]	1444	! 1 clear all the variables defined as common for the routines in slowproc
[8]	1445
	1446	IF (ALLOCATED (clayfraction)) DEALLOCATE (clayfraction)
[4808]	1447	IF (ALLOCATED (sandfraction)) DEALLOCATE (sandfraction)
	1448	IF (ALLOCATED (siltfraction)) DEALLOCATE (siltfraction)
[8]	1449	IF (ALLOCATED (laimap)) DEALLOCATE (laimap)
[3094]	1450	IF (ALLOCATED (veget_max_new)) DEALLOCATE (veget_max_new)
[4657]	1451	IF (ALLOCATED (woodharvest)) DEALLOCATE (woodharvest)
[3094]	1452	IF (ALLOCATED (frac_nobio_new)) DEALLOCATE (frac_nobio_new)
[947]	1453	IF ( ALLOCATED (soilclass_default)) DEALLOCATE (soilclass_default)
	1454
	1455	! 2. Clear all the variables in stomate
	1456
[8]	1457	CALL stomate_clear
	1458	!
	1459	END SUBROUTINE slowproc_clear
[947]	1460
	1461	!! ================================================================================================================================
	1462	!! SUBROUTINE : slowproc_derivvar
	1463	!!
	1464	!>\BRIEF Initializes variables related to the
	1465	!! parameters to be assimilated, the maximum water on vegetation, the vegetation height,
	1466	!! and the fraction of soil covered by dead leaves and the vegetation height
	1467	!!
	1468	!! DESCRIPTION : (definitions, functional, design, flags):
	1469	!! (1) Initialization of the variables relevant for the assimilation parameters
	1470	!! (2) Intialization of the fraction of soil covered by dead leaves
	1471	!! (3) Initialization of the Vegetation height per PFT
	1472	!! (3) Initialization the maximum water on vegetation for interception with a particular treatement of the PFT no.1
	1473	!!
	1474	!! RECENT CHANGE(S): None
	1475	!!
	1476	!! MAIN OUTPUT VARIABLE(S): ::qsintmax, ::deadleaf_cover, ::assim_param, ::height
	1477	!!
	1478	!! REFERENCE(S) : None
	1479	!!
	1480	!! FLOWCHART : None
	1481	!! \n
	1482	!_ ================================================================================================================================
	1483
[8]	1484	SUBROUTINE slowproc_derivvar (kjpindex, veget, lai, &
[2031]	1485	qsintmax, deadleaf_cover, assim_param, height, temp_growth)
[8]	1486
[947]	1487	!! INTERFACE DESCRIPTION
	1488
	1489	!! 0.1 Input scalar and fields
	1490	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only
	1491	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of pixel covered by PFT. Fraction accounts for none-biological land covers (unitless)
	1492	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! PFT leaf area index (m^{2} m^{-2})
	1493
	1494	!! 0.2. Output scalar and fields
	1495	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintmax !! Maximum water on vegetation for interception(mm)
	1496	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: deadleaf_cover !! fraction of soil covered by dead leaves (unitless)
	1497	REAL(r_std), DIMENSION (kjpindex,nvm,npco2), INTENT (out) :: assim_param !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
	1498	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! height of the vegetation or surface in general ??? (m)
[5364]	1499	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_growth !! growth temperature (Â°C)
[8]	1500	!
[947]	1501	!! 0.3 Local declaration
[3831]	1502	INTEGER(i_std) :: jv !! Local indices
[947]	1503	!_ ================================================================================================================================
[8]	1504
	1505	!
[947]	1506	! 1. Initialize (why here ??) the variables revelant for the assimilation parameters
[8]	1507	!
	1508	DO jv = 1, nvm
	1509	assim_param(:,jv,ivcmax) = vcmax_fix(jv)
	1510	ENDDO
	1511
	1512	!
[947]	1513	! 2. Intialize the fraction of soil covered by dead leaves
[8]	1514	!
	1515	deadleaf_cover(:) = zero
	1516
	1517	!
[947]	1518	! 3. Initialize the Vegetation height per PFT
[8]	1519	!
	1520	DO jv = 1, nvm
	1521	height(:,jv) = height_presc(jv)
	1522	ENDDO
	1523	!
[947]	1524	! 4. Initialize the maximum water on vegetation for interception
[8]	1525	!
	1526	qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:)
[947]	1527
	1528	! Added by Nathalie - July 2006
	1529	! Initialize the case of the PFT no.1 to zero
[8]	1530	qsintmax(:,1) = zero
	1531
[2031]	1532	temp_growth(:)=25.
	1533
[8]	1534	END SUBROUTINE slowproc_derivvar
	1535
	1536
[947]	1537	!! ================================================================================================================================
	1538	!! SUBROUTINE : slowproc_mean
	1539	!!
	1540	!>\BRIEF Accumulates field_in over a period of dt_tot.
	1541	!! Has to be called at every time step (dt).
	1542	!! Mean value is calculated if ldmean=.TRUE.
	1543	!! field_mean must be initialized outside of this routine!
	1544	!!
	1545	!! DESCRIPTION : (definitions, functional, design, flags):
	1546	!! (1) AcumAcuumlm
	1547	!!
	1548	!! RECENT CHANGE(S): None
	1549	!!
	1550	!! MAIN OUTPUT VARIABLE(S): ::field_main
	1551	!!
	1552	!! REFERENCE(S) : None
	1553	!!
	1554	!! FLOWCHART : None
	1555	!! \n
	1556	!_ ================================================================================================================================
[8]	1557
[947]	1558	SUBROUTINE slowproc_mean (npts, n_dim2, dt_tot, dt, ldmean, field_in, field_mean)
[8]	1559
[947]	1560	!
	1561	!! 0 declarations
[8]	1562
[947]	1563	!! 0.1 input scalar and variables
	1564	INTEGER(i_std), INTENT(in) :: npts !! Domain size- terrestrial pixels only
	1565	INTEGER(i_std), INTENT(in) :: n_dim2 !! Number of PFTs
	1566	REAL(r_std), INTENT(in) :: dt_tot !! Time step of stomate (in days). The period over which the accumulation or the mean is computed
	1567	REAL(r_std), INTENT(in) :: dt !! Time step in days
	1568	LOGICAL, INTENT(in) :: ldmean !! Flag to calculate the mean after the accumulation ???
	1569	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_in !! Daily field
[8]	1570
[947]	1571	!! 0.3 Modified field; The computed sum or mean field over dt_tot time period depending on the flag ldmean
	1572	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_mean !! Accumulated field at dt_tot time period or mean field over dt_tot
	1573
[8]	1574
[1082]	1575	!_ ================================================================================================================================
[8]	1576
	1577	!
[947]	1578	! 1. Accumulation the field over dt_tot period
[8]	1579	!
	1580	field_mean(:,:) = field_mean(:,:) + field_in(:,:) * dt
	1581
	1582	!
[947]	1583	! 2. If the flag ldmean set, the mean field is computed over dt_tot period
[8]	1584	!
	1585	IF (ldmean) THEN
	1586	field_mean(:,:) = field_mean(:,:) / dt_tot
	1587	ENDIF
	1588
	1589	END SUBROUTINE slowproc_mean
[947]	1590
	1591
	1592
	1593	!! ================================================================================================================================
	1594	!! SUBROUTINE : slowproc_long
	1595	!!
	1596	!>\BRIEF Calculates a temporally smoothed field (field_long) from
	1597	!! instantaneous input fields.Time constant tau determines the strength of the smoothing.
	1598	!! For tau -> infinity??, field_long becomes the true mean value of field_inst
	1599	!! (but the spinup becomes infinietly long, too).
	1600	!! field_long must be initialized outside of this routine!
	1601	!!
	1602	!! DESCRIPTION : (definitions, functional, design, flags):
	1603	!! (1) Testing the time coherence betwen the time step dt and the time tau over which
	1604	!! the rescaled of the mean is performed
	1605	!! (2) Computing the rescaled mean over tau period
	1606	!! MAIN OUTPUT VARIABLE(S): field_long
	1607	!!
	1608	!! RECENT CHANGE(S): None
	1609	!!
	1610	!! MAIN OUTPUT VARIABLE(S): ::field_long
	1611	!!
	1612	!! REFERENCE(S) : None
	1613	!!
	1614	!! FLOWCHART : None
	1615	!! \n
	1616	!_ ================================================================================================================================
	1617
[8]	1618	SUBROUTINE slowproc_long (npts, n_dim2, dt, tau, field_inst, field_long)
	1619
	1620	!
	1621	! 0 declarations
	1622	!
	1623
[947]	1624	! 0.1 input scalar and fields
[8]	1625
[947]	1626	INTEGER(i_std), INTENT(in) :: npts !! Domain size- terrestrial pixels only
	1627	INTEGER(i_std), INTENT(in) :: n_dim2 !! Second dimension of the fields, which represents the number of PFTs
	1628	REAL(r_std), INTENT(in) :: dt !! Time step in days
	1629	REAL(r_std), INTENT(in) :: tau !! Integration time constant (has to have same unit as dt!)
	1630	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_inst !! Instantaneous field
[8]	1631
[947]	1632
[8]	1633	! 0.2 modified field
	1634
	1635	! Long-term field
[947]	1636	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_long !! Mean value of the instantaneous field rescaled at tau time period
[8]	1637
[1082]	1638	!_ ================================================================================================================================
[8]	1639
	1640	!
[947]	1641	! 1 test coherence of the time
[8]	1642
[42]	1643	IF ( ( tau .LT. dt ) .OR. ( dt .LE. zero ) .OR. ( tau .LE. zero ) ) THEN
[8]	1644	WRITE(numout,*) 'slowproc_long: Problem with time steps'
	1645	WRITE(numout,*) 'dt=',dt
	1646	WRITE(numout,*) 'tau=',tau
	1647	ENDIF
	1648
	1649	!
[947]	1650	! 2 integration of the field over tau
[8]	1651
	1652	field_long(:,:) = ( field_inst(:,:)dt + field_long(:,:)(tau-dt) ) / tau
	1653
	1654	END SUBROUTINE slowproc_long
[947]	1655
	1656
	1657	!! ================================================================================================================================
[6145]	1658	!! SUBROUTINE : slowproc_veget_max_limit
	1659	!!
	1660	!>\BRIEF Set small fractions of veget_max to zero and normalize to keep the sum equal 1
	1661	!!
	1662	!! DESCRIPTION : Set small fractions of veget_max to zero and normalize to keep the sum equal 1
	1663	!!
	1664	!! RECENT CHANGE(S): The subroutine was previously a part of slowproc_veget,
	1665	!! but was separated to be called also from slowproc_readvegetmax in order
	1666	!! to have limited/normalized vegetation fractions right after its reading
	1667	!! from the file (added by V.Bastrikov, 15/06/2019)
	1668	!!
	1669	!! MAIN OUTPUT VARIABLE(S): :: frac_nobio, veget_max
	1670	!!
	1671	!! REFERENCE(S) : None
	1672	!!
	1673	!! FLOWCHART : None
	1674	!! \n
	1675	!_ ================================================================================================================================
	1676
	1677	SUBROUTINE slowproc_veget_max_limit (kjpindex, frac_nobio, veget_max)
	1678	!
	1679	! 0. Declarations
	1680	!
	1681	! 0.1 Input variables
	1682	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
	1683
	1684	! 0.2 Modified variables
	1685	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio !! Fraction of the mesh which is covered by ice, lakes, ...
	1686	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
	1687
	1688	! 0.4 Local scalar and varaiables
	1689	INTEGER(i_std) :: ji, jv !! indices
	1690	REAL(r_std) :: SUMveg !! Total vegetation summed across PFTs
	1691
	1692	!_ ================================================================================================================================
	1693	IF (printlev_loc >= 3) WRITE(numout,*) 'Entering slowproc_veget_max_limit'
	1694
	1695	!! Set to zero fractions of frac_nobio and veget_max smaller than min_vegfrac
	1696	DO ji = 1, kjpindex
	1697	IF ( SUM(frac_nobio(ji,:)) .LT. min_vegfrac ) THEN
	1698	frac_nobio(ji,:) = zero
	1699	ENDIF
	1700
	1701	IF (.NOT. ok_dgvm) THEN
	1702	DO jv = 1, nvm
	1703	IF ( veget_max(ji,jv) .LT. min_vegfrac ) THEN
	1704	veget_max(ji,jv) = zero
	1705	ENDIF
	1706	ENDDO
	1707	END IF
	1708
	1709	!! Normalize to keep the sum equal 1.
	1710	SUMveg = SUM(frac_nobio(ji,:))+SUM(veget_max(ji,:))
	1711	frac_nobio(ji,:) = frac_nobio(ji,:)/SUMveg
	1712	veget_max(ji,:) = veget_max(ji,:)/SUMveg
	1713	ENDDO
	1714
	1715	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_veget_max_limit ended'
	1716
	1717	END SUBROUTINE slowproc_veget_max_limit
	1718
	1719
	1720	!! ================================================================================================================================
[947]	1721	!! SUBROUTINE : slowproc_veget
	1722	!!
[2718]	1723	!>\BRIEF Set small fractions to zero and normalize to keep the sum equal 1. Calucate veget and soiltile.
[947]	1724	!!
[2718]	1725	!! DESCRIPTION : Set small fractions to zero and normalize to keep the sum equal 1. Calucate veget and soiltile.
	1726	!! (1) Set veget_max and frac_nobio for fraction smaller than min_vegfrac.
[4723]	1727	!! (2) Calculate veget
	1728	!! (3) Calculate totfrac_nobio
	1729	!! (4) Calculate soiltile
	1730	!! (5) Calculate fraclut
[947]	1731	!!
	1732	!! RECENT CHANGE(S): None
	1733	!!
[4723]	1734	!! MAIN OUTPUT VARIABLE(S): :: frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut
[947]	1735	!!
	1736	!! REFERENCE(S) : None
	1737	!!
	1738	!! FLOWCHART : None
	1739	!! \n
	1740	!_ ================================================================================================================================
	1741
[4825]	1742	SUBROUTINE slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
[8]	1743	!
	1744	! 0. Declarations
	1745	!
[947]	1746	! 0.1 Input variables
	1747	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
	1748	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! PFT leaf area index (m^{2} m^{-2})
[8]	1749
[947]	1750	! 0.2 Modified variables
	1751	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio !! Fraction of the mesh which is covered by ice, lakes, ...
	1752	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
[8]	1753
[947]	1754	! 0.3 Output variables
	1755	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: veget !! Fraction of pixel covered by PFT. Fraction accounts for none-biological land covers (unitless)
[2718]	1756	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio
[3969]	1757	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
[4723]	1758	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
[4825]	1759	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: nwdFraclut !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
[8]	1760
[947]	1761	! 0.4 Local scalar and varaiables
	1762	INTEGER(i_std) :: ji, jv, jst !! indices
	1763
	1764	!_ ================================================================================================================================
[2718]	1765	IF (printlev_loc > 8) WRITE(numout,*) 'Entering slowproc_veget'
[8]	1766
[2718]	1767	!! 1. Set to zero fractions of frac_nobio and veget_max smaller than min_vegfrac
	1768	!! Normalize to have the sum equal 1.
[6145]	1769	CALL slowproc_veget_max_limit(kjpindex, frac_nobio, veget_max)
[1118]	1770
[3094]	1771	!! 2. Calculate veget
[2718]	1772	!! If lai of a vegetation type (jv > 1) is small, increase soil part
	1773	!! stomate-like calculation
[8]	1774	DO ji = 1, kjpindex
[2718]	1775	veget(ji,1)=veget_max(ji,1)
[8]	1776	DO jv = 2, nvm
[3524]	1777	veget(ji,jv) = veget_max(ji,jv) * ( un - exp( - lai(ji,jv) * ext_coeff_vegetfrac(jv) ) )
[8]	1778	ENDDO
	1779	ENDDO
[947]	1780
[8]	1781
[3094]	1782	!! 3. Calculate totfrac_nobio
[947]	1783	totfrac_nobio(:) = zero
	1784	DO jv = 1, nnobio
	1785	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv)
	1786	ENDDO
	1787
[2718]	1788
[3094]	1789	!! 4. Calculate soiltiles
[2718]	1790	!! Soiltiles are only used in hydrol, but we fix them in here because some time it might depend
	1791	!! on a changing vegetation (but then some adaptation should be made to hydrol) and be also used
	1792	!! in the other modules to perform separated energy balances
[3969]	1793	! The sum of all soiltiles makes one, and corresponds to the bio fraction
	1794	! of the grid cell (called vegtot in hydrol)
[947]	1795	soiltile(:,:) = zero
	1796	DO jv = 1, nvm
	1797	jst = pref_soil_veg(jv)
	1798	DO ji = 1, kjpindex
	1799	soiltile(ji,jst) = soiltile(ji,jst) + veget_max(ji,jv)
	1800	ENDDO
	1801	ENDDO
[3588]	1802	DO ji = 1, kjpindex
	1803	IF (totfrac_nobio(ji) .LT. (1-min_sechiba)) THEN
	1804	soiltile(ji,:)=soiltile(ji,:)/(1.-totfrac_nobio(ji))
	1805	ENDIF
[3969]	1806	ENDDO
[3588]	1807
[4723]	1808	!! 5. Calculate fraction of landuse tiles to be used only for diagnostic variables
	1809	fraclut(:,:)=0
[4825]	1810	nwdFraclut(:,id_psl)=0
	1811	nwdFraclut(:,id_crp)=1.
[4872]	1812	nwdFraclut(:,id_urb)=xios_default_val
	1813	nwdFraclut(:,id_pst)=xios_default_val
[4723]	1814	DO jv=1,nvm
	1815	IF (natural(jv)) THEN
	1816	fraclut(:,id_psl) = fraclut(:,id_psl) + veget_max(:,jv)
[4825]	1817	IF(.NOT. is_tree(jv)) THEN
	1818	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl) + veget_max(:,jv)
	1819	ENDIF
[4723]	1820	ELSE
	1821	fraclut(:,id_crp) = fraclut(:,id_crp) + veget_max(:,jv)
	1822	ENDIF
	1823	END DO
[4857]	1824
	1825	WHERE (fraclut(:,id_psl) > min_sechiba)
	1826	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl)/fraclut(:,id_psl)
	1827	ELSEWHERE
[4872]	1828	nwdFraclut(:,id_psl) = xios_default_val
[4857]	1829	END WHERE
[4723]	1830
[8]	1831	END SUBROUTINE slowproc_veget
[947]	1832
	1833
	1834	!! ================================================================================================================================
	1835	!! SUBROUTINE : slowproc_lai
	1836	!!
	1837	!>\BRIEF Do the interpolation of lai for the PFTs in case the laimap is not read
	1838	!!
	1839	!! DESCRIPTION : (definitions, functional, design, flags):
	1840	!! (1) Interplation by using the mean value of laimin and laimax for the PFTs
	1841	!! (2) Interpolation between laimax and laimin values by using the temporal
	1842	!! variations
	1843	!! (3) If problem occurs during the interpolation, the routine stops
	1844	!!
	1845	!! RECENT CHANGE(S): None
	1846	!!
	1847	!! MAIN OUTPUT VARIABLE(S): ::lai
	1848	!!
	1849	!! REFERENCE(S) : None
	1850	!!
	1851	!! FLOWCHART : None
	1852	!! \n
	1853	!_ ================================================================================================================================
	1854
[4646]	1855	SUBROUTINE slowproc_lai (kjpindex,lcanop,stempdiag,lalo,resolution,lai,laimap)
[8]	1856	!
	1857	! 0. Declarations
	1858	!
[947]	1859	!! 0.1 Input variables
	1860	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
	1861	INTEGER(i_std), INTENT(in) :: lcanop !! soil level used for LAI
[4631]	1862	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: stempdiag !! Soil temperature (K) ???
[8]	1863	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
[947]	1864	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size in x an y of the grid (m) - surface area of the gridbox
	1865	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: laimap !! map of lai read
[8]	1866
[947]	1867	!! 0.2 Output
	1868	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: lai !! PFT leaf area index (m^{2} m^{-2})LAI
[8]	1869
[947]	1870	!! 0.4 Local
	1871	INTEGER(i_std) :: ji,jv !! Local indices
	1872	!_ ================================================================================================================================
[8]	1873
[947]	1874	!
	1875	IF ( .NOT. read_lai ) THEN
[8]	1876
[947]	1877	lai(: ,1) = zero
	1878	! On boucle sur 2,nvm au lieu de 1,nvm
	1879	DO jv = 2,nvm
	1880	SELECT CASE (type_of_lai(jv))
	1881
	1882	CASE ("mean ")
	1883	!
	1884	! 1. do the interpolation between laimax and laimin
	1885	!
[8]	1886	lai(:,jv) = undemi * (llaimax(jv) + llaimin(jv))
[947]	1887	!
	1888	CASE ("inter")
	1889	!
	1890	! 2. do the interpolation between laimax and laimin
	1891	!
[8]	1892	DO ji = 1,kjpindex
[947]	1893	lai(ji,jv) = llaimin(jv) + tempfunc(stempdiag(ji,lcanop)) * (llaimax(jv) - llaimin(jv))
[8]	1894	ENDDO
[947]	1895	!
	1896	CASE default
	1897	!
	1898	! 3. Problem
	1899	!
	1900	WRITE (numout,*) 'This kind of lai choice is not possible. '// &
	1901	' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
[2613]	1902	CALL ipslerr_p(3,'slowproc_lai','Bad value for type_of_lai','read_lai=false','')
[947]	1903	END SELECT
	1904
	1905	ENDDO
	1906	!
	1907	ELSE
	1908	lai(: ,1) = zero
	1909	! On boucle sur 2,nvm au lieu de 1,nvm
	1910	DO jv = 2,nvm
	1911
	1912	SELECT CASE (type_of_lai(jv))
	1913
	1914	CASE ("mean ")
	1915	!
	1916	! 1. force MAXVAL of laimap on lai on this PFT
	1917	!
[8]	1918	DO ji = 1,kjpindex
[947]	1919	lai(ji,jv) = MAXVAL(laimap(ji,jv,:))
[8]	1920	ENDDO
[947]	1921	!
	1922	CASE ("inter")
	1923	!
	1924	! 2. do the interpolation between laimax and laimin
	1925	!
	1926	!
	1927	! If January
	1928	!
[4646]	1929	IF (month_end .EQ. 1 ) THEN
	1930	IF (day_end .LE. 15) THEN
	1931	lai(:,jv) = laimap(:,jv,12)(1-(day_end+15)/30.) + laimap(:,jv,1)((day_end+15)/30.)
[8]	1932	ELSE
[4646]	1933	lai(:,jv) = laimap(:,jv,1)(1-(day_end-15)/30.) + laimap(:,jv,2)((day_end-15)/30.)
[8]	1934	ENDIF
[947]	1935	!
	1936	! If December
	1937	!
[4646]	1938	ELSE IF (month_end .EQ. 12) THEN
	1939	IF (day_end .LE. 15) THEN
	1940	lai(:,jv) = laimap(:,jv,11)(1-(day_end+15)/30.) + laimap(:,jv,12)((day_end+15)/30.)
[8]	1941	ELSE
[4646]	1942	lai(:,jv) = laimap(:,jv,12)(1-(day_end-15)/30.) + laimap(:,jv,1)((day_end-15)/30.)
[8]	1943	ENDIF
[947]	1944	!
	1945	! ELSE
	1946	!
[8]	1947	ELSE
[4646]	1948	IF (day_end .LE. 15) THEN
	1949	lai(:,jv) = laimap(:,jv,month_end-1)(1-(day_end+15)/30.) + laimap(:,jv,month_end)((day_end+15)/30.)
[8]	1950	ELSE
[4646]	1951	lai(:,jv) = laimap(:,jv,month_end)(1-(day_end-15)/30.) + laimap(:,jv,month_end+1)((day_end-15)/30.)
[8]	1952	ENDIF
	1953	ENDIF
[947]	1954	!
	1955	CASE default
	1956	!
	1957	! 3. Problem
	1958	!
	1959	WRITE (numout,*) 'This kind of lai choice is not possible. '// &
	1960	' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
[2613]	1961	CALL ipslerr_p(3,'slowproc_lai','Bad value for type_of_lai','read_lai=true','')
[947]	1962	END SELECT
	1963
	1964	ENDDO
	1965	ENDIF
[8]	1966
[947]	1967	END SUBROUTINE slowproc_lai
[8]	1968
[947]	1969	!! ================================================================================================================================
[2532]	1970	!! SUBROUTINE : slowproc_interlai
[947]	1971	!!
	1972	!>\BRIEF Interpolate the LAI map to the grid of the model
	1973	!!
	1974	!! DESCRIPTION : (definitions, functional, design, flags):
	1975	!!
	1976	!! RECENT CHANGE(S): None
	1977	!!
	1978	!! MAIN OUTPUT VARIABLE(S): ::laimap
	1979	!!
	1980	!! REFERENCE(S) : None
	1981	!!
	1982	!! FLOWCHART : None
	1983	!! \n
	1984	!_ ================================================================================================================================
	1985
[3831]	1986	SUBROUTINE slowproc_interlai(nbpt, lalo, resolution, neighbours, contfrac, laimap)
	1987
	1988	USE interpweight
	1989
	1990	IMPLICIT NONE
	1991
[8]	1992	!
	1993	!
	1994	!
	1995	! 0.1 INPUT
	1996	!
[947]	1997	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
	1998	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes
	1999	!! (beware of the order = 1 : latitude, 2 : longitude)
	2000	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
[3831]	2001	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
	2002	!! (1=North and then clockwise)
[947]	2003	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
[8]	2004	!
	2005	! 0.2 OUTPUT
	2006	!
[947]	2007	REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) !! lai read variable and re-dimensioned
[8]	2008	!
	2009	! 0.3 LOCAL
	2010	!
[947]	2011	CHARACTER(LEN=80) :: filename !! name of the LAI map read
[3831]	2012	INTEGER(i_std) :: ib, ip, jp, it, jv
	2013	REAL(r_std) :: lmax, lmin, ldelta
	2014	LOGICAL :: renormelize_lai ! flag to force LAI renormelization
	2015	INTEGER :: ier
[947]	2016
[3831]	2017	REAL(r_std), DIMENSION(nbpt) :: alaimap !! availability of the lai interpolation
	2018	INTEGER, DIMENSION(4) :: invardims
[5364]	2019	REAL(r_std), DIMENSION(nbpt,nvm,12) :: lairefrac !! lai fractions re-dimensioned
	2020	REAL(r_std), DIMENSION(nbpt,nvm,12) :: fraclaiinterp !! lai fractions re-dimensioned
[3831]	2021	REAL(r_std), DIMENSION(:), ALLOCATABLE :: vmin, vmax !! min/max values to use for the
	2022	!! renormalization
	2023	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
	2024	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
	2025	REAL(r_std), DIMENSION(nvm) :: variabletypevals !! Values for all the types of the variable
	2026	!! (variabletypevals(1) = -un, not used)
	2027	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
	2028	!! 'XYKindTime': Input values are kinds
	2029	!! of something with a temporal
	2030	!! evolution on the dx*dy matrix'
	2031	LOGICAL :: nonegative !! whether negative values should be removed
	2032	CHARACTER(LEN=50) :: maskingtype !! Type of masking
	2033	!! 'nomask': no-mask is applied
	2034	!! 'mbelow': take values below maskvals(1)
	2035	!! 'mabove': take values above maskvals(1)
	2036	!! 'msumrange': take values within 2 ranges;
	2037	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
	2038	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
	2039	!! (normalized by maskvals(3))
	2040	!! 'var': mask values are taken from a
	2041	!! variable inside the file (>0)
	2042	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
	2043	!! `maskingtype')
	2044	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
[947]	2045	!_ ================================================================================================================================
[3831]	2046
[8]	2047	!
[566]	2048	!Config Key = LAI_FILE
	2049	!Config Desc = Name of file from which the vegetation map is to be read
	2050	!Config If = LAI_MAP
[620]	2051	!Config Def = lai2D.nc
[566]	2052	!Config Help = The name of the file to be opened to read the LAI
	2053	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
	2054	!Config map which is derived from a Nicolas VIOVY one.
[843]	2055	!Config Units = [FILE]
[8]	2056	!
	2057	filename = 'lai2D.nc'
	2058	CALL getin_p('LAI_FILE',filename)
[3831]	2059	variablename = 'LAI'
[2613]	2060
[5364]	2061	IF (xios_interpolation) THEN
	2062	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_interlai: Use XIOS to read and interpolate " &
	2063	// TRIM(filename) //" for variable " //TRIM(variablename)
	2064
	2065	CALL xios_orchidee_recv_field('lai_interp',lairefrac)
	2066	CALL xios_orchidee_recv_field('frac_lai_interp',fraclaiinterp)
	2067	alaimap(:) = fraclaiinterp(:,1,1)
	2068	ELSE
[2613]	2069
[5364]	2070	IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_interlai: Start interpolate " &
	2071	// TRIM(filename) //" for variable " //TRIM(variablename)
[2613]	2072
[5364]	2073	! invardims: shape of variable in input file to interpolate
	2074	invardims = interpweight_get_var4dims_file(filename, variablename)
	2075	! Check coherence of dimensions read from the file
	2076	IF (invardims(4) /= 12) CALL ipslerr_p(3,'slowproc_interlai','Wrong dimension of time dimension in input file for lai','','')
	2077	IF (invardims(3) /= nvm) CALL ipslerr_p(3,'slowproc_interlai','Wrong dimension of PFT dimension in input file for lai','','')
[2613]	2078
[5364]	2079	ALLOCATE(vmin(nvm),stat=ier)
	2080	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_interlai','Problem in allocation of variable vmin','','')
[2613]	2081
[5364]	2082	ALLOCATE(vmax(nvm), STAT=ier)
	2083	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_interlai','Problem in allocation of variable vmax','','')
[2613]	2084
[5364]	2085
[3831]	2086	! Assigning values to vmin, vmax
[5364]	2087	vmin = un
	2088	vmax = nvm*un
[2613]	2089
[5364]	2090	variabletypevals = -un
[2613]	2091
[5364]	2092	!! Variables for interpweight
	2093	! Type of calculation of cell fractions
	2094	fractype = 'default'
	2095	! Name of the longitude and latitude in the input file
	2096	lonname = 'longitude'
	2097	latname = 'latitude'
	2098	! Should negative values be set to zero from input file?
	2099	nonegative = .TRUE.
	2100	! Type of mask to apply to the input data (see header for more details)
	2101	maskingtype = 'mbelow'
	2102	! Values to use for the masking
	2103	maskvals = (/ 20., undef_sechiba, undef_sechiba /)
	2104	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
	2105	namemaskvar = ''
[2613]	2106
[5364]	2107	CALL interpweight_4D(nbpt, nvm, variabletypevals, lalo, resolution, neighbours, &
	2108	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
	2109	maskvals, namemaskvar, nvm, invardims(4), -1, fractype, &
	2110	-1., -1., lairefrac, alaimap)
[3831]	2111
[5364]	2112	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_interlai after interpweight_4D'
[3831]	2113
[5364]	2114	ENDIF
	2115
	2116
	2117
[8]	2118	!
	2119	!
[3831]	2120	!Config Key = RENORM_LAI
	2121	!Config Desc = flag to force LAI renormelization
	2122	!Config If = LAI_MAP
	2123	!Config Def = n
	2124	!Config Help = If true, the laimap will be renormalize between llaimin and llaimax parameters.
	2125	!Config Units = [FLAG]
[8]	2126	!
[3831]	2127	renormelize_lai = .FALSE.
	2128	CALL getin_p('RENORM_LAI',renormelize_lai)
[2613]	2129
[8]	2130	!
	2131	laimap(:,:,:) = zero
	2132	!
[3831]	2133	IF (printlev_loc >= 5) THEN
	2134	WRITE(numout,*)' slowproc_interlai before starting loop nbpt:', nbpt
	2135	END IF
	2136
[3858]	2137	! Assigning the right values and giving a value where information was not found
[8]	2138	DO ib=1,nbpt
[5364]	2139	IF (alaimap(ib) < min_sechiba) THEN
[3858]	2140	DO jv=1,nvm
	2141	laimap(ib,jv,:) = (llaimax(jv)+llaimin(jv))/deux
[3831]	2142	ENDDO
[3858]	2143	ELSE
	2144	DO jv=1, nvm
[5364]	2145	DO it=1, 12
[3858]	2146	laimap(ib,jv,it) = lairefrac(ib,jv,it)
	2147	ENDDO
	2148	ENDDO
	2149	END IF
[8]	2150	ENDDO
	2151	!
[947]	2152	! Normelize the read LAI by the values SECHIBA is used to
	2153	!
	2154	IF ( renormelize_lai ) THEN
	2155	DO ib=1,nbpt
[3831]	2156	DO jv=1, nvm
[947]	2157	lmax = MAXVAL(laimap(ib,jv,:))
	2158	lmin = MINVAL(laimap(ib,jv,:))
	2159	ldelta = lmax-lmin
	2160	IF ( ldelta < min_sechiba) THEN
	2161	! LAI constante ... keep it constant
	2162	laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)+(llaimax(jv)+llaimin(jv))/deux
	2163	ELSE
	2164	laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)/(lmax-lmin)*(llaimax(jv)-llaimin(jv))+llaimin(jv)
	2165	ENDIF
	2166	ENDDO
	2167	ENDDO
	2168	ENDIF
[2581]	2169
[3831]	2170	! Write diagnostics
	2171	CALL xios_orchidee_send_field("alaimap",alaimap)
[5364]	2172	CALL xios_orchidee_send_field("interp_diag_lai",laimap)
	2173
[3831]	2174	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_interlai ended'
	2175
[2532]	2176	END SUBROUTINE slowproc_interlai
[947]	2177
	2178	!! ================================================================================================================================
[3094]	2179	!! SUBROUTINE : slowproc_readvegetmax
[947]	2180	!!
[3094]	2181	!>\BRIEF Read and interpolate a vegetation map (by pft)
[947]	2182	!!
	2183	!! DESCRIPTION : (definitions, functional, design, flags):
	2184	!!
[3094]	2185	!! RECENT CHANGE(S): The subroutine was previously called slowproc_update.
[947]	2186	!!
	2187	!! MAIN OUTPUT VARIABLE(S):
	2188	!!
	2189	!! REFERENCE(S) : None
	2190	!!
	2191	!! FLOWCHART : None
	2192	!! \n
	2193	!_ ================================================================================================================================
	2194
[3831]	2195	SUBROUTINE slowproc_readvegetmax(nbpt, lalo, neighbours, resolution, contfrac, veget_last, &
[2613]	2196	veget_next, frac_nobio_next, veget_year, init)
[3831]	2197
	2198	USE interpweight
	2199	IMPLICIT NONE
	2200
[8]	2201	!
	2202	!
	2203	!
	2204	! 0.1 INPUT
	2205	!
[947]	2206	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs
	2207	!! to be interpolated
	2208	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
[3447]	2209	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
	2210	!! (1=North and then clockwise)
[947]	2211	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
	2212	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of continent in the grid
[8]	2213	!
[947]	2214	REAL(r_std), DIMENSION(nbpt,nvm), INTENT(in) :: veget_last !! old max vegetfrac
	2215	INTEGER(i_std), INTENT(in) :: veget_year !! first year for landuse (0 == NO TIME AXIS)
[2613]	2216	LOGICAL, INTENT(in) :: init !! initialisation : in case of dgvm, it forces update of all PFTs
[8]	2217	!
	2218	! 0.2 OUTPUT
	2219	!
[947]	2220	REAL(r_std), DIMENSION(nbpt,nvm), INTENT(out) :: veget_next !! new max vegetfrac
	2221	REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(out) :: frac_nobio_next !! new fraction of the mesh which is
[2613]	2222	!! covered by ice, lakes, ...
[3831]	2223
[8]	2224	!
	2225	! 0.3 LOCAL
	2226	!
	2227	!
	2228	CHARACTER(LEN=80) :: filename
[3831]	2229	INTEGER(i_std) :: ib, inobio, jv
[8]	2230	REAL(r_std) :: sumf, err, norm
	2231	!
	2232	! for DGVM case :
	2233	REAL(r_std) :: sum_veg ! sum of vegets
	2234	REAL(r_std) :: sum_nobio ! sum of nobios
	2235	REAL(r_std) :: sumvAnthro_old, sumvAnthro ! last an new sum of antrhopic vegets
	2236	REAL(r_std) :: rapport ! (S-B) / (S-A)
[2613]	2237	LOGICAL :: partial_update ! if TRUE, partialy update PFT (only anthropic ones)
	2238	! e.g. in case of DGVM and not init (optional parameter)
[3831]	2239	REAL(r_std), DIMENSION(nbpt,nvm) :: vegetrefrac !! veget fractions re-dimensioned
	2240	REAL(r_std), DIMENSION(nbpt) :: aveget !! Availability of the soilcol interpolation
[5364]	2241	REAL(r_std), DIMENSION(nbpt,nvm) :: aveget_nvm !! Availability of the soilcol interpolation
[3831]	2242	REAL(r_std), DIMENSION(nvm) :: vmin, vmax !! min/max values to use for the renormalization
	2243	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
	2244	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
	2245	REAL(r_std), DIMENSION(nvm) :: variabletypevals !! Values for all the types of the variable
	2246	!! (variabletypevals(1) = -un, not used)
	2247	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
	2248	!! 'XYKindTime': Input values are kinds
	2249	!! of something with a temporal
	2250	!! evolution on the dx*dy matrix'
	2251	LOGICAL :: nonegative !! whether negative values should be removed
	2252	CHARACTER(LEN=50) :: maskingtype !! Type of masking
	2253	!! 'nomask': no-mask is applied
	2254	!! 'mbelow': take values below maskvals(1)
	2255	!! 'mabove': take values above maskvals(1)
	2256	!! 'msumrange': take values within 2 ranges;
	2257	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
	2258	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
	2259	!! (normalized by maskvals(3))
	2260	!! 'var': mask values are taken from a
	2261	!! variable inside the file (>0)
	2262	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
	2263	!! `maskingtype')
	2264	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
	2265	CHARACTER(LEN=250) :: msg
[1082]	2266
	2267	!_ ================================================================================================================================
[3831]	2268
	2269	IF (printlev_loc >= 5) PRINT *,' In slowproc_readvegetmax'
	2270
[8]	2271	!
[566]	2272	!Config Key = VEGETATION_FILE
	2273	!Config Desc = Name of file from which the vegetation map is to be read
[5518]	2274	!Config If =
[620]	2275	!Config Def = PFTmap.nc
[566]	2276	!Config Help = The name of the file to be opened to read a vegetation
	2277	!Config map (in pft) is to be given here.
[843]	2278	!Config Units = [FILE]
[8]	2279	!
[620]	2280	filename = 'PFTmap.nc'
[8]	2281	CALL getin_p('VEGETATION_FILE',filename)
[3831]	2282	variablename = 'maxvegetfrac'
[1850]	2283
[8]	2284
[5364]	2285	IF (xios_interpolation) THEN
	2286	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readvegetmax: Use XIOS to read and interpolate " &
	2287	// TRIM(filename) // " for variable " // TRIM(variablename)
[2613]	2288
[5364]	2289	CALL xios_orchidee_recv_field('frac_veget',vegetrefrac)
	2290	CALL xios_orchidee_recv_field('frac_veget_frac',aveget_nvm)
	2291	aveget(:)=aveget_nvm(:,1)
	2292
	2293	DO ib = 1, nbpt
	2294	IF (aveget(ib) > min_sechiba) THEN
	2295	vegetrefrac(ib,:) = vegetrefrac(ib,:)/aveget(ib) ! intersected area normalization
	2296	vegetrefrac(ib,:) = vegetrefrac(ib,:)/SUM(vegetrefrac(ib,:))
	2297	ENDIF
	2298	ENDDO
	2299
	2300	ELSE
[8]	2301
[5364]	2302	IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_readvegetmax: Start interpolate " &
	2303	// TRIM(filename) // " for variable " // TRIM(variablename)
[2613]	2304
[5364]	2305	! Assigning values to vmin, vmax
	2306	vmin = 1
	2307	vmax = nvm*1._r_std
	2308
	2309	variabletypevals = -un
	2310
	2311	!! Variables for interpweight
	2312	! Type of calculation of cell fractions
	2313	fractype = 'default'
	2314	! Name of the longitude and latitude in the input file
	2315	lonname = 'lon'
	2316	latname = 'lat'
	2317	! Should negative values be set to zero from input file?
	2318	nonegative = .FALSE.
	2319	! Type of mask to apply to the input data (see header for more details)
	2320	maskingtype = 'msumrange'
	2321	! Values to use for the masking
	2322	maskvals = (/ 1.-1.e-7, 0., 2. /)
	2323	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
	2324	namemaskvar = ''
	2325
	2326	CALL interpweight_3D(nbpt, nvm, variabletypevals, lalo, resolution, neighbours, &
	2327	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
	2328	maskvals, namemaskvar, nvm, 0, veget_year, fractype, &
	2329	-1., -1., vegetrefrac, aveget)
	2330	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_readvegetmax after interpeeight_3D'
	2331	ENDIF
[8]	2332	!
	2333	! Compute the logical for partial (only anthropic) PTFs update
[2613]	2334	IF (ok_dgvm .AND. .NOT. init) THEN
	2335	partial_update= .TRUE.
[8]	2336	ELSE
[2613]	2337	partial_update=.FALSE.
	2338	END IF
	2339
[3831]	2340	IF (printlev_loc >= 5) THEN
	2341	WRITE(numout,*)' slowproc_readvegetmax before updating loop nbpt:', nbpt
	2342	END IF
	2343
[8]	2344	IF ( .NOT. partial_update ) THEN
	2345	veget_next(:,:)=zero
[2613]	2346
[3831]	2347	IF (printlev_loc >=3 .AND. ANY(aveget < min_sechiba)) THEN
	2348	WRITE(numout,*) 'Some grid cells on the model grid did not have any points on the source grid.'
	2349	IF (init) THEN
	2350	WRITE(numout,*) 'Initialization with full fraction of bare soil are done for the below grid cells.'
[8]	2351	ELSE
[3831]	2352	WRITE(numout,*) 'Old values are kept for the below grid cells.'
	2353	ENDIF
	2354	WRITE(numout,*) 'List of grid cells (ib, lat, lon):'
	2355	END IF
	2356
	2357	DO ib = 1, nbpt
	2358	! vegetrefrac is already normalized to sum equal one for each grid cell
	2359	veget_next(ib,:) = vegetrefrac(ib,:)
	2360
	2361	IF (aveget(ib) < min_sechiba) THEN
	2362	IF (printlev_loc >=3) WRITE(numout,*) ib,lalo(ib,1),lalo(ib,2)
[2613]	2363	IF (init) THEN
	2364	veget_next(ib,1) = un
	2365	veget_next(ib,2:nvm) = zero
[8]	2366	ELSE
	2367	veget_next(ib,:) = veget_last(ib,:)
	2368	ENDIF
	2369	ENDIF
	2370	ENDDO
	2371	ELSE
[3831]	2372	! Partial update
[8]	2373	DO ib = 1, nbpt
[3831]	2374	IF (aveget(ib) > min_sechiba) THEN
	2375	! For the case with properly interpolated grid cells (aveget>0)
	2376
	2377	! last veget for this point
	2378	sum_veg=SUM(veget_last(ib,:))
	2379	!
	2380	! If the DGVM is activated, only anthropic PFTs are utpdated, the others are copied from previous time-step
	2381	veget_next(ib,:) = veget_last(ib,:)
	2382
[8]	2383	DO jv = 2, nvm
	2384	IF ( .NOT. natural(jv) ) THEN
[3831]	2385	veget_next(ib,jv) = vegetrefrac(ib,jv)
[8]	2386	ENDIF
	2387	ENDDO
[3831]	2388
[8]	2389	sumvAnthro_old = zero
	2390	sumvAnthro = zero
	2391	DO jv = 2, nvm
	2392	IF ( .NOT. natural(jv) ) THEN
[317]	2393	sumvAnthro = sumvAnthro + veget_next(ib,jv)
[8]	2394	sumvAnthro_old = sumvAnthro_old + veget_last(ib,jv)
	2395	ENDIF
	2396	ENDDO
[372]	2397
	2398	IF ( sumvAnthro_old < sumvAnthro ) THEN
[3831]	2399	! Increase of non natural vegetations (increase of agriculture)
	2400	! The proportion of natural PFT's must be preserved
	2401	! ie the sum of vegets is preserved
	2402	! and natural PFT / (sum of veget - sum of antropic veget)
	2403	! is preserved.
[372]	2404	rapport = ( sum_veg - sumvAnthro ) / ( sum_veg - sumvAnthro_old )
	2405	DO jv = 1, nvm
	2406	IF ( natural(jv) ) THEN
	2407	veget_next(ib,jv) = veget_last(ib,jv) * rapport
	2408	ENDIF
	2409	ENDDO
	2410	ELSE
[3831]	2411	! Increase of natural vegetations (decrease of agriculture)
	2412	! The decrease of agriculture is replaced by bare soil. The DGVM will
	2413	! re-introduce natural PFT's.
[372]	2414	DO jv = 1, nvm
	2415	IF ( natural(jv) ) THEN
	2416	veget_next(ib,jv) = veget_last(ib,jv)
	2417	ENDIF
	2418	ENDDO
	2419	veget_next(ib,1) = veget_next(ib,1) + sumvAnthro_old - sumvAnthro
	2420	ENDIF
	2421
[8]	2422	! test
[317]	2423	IF ( ABS( SUM(veget_next(ib,:)) - sum_veg ) > 10*EPSILON(un) ) THEN
[3831]	2424	WRITE(numout,*) 'slowproc_readvegetmax _______'
	2425	msg = " No conservation of sum of veget for point "
	2426	WRITE(numout,*) TRIM(msg), ib, ",(", lalo(ib,1),",", lalo(ib,2), ")"
	2427	WRITE(numout,*) " last sum of veget ", sum_veg, " new sum of veget ", &
	2428	SUM(veget_next(ib,:)), " error : ", SUM(veget_next(ib,:))-sum_veg
	2429	WRITE(numout,*) " Anthropic modifications : last ",sumvAnthro_old," new ",sumvAnthro
	2430	CALL ipslerr_p (3,'slowproc_readvegetmax', &
	2431	& 'No conservation of sum of veget_next', &
	2432	& "The sum of veget_next is different after reading Land Use map.", &
[511]	2433	& '(verify the dgvm case model.)')
[8]	2434	ENDIF
	2435	ELSE
[3831]	2436	! For the case when there was a propblem with the interpolation, aveget < min_sechiba
	2437	WRITE(numout,*) 'slowproc_readvegetmax _______'
	2438	WRITE(numout,*) " No land point in the map for point ", ib, ",(", lalo(ib,1), ",", &
	2439	lalo(ib,2),")"
	2440	CALL ipslerr_p (2,'slowproc_readvegetmax', &
	2441	& 'Problem with vegetation file for Land Use.', &
	2442	& "No land point in the map for point", &
[8]	2443	& '(verify your land use file.)')
	2444	veget_next(ib,:) = veget_last(ib,:)
	2445	ENDIF
	2446
[3831]	2447	ENDDO
[8]	2448	ENDIF
[3831]	2449
	2450	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_readvegetmax after updating'
[8]	2451	!
	2452	frac_nobio_next (:,:) = un
	2453	!
	2454	!MM
	2455	! Work only for one nnobio !! (ie ice)
	2456	DO inobio=1,nnobio
	2457	DO jv=1,nvm
	2458	DO ib = 1, nbpt
	2459	frac_nobio_next(ib,inobio) = frac_nobio_next(ib,inobio) - veget_next(ib,jv)
	2460	ENDDO
	2461	ENDDO
	2462	ENDDO
[2718]	2463
[8]	2464	DO ib = 1, nbpt
	2465	sum_veg = SUM(veget_next(ib,:))
	2466	sum_nobio = SUM(frac_nobio_next(ib,:))
	2467	IF (sum_nobio < 0.) THEN
	2468	frac_nobio_next(ib,:) = zero
[2718]	2469	veget_next(ib,1) = veget_next(ib,1) + sum_nobio
[8]	2470	sum_veg = SUM(veget_next(ib,:))
	2471	ENDIF
	2472	sumf = sum_veg + sum_nobio
	2473	IF (sumf > min_sechiba) THEN
	2474	veget_next(ib,:) = veget_next(ib,:) / sumf
	2475	frac_nobio_next(ib,:) = frac_nobio_next(ib,:) / sumf
	2476	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
	2477	err=norm-un
[3831]	2478	IF (printlev_loc >=5) WRITE(numout,*) " slowproc_readvegetmax: ib ",ib, &
	2479	" SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un, sumf",err,sumf
[42]	2480	IF (abs(err) > -EPSILON(un)) THEN
[8]	2481	IF ( SUM(frac_nobio_next(ib,:)) > min_sechiba ) THEN
	2482	frac_nobio_next(ib,1) = frac_nobio_next(ib,1) - err
	2483	ELSE
	2484	veget_next(ib,1) = veget_next(ib,1) - err
	2485	ENDIF
	2486	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
	2487	err=norm-un
[3831]	2488	IF (printlev_loc >=5) WRITE(numout,*) " slowproc_readvegetmax: ib ", ib, &
	2489	" SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un",err
[42]	2490	IF (abs(err) > EPSILON(un)) THEN
[3831]	2491	WRITE(numout,*) ' slowproc_readvegetmax _______'
[8]	2492	WRITE(numout,*) "update : Problem with point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")"
	2493	WRITE(numout,*) " err(sum-1.) = ",abs(err)
[3094]	2494	CALL ipslerr_p (2,'slowproc_readvegetmax', &
[3831]	2495	& 'Problem with sum vegetation + sum fracnobio for Land Use.', &
[8]	2496	& "sum not equal to 1.", &
	2497	& '(verify your land use file.)')
[3831]	2498	aveget(ib) = -0.6
[8]	2499	ENDIF
	2500	ENDIF
	2501	ELSE
[3831]	2502	! sumf < min_sechiba
	2503	WRITE(numout,*) ' slowproc_readvegetmax _______'
	2504	WRITE(numout,*)" No vegetation nor frac_nobio for point ", ib, ",(", lalo(ib,1), ",", &
	2505	lalo(ib,2),")"
	2506	WRITE(numout,*)" Replaced by bare_soil !! "
[8]	2507	veget_next(ib,1) = un
	2508	veget_next(ib,2:nvm) = zero
	2509	frac_nobio_next(ib,:) = zero
[3094]	2510	!!!$ CALL ipslerr_p (3,'slowproc_readvegetmax', &
[1078]	2511	!!!$ & 'Problem with vegetation file for Land Use.', &
	2512	!!!$ & "No vegetation nor frac_nobio for point ", &
	2513	!!!$ & '(verify your land use file.)')
[8]	2514	ENDIF
	2515	ENDDO
[3831]	2516
[6145]	2517	!! Set to zero fractions of frac_nobio and veget_max smaller than min_vegfrac
	2518	!! Normalize to have the sum equal 1.
	2519	CALL slowproc_veget_max_limit(nbpt, frac_nobio_next, veget_next)
	2520
[3831]	2521	! Write diagnostics
	2522	CALL xios_orchidee_send_field("aveget",aveget)
[5364]	2523	CALL xios_orchidee_send_field("interp_diag_aveget",aveget)
	2524	CALL xios_orchidee_send_field("interp_diag_vegetrefrac",vegetrefrac)
[3831]	2525
	2526	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_readvegetmax ended'
[2613]	2527
[3094]	2528	END SUBROUTINE slowproc_readvegetmax
[8]	2529
[947]	2530
	2531	!! ================================================================================================================================
	2532	!! SUBROUTINE : slowproc_nearest
	2533	!!
	2534	!>\BRIEF looks for nearest grid point on the fine map
	2535	!!
	2536	!! DESCRIPTION : (definitions, functional, design, flags):
	2537	!!
	2538	!! RECENT CHANGE(S): None
	2539	!!
	2540	!! MAIN OUTPUT VARIABLE(S): ::inear
	2541	!!
	2542	!! REFERENCE(S) : None
	2543	!!
	2544	!! FLOWCHART : None
	2545	!! \n
	2546	!_ ================================================================================================================================
	2547
[8]	2548	SUBROUTINE slowproc_nearest(iml, lon5, lat5, lonmod, latmod, inear)
	2549
[947]	2550	!! INTERFACE DESCRIPTION
	2551
	2552	!! 0.1 input variables
[8]	2553
[947]	2554	INTEGER(i_std), INTENT(in) :: iml !! size of the vector
	2555	REAL(r_std), DIMENSION(iml), INTENT(in) :: lon5, lat5 !! longitude and latitude vector, for the 5km vegmap
	2556	REAL(r_std), INTENT(in) :: lonmod, latmod !! longitude and latitude modelled
[8]	2557
[947]	2558	!! 0.2 output variables
	2559
	2560	INTEGER(i_std), INTENT(out) :: inear !! location of the grid point from the 5km vegmap grid
	2561	!! closest from the modelled grid point
	2562
	2563	!! 0.4 Local variables
	2564
[8]	2565	REAL(r_std) :: pa, p
	2566	REAL(r_std) :: coscolat, sincolat
	2567	REAL(r_std) :: cospa, sinpa
	2568	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: cosang
	2569	INTEGER(i_std) :: i
	2570	INTEGER(i_std), DIMENSION(1) :: ineartab
[947]	2571	INTEGER :: ALLOC_ERR
[8]	2572
[947]	2573	!_ ================================================================================================================================
	2574
[8]	2575	ALLOCATE(cosang(iml), STAT=ALLOC_ERR)
[2613]	2576	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_nearest','Error in allocation for cosang','','')
[8]	2577
[947]	2578	pa = pi/2.0 - latmod*pi/180.0 ! dist. between north pole and the point a
	2579	!! COLATITUDE, in radian
[8]	2580	cospa = COS(pa)
	2581	sinpa = SIN(pa)
	2582
	2583	DO i = 1, iml
	2584
[947]	2585	sincolat = SIN( pi/2.0 - lat5(i)*pi/180.0 ) !! sinus of the colatitude
	2586	coscolat = COS( pi/2.0 - lat5(i)*pi/180.0 ) !! cosinus of the colatitude
[8]	2587
[947]	2588	p = (lonmod-lon5(i))*pi/180.0 !! angle between a & b (between their meridian)in radians
[8]	2589
[947]	2590	!! dist(i) = ACOS( cospacoscolat + sinpasincolat*COS(p))
	2591	cosang(i) = cospacoscolat + sinpasincolat*COS(p) !! TL : cosang is maximum when angle is at minimal value
	2592	!! orthodromic distance between 2 points : cosang = cosinus (arc(AB)/R), with
	2593	!R = Earth radius, then max(cosang) = max(cos(arc(AB)/R)), reached when arc(AB)/R is minimal, when
	2594	! arc(AB) is minimal, thus when point B (corresponding grid point from LAI MAP) is the nearest from
	2595	! modelled A point
[8]	2596	ENDDO
	2597
	2598	ineartab = MAXLOC( cosang(:) )
	2599	inear = ineartab(1)
	2600
	2601	DEALLOCATE(cosang)
	2602	END SUBROUTINE slowproc_nearest
	2603
[947]	2604	!! ================================================================================================================================
	2605	!! SUBROUTINE : slowproc_soilt
	2606	!!
[2419]	2607	!>\BRIEF Interpolate the Zobler or Reynolds/USDA soil type map
[947]	2608	!!
	2609	!! DESCRIPTION : (definitions, functional, design, flags):
	2610	!!
[2419]	2611	!! RECENT CHANGE(S): Nov 2014, ADucharne
[947]	2612	!!
[4808]	2613	!! MAIN OUTPUT VARIABLE(S): ::soiltype, ::clayfraction, sandfraction, siltfraction
[947]	2614	!!
[2419]	2615	!! REFERENCE(S) : Reynold, Jackson, and Rawls (2000). Estimating soil water-holding capacities
	2616	!! by linking the Food and Agriculture Organization soil map of the world with global pedon
	2617	!! databases and continuous pedotransfer functions, WRR, 36, 3653-3662
[947]	2618	!!
	2619	!! FLOWCHART : None
	2620	!! \n
	2621	!_ ================================================================================================================================
[4808]	2622	SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction, sandfraction, siltfraction)
[947]	2623
[3831]	2624	USE interpweight
	2625
	2626	IMPLICIT NONE
[8]	2627	!
	2628	!
[2419]	2629	! This subroutine should read the Zobler/Reynolds map and interpolate to the model grid.
	2630	! The method is to get fraction of the three/12 main soiltypes for each grid box.
	2631	! For the Zobler case, also called FAO in the code, the soil fraction are going to be put
	2632	! into the array soiltype in the following order : coarse, medium and fine.
	2633	! For the Reynolds/USDA case, the soiltype array follows the order defined in constantes_soil_var.f90
[8]	2634	!
	2635	!
[947]	2636	!! 0.1 INPUT
[8]	2637	!
[947]	2638	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
	2639	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
[3447]	2640	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
	2641	!! (1=North and then clockwise)
[947]	2642	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
	2643	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
[8]	2644	!
	2645	! 0.2 OUTPUT
	2646	!
[2344]	2647	REAL(r_std), INTENT(out) :: soilclass(nbpt, nscm) !! Soil type map to be created from the Zobler map
[2419]	2648	!! or a map defining the 12 USDA classes (e.g. Reynolds)
	2649	!! Holds the area of each texture class in the ORCHIDEE grid cells
	2650	!! Final unit = fraction of ORCHIDEE grid-cell (unitless)
[947]	2651	REAL(r_std), INTENT(out) :: clayfraction(nbpt) !! The fraction of clay as used by STOMATE
[4808]	2652	REAL(r_std), INTENT(out) :: sandfraction(nbpt) !! The fraction of sand (for SP-MIP)
	2653	REAL(r_std), INTENT(out) :: siltfraction(nbpt) !! The fraction of silt (for SP-MIP)
[8]	2654	!
	2655	!
	2656	! 0.3 LOCAL
	2657	!
	2658	CHARACTER(LEN=80) :: filename
[3831]	2659	INTEGER(i_std) :: ib, ilf, nbexp, i
[2419]	2660	INTEGER(i_std) :: fopt !! Nb of pts from the texture map within one ORCHIDEE grid-cell
	2661	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: solt !! Texture the different points from the input texture map
	2662	!! in one ORCHIDEE grid cell (unitless)
[8]	2663	!
	2664	! Number of texture classes in Zobler
	2665	!
[2419]	2666	INTEGER(i_std), PARAMETER :: nzobler = 7 !! Nb of texture classes according in the Zobler map
	2667	REAL(r_std),ALLOCATABLE :: textfrac_table(:,:) !! conversion table between the texture index
	2668	!! and the granulometric composition
[8]	2669	!
	2670	INTEGER :: ALLOC_ERR
[3831]	2671	INTEGER :: ntextinfile !! number of soil textures in the in the file
	2672	REAL(r_std), DIMENSION(:,:), ALLOCATABLE :: textrefrac !! text fractions re-dimensioned
[4808]	2673	REAL(r_std), DIMENSION(nbpt) :: atext !! Availability of the texture interpolation
[3831]	2674	REAL(r_std) :: vmin, vmax !! min/max values to use for the
[1082]	2675
[3831]	2676	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
	2677	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in input file
	2678	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
	2679	!! (variabletypevals(1) = -un, not used)
	2680	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
	2681	!! 'XYKindTime': Input values are kinds
	2682	!! of something with a temporal
	2683	!! evolution on the dx*dy matrix'
	2684	LOGICAL :: nonegative !! whether negative values should be removed
	2685	CHARACTER(LEN=50) :: maskingtype !! Type of masking
	2686	!! 'nomask': no-mask is applied
	2687	!! 'mbelow': take values below maskvals(1)
	2688	!! 'mabove': take values above maskvals(1)
	2689	!! 'msumrange': take values within 2 ranges;
	2690	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
	2691	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
	2692	!! (normalized by maskvals(3))
	2693	!! 'var': mask values are taken from a
	2694	!! variable inside the file (>0)
	2695	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
	2696	!! `maskingtype')
	2697	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
	2698	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: vecpos
[4576]	2699	REAL(r_std) :: sgn !! sum of fractions excluding glaciers and ocean
[1082]	2700	!_ ================================================================================================================================
[3831]	2701
	2702	IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt'
[8]	2703	!
	2704	! Needs to be a configurable variable
	2705	!
	2706	!
[947]	2707	!Config Key = SOILCLASS_FILE
[566]	2708	!Config Desc = Name of file from which soil types are read
[620]	2709	!Config Def = soils_param.nc
[566]	2710	!Config If = NOT(IMPOSE_VEG)
	2711	!Config Help = The name of the file to be opened to read the soil types.
	2712	!Config The data from this file is then interpolated to the grid of
	2713	!Config of the model. The aim is to get fractions for sand loam and
	2714	!Config clay in each grid box. This information is used for soil hydrology
	2715	!Config and respiration.
[843]	2716	!Config Units = [FILE]
[8]	2717	!
[2419]	2718	! soils_param.nc file is 1deg soil texture file (Zobler)
	2719	! The USDA map from Reynolds is soils_param_usda.nc (1/12deg resolution)
	2720
[620]	2721	filename = 'soils_param.nc'
[947]	2722	CALL getin_p('SOILCLASS_FILE',filename)
[2613]	2723
[3831]	2724	variablename = 'soiltext'
[2613]	2725
[3831]	2726	!! Variables for interpweight
	2727	! Type of calculation of cell fractions
	2728	fractype = 'default'
[2613]	2729
[5364]	2730	IF (xios_interpolation) THEN
	2731	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Use XIOS to read and interpolate " &
	2732	// TRIM(filename) // " for variable " // TRIM(variablename)
	2733
	2734	SELECT CASE(soil_classif)
[2613]	2735
[5364]	2736	CASE('none')
	2737	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
	2738	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
	2739	DO ib=1, nbpt
	2740	soilclass(ib,:) = soilclass_default_fao
	2741	clayfraction(ib) = clayfraction_default
	2742	ENDDO
[2613]	2743
[5364]	2744
[3831]	2745	CASE('zobler')
[2613]	2746
[5364]	2747	!
	2748	soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
	2749	!
	2750	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
	2751	!
	2752	ALLOCATE(textrefrac(nbpt,nzobler))
	2753	ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
	2754	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
	2755	CALL get_soilcorr_zobler (nzobler, textfrac_table)
	2756
	2757	CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
	2758	CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
	2759	CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
	2760	CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
	2761	CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
	2762	CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
	2763	CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
[2613]	2764
[5364]	2765
[3831]	2766
[5364]	2767	CALL get_soilcorr_zobler (nzobler, textfrac_table)
	2768	!
	2769	!
	2770	DO ib =1, nbpt
	2771	soilclass(ib,1)=textrefrac(ib,1)
	2772	soilclass(ib,2)=textrefrac(ib,2)+textrefrac(ib,3)+textrefrac(ib,4)+textrefrac(ib,7)
	2773	soilclass(ib,3)=textrefrac(ib,5)
	2774
	2775	! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
	2776	! over the zobler pixels composing the ORCHIDEE grid-cell
	2777	clayfraction(ib) = textfrac_table(1,3) * textrefrac(ib,1)+textfrac_table(2,3) * textrefrac(ib,2) + &
	2778	textfrac_table(3,3) * textrefrac(ib,3)+textfrac_table(4,3) * textrefrac(ib,4) + &
	2779	textfrac_table(5,3) * textrefrac(ib,5)+textfrac_table(7,3) * textrefrac(ib,7)
[6144]	2780
[5364]	2781	sandfraction(ib) = textfrac_table(1,2) * textrefrac(ib,1)+textfrac_table(2,2) * textrefrac(ib,2) + &
	2782	textfrac_table(3,2) * textrefrac(ib,3)+textfrac_table(4,2) * textrefrac(ib,4) + &
	2783	textfrac_table(5,2) * textrefrac(ib,5)+textfrac_table(7,2) * textrefrac(ib,7)
	2784
[6144]	2785	siltfraction(ib) = textfrac_table(1,1) * textrefrac(ib,1)+textfrac_table(2,1) * textrefrac(ib,2) + &
	2786	textfrac_table(3,1) * textrefrac(ib,3)+textfrac_table(4,1) * textrefrac(ib,4) + &
	2787	textfrac_table(5,1) * textrefrac(ib,5)+textfrac_table(7,1) * textrefrac(ib,7)
[5364]	2788
	2789	sgn=SUM(soilclass(ib,1:3))
	2790
	2791	IF (sgn < min_sechiba) THEN
	2792	soilclass(ib,:) = soilclass_default(:)
	2793	clayfraction(ib) = clayfraction_default
	2794	sandfraction(ib) = sandfraction_default
	2795	siltfraction(ib) = siltfraction_default
	2796	atext(ib)=0.
	2797	ELSE
	2798	atext(ib)=sgn
	2799	clayfraction(ib) = clayfraction(ib) / sgn
	2800	sandfraction(ib) = sandfraction(ib) / sgn
	2801	siltfraction(ib) = siltfraction(ib) / sgn
	2802	soilclass(ib,1:3) = soilclass(ib,1:3) / sgn
	2803	ENDIF
	2804
	2805	ENDDO
	2806
	2807
[3831]	2808
[5364]	2809	CASE('usda')
	2810
	2811	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: start case usda'
	2812
	2813	soilclass_default=soilclass_default_usda
	2814	!
	2815	WRITE(numout,*) "Using a soilclass map with usda classification"
	2816	!
	2817	ALLOCATE(textrefrac(nbpt,nscm))
	2818	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
	2819	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
	2820
	2821	CALL get_soilcorr_usda (nscm, textfrac_table)
	2822
	2823	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
	2824
	2825	CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
	2826	CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
	2827	CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
	2828	CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
	2829	CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
	2830	CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
	2831	CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
	2832	CALL xios_orchidee_recv_field('soiltext8',textrefrac(:,8))
	2833	CALL xios_orchidee_recv_field('soiltext9',textrefrac(:,9))
	2834	CALL xios_orchidee_recv_field('soiltext10',textrefrac(:,10))
	2835	CALL xios_orchidee_recv_field('soiltext11',textrefrac(:,11))
	2836	CALL xios_orchidee_recv_field('soiltext12',textrefrac(:,12))
	2837
	2838
[3831]	2839
[5364]	2840	CALL get_soilcorr_usda (nscm, textfrac_table)
	2841	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
[2613]	2842
[5364]	2843	DO ib =1, nbpt
	2844	clayfraction(ib) = 0.0
	2845	DO ilf = 1,nscm
	2846	soilclass(ib,ilf)=textrefrac(ib,ilf)
	2847	clayfraction(ib) = clayfraction(ib) + textfrac_table(ilf,3)*textrefrac(ib,ilf)
	2848	sandfraction(ib) = sandfraction(ib) + textfrac_table(ilf,2)*textrefrac(ib,ilf)
	2849	siltfraction(ib) = siltfraction(ib) + textfrac_table(ilf,1)*textrefrac(ib,ilf)
	2850	ENDDO
	2851
	2852
	2853	sgn=SUM(soilclass(ib,:))
	2854
	2855	IF (sgn < min_sechiba) THEN
	2856	soilclass(ib,:) = soilclass_default(:)
	2857	clayfraction(ib) = clayfraction_default
	2858	sandfraction(ib) = sandfraction_default
	2859	siltfraction(ib) = siltfraction_default
	2860	atext(ib)=0
	2861	ELSE
	2862	soilclass(ib,:) = soilclass(ib,:) / sgn
	2863	clayfraction(ib) = clayfraction(ib) / sgn
	2864	sandfraction(ib) = sandfraction(ib) / sgn
	2865	siltfraction(ib) = siltfraction(ib) / sgn
	2866	atext(ib)=sgn
	2867	ENDIF
	2868	ENDDO
	2869
	2870	CASE DEFAULT
	2871	WRITE(numout,*) 'slowproc_soilt:'
	2872	WRITE(numout,*) ' A non supported soil type classification has been chosen'
	2873	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
	2874	END SELECT
	2875
	2876
	2877
	2878	ELSE ! xios_interpolation
	2879	! Read and interpolate using stardard method with IOIPSL and aggregate
	2880
	2881	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " &
	2882	// TRIM(filename) // " for variable " // TRIM(variablename)
	2883
	2884
	2885	! Name of the longitude and latitude in the input file
	2886	lonname = 'nav_lon'
	2887	latname = 'nav_lat'
	2888
	2889	IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_soilt: Start interpolate " &
	2890	// TRIM(filename) // " for variable " // TRIM(variablename)
	2891
	2892	IF ( TRIM(soil_classif) /= 'none' ) THEN
	2893
	2894	! Define a variable for the number of soil textures in the input file
	2895	SELECTCASE(soil_classif)
	2896	CASE('zobler')
	2897	ntextinfile=nzobler
	2898	CASE('usda')
	2899	ntextinfile=nscm
	2900	CASE DEFAULT
	2901	WRITE(numout,*) 'slowproc_soilt:'
	2902	WRITE(numout,*) ' A non supported soil type classification has been chosen'
	2903	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
	2904	ENDSELECT
	2905
	2906	ALLOCATE(textrefrac(nbpt,ntextinfile), STAT=ALLOC_ERR)
	2907	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable textrefrac',&
	2908	'','')
	2909
	2910	! Assigning values to vmin, vmax
	2911	vmin = un
	2912	vmax = ntextinfile*un
	2913
	2914	ALLOCATE(variabletypevals(ntextinfile), STAT=ALLOC_ERR)
	2915	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variabletypevals','','')
	2916	variabletypevals = -un
	2917
	2918	!! Variables for interpweight
	2919	! Should negative values be set to zero from input file?
	2920	nonegative = .FALSE.
	2921	! Type of mask to apply to the input data (see header for more details)
	2922	maskingtype = 'mabove'
	2923	! Values to use for the masking
	2924	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
	2925	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used)
	2926	namemaskvar = ''
	2927
	2928	CALL interpweight_2D(nbpt, ntextinfile, variabletypevals, lalo, resolution, neighbours, &
	2929	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
	2930	maskvals, namemaskvar, 0, 0, -1, fractype, -1., -1., textrefrac, atext)
	2931
	2932	ALLOCATE(vecpos(ntextinfile), STAT=ALLOC_ERR)
	2933	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable vecpos','','')
	2934	ALLOCATE(solt(ntextinfile), STAT=ALLOC_ERR)
	2935	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable solt','','')
	2936
	2937	IF (printlev_loc >= 5) THEN
	2938	WRITE(numout,*)' slowproc_soilt after interpweight_2D'
	2939	WRITE(numout,*)' slowproc_soilt before starting loop nbpt:', nbpt
	2940	WRITE(numout,*)" slowproc_soilt starting classification '" // TRIM(soil_classif) // "'..."
	2941	END IF
	2942	ELSE
	2943	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_soilt using default values all points are propertly ' // &
	2944	'interpolated atext = 1. everywhere!'
	2945	atext = 1.
[3831]	2946	END IF
[2613]	2947
[3831]	2948	nbexp = 0
[947]	2949	SELECTCASE(soil_classif)
	2950	CASE('none')
[2613]	2951	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
	2952	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
[947]	2953	DO ib=1, nbpt
	2954	soilclass(ib,:) = soilclass_default_fao
	2955	clayfraction(ib) = clayfraction_default
[4808]	2956	sandfraction(ib) = sandfraction_default
	2957	siltfraction(ib) = siltfraction_default
[947]	2958	ENDDO
	2959	CASE('zobler')
[8]	2960	!
[2419]	2961	soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
[8]	2962	!
[3831]	2963	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
[8]	2964	!
[2613]	2965	ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
	2966	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
[2419]	2967	CALL get_soilcorr_zobler (nzobler, textfrac_table)
[8]	2968	!
	2969	!
[3831]	2970	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_soilt after getting table of textures'
[947]	2971	DO ib =1, nbpt
	2972	soilclass(ib,:) = zero
	2973	clayfraction(ib) = zero
[4808]	2974	sandfraction(ib) = zero
	2975	siltfraction(ib) = zero
[947]	2976	!
[3831]	2977	! vecpos: List of positions where textures were not zero
	2978	! vecpos(1): number of not null textures found
	2979	vecpos = interpweight_ValVecR(textrefrac(ib,:),nzobler,zero,'neq')
	2980	fopt = vecpos(1)
[4576]	2981
[3831]	2982	IF ( fopt .EQ. 0 ) THEN
[2419]	2983	! No points were found for current grid box, use default values
[947]	2984	nbexp = nbexp + 1
	2985	soilclass(ib,:) = soilclass_default(:)
	2986	clayfraction(ib) = clayfraction_default
[4808]	2987	sandfraction(ib) = sandfraction_default
	2988	siltfraction(ib) = siltfraction_default
[947]	2989	ELSE
[3831]	2990	IF (fopt == nzobler) THEN
	2991	! All textures are not zero
	2992	solt=(/(i,i=1,nzobler)/)
	2993	ELSE
	2994	DO ilf = 1,fopt
	2995	solt(ilf) = vecpos(ilf+1)
	2996	END DO
	2997	END IF
[8]	2998	!
[2419]	2999	! Compute the fraction of each textural class
[947]	3000	!
[4576]	3001	sgn = 0.
[947]	3002	DO ilf = 1,fopt
[5364]	3003	!
	3004	! Here we make the correspondance between the 7 zobler textures and the 3 textures in ORCHIDEE
	3005	! and soilclass correspond to surfaces covered by the 3 textures of ORCHIDEE (coase,medium,fine)
	3006	! For type 6 = glacier, default values are set and it is also taken into account during the normalization
	3007	! of the fractions (done in interpweight_2D)
	3008	! Note that type 0 corresponds to ocean but it is already removed using the mask above.
	3009	!
[4576]	3010	IF ( (solt(ilf) .LE. nzobler) .AND. (solt(ilf) .GT. 0) .AND. &
	3011	(solt(ilf) .NE. 6) ) THEN
[3831]	3012	SELECT CASE(solt(ilf))
	3013	CASE(1)
	3014	soilclass(ib,1) = soilclass(ib,1) + textrefrac(ib,solt(ilf))
	3015	CASE(2)
	3016	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
	3017	CASE(3)
	3018	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
	3019	CASE(4)
	3020	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
	3021	CASE(5)
	3022	soilclass(ib,3) = soilclass(ib,3) + textrefrac(ib,solt(ilf))
	3023	CASE(7)
	3024	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
	3025	CASE DEFAULT
	3026	WRITE(numout,*) 'We should not be here, an impossible case appeared'
	3027	CALL ipslerr_p(3,'slowproc_soilt','Bad value for solt','','')
[947]	3028	END SELECT
[2419]	3029	! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
	3030	! over the zobler pixels composing the ORCHIDEE grid-cell
[947]	3031	clayfraction(ib) = clayfraction(ib) + &
[3831]	3032	& textfrac_table(solt(ilf),3) * textrefrac(ib,solt(ilf))
[4808]	3033	sandfraction(ib) = sandfraction(ib) + &
	3034	& textfrac_table(solt(ilf),2) * textrefrac(ib,solt(ilf))
	3035	siltfraction(ib) = siltfraction(ib) + &
	3036	& textfrac_table(solt(ilf),1) * textrefrac(ib,solt(ilf))
[4576]	3037	! Sum the fractions which are not glaciers nor ocean
	3038	sgn = sgn + textrefrac(ib,solt(ilf))
[947]	3039	ELSE
	3040	IF (solt(ilf) .GT. nzobler) THEN
	3041	WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
[2613]	3042	CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible','','')
[947]	3043	ENDIF
[5364]	3044	END IF
[947]	3045	ENDDO
[3831]	3046
[4576]	3047	IF ( sgn .LT. min_sechiba) THEN
	3048	! Set default values if grid cells were only covered by glaciers or ocean
	3049	! or if now information on the source grid was found.
[947]	3050	nbexp = nbexp + 1
	3051	soilclass(ib,:) = soilclass_default(:)
	3052	clayfraction(ib) = clayfraction_default
[4808]	3053	sandfraction(ib) = sandfraction_default
	3054	siltfraction(ib) = siltfraction_default
[4576]	3055	ELSE
	3056	! Normalize using the fraction of surface not including glaciers and ocean
	3057	soilclass(ib,:) = soilclass(ib,:)/sgn
	3058	clayfraction(ib) = clayfraction(ib)/sgn
[4808]	3059	sandfraction(ib) = sandfraction(ib)/sgn
	3060	siltfraction(ib) = siltfraction(ib)/sgn
[947]	3061	ENDIF
	3062	ENDIF
	3063	ENDDO
[5364]	3064
	3065	! The "USDA" case reads a map of the 12 USDA texture classes,
	3066	! such as to assign the corresponding soil properties
	3067	CASE("usda")
	3068	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with usda classification"
[2419]	3069
[5364]	3070	soilclass_default=soilclass_default_usda
[3831]	3071
[5364]	3072	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
	3073	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
[2613]	3074
[5364]	3075	CALL get_soilcorr_usda (nscm, textfrac_table)
[2419]	3076
[5364]	3077	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
[947]	3078	!
[5364]	3079	DO ib =1, nbpt
[947]	3080	! GO through the point we have found
	3081	!
	3082	!
[3831]	3083	! Provide which textures were found
	3084	! vecpos: List of positions where textures were not zero
	3085	! vecpos(1): number of not null textures found
	3086	vecpos = interpweight_ValVecR(textrefrac(ib,:),ntextinfile,zero,'neq')
	3087	fopt = vecpos(1)
	3088
[947]	3089	!
	3090	! Check that we found some points
	3091	!
	3092	soilclass(ib,:) = 0.0
	3093	clayfraction(ib) = 0.0
[3831]	3094
[947]	3095	IF ( fopt .EQ. 0) THEN
[3831]	3096	! No points were found for current grid box, use default values
	3097	IF (printlev_loc>=3) WRITE(numout,*)'slowproc_soilt: no soil class in input file found for point=', ib
[947]	3098	nbexp = nbexp + 1
	3099	soilclass(ib,:) = soilclass_default
	3100	clayfraction(ib) = clayfraction_default
[4808]	3101	sandfraction(ib) = sandfraction_default
	3102	siltfraction(ib) = siltfraction_default
[947]	3103	ELSE
[3831]	3104	IF (fopt == nscm) THEN
	3105	! All textures are not zero
	3106	solt(:) = (/(i,i=1,nscm)/)
	3107	ELSE
	3108	DO ilf = 1,fopt
	3109	solt(ilf) = vecpos(ilf+1)
	3110	END DO
	3111	END IF
[5364]	3112
	3113	!
	3114	!
	3115	! Compute the fraction of each textural class
	3116	!
	3117	!
	3118	DO ilf = 1,fopt
	3119	IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
	3120	soilclass(ib,solt(ilf)) = textrefrac(ib,solt(ilf))
	3121	clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) * &
	3122	textrefrac(ib,solt(ilf))
	3123	sandfraction(ib) = sandfraction(ib) + textfrac_table(solt(ilf),2) * &
	3124	textrefrac(ib,solt(ilf))
	3125	siltfraction(ib) = siltfraction(ib) + textfrac_table(solt(ilf),1) * &
[3831]	3126	textrefrac(ib,solt(ilf))
[5364]	3127	ELSE
	3128	IF (solt(ilf) .GT. nscm) THEN
	3129	WRITE(,) 'The file contains a soil color class which is incompatible with this program'
	3130	CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible 2','','')
	3131	ENDIF
[947]	3132	ENDIF
[5364]	3133	!
	3134	ENDDO
	3135
	3136	! Set default values if the surface in source file is too small
	3137	IF ( atext(ib) .LT. min_sechiba) THEN
	3138	nbexp = nbexp + 1
	3139	soilclass(ib,:) = soilclass_default(:)
	3140	clayfraction(ib) = clayfraction_default
	3141	sandfraction(ib) = sandfraction_default
	3142	siltfraction(ib) = siltfraction_default
[947]	3143	ENDIF
	3144	ENDIF
[3831]	3145
[5364]	3146	ENDDO
	3147
	3148	IF (printlev_loc>=4) WRITE (numout,*) ' slowproc_soilt: End case usda'
	3149
	3150	CASE DEFAULT
	3151	WRITE(numout,*) 'slowproc_soilt _______'
	3152	WRITE(numout,*) ' A non supported soil type classification has been chosen'
	3153	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
	3154	ENDSELECT
	3155	IF (printlev_loc >= 5 ) WRITE(numout,*)' slowproc_soilt end of type classification'
[3831]	3156
[5364]	3157	IF ( nbexp .GT. 0 ) THEN
	3158	WRITE(numout,*) 'slowproc_soilt:'
	3159	WRITE(numout,*) ' The interpolation of the bare soil albedo had ', nbexp
	3160	WRITE(numout,*) ' points without data. This are either coastal points or ice covered land.'
	3161	WRITE(numout,*) ' The problem was solved by using the default soil types.'
	3162	ENDIF
[2581]	3163
[5364]	3164	IF (ALLOCATED(variabletypevals)) DEALLOCATE (variabletypevals)
	3165	IF (ALLOCATED(textrefrac)) DEALLOCATE (textrefrac)
	3166	IF (ALLOCATED(solt)) DEALLOCATE (solt)
	3167	IF (ALLOCATED(textfrac_table)) DEALLOCATE (textfrac_table)
	3168
	3169	ENDIF ! xios_interpolation
	3170
[3831]	3171	! Write diagnostics
	3172	CALL xios_orchidee_send_field("atext",atext)
[5364]	3173
	3174	CALL xios_orchidee_send_field("interp_diag_atext",atext)
	3175	CALL xios_orchidee_send_field("interp_diag_soilclass",soilclass)
	3176	CALL xios_orchidee_send_field("interp_diag_clayfraction",clayfraction)
[3831]	3177
	3178	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_soilt ended'
	3179
[8]	3180	END SUBROUTINE slowproc_soilt
[511]	3181
[721]	3182	!! ================================================================================================================================
[947]	3183	!! SUBROUTINE : slowproc_slope
	3184	!!
	3185	!>\BRIEF Calculate mean slope coef in each model grid box from the slope map
	3186	!!
	3187	!! DESCRIPTION : (definitions, functional, design, flags):
	3188	!!
	3189	!! RECENT CHANGE(S): None
	3190	!!
	3191	!! MAIN OUTPUT VARIABLE(S): ::reinf_slope
	3192	!!
	3193	!! REFERENCE(S) : None
	3194	!!
	3195	!! FLOWCHART : None
	3196	!! \n
	3197	!_ ================================================================================================================================
	3198
	3199	SUBROUTINE slowproc_slope(nbpt, lalo, neighbours, resolution, contfrac, reinf_slope)
[3831]	3200
	3201	USE interpweight
	3202
	3203	IMPLICIT NONE
	3204
[947]	3205	!
	3206	!
	3207	!
	3208	! 0.1 INPUT
	3209	!
	3210	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
	3211	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
[3447]	3212	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)! Vector of neighbours for each grid point
	3213	! (1=North and then clockwise)
[947]	3214	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
	3215	REAL(r_std), INTENT (in) :: contfrac(nbpt) !! Fraction of continent in the grid
	3216	!
	3217	! 0.2 OUTPUT
	3218	!
	3219	REAL(r_std), INTENT(out) :: reinf_slope(nbpt) ! slope coef
	3220	!
	3221	! 0.3 LOCAL
	3222	!
	3223	!
	3224	REAL(r_std) :: slope_noreinf ! Slope above which runoff is maximum
	3225	CHARACTER(LEN=80) :: filename
[3831]	3226	REAL(r_std) :: vmin, vmax !! min/max values to use for the
	3227	!! renormalization
	3228	REAL(r_std), DIMENSION(nbpt) :: aslope !! slope availability
	3229
	3230	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
	3231	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in the input file
	3232	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
	3233	!! 'XYKindTime': Input values are kinds
	3234	!! of something with a temporal
	3235	!! evolution on the dx*dy matrix'
	3236	LOGICAL :: nonegative !! whether negative values should be removed
	3237	CHARACTER(LEN=50) :: maskingtype !! Type of masking
	3238	!! 'nomask': no-mask is applied
	3239	!! 'mbelow': take values below maskvals(1)
	3240	!! 'mabove': take values above maskvals(1)
	3241	!! 'msumrange': take values within 2 ranges;
	3242	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
	3243	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
	3244	!! (normalized by maskvals(3))
	3245	!! 'var': mask values are taken from a
	3246	!! variable inside the file (>0)
	3247	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
	3248	!! `maskingtype')
	3249	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
	3250
[947]	3251	!_ ================================================================================================================================
[3831]	3252
[947]	3253	!
	3254	!Config Key = SLOPE_NOREINF
	3255	!Config Desc = See slope_noreinf above
	3256	!Config If =
	3257	!Config Def = 0.5
	3258	!Config Help = The slope above which there is no reinfiltration
	3259	!Config Units = [-]
	3260	!
	3261	slope_noreinf = 0.5
	3262	!
	3263	CALL getin_p('SLOPE_NOREINF',slope_noreinf)
	3264	!
	3265	!Config Key = TOPOGRAPHY_FILE
	3266	!Config Desc = Name of file from which the topography map is to be read
	3267	!Config If =
	3268	!Config Def = cartepente2d_15min.nc
	3269	!Config Help = The name of the file to be opened to read the orography
	3270	!Config map is to be given here. Usualy SECHIBA runs with a 2'
	3271	!Config map which is derived from the NGDC one.
	3272	!Config Units = [FILE]
	3273	!
	3274	filename = 'cartepente2d_15min.nc'
	3275	CALL getin_p('TOPOGRAPHY_FILE',filename)
[2613]	3276
[5364]	3277	IF (xios_interpolation) THEN
	3278
	3279	CALL xios_orchidee_recv_field('reinf_slope_interp',reinf_slope)
	3280	CALL xios_orchidee_recv_field('frac_slope_interp',aslope)
[2613]	3281
	3282
[5364]	3283	ELSE
	3284
	3285	variablename = 'pente'
	3286	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_slope: Read and interpolate " &
	3287	// TRIM(filename) // " for variable " // TRIM(variablename)
[2613]	3288
[5364]	3289	! For this case there are not types/categories. We have 'only' a continuos field
	3290	! Assigning values to vmin, vmax
	3291	vmin = 0.
	3292	vmax = 9999.
[2613]	3293
[5364]	3294	!! Variables for interpweight
	3295	! Type of calculation of cell fractions
	3296	fractype = 'slopecalc'
	3297	! Name of the longitude and latitude in the input file
	3298	lonname = 'longitude'
	3299	latname = 'latitude'
	3300	! Should negative values be set to zero from input file?
	3301	nonegative = .FALSE.
	3302	! Type of mask to apply to the input data (see header for more details)
	3303	maskingtype = 'mabove'
	3304	! Values to use for the masking
	3305	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
	3306	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
	3307	namemaskvar = ''
	3308
	3309	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
	3310	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
	3311	maskvals, namemaskvar, -1, fractype, slope_default, slope_noreinf, &
	3312	reinf_slope, aslope)
	3313	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_slope after interpweight_2Dcont'
	3314
	3315	ENDIF
	3316
	3317	! Write diagnostics
[3831]	3318	CALL xios_orchidee_send_field("aslope",aslope)
[5364]	3319	CALL xios_orchidee_send_field("interp_diag_aslope",aslope)
[2613]	3320
[5364]	3321	CALL xios_orchidee_send_field("interp_diag_reinf_slope",reinf_slope)
	3322
[3831]	3323	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_slope ended'
[2613]	3324
[947]	3325	END SUBROUTINE slowproc_slope
	3326
[4657]	3327
[947]	3328	!! ================================================================================================================================
[4657]	3329	!! SUBROUTINE : slowproc_woodharvest
	3330	!!
	3331	!>\BRIEF
	3332	!!
	3333	!! DESCRIPTION :
	3334	!!
	3335	!! RECENT CHANGE(S): None
	3336	!!
	3337	!! MAIN OUTPUT VARIABLE(S): ::
	3338	!!
	3339	!! REFERENCE(S) : None
	3340	!!
	3341	!! FLOWCHART : None
	3342	!! \n
	3343	!_ ================================================================================================================================
	3344
	3345	SUBROUTINE slowproc_woodharvest(nbpt, lalo, neighbours, resolution, contfrac, woodharvest)
	3346
	3347	USE interpweight
	3348
	3349	IMPLICIT NONE
	3350
	3351	!
	3352	!
	3353	!
	3354	! 0.1 INPUT
	3355	!
	3356	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
	3357	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
	3358	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
	3359	!! (1=North and then clockwise)
	3360	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
	3361	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of continent in the grid
	3362	!
	3363	! 0.2 OUTPUT
	3364	!
	3365	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: woodharvest !! Wood harvest
	3366	!
	3367	! 0.3 LOCAL
	3368	!
	3369	CHARACTER(LEN=80) :: filename
	3370	REAL(r_std) :: vmin, vmax
	3371	REAL(r_std), DIMENSION(nbpt) :: aoutvar !! availability of input data to
	3372	!! interpolate output variable
	3373	!! (on the nbpt space)
	3374	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
	3375	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in the input file
	3376	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
	3377	!! 'XYKindTime': Input values are kinds
	3378	!! of something with a temporal
	3379	!! evolution on the dx*dy matrix'
	3380	LOGICAL :: nonegative !! whether negative values should be removed
	3381	CHARACTER(LEN=50) :: maskingtype !! Type of masking
	3382	!! 'nomask': no-mask is applied
	3383	!! 'mbelow': take values below maskvals(1)
	3384	!! 'mabove': take values above maskvals(1)
	3385	!! 'msumrange': take values within 2 ranges;
	3386	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
	3387	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
	3388	!! (normalized by maskvals(3))
	3389	!! 'var': mask values are taken from a
	3390	!! variable inside the file (>0)
	3391	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
	3392	!! `maskingtype')
	3393	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
	3394	REAL(r_std), DIMENSION(1) :: variabletypevals !!
[5364]	3395	! REAL(r_std), DIMENSION(nbp_mpi) :: woodharvest_mpi !! Wood harvest where all thredds OMP are gatherd
[4657]	3396	!_ ================================================================================================================================
	3397
	3398
	3399	!Config Key = WOODHARVEST_FILE
[4677]	3400	!Config Desc = Name of file from which the wood harvest will be read
	3401	!Config If = DO_WOOD_HARVEST
	3402	!Config Def = woodharvest.nc
	3403	!Config Help =
[4657]	3404	!Config Units = [FILE]
	3405	filename = 'woodharvest.nc'
	3406	CALL getin_p('WOODHARVEST_FILE',filename)
	3407	variablename = 'woodharvest'
	3408
	3409
[5364]	3410	IF (xios_interpolation) THEN
	3411	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readwoodharvest: Use XIOS to read and interpolate " &
	3412	// TRIM(filename) // " for variable " // TRIM(variablename)
[4657]	3413
[5364]	3414	CALL xios_orchidee_recv_field('woodharvest_interp',woodharvest)
[4657]	3415
[5364]	3416	aoutvar = 1.0
	3417	ELSE
[4657]	3418
[5364]	3419	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readwoodharvest: Read and interpolate " &
	3420	// TRIM(filename) // " for variable " // TRIM(variablename)
[4657]	3421
[5364]	3422	! For this case there are not types/categories. We have 'only' a continuos field
	3423	! Assigning values to vmin, vmax
	3424	vmin = 0.
	3425	vmax = 9999.
	3426
	3427	!! Variables for interpweight
	3428	! Type of calculation of cell fractions
	3429	fractype = 'default'
	3430	! Name of the longitude and latitude in the input file
	3431	lonname = 'longitude'
	3432	latname = 'latitude'
	3433	! Should negative values be set to zero from input file?
	3434	nonegative = .TRUE.
	3435	! Type of mask to apply to the input data (see header for more details)
	3436	maskingtype = 'nomask'
	3437	! Values to use for the masking
	3438	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
	3439	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
	3440	namemaskvar = ''
	3441
	3442	variabletypevals=-un
	3443	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
	3444	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
	3445	maskvals, namemaskvar, -1, fractype, 0., 0., woodharvest, aoutvar)
	3446	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_wodharvest after interpweight_2Dcont'
	3447
	3448	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_woodharvest ended'
	3449	END IF
[4657]	3450	END SUBROUTINE slowproc_woodharvest
	3451
[721]	3452
	3453	!! ================================================================================================================================
[2419]	3454	!! SUBROUTINE : get_soilcorr_zobler
[721]	3455	!!
	3456	!>\BRIEF The "get_soilcorr" routine defines the table of correspondence
	3457	!! between the Zobler types and the three texture types known by SECHIBA and STOMATE :
	3458	!! silt, sand and clay.
	3459	!!
	3460	!! DESCRIPTION : get_soilcorr is needed if you use soils_param.nc .\n
	3461	!! The data from this file is then interpolated to the grid of the model. \n
	3462	!! The aim is to get fractions for sand loam and clay in each grid box.\n
	3463	!! This information is used for soil hydrology and respiration.
	3464	!!
	3465	!!
	3466	!! RECENT CHANGE(S): None
	3467	!!
	3468	!! MAIN OUTPUT VARIABLE(S) : ::texfrac_table
	3469	!!
	3470	!! REFERENCE(S) :
	3471	!! - Zobler L., 1986, A World Soil File for global climate modelling. NASA Technical memorandum 87802. NASA
	3472	!! Goddard Institute for Space Studies, New York, U.S.A.
	3473	!!
	3474	!! FLOWCHART : None
	3475	!! \n
	3476	!_ ================================================================================================================================
	3477
[2419]	3478	SUBROUTINE get_soilcorr_zobler (nzobler,textfrac_table)
[721]	3479
[733]	3480	IMPLICIT NONE
	3481
[721]	3482	!! 0. Variables and parameters declaration
	3483
	3484	INTEGER(i_std),PARAMETER :: nbtypes_zobler = 7 !! Number of Zobler types (unitless)
	3485
	3486	!! 0.1 Input variables
	3487
	3488	INTEGER(i_std),INTENT(in) :: nzobler !! Size of the array (unitless)
	3489
	3490	!! 0.2 Output variables
	3491
[2419]	3492	REAL(r_std),DIMENSION(nzobler,ntext),INTENT(out) :: textfrac_table !! Table of correspondence between soil texture class
	3493	!! and granulometric composition (0-1, unitless)
[721]	3494
	3495	!! 0.4 Local variables
	3496
	3497	INTEGER(i_std) :: ib !! Indice (unitless)
	3498
	3499	!_ ================================================================================================================================
	3500
[511]	3501	!-
[721]	3502	! 0. Check consistency
	3503	!-
	3504	IF (nzobler /= nbtypes_zobler) THEN
[1078]	3505	CALL ipslerr_p(3,'get_soilcorr', 'nzobler /= nbtypes_zobler',&
[721]	3506	& 'We do not have the correct number of classes', &
	3507	& ' in the code for the file.') ! Fatal error
[511]	3508	ENDIF
[721]	3509
[511]	3510	!-
[721]	3511	! 1. Textural fraction for : silt sand clay
[511]	3512	!-
	3513	textfrac_table(1,:) = (/ 0.12, 0.82, 0.06 /)
	3514	textfrac_table(2,:) = (/ 0.32, 0.58, 0.10 /)
	3515	textfrac_table(3,:) = (/ 0.39, 0.43, 0.18 /)
	3516	textfrac_table(4,:) = (/ 0.15, 0.58, 0.27 /)
	3517	textfrac_table(5,:) = (/ 0.34, 0.32, 0.34 /)
	3518	textfrac_table(6,:) = (/ 0.00, 1.00, 0.00 /)
	3519	textfrac_table(7,:) = (/ 0.39, 0.43, 0.18 /)
	3520
[721]	3521
	3522	!-
	3523	! 2. Check the mapping for the Zobler types which are going into the ORCHIDEE textures classes
	3524	!-
	3525	DO ib=1,nzobler ! Loop over # classes soil
	3526
	3527	IF (ABS(SUM(textfrac_table(ib,:))-1.0) > EPSILON(1.0)) THEN ! The sum of the textural fractions should not exceed 1 !
[511]	3528	WRITE(numout,*) &
	3529	& 'Error in the correspondence table', &
	3530	& ' sum is not equal to 1 in', ib
	3531	WRITE(numout,*) textfrac_table(ib,:)
[1078]	3532	CALL ipslerr_p(3,'get_soilcorr', 'SUM(textfrac_table(ib,:)) /= 1.0',&
[721]	3533	& '', 'Error in the correspondence table') ! Fatal error
[511]	3534	ENDIF
[721]	3535
[733]	3536	ENDDO ! Loop over # classes soil
	3537
	3538
[2419]	3539	END SUBROUTINE get_soilcorr_zobler
[511]	3540
[721]	3541	!! ================================================================================================================================
[2419]	3542	!! SUBROUTINE : get_soilcorr_usda
	3543	!!
	3544	!>\BRIEF The "get_soilcorr_usda" routine defines the table of correspondence
	3545	!! between the 12 USDA textural classes and their granulometric composition,
	3546	!! as % of silt, sand and clay. This is used to further defien clayfraction.
	3547	!!
	3548	!! DESCRIPTION : get_soilcorr is needed if you use soils_param.nc .\n
	3549	!! The data from this file is then interpolated to the grid of the model. \n
	3550	!! The aim is to get fractions for sand loam and clay in each grid box.\n
	3551	!! This information is used for soil hydrology and respiration.
	3552	!! The default map in this case is derived from Reynolds et al 2000, \n
	3553	!! at the 1/12deg resolution, with indices that are consistent with the \n
	3554	!! textures tabulated below
	3555	!!
	3556	!! RECENT CHANGE(S): Created by A. Ducharne on July 02, 2014
	3557	!!
	3558	!! MAIN OUTPUT VARIABLE(S) : ::texfrac_table
	3559	!!
	3560	!! REFERENCE(S) :
	3561	!!
	3562	!! FLOWCHART : None
	3563	!! \n
	3564	!_ ================================================================================================================================
	3565
	3566	SUBROUTINE get_soilcorr_usda (nusda,textfrac_table)
	3567
	3568	IMPLICIT NONE
	3569
	3570	!! 0. Variables and parameters declaration
	3571
	3572	!! 0.1 Input variables
	3573
	3574	INTEGER(i_std),INTENT(in) :: nusda !! Size of the array (unitless)
	3575
	3576	!! 0.2 Output variables
	3577
	3578	REAL(r_std),DIMENSION(nusda,ntext),INTENT(out) :: textfrac_table !! Table of correspondence between soil texture class
	3579	!! and granulometric composition (0-1, unitless)
	3580
	3581	!! 0.4 Local variables
	3582
	3583	INTEGER(i_std),PARAMETER :: nbtypes_usda = 12 !! Number of USDA texture classes (unitless)
	3584	INTEGER(i_std) :: n !! Index (unitless)
	3585
	3586	!_ ================================================================================================================================
	3587
	3588	!-
	3589	! 0. Check consistency
	3590	!-
	3591	IF (nusda /= nbtypes_usda) THEN
	3592	CALL ipslerr_p(3,'get_soilcorr', 'nusda /= nbtypes_usda',&
	3593	& 'We do not have the correct number of classes', &
	3594	& ' in the code for the file.') ! Fatal error
	3595	ENDIF
	3596
	3597	!! Parameters for soil type distribution :
	3598	!! Sand, Loamy Sand, Sandy Loam, Silt Loam, Silt, Loam, Sandy Clay Loam, Silty Clay Loam, Clay Loam, Sandy Clay, Silty Clay, Clay
	3599	! The order comes from constantes_soil.f90
	3600	! The corresponding granulometric composition comes from Carsel & Parrish, 1988
	3601
	3602	!-
	3603	! 1. Textural fractions for : sand, clay
	3604	!-
	3605	textfrac_table(1,2:3) = (/ 0.93, 0.03 /) ! Sand
	3606	textfrac_table(2,2:3) = (/ 0.81, 0.06 /) ! Loamy Sand
	3607	textfrac_table(3,2:3) = (/ 0.63, 0.11 /) ! Sandy Loam
	3608	textfrac_table(4,2:3) = (/ 0.17, 0.19 /) ! Silt Loam
	3609	textfrac_table(5,2:3) = (/ 0.06, 0.10 /) ! Silt
	3610	textfrac_table(6,2:3) = (/ 0.40, 0.20 /) ! Loam
	3611	textfrac_table(7,2:3) = (/ 0.54, 0.27 /) ! Sandy Clay Loam
	3612	textfrac_table(8,2:3) = (/ 0.08, 0.33 /) ! Silty Clay Loam
	3613	textfrac_table(9,2:3) = (/ 0.30, 0.33 /) ! Clay Loam
	3614	textfrac_table(10,2:3) = (/ 0.48, 0.41 /) ! Sandy Clay
	3615	textfrac_table(11,2:3) = (/ 0.06, 0.46 /) ! Silty Clay
	3616	textfrac_table(12,2:3) = (/ 0.15, 0.55 /) ! Clay
	3617
	3618	! Fraction of silt
	3619
	3620	DO n=1,nusda
	3621	textfrac_table(n,1) = 1. - textfrac_table(n,2) - textfrac_table(n,3)
	3622	END DO
	3623
	3624	END SUBROUTINE get_soilcorr_usda
	3625
	3626	!! ================================================================================================================================
[721]	3627	!! FUNCTION : tempfunc
	3628	!!
	3629	!>\BRIEF ! This function interpolates value between ztempmin and ztempmax
	3630	!! used for lai detection.
	3631	!!
	3632	!! DESCRIPTION : This subroutine calculates a scalar between 0 and 1 with the following equation :\n
	3633	!! \latexonly
	3634	!! \input{constantes_veg_tempfunc.tex}
	3635	!! \endlatexonly
	3636	!!
	3637	!! RECENT CHANGE(S): None
	3638	!!
	3639	!! RETURN VALUE : tempfunc_result
	3640	!!
	3641	!! REFERENCE(S) : None
	3642	!!
	3643	!! FLOWCHART : None
	3644	!! \n
	3645	!_ ================================================================================================================================
[511]	3646
	3647	FUNCTION tempfunc (temp_in) RESULT (tempfunc_result)
[721]	3648
	3649
	3650	!! 0. Variables and parameters declaration
	3651
	3652	REAL(r_std),PARAMETER :: ztempmin=273._r_std !! Temperature for laimin (K)
	3653	REAL(r_std),PARAMETER :: ztempmax=293._r_std !! Temperature for laimax (K)
	3654	REAL(r_std) :: zfacteur !! Interpolation factor (K^{-2})
	3655
	3656	!! 0.1 Input variables
	3657
	3658	REAL(r_std),INTENT(in) :: temp_in !! Temperature (K)
	3659
[733]	3660	!! 0.2 Result
[721]	3661
	3662	REAL(r_std) :: tempfunc_result !! (unitless)
	3663
	3664	!_ ================================================================================================================================
	3665
[733]	3666	!! 1. Define a coefficient
[511]	3667	zfacteur = un/(ztempmax-ztempmin)**2
[721]	3668
[733]	3669	!! 2. Computes tempfunc
[511]	3670	IF (temp_in > ztempmax) THEN
	3671	tempfunc_result = un
	3672	ELSEIF (temp_in < ztempmin) THEN
	3673	tempfunc_result = zero
	3674	ELSE
	3675	tempfunc_result = un-zfacteur(ztempmax-temp_in)*2
[721]	3676	ENDIF !(temp_in > ztempmax)
[733]	3677
	3678
[511]	3679	END FUNCTION tempfunc
	3680
	3681
[1136]	3682	!! ================================================================================================================================
	3683	!! SUBROUTINE : slowproc_checkveget
	3684	!!
	3685	!>\BRIEF To verify the consistency of the various fractions defined within the grid box after having been
	3686	!! been updated by STOMATE or the standard procedures.
	3687	!!
	3688	!! DESCRIPTION : (definitions, functional, design, flags):
	3689	!!
	3690	!! RECENT CHANGE(S): None
	3691	!!
	3692	!! MAIN OUTPUT VARIABLE(S): :: none
	3693	!!
	3694	!! REFERENCE(S) : None
	3695	!!
	3696	!! FLOWCHART : None
	3697	!! \n
	3698	!_ ================================================================================================================================
	3699	!
[2718]	3700	SUBROUTINE slowproc_checkveget(nbpt, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
	3701
[1136]	3702	! 0.1 INPUT
	3703	!
	3704	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
	3705	REAL(r_std),DIMENSION (nbpt,nnobio), INTENT(in) :: frac_nobio ! Fraction of ice,lakes,cities, ... (unitless)
	3706	REAL(r_std),DIMENSION (nbpt,nvm), INTENT(in) :: veget_max ! Maximum fraction of vegetation type including none biological fraction (unitless)
	3707	REAL(r_std),DIMENSION (nbpt,nvm), INTENT(in) :: veget ! Vegetation fractions
[3969]	3708	REAL(r_std),DIMENSION (nbpt), INTENT(in) :: tot_bare_soil ! Total evaporating bare soil fraction within the mesh
[1136]	3709	REAL(r_std),DIMENSION (nbpt,nstm), INTENT(in) :: soiltile ! Fraction of soil tiles in the gridbox (unitless)
	3710
	3711	! 0.3 LOCAL
	3712	!
[3831]	3713	INTEGER(i_std) :: ji, jn, jv
[2718]	3714	REAL(r_std) :: epsilocal !! A very small value
[1136]	3715	REAL(r_std) :: totfrac
	3716	CHARACTER(len=80) :: str1, str2
[2718]	3717
[1136]	3718	!_ ================================================================================================================================
[2718]	3719
[1136]	3720	!
	3721	! There is some margin added as the computing errors might bring us above EPSILON(un)
	3722	!
	3723	epsilocal = EPSILON(un)*1000.
[2718]	3724
	3725	!! 1.0 Verify that none of the fractions are smaller than min_vegfrac, without beeing zero.
	3726	!!
[1136]	3727	DO ji=1,nbpt
	3728	DO jn=1,nnobio
	3729	IF ( frac_nobio(ji,jn) > epsilocal .AND. frac_nobio(ji,jn) < min_vegfrac ) THEN
	3730	WRITE(str1,'("Occurs on grid box", I8," and nobio type ",I3 )') ji, jn
	3731	WRITE(str2,'("The small value obtained is ", E14.4)') frac_nobio(ji,jn)
	3732	CALL ipslerr_p (3,'slowproc_checkveget', &
[2668]	3733	"frac_nobio is larger than zero but smaller than min_vegfrac.", str1, str2)
[1136]	3734	ENDIF
	3735	ENDDO
[2668]	3736	END DO
	3737
	3738	IF (.NOT. ok_dgvm) THEN
	3739	DO ji=1,nbpt
	3740	DO jv=1,nvm
	3741	IF ( veget_max(ji,jv) > epsilocal .AND. veget_max(ji,jv) < min_vegfrac ) THEN
	3742	WRITE(str1,'("Occurs on grid box", I8," and nobio type ",I3 )') ji, jn
	3743	WRITE(str2,'("The small value obtained is ", E14.4)') veget_max(ji,jv)
	3744	CALL ipslerr_p (3,'slowproc_checkveget', &
	3745	"veget_max is larger than zero but smaller than min_vegfrac.", str1, str2)
	3746	ENDIF
	3747	ENDDO
[1136]	3748	ENDDO
[2668]	3749	END IF
[2718]	3750
	3751	!! 2.0 verify that with all the fractions we cover the entire grid box
	3752	!!
[1136]	3753	DO ji=1,nbpt
	3754	totfrac = zero
	3755	DO jn=1,nnobio
	3756	totfrac = totfrac + frac_nobio(ji,jn)
	3757	ENDDO
	3758	DO jv=1,nvm
	3759	totfrac = totfrac + veget_max(ji,jv)
	3760	ENDDO
	3761	IF ( ABS(totfrac - un) > epsilocal) THEN
	3762	WRITE(str1,'("This occurs on grid box", I8)') ji
	3763	WRITE(str2,'("The sum over all fraction and error are ", E14.4, E14.4)') totfrac, ABS(totfrac - un)
	3764	CALL ipslerr_p (3,'slowproc_checkveget', &
	3765	"veget_max + frac_nobio is not equal to 1.", str1, str2)
	3766	WRITE(,) "EPSILON =", epsilocal
	3767	ENDIF
	3768	ENDDO
[2718]	3769
	3770	!! 3.0 Verify that veget is smaller or equal to veget_max
	3771	!!
[1136]	3772	DO ji=1,nbpt
	3773	DO jv=1,nvm
[2962]	3774	IF ( jv == ibare_sechiba ) THEN
[1136]	3775	IF ( ABS(veget(ji,jv) - veget_max(ji,jv)) > epsilocal ) THEN
	3776	WRITE(str1,'("This occurs on grid box", I8)') ji
	3777	WRITE(str2,'("The difference is ", E14.4)') veget(ji,jv) - veget_max(ji,jv)
	3778	CALL ipslerr_p (3,'slowproc_checkveget', &
	3779	"veget is not equal to veget_max on bare soil.", str1, str2)
	3780	ENDIF
	3781	ELSE
	3782	IF ( veget(ji,jv) > veget_max(ji,jv) ) THEN
	3783	WRITE(str1,'("This occurs on grid box", I8)') ji
	3784	WRITE(str2,'("The values for veget and veget_max :", F8.4, F8.4)') veget(ji,jv), veget_max(ji,jv)
	3785	CALL ipslerr_p (3,'slowproc_checkveget', &
	3786	"veget is greater than veget_max.", str1, str2)
	3787	ENDIF
	3788	ENDIF
	3789	ENDDO
	3790	ENDDO
[2718]	3791
	3792	!! 4.0 Test tot_bare_soil in relation to the other variables
	3793	!!
[1136]	3794	DO ji=1,nbpt
	3795	totfrac = zero
	3796	DO jv=1,nvm
	3797	totfrac = totfrac + (veget_max(ji,jv) - veget(ji,jv))
	3798	ENDDO
	3799	! add the bare soil fraction to totfrac
[2962]	3800	totfrac = totfrac + veget(ji,ibare_sechiba)
[1136]	3801	! do the test
	3802	IF ( ABS(totfrac - tot_bare_soil(ji)) > epsilocal ) THEN
	3803	WRITE(str1,'("This occurs on grid box", I8)') ji
	3804	WRITE(str2,'("The values for tot_bare_soil, tot frac and error :", F8.4, F8.4, E14.4)') &
	3805	& tot_bare_soil(ji), totfrac, ABS(totfrac - tot_bare_soil(ji))
	3806	CALL ipslerr_p (3,'slowproc_checkveget', &
	3807	"tot_bare_soil does not correspond to the total bare soil fraction.", str1, str2)
	3808	ENDIF
	3809	ENDDO
[2718]	3810
	3811	!! 5.0 Test that soiltile has the right sum
	3812	!!
[1136]	3813	DO ji=1,nbpt
[3593]	3814	totfrac = SUM(soiltile(ji,:))
[1136]	3815	IF ( ABS(totfrac - un) > epsilocal ) THEN
[3593]	3816	WRITE(numout,*) "soiltile does not sum-up to one. This occurs on grid box", ji
	3817	WRITE(numout,*) "The soiltile for ji are :", soiltile(ji,:)
	3818	CALL ipslerr_p (2,'slowproc_checkveget', &
	3819	"soiltile does not sum-up to one.", "", "")
[1136]	3820	ENDIF
	3821	ENDDO
[2718]	3822
[1136]	3823	END SUBROUTINE slowproc_checkveget
[3094]	3824
	3825
	3826	!! ================================================================================================================================
	3827	!! SUBROUTINE : slowproc_change_frac
	3828	!!
	3829	!>\BRIEF Update the vegetation fractions
	3830	!!
	3831	!! DESCRIPTION : Update the vegetation fractions. This subroutine is called in the same time step as lcchange in stomatelpj has
	3832	!! has been done. This subroutine is called after the diagnostics have been written in sechiba_main.
	3833	!!
	3834	!! RECENT CHANGE(S): None
	3835	!!
	3836	!! MAIN OUTPUT VARIABLE(S): :: veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile
	3837	!!
	3838	!! REFERENCE(S) : None
	3839	!!
	3840	!! FLOWCHART : None
	3841	!! \n
	3842	!_ ================================================================================================================================
[3831]	3843
[3094]	3844	SUBROUTINE slowproc_change_frac(kjpindex, lai, &
[4825]	3845	veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile, fraclut, nwdFraclut)
[3094]	3846	!
	3847	! 0. Declarations
	3848	!
[3831]	3849	! 0.1 Input variables
	3850	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
	3851	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m^2 m^{-2})
	3852
	3853	! 0.2 Output variables
[3969]	3854	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: veget_max !! Maximum fraction of vegetation type in the mesh (unitless)
	3855	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: veget !! Fraction of vegetation type in the mesh (unitless)
[3094]	3856	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
	3857	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh
[3969]	3858	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_bare_soil !! Total evaporating bare soil fraction in the mesh
	3859	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
[4723]	3860	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
[4825]	3861	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: nwdfraclut !! Fraction of non woody vegetation in each landuse tile (0-1, unitless)
[3831]	3862
[3094]	3863	! 0.3 Local variables
	3864	INTEGER(i_std) :: ji, jv !! Loop index
[3165]	3865
[3831]	3866
[3165]	3867	!! Update vegetation fractions with the values coming from the vegetation file read in slowproc_readvegetmax.
[3831]	3868	!! Partial update has been taken into account for the case with DGVM and AGRICULTURE in slowproc_readvegetmax.
[3165]	3869	veget_max = veget_max_new
	3870	frac_nobio = frac_nobio_new
[3831]	3871
[3094]	3872	!! Verification and correction on veget_max, calculation of veget and soiltile.
[4825]	3873	CALL slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
[3831]	3874
[3969]	3875	!! Calculate tot_bare_soil needed in hydrol, diffuco and condveg (fraction of bare soil in the mesh)
[3094]	3876	tot_bare_soil(:) = veget_max(:,1)
	3877	DO jv = 2, nvm
	3878	DO ji =1, kjpindex
	3879	tot_bare_soil(ji) = tot_bare_soil(ji) + (veget_max(ji,jv) - veget(ji,jv))
	3880	ENDDO
	3881	END DO
[3831]	3882
	3883	!! Do some basic tests on the surface fractions updated above
[3094]	3884	CALL slowproc_checkveget(kjpindex, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
[3831]	3885
	3886	END SUBROUTINE slowproc_change_frac
[1136]	3887
[8]	3888	END MODULE slowproc

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

Download in other formats: