Context Navigation

slowproc.f90 @ 8579

Last change on this file since 8579 was 8579, checked in by antoine.bierjon, 3 weeks ago
Adapted read and interpolation of irrigation input files with XIOS_interpolation. Added a default value for AEI_SW interpolation (in case of missing value, with and without XIOS), which is chosen to be the average of the GMIA5 dataset from FAO. Related to ticket #857
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 248.9 KB

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

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

Context Navigation

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

Download in other formats: