Context Navigation

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

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

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

Download in other formats: