Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90 @ 7265

Last change on this file since 7265 was 7207, checked in by agnes.ducharne, 3 years ago
Bugfix of r7206 (the checks were wrong as compilation has been forgotten! Now, with bugfix, compilation is ok, and the 5d simulation is identical to the one with r7200).
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 100.8 KB

Rev	Line
[947]	1	! =================================================================================================================================
	2	! MODULE : thermosoil
	3	!
[4470]	4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
[947]	5	!
	6	! LICENCE : IPSL (2006)
	7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
	8	!
	9	!>\BRIEF Calculates the soil temperatures by solving the heat
[2892]	10	!! diffusion equation within the soil. This module is only used with CWRR hydrology.
[8]	11	!!
[947]	12	!!\n DESCRIPTION : General important informations about the numerical scheme and
	13	!! the soil vertical discretization:\n
[2956]	14	!! - the soil is zmaxt deep (by default 10m) and divided into "ngrnd" layers.
	15	!! From 0-zmaxh(default 2m), the discretization is the same as for hydrology.
	16	!! From zmaxh(2m) and below, the depth increase linearly (by default) or geometrically. \n
[947]	17	!! - "jg" is usually used as the index going from 1 to ngrnd to describe the
	18	!! layers, from top (jg=1) to bottom (jg=ngrnd)\n
	19	!! - the thermal numerical scheme is implicit finite differences.\n
	20	!! -- When it is resolved in thermosoil_profile at the present timestep t, the
	21	!! dependancy from the previous timestep (t-1) is hidden in the
	22	!! integration coefficients cgrnd and dgrnd, which are therefore
	23	!! calculated at the very end of thermosoil_main (call to
	24	!! thermosoil_coef) for use in the next timestep.\n
	25	!! -- At timestep t, the system becomes :\n
[8]	26	!!
[947]	27	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n
	28	!! -- EQ1 -- \n
	29	!!
	30	!! (the bottom boundary condition has been used to obtained this equation).\n
	31	!! To solve it, the uppermost soil temperature T(1) is required.
	32	!! It is obtained from the surface temperature Ts, which is
	33	!! considered a linear extrapolation of T(1) and T(2)\n
	34	!!
[2956]	35	!! Ts=(1+lambda)T(1) -lambdaT(2) \n
[947]	36	!! -- EQ2--\n
	37	!!
	38	!! -- caveat 1 : Ts is called 'temp_soil_new' in this routine,
	39	!! don' t act.\n
	40	!! -- caveat 2 : actually, the surface temperature at time t Ts
	41	!! depends on the soil temperature at time t through the
	42	!! ground heat flux. This is again implicitly solved, with Ts(t)
	43	!! expressed as :\n
	44	!!
[2956]	45	!! soilcap*(Ts(t)-Ts(t-1))/dt=soilflx+otherfluxes(Ts(t))\n
[947]	46	!! -- EQ3 --\n
	47	!!
	48	!! and the dependency from the previous timestep is hidden in
[2956]	49	!! soilcap and soilflx (apparent surface heat capacity and heat
	50	!! flux respectively). Soilcap and soilflx are therefore
	51	!! calculated at the previous timestep, at the very end of thermosoil
[947]	52	!! (final call to thermosoil_coef) and stored to be used at the next time step.
	53	!! At timestep t, EQ3 is solved for Ts in enerbil, and Ts
	54	!! is used in thermosoil to get T(1) and solve EQ1.\n
	55	!!
	56	!! - lambda is the @tex $\mu$ @endtex of F. Hourdin' s PhD thesis, equation (A28); ie the
	57	!! coefficient of the linear extrapolation of Ts (surface temperature) from T1 and T2 (ptn(jg=1) and ptn(jg=2)), so that:\n
	58	!! Ts= (1+lambda)T(1)-lambdaT(2) --EQ2-- \n
[2942]	59	!! lambda = (zlt(1))/((zlt(2)-zlt(1))) \n
[947]	60	!!
[2956]	61	!! RECENT CHANGE(S) : - Change soil thermal properties to consider also soil texture, rev 2922.
	62	!! - Change vertical discretization, rev 2917. Note: In the revised thermosoil,
[2942]	63	!! cstgrnd and lskin are not needed any more. The depth znt, zlt and dlt
[2956]	64	!! are computed in vertical_soil and are in meter
[947]	65	!!
	66	!! REFERENCE(S) : None
	67	!!
	68	!! SVN :
	69	!! $HeadURL$
	70	!! $Date$
	71	!! $Revision$
	72	!! \n
	73	!_ ================================================================================================================================
	74
[8]	75	MODULE thermosoil
	76
[1392]	77	! modules used :
[8]	78	USE ioipsl
[1392]	79	USE ioipsl_para
[1788]	80	USE xios_orchidee
[8]	81	USE constantes
[4646]	82	USE time, ONLY : one_day, dt_sechiba
[947]	83	USE constantes_soil
[4281]	84	USE sechiba_io_p
[8]	85	USE grid
	86
	87	IMPLICIT NONE
	88
[947]	89	!private and public routines :
[8]	90	PRIVATE
[5364]	91	PUBLIC :: thermosoil_main, thermosoil_clear, thermosoil_initialize, thermosoil_finalize, thermosoil_xios_initialize
[8]	92
[2917]	93	REAL(r_std), SAVE :: lambda !! See Module description
	94	!$OMP THREADPRIVATE(lambda)
[947]	95	REAL(r_std), SAVE :: fz1, zalph !! usefull constants for diverse use
[1078]	96	!$OMP THREADPRIVATE(fz1, zalph)
[947]	97	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn !! vertically discretized
	98	!! soil temperatures @tex ($K$) @endtex.
[1078]	99	!$OMP THREADPRIVATE(ptn)
[2001]	100	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz1 !! numerical constant used in the thermal numerical
[947]	101	!! scheme @tex ($m^{-1}$) @endtex. ; it corresponds
	102	!! to the coefficient @tex $d_k$ @endtex of equation
	103	!! (A.12) in F. Hourdin PhD thesis.
[1078]	104	!$OMP THREADPRIVATE(dz1)
[947]	105	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !! integration coefficient for the numerical scheme,
	106	!! see eq.1
[1078]	107	!$OMP THREADPRIVATE(cgrnd)
[947]	108	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !! integration coefficient for the numerical scheme,
	109	!! see eq.1
[1078]	110	!$OMP THREADPRIVATE(dgrnd)
[3059]	111
[947]	112	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !! volumetric vertically discretized soil heat
	113	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
	114	!! It depends on the soil
[2222]	115	!! moisture content (shum_ngrnd_perma) and is calculated at
[947]	116	!! each time step in thermosoil_coef.
[1078]	117	!$OMP THREADPRIVATE(pcapa)
[947]	118	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !! vertically discretized soil thermal conductivity
	119	!! @tex ($W K^{-1} m^{-1}$) @endtex. Same as pcapa.
[1078]	120	!$OMP THREADPRIVATE(pkappa)
[3059]	121	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_snow !! volumetric vertically discretized snow heat
	122	!! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
	123	!$OMP THREADPRIVATE(pcapa_snow)
	124	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa_snow !! vertically discretized snow thermal conductivity
	125	!! @tex ($W K^{-1} m^{-1}$) @endtex.
	126	!$OMP THREADPRIVATE(pkappa_snow)
	127
[947]	128	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! heat capacity used for surfheat_incr and
	129	!! coldcont_incr
[1078]	130	!$OMP THREADPRIVATE(pcapa_en)
[4682]	131	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg !! ptn as it is after thermosoil_profile but before thermosoil_coef,
	132	!! used in thermosoil_readjust
[1078]	133	!$OMP THREADPRIVATE(ptn_beg)
[4682]	134	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at previous timestep (K)
[1078]	135	!$OMP THREADPRIVATE(temp_sol_beg)
[947]	136	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat content during the timestep
	137	!! @tex ($J$) @endtex.
[1078]	138	!$OMP THREADPRIVATE(surfheat_incr)
[947]	139	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow heat content @tex ($J$) @endtex.
[1078]	140	!$OMP THREADPRIVATE(coldcont_incr)
[2222]	141	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shum_ngrnd_perma !! Saturation degree on the thermal axes (0-1, dimensionless)
	142	!$OMP THREADPRIVATE(shum_ngrnd_perma)
	143
	144	! Variables related to soil freezing
	145	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: profil_froz !! Frozen fraction of the soil on hydrological levels (-)
	146	!$OMP THREADPRIVATE(profil_froz)
[7199]	147	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat !! Latent heat released or consumed in the freezing/thawing processes summed vertically
	148	!! for the whole soil (J/m2) and on the whole simulation to check/correct energy conservation
[2222]	149	!$OMP THREADPRIVATE(e_soil_lat)
[7199]	150	REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp !! Additional surfacic heat capacity due to soil freezing for each soil layer (J/K/m2)
[2222]	151	!$OMP THREADPRIVATE(pcappa_supp)
[2922]	152	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz5 !! Used for numerical calculation [-]
	153	!$OMP THREADPRIVATE(dz5)
	154	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: QZ !! quartz content [-]
	155	!$OMP THREADPRIVATE(QZ)
	156	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_dry_ns !! Dry soil Heat capacity of soils,J.m^{-3}.K^{-1}
	157	!$OMP THREADPRIVATE(so_capa_dry_ns)
[7199]	158	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_ice !! Heat capacity of saturated frozen soil (J/K/m3)
	159	!$OMP THREADPRIVATE(so_capa_ice)
[2922]	160	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_layt !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface
	161	!$OMP THREADPRIVATE(mc_layt)
	162	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mcl_layt !! Volumetric soil moisture (liquid) (m3/m3) on the thermodynamical levels at interface
	163	!$OMP THREADPRIVATE(mcl_layt)
	164	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_layt !! Total soil moisture content for each layer (liquid+ice) (mm) on the thermodynamical levels
	165	!$OMP THREADPRIVATE(tmc_layt)
	166	INTEGER(i_std), SAVE :: brk_flag = 0 !! Flag to consider bedrock: 0.no; 1.yes
	167	!$OMP THREADPRIVATE(brk_flag)
[2222]	168
	169
[947]	170	CONTAINS
[5364]	171
	172
	173	!! =============================================================================================================================
	174	!! SUBROUTINE: thermosoil_xios_initialize
	175	!!
	176	!>\BRIEF Initialize xios dependant defintion before closing context defintion
	177	!!
	178	!! DESCRIPTION: Initialize xios dependant defintion before closing context defintion
	179	!! Reading is deactivated if the sechiba restart file exists because the
	180	!! variable should be in the restart file already.
	181	!!
	182	!! \n
	183	!_ ==============================================================================================================================
	184
	185	SUBROUTINE thermosoil_xios_initialize
	186
	187	CHARACTER(LEN=255) :: filename, name
	188
	189	filename = 'reftemp.nc'
	190	CALL getin_p('REFTEMP_FILE',filename)
	191
	192	name = filename(1:LEN_TRIM(FILENAME)-3)
	193	CALL xios_orchidee_set_file_attr("reftemp_file",name=name)
	194
	195	! Check if the reftemp file will be read by XIOS, by IOIPSL or not at all
	196	IF (xios_interpolation .AND. read_reftemp .AND. restname_in=='NONE') THEN
	197	! The reftemp file will be read using XIOS
	198	IF (printlev>=2) WRITE(numout,*) 'Reading of reftemp file will be done later using XIOS. The filename is ', filename
	199	ELSE
	200	IF (.NOT. read_reftemp) THEN
	201	IF (printlev>=2) WRITE (numout,*) 'No reading of reftemp will be done because read_reftemp=FALSE'
	202	ELSE IF (restname_in=='NONE') THEN
	203	IF (printlev>=2) WRITE (numout,*) 'The reftemp file will be read later by IOIPSL'
	204	ELSE
	205	IF (printlev>=2) WRITE (numout,*) 'The reftemp file will not be read because the restart file exists.'
	206	END IF
	207
	208	! The reftemp file will not be read by XIOS. Now deactivate albedo for XIOS.
	209	CALL xios_orchidee_set_file_attr("reftemp_file",enabled=.FALSE.)
	210	CALL xios_orchidee_set_field_attr("reftemp_interp",enabled=.FALSE.)
	211	ENDIF
	212
	213	END SUBROUTINE thermosoil_xios_initialize
	214
[2581]	215	!! =============================================================================================================================
	216	!! SUBROUTINE : thermosoil_initialize
	217	!!
	218	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
	219	!!
	220	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values.
	221	!! Call thermosoil_var_init to calculate physical constants.
	222	!! Call thermosoil_coef to calculate thermal soil properties.
	223	!!
	224	!! RECENT CHANGE(S) : None
	225	!!
	226	!! REFERENCE(S) : None
	227	!!
	228	!! FLOWCHART : None
	229	!! \n
	230	!_ ==============================================================================================================================
[7206]	231	SUBROUTINE thermosoil_initialize (kjit, kjpindex, rest_id, mcs, &
[2650]	232	temp_sol_new, snow, shumdiag_perma, &
	233	soilcap, soilflx, stempdiag, &
[3059]	234	gtemp, &
[3020]	235	mc_layh, mcl_layh, tmc_layh, njsc, &
[3269]	236	frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
	237	snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb)
[8]	238
[2581]	239	!! 0. Variable and parameter declaration
	240	!! 0.1 Input variables
	241	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
	242	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
	243	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (unitless)
[7206]	244	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3)
[2581]	245	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step,
	246	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
[2942]	247	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree (0-1, unitless)
	248	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content (liquid+ice) for hydrological layers, at node (m3/m3)
	249	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content (liquid) for hydrological layers, at node (m3/m3)
	250	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content(liquid+ice) for hydrological layers (mm)
[2922]	251	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
[3020]	252	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
	253	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
	254	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
	255	!! (unitless,0-1)
[3059]	256	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth
[3269]	257	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
[3410]	258	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
[3269]	259	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
[2581]	260
	261	!! 0.2 Output variables
[3020]	262	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilcap !! apparent surface heat capacity considering snow and soil surface (J m-2 K-1)
	263	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilflx !! apparent soil heat flux considering snow and soil surface (W m-2)
[2942]	264	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile on the levels in hydrol(K)
[2650]	265	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
[2581]	266
[3269]	267	!! 0.3 Modified variables
	268	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
	269	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
	270	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
[2650]	271
[2581]	272	!! 0.4 Local variables
[7199]	273	INTEGER(i_std) :: ier, i, jg, jsc
[2868]	274	LOGICAL :: calculate_coef !! Local flag to initialize variables by call to thermosoil_coef
[2581]	275	!_ ================================================================================================================================
	276
[2922]	277
	278	!
	279	! !! Flag to consider bedrock at deeper layers
	280	! !! It affects heat capacity and thermal conductivity (energy balance).
	281	!
	282	!Config Key = BEDROCK_FLAG
[2956]	283	!Config Desc = Flag to consider bedrock at deeper layers.
[2922]	284	!Config If =
	285	!Config Def = 0
	286	!Config Help = 0, no, 1, yes.
	287	!Config Units = [FLAG]
	288	brk_flag = 0
	289	CALL getin_p('BEDROCK_FLAG', brk_flag)
	290
[2581]	291	IF (printlev >= 3) WRITE (numout,*) 'Start thermosoil_initialize '
	292
	293	!! 1. Allocate soil temperatures variables
	294	ALLOCATE (ptn(kjpindex,ngrnd),stat=ier)
	295	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn','','')
	296
	297	ALLOCATE (dz1(ngrnd),stat=ier)
	298	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz1','','')
	299
	300	ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier)
	301	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of cgrnd','','')
	302
	303	ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier)
	304	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dgrnd','','')
	305
	306	ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier)
	307	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa','','')
	308
	309	ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier)
	310	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa','','')
	311
[3059]	312	ALLOCATE (pcapa_snow(kjpindex,nsnow),stat=ier)
	313	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_snow','','')
	314
	315	ALLOCATE (pkappa_snow(kjpindex,nsnow),stat=ier)
	316	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa_snow','','')
	317
[3411]	318	! Temporary fix: Initialize following variable because they are output to xios before the first calculation
	319	pcapa = 0
	320	pkappa = 0
	321	pcapa_snow = 0
	322	pkappa_snow = 0
	323
[2581]	324	ALLOCATE (surfheat_incr(kjpindex),stat=ier)
	325	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of surfheat_incr','','')
	326
	327	ALLOCATE (coldcont_incr(kjpindex),stat=ier)
	328	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of coldcont_incr','','')
	329
	330	ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier)
	331	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_en','','')
[5293]	332	! Initialization to zero used at first time step in thermosoil_energy_diag, only for diagnostic variables coldcont_incr and surfheat_incr
	333	pcapa_en(:,:) = 0.
[2581]	334
	335	ALLOCATE (ptn_beg(kjpindex,ngrnd),stat=ier)
	336	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn_beg','','')
	337
	338	ALLOCATE (temp_sol_beg(kjpindex),stat=ier)
	339	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of temp_sol_beg','','')
	340
	341	ALLOCATE (shum_ngrnd_perma(kjpindex,ngrnd),stat=ier)
	342	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of shum_ngrnd_perma','','')
	343
	344	ALLOCATE (profil_froz(kjpindex,ngrnd),stat=ier)
	345	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of profil_froz','','')
	346
	347	IF (ok_freeze_thermix) THEN
	348	ALLOCATE (pcappa_supp(kjpindex,ngrnd),stat=ier)
	349	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ok_freeze_termix','','')
[5433]	350	! Initialization to zero used at first time step only for diagnostic output.
	351	! This variable is only used in thermosoil_readajust and always calculated before in thermosoil_getdiff.
	352	pcappa_supp(:,:) = 0.
[2581]	353	END IF
	354	IF (ok_Ecorr) THEN
	355	ALLOCATE (e_soil_lat(kjpindex),stat=ier)
	356	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of e_soil_lat','','')
	357	END IF
	358
[2922]	359	ALLOCATE (dz5(ngrnd),stat=ier)
	360	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz5','','')
	361
	362	ALLOCATE (mc_layt(kjpindex,ngrnd),stat=ier)
	363	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mc_layt','','')
	364
	365	ALLOCATE (mcl_layt(kjpindex,ngrnd),stat=ier)
	366	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcl_layt','','')
	367
	368	ALLOCATE (tmc_layt(kjpindex,ngrnd),stat=ier)
	369	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of tmc_layt','','')
	370
	371	ALLOCATE (QZ(nscm),stat=ier)
	372	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of QZ','','')
	373
	374	ALLOCATE (so_capa_dry_ns(nscm),stat=ier)
	375	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_dry_ns','','')
[7199]	376
[7206]	377	ALLOCATE (so_capa_ice(kjpindex),stat=ier)
[7199]	378	IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_ice','','')
[2922]	379
[3059]	380
[2922]	381	!! Soil texture choose
	382	SELECTCASE (nscm)
	383	CASE (3)
	384	QZ(:) = QZ_fao(:)
	385	so_capa_dry_ns(:) = so_capa_dry_ns_fao(:)
[6954]	386	CASE (13) !Salma changed from 12 to 13 for the new class Oxisols
[7199]	387	QZ(:) = QZ_usda(:)
[2922]	388	so_capa_dry_ns(:) = so_capa_dry_ns_usda(:)
	389	CASE DEFAULT
	390	WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler, fao and usda according to the map'
	391	STOP 'thermosoil_initialize'
	392	ENDSELECT
	393
[2581]	394
	395	!! 2. Initialize variable from restart file or with default values
	396
	397	!! Reads restart files for soil temperatures only. If no restart file is
	398	!! found, the initial soil temperature is by default set to 280K at all depths. The user
	399	!! can decide to initialize soil temperatures at an other value, in which case he should set the flag THERMOSOIL_TPRO
	400	!! to this specific value in the run.def.
	401	IF (printlev>=3) WRITE (numout,*) 'Read restart file for THERMOSOIL variables'
	402
	403	CALL ioconf_setatt_p('UNITS', 'K')
	404	CALL ioconf_setatt_p('LONG_NAME','Soil Temperature profile')
	405	CALL restget_p (rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, .TRUE., ptn, "gather", nbp_glo, index_g)
[3345]	406
	407	! Initialize ptn if it was not found in restart file
	408	IF (ALL(ptn(:,:)==val_exp)) THEN
	409	! ptn was not found in restart file
	410
	411	IF (read_reftemp) THEN
	412	! Read variable ptn from file
[3850]	413	CALL thermosoil_read_reftempfile(kjpindex,lalo,ptn)
[3345]	414	ELSE
	415	! Initialize ptn with a constant value which can be set in run.def
	416
	417	!Config Key = THERMOSOIL_TPRO
	418	!Config Desc = Initial soil temperature profile if not found in restart
	419	!Config Def = 280.
	420	!Config If = OK_SECHIBA
	421	!Config Help = The initial value of the temperature profile in the soil if
	422	!Config its value is not found in the restart file. Here
	423	!Config we only require one value as we will assume a constant
	424	!Config throughout the column.
	425	!Config Units = Kelvin [K]
	426	CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std)
	427	END IF
[2581]	428	END IF
	429
[4682]	430	! Initialize ptn_beg (variable needed in thermosoil_readadjust called from thermosoil_coef)
[2581]	431	ptn_beg(:,:) = ptn(:,:)
	432
	433	! Initialize temp_sol_beg with values from previous time-step
	434	temp_sol_beg(:) = temp_sol_new(:)
	435
	436	! Read e_soil_lat from restart file or initialize
	437	IF (ok_Ecorr) THEN
	438	CALL restget_p (rest_id, 'e_soil_lat', nbp_glo, 1, 1, kjit, .TRUE., &
	439	e_soil_lat, "gather", nbp_glo, index_g)
	440	CALL setvar_p (e_soil_lat, val_exp,'NO_KEYWORD',zero)
	441	END IF
	442
[2650]	443	! Read gtemp from restart file
	444	CALL restget_p (rest_id, 'gtemp', nbp_glo, 1, 1, kjit, .TRUE., &
	445	gtemp, "gather", nbp_glo, index_g)
	446	CALL setvar_p (gtemp, val_exp,'NO_KEYWORD',zero)
	447
	448
[2868]	449	! Read variables calculated in thermosoil_coef from restart file
	450	! If the variables were not found in the restart file, the logical
	451	! calculate_coef will be true and thermosoil_coef will be called further below.
	452	! These variables need to be in the restart file to avoid a time shift that
	453	! would be done using thermosoil_coef at this stage.
	454	calculate_coef=.FALSE.
	455	CALL ioconf_setatt_p('UNITS', 'J m-2 K-1')
	456	CALL ioconf_setatt_p('LONG_NAME','Apparent surface heat capacity')
	457	CALL restget_p (rest_id, 'soilcap', nbp_glo, 1, 1, kjit, .TRUE., &
	458	soilcap, "gather", nbp_glo, index_g)
	459	IF (ALL(soilcap(:)==val_exp)) calculate_coef=.TRUE.
	460
	461	CALL ioconf_setatt_p('UNITS', 'W m-2')
	462	CALL ioconf_setatt_p('LONG_NAME','Apparent soil heat flux')
	463	CALL restget_p (rest_id, 'soilflx', nbp_glo, 1, 1, kjit, .TRUE., &
	464	soilflx, "gather", nbp_glo, index_g)
	465	IF (ALL(soilflx(:)==val_exp)) calculate_coef=.TRUE.
	466
	467	CALL ioconf_setatt_p('UNITS', '')
	468	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
	469	CALL restget_p (rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
	470	cgrnd, "gather", nbp_glo, index_g)
	471	IF (ALL(cgrnd(:,:)==val_exp)) calculate_coef=.TRUE.
	472
	473	CALL ioconf_setatt_p('UNITS', '')
	474	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
	475	CALL restget_p (rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
	476	dgrnd, "gather", nbp_glo, index_g)
	477	IF (ALL(dgrnd(:,:)==val_exp)) calculate_coef=.TRUE.
	478
	479	CALL ioconf_setatt_p('UNITS', '')
[3059]	480	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
	481	CALL restget_p (rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
	482	cgrnd_snow, "gather", nbp_glo, index_g)
	483	IF (ALL(cgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.
[2868]	484
	485	CALL ioconf_setatt_p('UNITS', '')
[3059]	486	CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
	487	CALL restget_p (rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
	488	dgrnd_snow, "gather", nbp_glo, index_g)
	489	IF (ALL(dgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.
[2868]	490
[3059]	491	CALL ioconf_setatt_p('UNITS', '')
	492	CALL ioconf_setatt_p('LONG_NAME','Coefficient of the linear extrapolation of surface temperature')
	493	CALL restget_p (rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, .TRUE., &
	494	lambda_snow, "gather", nbp_glo, index_g)
	495	IF (ALL(lambda_snow(:)==val_exp)) calculate_coef=.TRUE.
	496
[2942]	497	!! 2.2 Computes some physical constants and arrays depending on the soil vertical discretization
	498
	499	! Calculate so_capa_ice
[7206]	500	so_capa_ice(:) = so_capa_dry + mcs(:)capa_icerho_ice
[7199]	501	IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice(:)=', so_capa_ice(:),' and capa_ice=',capa_ice
[2581]	502
[2942]	503	! Computing some usefull constants for the numerical scheme
	504	! Use znt(depth of nodes) and zlt(depth of deeper layer interface) from vertical_soil module.
	505	DO jg=1,ngrnd-1
	506	dz1(jg) = un / (znt(jg+1) - znt(jg))
	507	dz5(jg) = (zlt(jg) - znt(jg)) * dz1(jg)
	508	ENDDO
	509	dz5(ngrnd) = 0.0
	510	lambda = znt(1) * dz1(1)
	511
[2956]	512	! Send out the temperature profile on the first nslm levels(the levels treated in hydrol)
[2942]	513	stempdiag(:,:) = ptn(:,1:nslm)
	514
	515
[2956]	516	!! 2.3. Computes cgrnd, dgrnd, soilflx and soilcap coefficients only if they were not found in restart file.
[2868]	517	IF (calculate_coef) THEN
[2922]	518	! Interpolate variables needed by thermosoil_coef to the thermal levels
	519	CALL thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh)
	520
[2868]	521	IF (printlev>=3) WRITE (numout,*) 'thermosoil_coef will be called in the intialization phase'
	522	CALL thermosoil_coef (&
[3059]	523	kjpindex, temp_sol_new, snow, njsc, &
[7206]	524	mcs, frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
[3059]	525	snowdz, snowrho, snowtemp, pb, &
	526	ptn, &
[3269]	527	soilcap, soilflx, cgrnd, dgrnd,&
	528	lambda_snow, cgrnd_snow, dgrnd_snow)
[2868]	529	END IF
[2650]	530
[2581]	531	END SUBROUTINE thermosoil_initialize
	532
	533
[947]	534	!! ================================================================================================================================
	535	!! SUBROUTINE : thermosoil_main
	536	!!
	537	!>\BRIEF Thermosoil_main computes the soil thermal properties and dynamics, ie solves
[2942]	538	!! the heat diffusion equation within the soil.
[947]	539	!!
	540	!! DESCRIPTION : The resolution of the soil heat diffusion equation
	541	!! relies on a numerical finite-difference implicit scheme
	542	!! fully described in the reference and in the header of the thermosoil module.
	543	!! - The dependency of the previous timestep hidden in the
	544	!! integration coefficients cgrnd and dgrnd (EQ1), calculated in thermosoil_coef, and
	545	!! called at the end of the routine to prepare for the next timestep.
	546	!! - The effective computation of the new soil temperatures is performed in thermosoil_profile.
	547	!!
	548	!! - thermosoil_coef calculates the coefficients for the numerical scheme for the very first iteration of thermosoil;
	549	!! after that, thermosoil_coef is called only at the end of the module to calculate the coefficients for the next timestep.
	550	!! - thermosoil_profile solves the numerical scheme.\n
	551	!!
	552	!! - Flags : one unique flag : THERMOSOIL_TPRO (to be set to the desired initial soil in-depth temperature in K; by default 280K)
	553	!!
[2956]	554	!! RECENT CHANGE(S) : Change vertical discretization (consistent with hydrology layers) and soil thermal properties (taking into account soil texture effects).
[947]	555	!!
	556	!! MAIN OUTPUT VARIABLE(S): vertically discretized soil temperatures ptn, soil
	557	!! thermal properties (pcapa, pkappa), apparent surface heat capacity (soilcap)
[2956]	558	!! and heat flux (soilflx) to be used in enerbil at the next timestep to solve
[947]	559	!! the surface energy balance.
	560	!!
	561	!! REFERENCE(S) :
	562	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
	563	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin' s PhD thesis relative to the thermal
	564	!! integration scheme has been scanned and is provided along with the documentation, with name :
	565	!! Hourdin_1992_PhD_thermal_scheme.pdf
	566	!!
	567	!! FLOWCHART :
	568	!! \latexonly
	569	!! \includegraphics[scale = 1]{thermosoil_flowchart.png}
	570	!! \endlatexonly
	571	!!
	572	!! \n
	573	!_ ================================================================================================================================
[8]	574
[2591]	575	SUBROUTINE thermosoil_main (kjit, kjpindex, &
[7206]	576	index, indexgrnd, mcs, &
[2222]	577	temp_sol_new, snow, soilcap, soilflx, &
	578	shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
[3059]	579	snowdz,snowrho,snowtemp,gtemp,pb,&
[3269]	580	mc_layh, mcl_layh, tmc_layh, njsc, frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
	581	lambda_snow, cgrnd_snow, dgrnd_snow)
[8]	582
[2222]	583	!! 0. Variable and parameter declaration
[8]	584
[947]	585	!! 0.1 Input variables
[8]	586
[947]	587	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
	588	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
	589	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! Restart_ file and history file identifier
	590	!! (unitless)
	591	INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier (unitless)
	592	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless)
	593	INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd !! Indeces of the points on the 3D map (vertical
	594	!! dimension towards the ground) (unitless)
[7207]	595	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3)
[3313]	596	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! Surface temperature at the present time-step,
[947]	597	!! Ts @tex ($K$) @endtex
	598	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass @tex ($kg$) @endtex.
	599	!! Caveat: when there is snow on the
	600	!! ground, the snow is integrated into the soil for
	601	!! the calculation of the thermal dynamics. It means
	602	!! that the uppermost soil layers can completely or
	603	!! partially consist in snow. In the second case, zx1
	604	!! and zx2 are the fraction of the soil layer
	605	!! consisting in snow and 'normal' soil, respectively
	606	!! This is calculated in thermosoil_coef.
[2942]	607	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree (0-1, unitless)
[2222]	608	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth
[3059]	609	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowrho !! Snow density
[3410]	610	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (inout) :: snowtemp !! Snow temperature (K)
	611	REAL(r_std), DIMENSION (kjpindex),INTENT (in) :: pb !! Surface presure (hPa)
[2942]	612	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid + ice) (m3/m3)
	613	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) (m3/m3)
	614	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
[2922]	615	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
[3020]	616	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
	617	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
	618	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
	619	!!(unitless,0-1)
[3059]	620	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_add !! additional surface temperature due to the melt of first layer
	621	!! at the present time-step @tex ($K$) @endtex
[2222]	622
[947]	623	!! 0.2 Output variables
[8]	624
[2222]	625	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ptnlev1 !! 1st level soil temperature
	626	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature
	627
	628
	629	!! 0.3 Modified variables
	630
[3020]	631	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! apparent surface heat capacity considering snow and soil surface
[947]	632	!! @tex ($J m^{-2} K^{-1}$) @endtex
[3020]	633	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx !! apparent soil heat flux considering snow and soil surface
	634	!! @tex ($W m^{-2}$) @endtex
[947]	635	!! , positive
	636	!! towards the soil, writen as Qg (ground heat flux)
	637	!! in the history files, and computed at the end of
	638	!! thermosoil for the calculation of Ts in enerbil,
	639	!! see EQ3.
[2942]	640	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile @tex ($K$) @endtex
[3269]	641	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
	642	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
	643	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
[8]	644
[947]	645	!! 0.4 Local variables
[8]	646
[2222]	647	INTEGER(i_std) :: jv,ji,ii
[4881]	648	REAL(r_std),DIMENSION (kjpindex) :: snowtemp_weighted!! Snow temperature weighted by snow density, only for diag (K)
	649	REAL(r_std),DIMENSION (kjpindex, nsnow) :: pkappa_snow_diag !! Only for diag, containing xios_default_val
	650	REAL(r_std),DIMENSION (kjpindex, nsnow) :: pcapa_snow_diag !! Only for diag, containing xios_default_val
	651	REAL(r_std),DIMENSION (kjpindex, nsnow) :: snowtemp_diag !! Only for diag, containing xios_default_val
[8]	652
[947]	653	!_ ================================================================================================================================
	654
	655	!! 3. Put the soil wetness diagnostic on the levels of the soil temperature
[8]	656
[947]	657	!!?? this could logically be put just before the last call to
	658	!!thermosoil_coef, as the results are used there...
[2922]	659	CALL thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh)
[947]	660
	661
[4594]	662	!! 4. Effective computation of the soil temperatures profile.
	663	!! cgrnd and dgrnd have been calculated in thermosoil_coef at the previous time step
	664	!! but they are correct for the actual time-step.
[3269]	665	CALL thermosoil_profile (kjpindex, temp_sol_new, &
[3059]	666	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
[3269]	667	ptn, stempdiag, snowtemp, &
	668	cgrnd_snow, dgrnd_snow)
[947]	669
	670
[4682]	671	!! 5. Call to thermosoil_energy_diag for calculation of diagnostic variables
	672	CALL thermosoil_energy_diag(kjpindex, temp_sol_new, soilcap)
[2222]	673
[4682]	674	!! Save ptn at current stage, to be used in thermosoil_readjust
	675	ptn_beg(:,:) = ptn(:,:)
	676
[947]	677	!! 6. Writing the history files according to the ALMA standards (or not..)
	678
[4881]	679	! Add XIOS default value where no snow
	680	DO ji=1,kjpindex
	681	IF (snow(ji) .GT. zero) THEN
	682	pkappa_snow_diag(ji,:) = pkappa_snow(ji,:)
	683	pcapa_snow_diag(ji,:) = pcapa_snow(ji,:)
	684	snowtemp_diag(ji,:) = snowtemp(ji,:)
	685	ELSE
	686	pkappa_snow_diag(ji,:) = xios_default_val
	687	pcapa_snow_diag(ji,:) = xios_default_val
	688	snowtemp_diag(ji,:) = xios_default_val
	689	END IF
[5096]	690	END DO
[4881]	691
[5470]	692	DO ji=1,kjpindex
	693	! Use min_sechiba instead of zero to avoid problem with division by zero
	694	IF (snow(ji) .GT. min_sechiba) THEN
	695	snowtemp_weighted(ji) = SUM(snowtemp(ji,:)*snowrho(ji,:))/SUM(snowrho(ji,:))
	696	ELSE
	697	snowtemp_weighted(ji) = xios_default_val
	698	END IF
	699	END DO
	700	CALL xios_orchidee_send_field("snowtemp_weighted",snowtemp_weighted)
	701
[1788]	702	CALL xios_orchidee_send_field("ptn",ptn)
[3408]	703	CALL xios_orchidee_send_field("soilflx",soilflx)
	704	CALL xios_orchidee_send_field("surfheat_incr",surfheat_incr)
	705	CALL xios_orchidee_send_field("coldcont_incr",coldcont_incr)
[3410]	706	CALL xios_orchidee_send_field("pkappa",pkappa)
[4881]	707	CALL xios_orchidee_send_field("pkappa_snow",pkappa_snow_diag)
[3410]	708	CALL xios_orchidee_send_field("pcapa",pcapa)
[4881]	709	CALL xios_orchidee_send_field("pcapa_snow",pcapa_snow_diag)
	710	CALL xios_orchidee_send_field("snowtemp",snowtemp_diag)
[4725]	711
[5096]	712
[8]	713	IF ( .NOT. almaoutput ) THEN
[1078]	714	CALL histwrite_p(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
[2222]	715	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
[5454]	716	CALL histwrite_p(hist_id, 'ptn_beg', kjit, ptn_beg, kjpindex*ngrnd, indexgrnd)
	717	CALL histwrite_p(hist_id, 'pkappa', kjit, pkappa, kjpindex*ngrnd, indexgrnd)
	718	CALL histwrite_p(hist_id, 'pcapa', kjit, pcapa, kjpindex*ngrnd, indexgrnd)
	719
	720	IF (ok_freeze_thermix) THEN
	721	CALL histwrite_p(hist_id, 'profil_froz', kjit, profil_froz, kjpindex*ngrnd, indexgrnd)
	722	CALL histwrite_p(hist_id, 'pcappa_supp', kjit, pcappa_supp, kjpindex*ngrnd, indexgrnd)
[2222]	723	END IF
[5454]	724	CALL histwrite_p(hist_id, 'shum_ngrnd_perma', kjit, shum_ngrnd_perma(:,:), kjpindex*ngrnd, indexgrnd)
	725
[8]	726	ELSE
[1078]	727	CALL histwrite_p(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
	728	CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
	729	CALL histwrite_p(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
	730	CALL histwrite_p(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
[8]	731	ENDIF
	732	IF ( hist2_id > 0 ) THEN
	733	IF ( .NOT. almaoutput ) THEN
[1078]	734	CALL histwrite_p(hist2_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
[8]	735	ELSE
[1078]	736	CALL histwrite_p(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
	737	CALL histwrite_p(hist2_id, 'Qg', kjit, soilflx, kjpindex, index)
	738	CALL histwrite_p(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
	739	CALL histwrite_p(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
[8]	740	ENDIF
	741	ENDIF
[947]	742
	743	!! 7. A last final call to thermosoil_coef
	744
	745	!! A last final call to thermosoil_coef, which calculates the different
[3059]	746	!!coefficients (cgrnd, dgrnd, soilcap, soilflx) from this time step to be
[947]	747	!!used at the next time step, either in the surface temperature calculation
	748	!!(soilcap, soilflx) or in the soil thermal numerical scheme.
[2222]	749	CALL thermosoil_coef (&
[3059]	750	kjpindex, temp_sol_new, snow, njsc, &
[7207]	751	mcs, frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
[3059]	752	snowdz, snowrho, snowtemp, pb, &
	753	ptn, &
[3269]	754	soilcap, soilflx, cgrnd, dgrnd,&
	755	lambda_snow, cgrnd_snow, dgrnd_snow)
[3059]	756
[8]	757
[2222]	758	! Save variables for explicit snow model
	759	gtemp(:) = ptn(:,1)
[2650]	760
[2222]	761	!! Initialize output arguments to be used in sechiba
[890]	762	ptnlev1(:) = ptn(:,1)
	763
[5470]	764	!! Surface temperature is forced to zero celcius if its value is larger than melting point
	765	DO ji=1,kjpindex
	766	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
	767	IF (temp_sol_new(ji) .GE. tp_00) THEN
	768	temp_sol_new(ji) = tp_00
	769	ENDIF
	770	END IF
	771	END DO
	772
[2348]	773	IF (printlev>=3) WRITE (numout,*) ' thermosoil_main done '
[8]	774
	775	END SUBROUTINE thermosoil_main
	776
[2581]	777	!! =============================================================================================================================
	778	!! SUBROUTINE : thermosoil_finalize
	779	!!
	780	!>\BRIEF Write to restart file
	781	!!
	782	!! DESCRIPTION : This subroutine writes the module variables and variables calculated in thermosoil
	783	!! to restart file
	784	!! \n
	785	!_ ==============================================================================================================================
[3059]	786	SUBROUTINE thermosoil_finalize (kjit, kjpindex, rest_id, gtemp, &
[3269]	787	soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
[947]	788
[2222]	789	!! 0. Variable and parameter declaration
[947]	790	!! 0.1 Input variables
[2581]	791	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
	792	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
	793	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier(unitless)
[2650]	794	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: gtemp !! First soil layer temperature
[2868]	795	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap
	796	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilflx
[3269]	797	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
	798	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
	799	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
[2650]	800
[947]	801	!_ ================================================================================================================================
[2222]	802
[2581]	803	!! 1. Write variables to restart file to be used for the next simulation
	804	IF (printlev>=3) WRITE (numout,*) 'Write restart file with THERMOSOIL variables'
	805
	806	CALL restput_p(rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, ptn, 'scatter', nbp_glo, index_g)
	807
[2222]	808	IF (ok_Ecorr) THEN
[2581]	809	CALL restput_p(rest_id, 'e_soil_lat', nbp_glo, 1 , 1, kjit, e_soil_lat, 'scatter', nbp_glo, index_g)
[2222]	810	END IF
	811
[2650]	812	CALL restput_p(rest_id, 'gtemp', nbp_glo, 1, 1, kjit, gtemp, 'scatter', nbp_glo, index_g)
[2868]	813	CALL restput_p(rest_id, 'soilcap', nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g)
	814	CALL restput_p(rest_id, 'soilflx', nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g)
	815	CALL restput_p(rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g)
	816	CALL restput_p(rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g)
[3059]	817	CALL restput_p(rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, cgrnd_snow, 'scatter', nbp_glo, index_g)
	818	CALL restput_p(rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, dgrnd_snow, 'scatter', nbp_glo, index_g)
	819	CALL restput_p(rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, lambda_snow, 'scatter', nbp_glo, index_g)
	820
[2581]	821	END SUBROUTINE thermosoil_finalize
[8]	822
[2436]	823
[947]	824	!! ================================================================================================================================
	825	!! SUBROUTINE : thermosoil_clear
	826	!!
[2956]	827	!>\BRIEF Deallocates the allocated arrays.
[2581]	828	!! The call of thermosoil_clear originates from sechiba_clear but the calling sequence and
[947]	829	!! its purpose require further investigation.
	830	!!
	831	!! DESCRIPTION : None
	832	!!
	833	!! RECENT CHANGE(S) : None
	834	!!
	835	!! MAIN OUTPUT VARIABLE(S): None
	836	!!
	837	!! REFERENCE(S) : None
	838	!!
	839	!! FLOWCHART : None
	840	!! \n
	841	!_ ================================================================================================================================
	842
[8]	843	SUBROUTINE thermosoil_clear()
	844
	845	IF ( ALLOCATED (ptn)) DEALLOCATE (ptn)
	846	IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd)
	847	IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd)
	848	IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa)
	849	IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa)
[3059]	850	IF ( ALLOCATED (pcapa_snow)) DEALLOCATE (pcapa_snow)
	851	IF ( ALLOCATED (pkappa_snow)) DEALLOCATE (pkappa_snow)
[8]	852	IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en)
	853	IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg)
	854	IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg)
	855	IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr)
	856	IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr)
[2222]	857	IF ( ALLOCATED (shum_ngrnd_perma)) DEALLOCATE (shum_ngrnd_perma)
	858	IF ( ALLOCATED (profil_froz)) DEALLOCATE (profil_froz)
[2922]	859	IF ( ALLOCATED (mc_layt)) DEALLOCATE (mc_layt)
	860	IF ( ALLOCATED (mcl_layt)) DEALLOCATE (mcl_layt)
	861	IF ( ALLOCATED (tmc_layt)) DEALLOCATE (tmc_layt)
[8]	862	END SUBROUTINE thermosoil_clear
[947]	863
	864
	865	!! ================================================================================================================================
	866	!! SUBROUTINE : thermosoil_coef
	867	!!
	868	!>\BRIEF Calculate soil thermal properties, integration coefficients, apparent heat flux,
	869	!! surface heat capacity,
	870	!!
	871	!! DESCRIPTION : This routine computes : \n
	872	!! 1. the soil thermal properties. \n
	873	!! 2. the integration coefficients of the thermal numerical scheme, cgrnd and dgrnd,
	874	!! which depend on the vertical grid and on soil properties, and are used at the next
	875	!! timestep.\n
[2956]	876	!! 3. the soil apparent heat flux and surface heat capacity (soilflx
	877	!! and soilcap), used by enerbil to compute the surface temperature at the next
[947]	878	!! timestep.\n
[2956]	879	!! - The soil thermal properties depend on water content (shum_ngrnd_perma, shumdiag_perma,
	880	!! mc_layt, mcl_layt, tmc_layt), dominant soil texture(njsc), and on the presence
[947]	881	!! of snow : snow is integrated into the soil for the thermal calculations, ie if there
	882	!! is snow on the ground, the first thermal layer(s) consist in snow, depending on the
	883	!! snow-depth. If a layer consists out of snow and soil, wheighed soil properties are
	884	!! calculated\n
	885	!! - The coefficients cgrnd and dgrnd are the integration
	886	!! coefficients for the thermal scheme \n
	887	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n
	888	!! -- EQ1 -- \n
	889	!! They correspond respectively to $\beta$ and $\alpha$ from F. Hourdin\'s thesis and
	890	!! their expression can be found in this document (eq A19 and A20)
[2956]	891	!! - soilcap and soilflx are the apparent surface heat capacity and flux
[947]	892	!! used in enerbil at the next timestep to solve the surface
	893	!! balance for Ts (EQ3); they correspond to $C_s$ and $F_s$ in F.
	894	!! Hourdin\'s PhD thesis and are expressed in eq. A30 and A31. \n
[2956]	895	!! soilcap*(Ts(t)-Ts(t-1))/dt=soilflx+otherfluxes(Ts(t)) \n
[947]	896	!! -- EQ3 --\n
	897	!!
	898	!! RECENT CHANGE(S) : None
	899	!!
[3059]	900	!! MAIN OUTPUT VARIABLE(S): cgrnd, dgrnd, pcapa, pkappa, soilcap, soilflx
[947]	901	!!
	902	!! REFERENCE(S) :
	903	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
	904	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
	905	!! integration scheme has been scanned and is provided along with the documentation, with name :
	906	!! Hourdin_1992_PhD_thermal_scheme.pdf
	907	!!
	908	!! FLOWCHART : None
	909	!! \n
	910	!_ ================================================================================================================================
	911
[3059]	912	SUBROUTINE thermosoil_coef (kjpindex, temp_sol_new, snow, njsc, &
[7206]	913	mcs, frac_snow_veg, frac_snow_nobio,totfrac_nobio, &
[3059]	914	snowdz, snowrho, snowtemp, pb, &
	915	ptn, &
[3269]	916	soilcap, soilflx, cgrnd, dgrnd,&
	917	lambda_snow, cgrnd_snow, dgrnd_snow)
[8]	918
[2222]	919	!! 0. Variables and parameter declaration
[8]	920
[947]	921	!! 0.1 Input variables
[8]	922
[947]	923	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
	924	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! soil surface temperature @tex ($K$) @endtex
	925	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! snow mass @tex ($Kg$) @endtex
[7207]	926	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class
	927	!! in the grid cell (1-nscm, unitless)
	928	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3)
[3020]	929	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
	930	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
	931	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
	932	!!(unitless,0-1)
[3410]	933	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth (m)
	934	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
	935	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
	936	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa)
[3020]	937
[947]	938	!! 0.2 Output variables
	939
[3020]	940	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !! surface heat capacity considering snow and soil surface
[947]	941	!! @tex ($J m^{-2} K^{-1}$) @endtex
[3020]	942	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !! surface heat flux considering snow and soil surface @tex ($W m^{-2}$) @endtex,
[947]	943	!! positive towards the
	944	!! soil, writen as Qg (ground heat flux) in the history
	945	!! files.
	946	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !! matrix coefficient for the computation of soil
	947	!! temperatures (beta in F. Hourdin thesis)
	948	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !! matrix coefficient for the computation of soil
	949	!! temperatures (alpha in F. Hourdin thesis)
	950
	951
	952	!! 0.3 Modified variable
	953
[3059]	954	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (inout):: ptn !! vertically discretized soil temperatures. ptn is only modified if ok_Ecorr.
[3269]	955	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
	956	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme
	957	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme
[947]	958
	959	!! 0.4 Local variables
	960
	961	INTEGER(i_std) :: ji, jg
[2436]	962	REAL(r_std), DIMENSION (kjpindex,ngrnd-1) :: zdz1 !! numerical (buffer) constant
	963	!! @tex ($W m^{-1} K^{-1}$) @endtex
	964	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: zdz2 !! numerical (buffer) constant
	965	!! @tex ($W m^{-1} K^{-1}$) @endtex
	966	REAL(r_std), DIMENSION (kjpindex) :: z1 !! numerical constant @tex ($W m^{-1} K^{-1}$) @endtex
[3020]	967	REAL(r_std), DIMENSION (kjpindex) :: soilcap_nosnow !! surface heat capacity
	968	!! @tex ($J m^{-2} K^{-1}$)
	969	!! @endtex
	970	REAL(r_std), DIMENSION (kjpindex) :: soilflx_nosnow !! surface heat flux @tex ($W m^{-2}$) @endtex,
	971	!! positive towards the soil, written as Qg
	972	!!(ground heat flux in the history files).
[3059]	973	REAL(r_std), DIMENSION (kjpindex) :: snowcap !! apparent snow heat capacity @tex ($J m^{-2} K^{-1}$)
	974	REAL(r_std), DIMENSION (kjpindex) :: snowflx !! apparent snow-atmosphere heat flux @tex ($W m^{-2}$) @endtex
	975	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz1_snow
	976	REAL(r_std), DIMENSION (kjpindex,nsnow) :: ZSNOWDZM
	977	REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz2_snow
	978	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz1_snow
	979	REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz2_snow
	980	REAL(r_std), DIMENSION (kjpindex) :: z1_snow
[4544]	981	REAL(r_std), DIMENSION (kjpindex) :: snowflxtot !! Total snow flux (including snow on vegetated and bare soil and nobio areas)
	982	!! @tex ($W m^{-2}$) @endtex
	983	!! positive towards the soil
[3020]	984
[947]	985	!_ ================================================================================================================================
	986
	987	!! 1. Computation of the soil thermal properties
	988
	989	! Computation of the soil thermal properties; snow properties are also accounted for
[2222]	990	IF (ok_freeze_thermix) THEN
[7206]	991	CALL thermosoil_getdiff( kjpindex, snow, ptn, mcs, njsc, snowrho, snowtemp, pb )
[5470]	992	ELSE
	993	! Special case without soil freezing
[7206]	994	CALL thermosoil_getdiff_old_thermix_without_snow( kjpindex, mcs, njsc, snowrho, snowtemp, pb )
[2222]	995	ENDIF
[947]	996
[2222]	997	! Energy conservation : Correction to make sure that the same latent heat is released and
	998	! consumed during freezing and thawing
	999	IF (ok_Ecorr) THEN
	1000	CALL thermosoil_readjust(kjpindex, ptn)
	1001	ENDIF
[947]	1002
	1003
[3059]	1004	!! 2. Computation of the coefficients of the numerical integration scheme for the soil layers
[947]	1005
[3059]	1006	!! 2.1 Calculate numerical coefficients zdz1 and zdz2
[8]	1007	DO jg=1,ngrnd
	1008	DO ji=1,kjpindex
[2942]	1009	zdz2(ji,jg)=pcapa(ji,jg) * dlt(jg)/dt_sechiba
[8]	1010	ENDDO
	1011	ENDDO
[947]	1012
[8]	1013	DO jg=1,ngrnd-1
	1014	DO ji=1,kjpindex
	1015	zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg)
	1016	ENDDO
	1017	ENDDO
[947]	1018
[3059]	1019	!! 2.2 Calculate coefficients cgrnd and dgrnd for soil
[8]	1020	DO ji = 1,kjpindex
	1021	z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1)
	1022	cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji)
	1023	dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji)
	1024	ENDDO
	1025
	1026	DO jg = ngrnd-1,2,-1
	1027	DO ji = 1,kjpindex
	1028	z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg)))
	1029	cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji)
	1030	dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji)
	1031	ENDDO
	1032	ENDDO
	1033
[2650]	1034
[3059]	1035	!! 3. Computation of the coefficients of the numerical integration scheme for the snow layers
	1036
	1037	!! 3.1 Calculate numerical coefficients zdz1_snow, zdz2_snow and lambda_snow
	1038	DO ji = 1, kjpindex
	1039
[5470]	1040	! Calculate internal values
	1041	DO jg = 1, nsnow
	1042	ZSNOWDZM(ji,jg) = MAX(snowdz(ji,jg),psnowdzmin)
	1043	ENDDO
	1044	dz2_snow(ji,:)=ZSNOWDZM(ji,:)
	1045
	1046	DO jg = 1, nsnow-1
	1047	dz1_snow(ji,jg) = 2.0 / (dz2_snow(ji,jg+1)+dz2_snow(ji,jg))
	1048	ENDDO
	1049
	1050	lambda_snow(ji) = dz2_snow(ji,1)/2.0 * dz1_snow(ji,1)
	1051
	1052	DO jg=1,nsnow
	1053	zdz2_snow(ji,jg)=pcapa_snow(ji,jg) * dz2_snow(ji,jg)/dt_sechiba
	1054	ENDDO
	1055
	1056	DO jg=1,nsnow-1
	1057	zdz1_snow(ji,jg) = dz1_snow(ji,jg) * pkappa_snow(ji,jg)
	1058	ENDDO
	1059
	1060	! the bottom snow
	1061	zdz1_snow(ji,nsnow) = pkappa_snow(ji,nsnow) / ( zlt(1) + dz2_snow(ji,nsnow)/2 )
	1062
[3059]	1063	ENDDO
	1064
	1065	!! 3.2 Calculate coefficients cgrnd_snow and dgrnd_snow for snow
	1066	DO ji = 1,kjpindex
[5470]	1067	! bottom level
	1068	z1_snow(ji) = zdz2(ji,1)+(un-dgrnd(ji,1))*zdz1(ji,1)+zdz1_snow(ji,nsnow)
	1069	cgrnd_snow(ji,nsnow) = (zdz2(ji,1) * ptn(ji,1) + zdz1(ji,1) * cgrnd(ji,1) ) / z1_snow(ji)
	1070	dgrnd_snow(ji,nsnow) = zdz1_snow(ji,nsnow) / z1_snow(ji)
	1071
	1072	! next-to-bottom level
	1073	z1_snow(ji) = zdz2_snow(ji,nsnow)+(un-dgrnd_snow(ji,nsnow))*zdz1_snow(ji,nsnow)+zdz1_snow(ji,nsnow-1)
	1074	cgrnd_snow(ji,nsnow-1) = (zdz2_snow(ji,nsnow)*snowtemp(ji,nsnow)+&
	1075	zdz1_snow(ji,nsnow)*cgrnd_snow(ji,nsnow))/z1_snow(ji)
	1076	dgrnd_snow(ji,nsnow-1) = zdz1_snow(ji,nsnow-1) / z1_snow(ji)
	1077
	1078	DO jg = nsnow-1,2,-1
	1079	z1_snow(ji) = un / (zdz2_snow(ji,jg) + zdz1_snow(ji,jg-1) + zdz1_snow(ji,jg) * (un - dgrnd_snow(ji,jg)))
	1080	cgrnd_snow(ji,jg-1) = (snowtemp(ji,jg) * zdz2_snow(ji,jg) + zdz1_snow(ji,jg) * cgrnd_snow(ji,jg)) * z1_snow(ji)
	1081	dgrnd_snow(ji,jg-1) = zdz1_snow(ji,jg-1) * z1_snow(ji)
	1082	ENDDO
[3059]	1083	ENDDO
	1084
	1085
	1086
	1087	!! 4. Computation of the apparent ground heat flux
	1088	!! Computation of apparent snow-atmosphere flux
	1089	DO ji = 1,kjpindex
[5470]	1090	snowflx(ji) = zdz1_snow(ji,1) * (cgrnd_snow(ji,1) + (dgrnd_snow(ji,1)-1.) * snowtemp(ji,1))
	1091	snowcap(ji) = (zdz2_snow(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd_snow(ji,1)) * zdz1_snow(ji,1))
	1092	z1_snow(ji) = lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un
	1093	snowcap(ji) = snowcap(ji) / z1_snow(ji)
	1094	snowflx(ji) = snowflx(ji) + &
	1095	snowcap(ji) * (snowtemp(ji,1) * z1_snow(ji) - lambda_snow(ji) * cgrnd_snow(ji,1) - temp_sol_new(ji)) / dt_sechiba
[3059]	1096	ENDDO
	1097
	1098
[947]	1099	!! Computation of the apparent ground heat flux (> towards the soil) and
	1100	!! apparent surface heat capacity, used at the next timestep by enerbil to
	1101	!! compute the surface temperature.
[8]	1102	DO ji = 1,kjpindex
[3020]	1103	soilflx_nosnow(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1))
	1104	soilcap_nosnow(ji) = (zdz2(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd(ji,1)) * zdz1(ji,1))
[8]	1105	z1(ji) = lambda * (un - dgrnd(ji,1)) + un
[3020]	1106	soilcap_nosnow(ji) = soilcap_nosnow(ji) / z1(ji)
	1107	soilflx_nosnow(ji) = soilflx_nosnow(ji) + &
	1108	& soilcap_nosnow(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dt_sechiba
[8]	1109	ENDDO
	1110
[3020]	1111	!! Add snow fraction
[5470]	1112	! Using an effective heat capacity and heat flux by a simple pondering of snow and soil fraction
	1113	DO ji = 1, kjpindex
	1114	soilcap(ji) = snowcap(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
	1115	soilcap_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
	1116	soilcap_nosnow(ji)(1-(frac_snow_veg(ji)(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
	1117	soilflx(ji) = snowflx(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
	1118	soilflx_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
	1119	soilflx_nosnow(ji)(1-(frac_snow_veg(ji)(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
	1120	ENDDO
[3020]	1121
[4544]	1122	! Total snow flux (including snow on vegetated and bare soil and nobio areas)
	1123	DO ji = 1, kjpindex
[4546]	1124	snowflxtot(ji) = snowflx(ji)frac_snow_veg(ji)(1-totfrac_nobio(ji)) + &
	1125	soilflx_nosnow(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)
[4544]	1126	ENDDO
	1127	CALL xios_orchidee_send_field("snowflxtot",snowflxtot(:))
	1128
[2348]	1129	IF (printlev>=3) WRITE (numout,*) ' thermosoil_coef done '
[8]	1130
	1131	END SUBROUTINE thermosoil_coef
[947]	1132
	1133
	1134	!! ================================================================================================================================
	1135	!! SUBROUTINE : thermosoil_profile
	1136	!!
	1137	!>\BRIEF In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile;
[2956]	1138	!!
[947]	1139	!!
	1140	!! DESCRIPTION : The calculation of the new soil temperature profile is based on
	1141	!! the cgrnd and dgrnd values from the previous timestep and the surface temperature Ts aka temp_sol_new. (see detailed
	1142	!! explanation in the header of the thermosoil module or in the reference).\n
	1143	!! T(k+1)=cgrnd(k)+dgrnd(k)*T(k)\n
	1144	!! -- EQ1 --\n
[2956]	1145	!! Ts=(1+lambda)T(1) -lambdaT(2)\n
[947]	1146	!! -- EQ2--\n
	1147	!!
	1148	!! RECENT CHANGE(S) : None
	1149	!!
	1150	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis),
	1151	!! stempdiag (soil temperature profile on the diagnostic axis)
	1152	!!
	1153	!! REFERENCE(S) :
	1154	!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
	1155	!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
	1156	!! integration scheme has been scanned and is provided along with the documentation, with name :
	1157	!! Hourdin_1992_PhD_thermal_scheme.pdf
	1158	!!
	1159	!! FLOWCHART : None
	1160	!! \n
	1161	!_ ================================================================================================================================
[8]	1162
[3269]	1163	SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, &
[3059]	1164	frac_snow_veg, frac_snow_nobio, totfrac_nobio, &
[3269]	1165	ptn, stempdiag, snowtemp, &
	1166	cgrnd_snow, dgrnd_snow)
[8]	1167
[947]	1168	!! 0. Variables and parameter declaration
[8]	1169
[947]	1170	!! 0.1 Input variables
	1171
	1172	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
	1173	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step
	1174	!! @tex ($K$) @endtex
[3020]	1175	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
	1176	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
	1177	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
	1178	!! (unitless,0-1)
[3410]	1179	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
[3269]	1180	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
	1181	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
[3020]	1182
[3059]	1183	!! 0.3 Modified variables
[3269]	1184
[3020]	1185
[947]	1186	!! 0.2 Output variables
[3059]	1187	REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (out) :: ptn !! vertically discretized soil temperatures
	1188	!! @tex ($K$) @endtex
[2942]	1189	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! diagnostic temperature profile
[947]	1190	!! @tex ($K$) @endtex
	1191
	1192	!! 0.4 Local variables
	1193
	1194	INTEGER(i_std) :: ji, jg
[3020]	1195	REAL(r_std) :: temp_sol_eff !! effective surface temperature including snow and soil
[2222]	1196
[947]	1197	!_ ================================================================================================================================
	1198
[3269]	1199	!! 1. Computes the soil temperatures ptn.
[2222]	1200
[3269]	1201	!! 1.1. ptn(jg=1) using EQ1 and EQ2
[3059]	1202	DO ji = 1,kjpindex
	1203
[5470]	1204	! Using an effective surface temperature by a simple pondering
	1205	temp_sol_eff=snowtemp(ji,nsnow)frac_snow_veg(ji)(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation
	1206	temp_sol_new(ji)SUM(frac_snow_nobio(ji,:))totfrac_nobio(ji)+ & ! weights related to SCF on nobio
	1207	temp_sol_new(ji)(1-(frac_snow_veg(ji)(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
	1208	! Soil temperature calculation with explicit snow if there is snow on the ground
	1209	ptn(ji,1) = cgrnd_snow(ji,nsnow) + dgrnd_snow(ji,nsnow) * temp_sol_eff
[8]	1210	ENDDO
	1211
[947]	1212	!! 1.2. ptn(jg=2:ngrnd) using EQ1.
[8]	1213	DO jg = 1,ngrnd-1
	1214	DO ji = 1,kjpindex
	1215	ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg)
	1216	ENDDO
	1217	ENDDO
	1218
[3269]	1219	!! 2. Assigne the soil temperature to the output variable. It is already on the right axis.
[2942]	1220	stempdiag(:,:) = ptn(:,1:nslm)
[8]	1221
[2348]	1222	IF (printlev>=3) WRITE (numout,*) ' thermosoil_profile done '
[8]	1223
	1224	END SUBROUTINE thermosoil_profile
[947]	1225
[2922]	1226	!================================================================================================================================
	1227	!! SUBROUTINE : thermosoil_cond
	1228	!!
	1229	!>\BRIEF Calculate soil thermal conductivity.
	1230	!!
	1231	!! DESCRIPTION : This routine computes soil thermal conductivity
	1232	!! Code introduced from NOAH LSM.
	1233	!!
	1234	!! RECENT CHANGE(S) : None
	1235	!!
	1236	!! MAIN OUTPUT VARIABLE(S): cnd
	1237	!!
	1238	!! REFERENCE(S) :
	1239	!! Farouki, O.T.,1986: Thermal Properties of Soils. Series on Rock
	1240	!! and Soil Mechanics, Vol. 11, Trans Tech, 136 PP.
	1241	!! Johansen, O., 1975: Thermal Conductivity of Soils. Ph.D. Thesis,
	1242	!! University of Trondheim,
	1243	!! Peters-Lidard, C. D., Blackburn, E., Liang, X., & Wood, E. F.,
	1244	!! 1998: The effect of soil thermal conductivity
	1245	!! Parameterization on Surface Energy fluxes
	1246	!! and Temperatures. J. of The Atmospheric Sciences,
	1247	!! Vol. 55, pp. 1209-1224.
	1248	!! Modify histroy:
	1249	!!
	1250	!! FLOWCHART : None
	1251	!! \n
	1252	!_
	1253	!================================================================================================================================
[947]	1254
[7206]	1255	SUBROUTINE thermosoil_cond (kjpindex, njsc, mcs, smc, qz, sh2o, cnd)
[2922]	1256
	1257	!! 0. Variables and parameter declaration
	1258
	1259	!! 0.1 Input variables
[2956]	1260	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
	1261	INTEGER(i_std), DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
[7206]	1262	REAL(r_std), DIMENSION (kjpindex), INTENT(IN) :: mcs !! Saturated moisture content (m3/m3)
[2956]	1263	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3)
[7199]	1264	REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1)
[2956]	1265	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o
[2922]	1266
	1267	!! 0.2 Output variables
[2956]	1268	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(OUT) :: cnd !! Soil Thermal Conductivity (W/m/k)
[2922]	1269
	1270	!! 0.3 Modified variables
	1271
	1272	!! 0.4 Local variables
[4767]	1273	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: ake !! Kersten Number (unitless)
[2956]	1274	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: thksat !! Saturated Thermal Conductivity (W/m/k)
	1275	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: satratio !! Degree of Saturation (0-1)
	1276	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xu !! Unfrozen Volume For Saturation (0-1)
	1277	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xunfroz !! Unfrozon Volume Fraction (0-1)
	1278	REAL(r_std) :: thko !! Thermal Conductivity for Other Ssoil Components (W/m/k)
	1279	REAL(r_std) :: gammd !! Dry Dendity (kg/m3)
	1280	REAL(r_std) :: thkdry !! Dry Thermal Conductivity (W/m/k)
	1281	REAL(r_std) :: thks !! Thermal Conductivity for the Solids Combined (Quartz + Other) (W/m/k)
	1282	REAL(r_std), PARAMETER :: THKICE = 2.2 !! Ice Thermal Conductivity (W/m/k)
	1283	REAL(r_std), PARAMETER :: THKQTZ = 7.7 !! Thermal Conductivity for Quartz (W/m/k)
	1284	REAL(r_std), PARAMETER :: THKW = 0.57 !! Water Thermal Conductivity (W/m/k)
	1285	INTEGER(i_std) :: ji, jg, jst
[2922]	1286
	1287	!_================================================================================================================================
	1288
[2956]	1289	!! 1. Dry and Saturated Thermal Conductivity.
	1290
	1291	DO ji = 1,kjpindex
	1292	jst = njsc(ji)
	1293
	1294	!! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1)
[7206]	1295	gammd = (1. - mcs(ji))*2700.
[2956]	1296	thkdry = (0.135* gammd+ 64.7)/ (2700. - 0.947* gammd)
	1297
	1298	!! 1.2. thermal conductivity of "other" soil components
	1299	IF (qz(jst) > 0.2) THEN
	1300	thko = 2.0
	1301	ELSEIF (qz(jst) <= 0.2) THEN
	1302	thko = 3.0
	1303	ENDIF
	1304
	1305	!! 1.3. Thermal conductivity of solids
	1306	thks = (THKQTZ ** qz(jst))* (thko ** (1. - qz(jst)))
	1307
	1308	DO jg = 1,ngrnd
	1309	!! 1.4. saturation ratio
[7206]	1310	satratio(ji,jg) = smc(ji,jg) / mcs(ji)
[2922]	1311
[2956]	1312	!! 1.5. Saturated Thermal Conductivity (thksat)
[3593]	1313	IF ( smc(ji,jg) > min_sechiba ) THEN
	1314	xunfroz(ji,jg) = sh2o(ji,jg) / smc(ji,jg) ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen))
[7206]	1315	xu(ji,jg) = xunfroz(ji,jg) * mcs(ji) ! Unfrozen volume for saturation (porosity*xunfroz)
	1316	thksat(ji,jg) = thks ** (1. - mcs(ji))* THKICE ** (mcs(ji) - xu(ji,jg))* THKW ** (xu(ji,jg))
[3593]	1317	ELSE
	1318	! this value will not be used since ake=0 for this case
	1319	thksat(ji,jg)=0
	1320	END IF
[2956]	1321	END DO ! DO jg = 1,ngrnd
	1322
[4767]	1323	!! 2. Kersten Number (ake)
[2956]	1324	DO jg = 1,ngrnd
	1325	IF ( (sh2o(ji,jg) + 0.0005) < smc(ji,jg) ) THEN
	1326	! Frozen
	1327	ake(ji,jg) = satratio(ji,jg)
	1328	ELSE
	1329	! Unfrozen
[4767]	1330	! Eq 11 in Peters-Lidard et al., 1998
[2956]	1331	IF ( satratio(ji,jg) > 0.1 ) THEN
[4767]	1332	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
	1333	! Coarse
	1334	ake(ji,jg) = 0.7 * LOG10 (SATRATIO(ji,jg)) + 1.0
	1335	ELSE
	1336	! Fine
	1337	ake(ji,jg) = LOG10 (satratio(ji,jg)) + 1.0
	1338	ENDIF
[2956]	1339	ELSEIF ( satratio(ji,jg) > 0.05 .AND. satratio(ji,jg) <= 0.1 ) THEN
[4767]	1340	IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN
	1341	! Coarse
	1342	ake(ji,jg) = 0.7 * LOG10 (satratio(ji,jg)) + 1.0
	1343	ELSE
	1344	! Fine
	1345	ake(ji,jg) = 0.0
	1346	ENDIF
[2956]	1347	ELSE
	1348	ake(ji,jg) = 0.0 ! use k = kdry
	1349	END IF
	1350	END IF
	1351	END DO ! DO jg = 1,ngrnd
	1352
	1353	!! 3. Thermal conductivity (cnd)
	1354	DO jg = 1,ngrnd
	1355	cnd(ji,jg) = ake(ji,jg) * (thksat(ji,jg) - thkdry) + thkdry
	1356	END DO ! DO jg = 1,ngrnd
	1357
	1358	END DO !DO ji = 1,kjpindex
[2922]	1359
	1360	END SUBROUTINE thermosoil_cond
	1361
	1362
[947]	1363	!! ================================================================================================================================
	1364	!! SUBROUTINE : thermosoil_humlev
[8]	1365	!!
[2922]	1366	!>\BRIEF Interpolate variables from the hydrology layers to the thermodynamic layers
[8]	1367	!!
[2922]	1368	!! DESCRIPTION : Interpolate the volumetric soil moisture content from the node to the interface of the layer.
	1369	!! The values for the deep layers in thermosoil where hydrology is not existing are constant.
	1370	!! No interpolation is needed for the total soil moisture content and for the soil saturation degree.
[8]	1371	!!
[947]	1372	!! RECENT CHANGE(S) : None
[8]	1373	!!
[2922]	1374	!! MAIN OUTPUT VARIABLE(S): mc_layt, mcl_layt, tmc_layt, shum_ngrnd_perma
[947]	1375	!!
	1376	!! REFERENCE(S) : None
	1377	!!
	1378	!! FLOWCHART : None
	1379	!! \n
	1380	!_ ================================================================================================================================
[2922]	1381	SUBROUTINE thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh)
[947]	1382
	1383	!! 0. Variables and parameter declaration
	1384
	1385	!! 0.1 Input variables
	1386
[2922]	1387	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
[2942]	1388	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree on the diagnostic axis (0-1, unitless)
	1389	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid+ice) [m/s]
	1390	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) [m/s]
	1391	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm]
[947]	1392
	1393	!! 0.2 Output variables
	1394
	1395	!! 0.3 Modified variables
	1396
	1397	!! 0.4 Local variables
[2917]	1398	INTEGER(i_std) :: ji, jd
[947]	1399
	1400	!_ ================================================================================================================================
	1401
[2917]	1402	IF (printlev >= 4) WRITE(numout,*) 'Start thermosoil_humlev'
	1403
[2922]	1404	! The values for the deep layers in thermosoil where hydrology is not existing are constant.
[2917]	1405	! For exemple if thermosoil uses 8m, and hydrol uses 2m vertical discretization,
	1406	! the values between 2m and 8m are constant.
[3480]	1407	! The moisture computed in hydrol is at the nodes (except for the
	1408	! top and bottom layer which are at interfaces)
	1409	! A linear interpolation is applied to obtain the moisture values at
	1410	! the interfaces (mc_layt), from the mc_layh at the nodes
	1411
[947]	1412	DO ji=1,kjpindex
[2942]	1413	DO jd = 1, nslm
[3480]	1414	IF(jd == 1) THEN ! the moisture at the 1st interface mc_layh(1) is at the surface, no interpolation
[3016]	1415	mc_layt(ji,jd) = mc_layh(ji,jd)
	1416	mcl_layt(ji,jd) = mcl_layh(ji,jd)
[3480]	1417	ELSEIF(jd == 2) THEN !! the mc_layt at the 2nd interface is interpolated using mc_layh(1) at surface and mc_layh(2) at the node
[3016]	1418	mc_layt(ji, jd) = mc_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
	1419	mc_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0)
	1420	mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + &
	1421	mcl_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0)
[3480]	1422	ELSEIF(jd == nslm) THEN ! the mc_layt at the nslm interface is interpolated using mc_layh(nslm) and mc_layh(nslm-1)
[3016]	1423	mc_layt(ji, jd) = mc_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
	1424	mc_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1))
	1425	mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + &
	1426	mcl_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1))
[3480]	1427	ELSE ! the mc_layt at the other interfaces are interpolated using mc_layh at adjacent nodes.
[3016]	1428	mc_layt(ji, jd) = mc_layh(ji, jd-1)(1-dz5(jd-1)) + mc_layh(ji,jd)dz5(jd-1)
	1429	mcl_layt(ji, jd) = mcl_layh(ji, jd-1)(1-dz5(jd-1)) + mcl_layh(ji,jd)dz5(jd-1)
[2922]	1430	ENDIF
[3016]	1431
[2917]	1432	shum_ngrnd_perma(ji,jd) = shumdiag_perma(ji,jd)
[2922]	1433	tmc_layt(ji,jd) = tmc_layh(ji,jd)
[947]	1434	ENDDO
[2922]	1435
	1436	! The deep layers in thermosoil where hydro is not existing
[2942]	1437	DO jd = nslm+1, ngrnd
	1438	shum_ngrnd_perma(ji,jd) = shumdiag_perma(ji,nslm)
	1439	mc_layt(ji,jd) = mc_layh(ji,nslm)
	1440	mcl_layt(ji,jd) = mcl_layh(ji,nslm)
[4619]	1441	tmc_layt(ji,jd) = tmc_layh(ji,nslm)/dlt(nslm) *dlt(jd)
[947]	1442	ENDDO
[8]	1443	ENDDO
[2922]	1444
[2917]	1445	IF (printlev >= 4) WRITE(numout,*) 'thermosoil_humlev done'
[2436]	1446
[8]	1447	END SUBROUTINE thermosoil_humlev
	1448
[947]	1449
	1450	!! ================================================================================================================================
[4682]	1451	!! SUBROUTINE : thermosoil_energy_diag
[947]	1452	!!
[4682]	1453	!>\BRIEF Calculate diagnostics
[947]	1454	!!
[4682]	1455	!! DESCRIPTION : Calculate diagnostic variables coldcont_incr and coldcont_incr
[947]	1456	!!
	1457	!! RECENT CHANGE(S) : None
	1458	!!
[4682]	1459	!! MAIN OUTPUT VARIABLE(S) :
[947]	1460	!!
	1461	!! REFERENCE(S) : None
	1462	!!
	1463	!! FLOWCHART : None
	1464	!! \n
	1465	!_ ================================================================================================================================
	1466
[4682]	1467	SUBROUTINE thermosoil_energy_diag(kjpindex, temp_sol_new, soilcap)
[947]	1468
	1469	!! 0. Variables and parameter declaration
	1470
	1471	!! 0.1 Input variables
	1472
	1473	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
	1474	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step, Ts
	1475	!! @tex ($K$) @endtex
	1476	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Apparent surface heat capacity
	1477	!! @tex ($J m^{-2} K^{-1}$) @endtex,
	1478	!! see eq. A29 of F. Hourdin\'s PhD thesis.
	1479
	1480	!! 0.2 Output variables
	1481
	1482	!! 0.3 Modified variables
	1483
	1484	!! 0.4 Local variables
	1485
	1486	INTEGER(i_std) :: ji, jg
	1487	!_ ================================================================================================================================
[5293]	1488
[8]	1489	! Sum up the energy content of all layers in the soil.
	1490	DO ji = 1, kjpindex
[947]	1491
[8]	1492	IF (pcapa_en(ji,1) .LE. sn_capa) THEN
[947]	1493
[8]	1494	! Verify the energy conservation in the surface layer
	1495	coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
	1496	surfheat_incr(ji) = zero
	1497	ELSE
[947]	1498
[8]	1499	! Verify the energy conservation in the surface layer
	1500	surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
	1501	coldcont_incr(ji) = zero
	1502	ENDIF
	1503	ENDDO
	1504
[4682]	1505	! Save temp_sol_new to be used at next timestep
[8]	1506	temp_sol_beg(:) = temp_sol_new(:)
	1507
[4682]	1508	END SUBROUTINE thermosoil_energy_diag
[8]	1509
[2222]	1510
	1511
	1512	!! ================================================================================================================================
	1513	!! SUBROUTINE : thermosoil_readjust
	1514	!!
	1515	!>\BRIEF
	1516	!!
	1517	!! DESCRIPTION : Energy conservation : Correction to make sure that the same latent heat is released and
	1518	!! consumed during freezing and thawing
	1519	!!
	1520	!! RECENT CHANGE(S) : None
	1521	!!
	1522	!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis),
	1523	!!
	1524	!! REFERENCE(S) :
	1525	!!
	1526	!! FLOWCHART : None
	1527	!! \n
	1528	!_ ================================================================================================================================
	1529
	1530	SUBROUTINE thermosoil_readjust(kjpindex, ptn)
	1531
	1532	!! 0. Variables and parameter declaration
	1533
	1534	!! 0.1 Input variables
	1535	INTEGER(i_std), INTENT(in) :: kjpindex
[2956]	1536
[2222]	1537	!! 0.2 Modified variables
	1538	REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(inout) :: ptn
	1539
	1540	!! 0.3 Local variables
	1541	INTEGER(i_std) :: ji, jg
[2957]	1542	INTEGER(i_std) :: lev3m !! Closest interface level to 3m
[2222]	1543	REAL(r_std) :: ptn_tmp
	1544
[2957]	1545	! The energy is spread over the layers down to approximatly 3m
	1546	! Find the closest level to 3m. It can be below or above 3m.
	1547	lev3m=MINLOC(ABS(zlt(:)-3.0),dim=1)
	1548	IF (printlev >= 3) WRITE(numout,*) 'In thermosoil_adjust: lev3m=',lev3m, ' zlt(lev3m)=', zlt(lev3m)
	1549
[2222]	1550	DO jg=1, ngrnd
	1551	DO ji=1, kjpindex
	1552	! All soil latent energy is put into e_soil_lat(ji)
	1553	! because the variable soil layers make it difficult to keep track of all
	1554	! layers in this version
	1555	! NOTE : pcapa has unit J/K/m3 and pcappa_supp has J/K
	1556	e_soil_lat(ji)=e_soil_lat(ji)+pcappa_supp(ji,jg)*(ptn(ji,jg)-ptn_beg(ji,jg))
	1557	END DO
[2957]	1558	END DO
[2222]	1559
	1560	DO ji=1, kjpindex
	1561	IF (e_soil_lat(ji).GT.min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
[2957]	1562	! The soil is thawed: we spread the excess of energy over the uppermost lev3m levels
[2222]	1563	! Here we increase the temperatures
[2957]	1564	DO jg=1, lev3m
[2222]	1565	ptn_tmp=ptn(ji,jg)
	1566
[2957]	1567	ptn(ji,jg)=ptn(ji,jg)+MIN(e_soil_lat(ji)/pcapa(ji,jg)/zlt(lev3m), 0.5)
[2942]	1568	e_soil_lat(ji)=e_soil_lat(ji)-(ptn(ji,jg)-ptn_tmp)pcapa(ji,jg)dlt(jg)
[2222]	1569	ENDDO
	1570	ELSE IF (e_soil_lat(ji).LT.-min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
	1571	! The soil is thawed
	1572	! Here we decrease the temperatures
[2957]	1573	DO jg=1, lev3m
[2222]	1574	ptn_tmp=ptn(ji,jg)
[2957]	1575	ptn(ji,jg)=MAX(ZeroCelsius+fr_dT/2., ptn_tmp+e_soil_lat(ji)/pcapa(ji,jg)/zlt(lev3m))
[2942]	1576	e_soil_lat(ji)=e_soil_lat(ji)+(ptn_tmp-ptn(ji,jg))pcapa(ji,jg)dlt(jg)
[2222]	1577	END DO
	1578	END IF
	1579	END DO
	1580
	1581	END SUBROUTINE thermosoil_readjust
	1582
	1583	!-------------------------------------------------------------------
	1584
	1585
	1586
	1587	!! ================================================================================================================================
	1588	!! SUBROUTINE : thermosoil_getdiff
	1589	!!
[3059]	1590	!>\BRIEF Computes soil and snow heat capacity and conductivity
[2222]	1591	!!
	1592	!! DESCRIPTION : Computation of the soil thermal properties; snow properties are also accounted for
	1593	!!
	1594	!! RECENT CHANGE(S) : None
	1595	!!
[3059]	1596	!! MAIN OUTPUT VARIABLE(S):
[2222]	1597	!!
	1598	!! REFERENCE(S) :
	1599	!!
	1600	!! FLOWCHART : None
	1601	!! \n
	1602	!_ ================================================================================================================================
	1603
[7206]	1604	SUBROUTINE thermosoil_getdiff( kjpindex, snow, ptn, mcs, njsc, snowrho, snowtemp, pb )
[2222]	1605
	1606	!! 0. Variables and parameter declaration
	1607
	1608	!! 0.1 Input variables
	1609	INTEGER(i_std),INTENT(in) :: kjpindex
	1610	REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow !! Snow mass
[2922]	1611	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
[7206]	1612	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3)
[3410]	1613	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
	1614	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
[7206]	1615	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
[3059]	1616	REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(in) :: ptn !! Soil temperature profile
[2222]	1617
	1618	!! 0.3 Local variables
	1619	REAL :: xx !! Unfrozen fraction of the soil
	1620	REAL(r_std), DIMENSION(kjpindex) :: snow_h
[7199]	1621	REAL(r_std), DIMENSION(kjpindex,ngrnd) :: zx1, zx2
[2222]	1622	INTEGER :: ji,jg
[2922]	1623	INTEGER :: jst
[4820]	1624	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! soil heat capacity (J/m3/K)
	1625	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_spec !! SPECIFIC soil heat capacity (J/kg/K)
	1626	REAL(r_std) :: rho_tot !! Soil density (kg/m3)
[2222]	1627
[2922]	1628	pcapa_tmp(:,:) = 0.0
[2222]	1629
[3059]	1630	!! Computes soil heat capacity and conductivity
[2222]	1631	DO ji = 1,kjpindex
	1632
[5470]	1633	! Since explicitsnow module is implemented set zx1=0 and zx2=1
	1634	zx1(ji,:) = 0.
	1635	zx2(ji,:) = 1.
	1636
[2222]	1637	DO jg = 1, ngrnd
[2922]	1638	jst = njsc(ji)
[7206]	1639	pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
[2222]	1640	!
[4820]	1641	! 2. Calculate volumetric heat capacity with allowance for permafrost
[2222]	1642	! 2.1. soil heat capacity depending on temperature and humidity
[4820]	1643	! For SP6MIP we also diagnose a specific heat capacity (pcapa_spec),
	1644	! which requires to calculate the total density of the soil (rho_tot), always >> 0
[2222]	1645
	1646	IF (ptn(ji,jg) .LT. ZeroCelsius-fr_dT/2.) THEN
	1647	! frozen soil
	1648	profil_froz(ji,jg) = 1.
	1649	pcappa_supp(ji,jg)= 0.
[7206]	1650	pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + so_capa_ice(ji) * tmc_layt(ji,jg) / mille / dlt(jg)
	1651	rho_tot = rho_soil * (1-mcs(ji)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg)
[4820]	1652	pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
[2222]	1653
	1654	ELSEIF (ptn(ji,jg) .GT. ZeroCelsius+fr_dT/2.) THEN
	1655	! unfrozen soil
[2922]	1656	pcapa(ji, jg) = pcapa_tmp(ji, jg)
[2222]	1657	profil_froz(ji,jg) = 0.
	1658	pcappa_supp(ji,jg)= 0.
[7206]	1659	rho_tot = rho_soil * (1-mcs(ji)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg)
[4820]	1660	pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
[2222]	1661	ELSE
	1662	! xx is the unfrozen fraction of soil water
	1663	xx = (ptn(ji,jg)-(ZeroCelsius-fr_dT/2.)) / fr_dT
	1664	profil_froz(ji,jg) = (1. - xx)
	1665
[5150]	1666	IF (ok_freeze_thaw_latent_heat) THEN
[7206]	1667	pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + &
[3016]	1668	water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
[7206]	1669	so_capa_ice(ji) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + &
	1670	shum_ngrnd_perma(ji,jg)mcs(ji)lhf*rho_water/fr_dT
[5150]	1671	ELSE
[7206]	1672	pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + &
[5150]	1673	water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
[7206]	1674	so_capa_ice(ji) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
[5150]	1675	ENDIF
	1676
[7206]	1677	rho_tot = rho_soil* (1-mcs(ji)) + &
[4820]	1678	rho_water * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
	1679	rho_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
	1680	pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
	1681
[7206]	1682	pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)mcs(ji)lhfrho_water/fr_dTzx2(ji,jg)*dlt(jg)
[2436]	1683
[2222]	1684	ENDIF
	1685
[4767]	1686	!
[2222]	1687	! 2.2. Take into account the snow and soil fractions in the layer
	1688	!
[2953]	1689	pcapa(ji,jg) = zx1(ji,jg) * sn_capa + zx2(ji,jg) * pcapa(ji,jg)
[2222]	1690
	1691	!
	1692	! 2.3. Calculate the heat capacity for energy conservation check
	1693	IF ( zx1(ji,jg).GT.0. ) THEN
	1694	pcapa_en(ji,jg) = sn_capa
	1695	ELSE
	1696	pcapa_en(ji,jg) = pcapa(ji,jg)
	1697	ENDIF
[3059]	1698
[2222]	1699	END DO
[4820]	1700	ENDDO
[2222]	1701
[4820]	1702	! Output the specific heat capcaity for SP-MIP
	1703	CALL xios_orchidee_send_field("pcapa_spec",pcapa_spec)
	1704
[2956]	1705	!
	1706	! 3. Calculate the heat conductivity with allowance for permafrost
	1707	!
[5150]	1708	IF (ok_freeze_thaw_latent_heat) THEN
[7206]	1709	CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt*(1-profil_froz), pkappa)
[5150]	1710	ELSE
[7206]	1711	CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt, pkappa)
[5150]	1712	ENDIF
[3059]	1713
	1714	!! Computes snow heat capacity and conductivity
	1715	DO ji = 1,kjpindex
	1716	pcapa_snow(ji,:) = snowrho(ji,:) * xci
	1717	pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
	1718	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
	1719	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
	1720	END DO
[2222]	1721
	1722	END SUBROUTINE thermosoil_getdiff
	1723
	1724
	1725	!! ================================================================================================================================
	1726	!! SUBROUTINE : thermosoil_getdiff_old_thermix_without_snow
	1727	!!
[3059]	1728	!>\BRIEF Computes soil and snow heat capacity and conductivity
[2222]	1729	!!
[5470]	1730	!! DESCRIPTION : Calculations of soil and snow thermal properties without effect of freezing.
[2222]	1731	!!
	1732	!!
	1733	!! RECENT CHANGE(S) : None
	1734	!!
	1735	!! MAIN OUTPUT VARIABLE(S):
	1736	!!
	1737	!! REFERENCE(S) :
	1738	!!
	1739	!! FLOWCHART : None
	1740	!! \n
	1741	!_ ================================================================================================================================
	1742
[7206]	1743	SUBROUTINE thermosoil_getdiff_old_thermix_without_snow( kjpindex, mcs, njsc, snowrho, snowtemp, pb )
[2222]	1744
	1745	!! 0. Variables and parameter declaration
	1746
	1747	!! 0.1 Input variables
	1748	INTEGER(i_std), INTENT(in) :: kjpindex
[2922]	1749	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
[7206]	1750	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3)
[3059]	1751	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
[3410]	1752	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
[7206]	1753	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa)
[2222]	1754
[3059]	1755
[2222]	1756	!! 0.1 Local variables
[2922]	1757	INTEGER(i_std) :: ji,jg, jst !! Index
	1758	REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! Soil heat capacity (J/m3/K)
[2222]	1759
[3059]	1760	!! Computes soil heat capacity and conductivity
[2222]	1761	DO jg = 1,ngrnd
	1762	DO ji = 1,kjpindex
[2922]	1763	jst = njsc(ji)
[7206]	1764	pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
[2922]	1765	pcapa(ji,jg) = pcapa_tmp(ji, jg)
	1766	pcapa_en(ji,jg) = pcapa_tmp(ji, jg)
[2222]	1767	ENDDO
	1768	ENDDO
	1769
[7206]	1770	CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt, pkappa)
[2956]	1771
[2922]	1772	IF (brk_flag == 1) THEN
	1773	! Bedrock flag is activated
	1774	DO jg = ngrnd-1,ngrnd
	1775	DO ji = 1,kjpindex
	1776	pcapa(ji,jg) = brk_capa
	1777	pcapa_en(ji,jg) = brk_capa
	1778	pkappa(ji,jg) = brk_cond
	1779	ENDDO
	1780	ENDDO
	1781	ENDIF
	1782
[3059]	1783	!! Computes snow heat capacity and conductivity
	1784	DO ji = 1,kjpindex
	1785	pcapa_snow(ji,:) = snowrho(ji,:) * xci
	1786	pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2snowrho(ji,:)snowrho(ji,:)) + &
	1787	MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
	1788	ZSNOWTHRMCOND_CVAP)))(XP00/(pb(ji)100.)))
	1789	END DO
	1790
[3345]	1791	END SUBROUTINE thermosoil_getdiff_old_thermix_without_snow
[3059]	1792
[2222]	1793
	1794	!! ================================================================================================================================
[3850]	1795	!! SUBROUTINE : thermosoil_read_reftempfile
[2222]	1796	!!
	1797	!>\BRIEF
	1798	!!
[3345]	1799	!! DESCRIPTION : Read file with longterm soil temperature
[2222]	1800	!!
	1801	!!
	1802	!! RECENT CHANGE(S) : None
	1803	!!
	1804	!! MAIN OUTPUT VARIABLE(S): reftemp : Reference temerature
	1805	!!
	1806	!! REFERENCE(S) :
	1807	!!
	1808	!! FLOWCHART : None
	1809	!! \n
	1810	!_ ================================================================================================================================
[3850]	1811	SUBROUTINE thermosoil_read_reftempfile(kjpindex,lalo,reftemp)
[3345]	1812
[3831]	1813	USE interpweight
	1814
	1815	IMPLICIT NONE
	1816
[2222]	1817	!! 0. Variables and parameter declaration
	1818
	1819	!! 0.1 Input variables
	1820	INTEGER(i_std), INTENT(in) :: kjpindex
	1821	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo
	1822
[3345]	1823	!! 0.2 Output variables
	1824	REAL(r_std), DIMENSION(kjpindex, ngrnd), INTENT(out) :: reftemp
[2222]	1825
	1826	!! 0.3 Local variables
[3831]	1827	INTEGER(i_std) :: ib
[2222]	1828	CHARACTER(LEN=80) :: filename
[3831]	1829	REAL(r_std),DIMENSION(kjpindex) :: reftemp_file !! Horizontal temperature field interpolated from file [C]
	1830	INTEGER(i_std),DIMENSION(kjpindex,8) :: neighbours
	1831	REAL(r_std) :: vmin, vmax !! min/max values to use for the
	1832	!! renormalization
	1833	REAL(r_std), DIMENSION(kjpindex) :: areftemp !! Availability of data for the interpolation
	1834	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
	1835	!! the file
	1836	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
	1837	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
	1838	!! (variabletypevals(1) = -un, not used)
	1839	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
	1840	!! 'XYKindTime': Input values are kinds
	1841	!! of something with a temporal
	1842	!! evolution on the dx*dy matrix'
	1843	LOGICAL :: nonegative !! whether negative values should be removed
	1844	CHARACTER(LEN=50) :: maskingtype !! Type of masking
	1845	!! 'nomask': no-mask is applied
	1846	!! 'mbelow': take values below maskvals(1)
	1847	!! 'mabove': take values above maskvals(1)
	1848	!! 'msumrange': take values within 2 ranges;
	1849	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
	1850	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
	1851	!! (normalized by maskvals(3))
	1852	!! 'var': mask values are taken from a
	1853	!! variable inside the file (>0)
	1854	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
	1855	!! `maskingtype')
	1856	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
	1857	REAL(r_std) :: reftemp_norefinf
	1858	REAL(r_std) :: reftemp_default !! Default value
[3345]	1859
	1860	!Config Key = SOIL_REFTEMP_FILE
	1861	!Config Desc = File with climatological soil temperature
	1862	!Config If = READ_REFTEMP
	1863	!Config Def = reftemp.nc
	1864	!Config Help =
	1865	!Config Units = [FILE]
[2222]	1866	filename = 'reftemp.nc'
[3831]	1867	CALL getin_p('REFTEMP_FILE',filename)
[2222]	1868
[3831]	1869	variablename = 'temperature'
[2222]	1870
[4693]	1871	IF (printlev >= 1) WRITE(numout,*) "thermosoil_read_reftempfile: Read and interpolate file " &
	1872	// TRIM(filename) //" for variable " //TRIM(variablename)
[3831]	1873
[5364]	1874	IF (xios_interpolation) THEN
[3831]	1875
[5364]	1876	CALL xios_orchidee_recv_field('reftemp_interp',reftemp(:,1))
[3831]	1877
[5364]	1878	DO ib=1, kjpindex
	1879	reftemp(ib,:) = reftemp(ib,1) + ZeroCelsius
	1880	END DO
	1881	areftemp = 1.0
	1882	ELSE
[3831]	1883
	1884
[5364]	1885	! For this case there are not types/categories. We have 'only' a continuos field
	1886	! Assigning values to vmin, vmax
	1887
	1888	vmin = 0.
	1889	vmax = 9999.
[3831]	1890
[5364]	1891	! For this file we do not need neightbours!
	1892	neighbours = 0
	1893
	1894	!! Variables for interpweight
	1895	! Type of calculation of cell fractions
	1896	fractype = 'default'
	1897	! Name of the longitude and latitude in the input file
	1898	lonname = 'nav_lon'
	1899	latname = 'nav_lat'
	1900	! Default value when no value is get from input file
	1901	reftemp_default = 1.
	1902	! Reference value when no value is get from input file
	1903	reftemp_norefinf = 1.
	1904	! Should negative values be set to zero from input file?
	1905	nonegative = .FALSE.
	1906	! Type of mask to apply to the input data (see header for more details)
	1907	maskingtype = 'nomask'
	1908	! Values to use for the masking (here not used)
	1909	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
	1910	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
	1911	namemaskvar = ''
	1912
	1913	CALL interpweight_2Dcont(kjpindex, 0, 0, lalo, resolution, neighbours, &
	1914	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
	1915	maskvals, namemaskvar, -1, fractype, reftemp_default, reftemp_norefinf, &
	1916	reftemp_file, areftemp)
	1917	IF (printlev >= 5) WRITE(numout,*)' thermosoil_read_reftempfile after interpweight_2Dcont'
	1918
	1919	! Copy reftemp_file temperature to all ground levels and transform into Kelvin
	1920	DO ib=1, kjpindex
	1921	reftemp(ib, :) = reftemp_file(ib)+ZeroCelsius
	1922	END DO
[3831]	1923
[5364]	1924	END IF
	1925
[3831]	1926	! Write diagnostics
	1927	CALL xios_orchidee_send_field("areftemp",areftemp)
[5364]	1928
[3850]	1929	END SUBROUTINE thermosoil_read_reftempfile
[2222]	1930
[8]	1931	END MODULE thermosoil

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

Download in other formats: