Context Navigation

enerbil.f90 @ 7325

Last change on this file since 7325 was 7265, checked in by agnes.ducharne, 3 years ago
Integrated r5705 (solay now in meters in output files), and small changes with no impact on code (r6220, r6565, r6567) from the trunk. Checked with a 5d run.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 82.9 KB

Rev	Line
[947]	1	! ==============================================================================================================================
	2	! MODULE : enerbil
	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 This module computes the energy balance on
[2581]	10	!! continental surfaces. The module contains the following subroutines: enerbil_initialize, enerbil_main,
[5470]	11	!! enerbil_clear, enerbil_begin, enerbil_surftemp, enerbil_flux and enerbil_evapveg
[8]	12	!!
[947]	13	!!\n DESCRIPTION :
	14	!! \n
	15	!! \latexonly
	16	!! \input{enerbil_intro2.tex}
	17	!! \endlatexonly
[8]	18	!!
[947]	19	!! IMPORTANT NOTE: The coefficients A and B are defined differently than those in the referenced
	20	!! literature and from those in the code and documentation of the atmospheric model LMDZ. For the
	21	!! avoidance of doubt, the coefficients as described here always refer to the ORCHIDEE coefficients,
	22	!! and are denoted as such (with the marker: ORC). The re-definition of the coefficients takes place
	23	!! within LMDZ before they are passed to ORCHIDEE. The following sequence of expressions is to be
	24	!! found within the LMDZ module 'surf_land_orchidee':\n
	25	!!
	26	!! \latexonly
	27	!! \input{surflandLMDZ1.tex}
	28	!! \input{surflandLMDZ2.tex}
	29	!! \input{surflandLMDZ3.tex}
	30	!! \input{surflandLMDZ4.tex}
	31	!! \endlatexonly
	32	!! \n
	33	!!
	34	!! \latexonly
	35	!! \input{enerbil_symbols.tex}
	36	!! \endlatexonly
	37	!!
	38	!! RECENT CHANGE(S) : None
	39	!!
	40	!! REFERENCE(S) : None
	41	!!
	42	!! SVN :
	43	!! $HeadURL$
	44	!! $Date$
	45	!! $Revision$
	46	!! \n
	47	!_ ================================================================================================================================
	48
[8]	49	MODULE enerbil
	50
	51	! routines called : restput, restget
	52	!
[1392]	53	! modules used
[8]	54	USE ioipsl
[1788]	55	USE xios_orchidee
[1392]	56	USE ioipsl_para
[8]	57	USE constantes
[4646]	58	USE time, ONLY : one_day, dt_sechiba
[511]	59	USE pft_parameters
	60	USE qsat_moisture
[4281]	61	USE sechiba_io_p
[2425]	62	USE constantes_soil
[2222]	63	USE explicitsnow
[1078]	64
[8]	65	IMPLICIT NONE
	66
	67	PRIVATE
[5470]	68	PUBLIC :: enerbil_main, enerbil_initialize, enerbil_finalize, enerbil_clear
[8]	69
	70	! variables used inside enerbil module : declaration and initialisation
[947]	71
[8]	72	! one dimension array allocated, computed and used in enerbil module exclusively
[947]	73	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: psold !! Old surface dry static energy (J kg^{-1})
[1078]	74	!$OMP THREADPRIVATE(psold)
[947]	75	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsol_sat !! Saturated specific humudity for old temperature (kg kg^{-1})
[1078]	76	!$OMP THREADPRIVATE(qsol_sat)
[947]	77	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: pdqsold !! Deriv. of saturated specific humidity at old temp
	78	!! (kg (kg s)^{-1})
[1078]	79	!$OMP THREADPRIVATE(pdqsold)
[947]	80	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: psnew !! New surface static energy (J kg^{-1})
[1078]	81	!$OMP THREADPRIVATE(psnew)
[947]	82	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsol_sat_new !! New saturated surface air moisture (kg kg^{-1})
[1078]	83	!$OMP THREADPRIVATE(qsol_sat_new)
[947]	84	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: netrad !! Net radiation (W m^{-2})
[1078]	85	!$OMP THREADPRIVATE(netrad)
[947]	86	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwabs !! LW radiation absorbed by the surface (W m^{-2})
[1078]	87	!$OMP THREADPRIVATE(lwabs)
[947]	88	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwup !! Long-wave up radiation (W m^{-2})
[1078]	89	!$OMP THREADPRIVATE(lwup)
[947]	90	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwnet !! Net Long-wave radiation (W m^{-2})
[1078]	91	!$OMP THREADPRIVATE(lwnet)
[947]	92	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: fluxsubli !! Energy of sublimation (mm day^{-1})
[1078]	93	!$OMP THREADPRIVATE(fluxsubli)
[947]	94	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsat_air !! Air saturated specific humidity (kg kg^{-1})
[1078]	95	!$OMP THREADPRIVATE(qsat_air)
[947]	96	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tair !! Air temperature (K)
[1078]	97	!$OMP THREADPRIVATE(tair)
[947]	98	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: q_sol_pot !! Potential surface humidity
[1078]	99	!$OMP THREADPRIVATE(q_sol_pot)
[2650]	100
[947]	101
[2581]	102	CONTAINS
	103	!! =============================================================================================================================
	104	!! SUBROUTINE: enerbil_initialize
	105	!!
	106	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
	107	!!
	108	!! DESCRIPTION: The purpose of this module is, firstly, to allocate space
	109	!! in memory for key variables within the 'enerbil' module. The second task is to retrieve previous data
	110	!! from the restart file. If the variables are not in restart file, default initialization is done.
	111	!!
	112	!! RECENT CHANGE(S): None
	113	!!
	114	!! MAIN OUTPUT VARIABLE(S):
	115	!!
	116	!! REFERENCE(S): None
	117	!!
	118	!! FLOWCHART: None
	119	!! \n
	120	!_ ==============================================================================================================================
	121
	122	SUBROUTINE enerbil_initialize (kjit, kjpindex, index, rest_id, &
[4384]	123	qair, &
	124	temp_sol, temp_sol_new, tsol_rad, &
[2581]	125	evapot, evapot_corr, qsurf, fluxsens, &
	126	fluxlat, vevapp )
	127
	128	!! 0 Variable and parameter description
	129	!! 0.1 Input variables
	130	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
	131	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
	132	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
	133	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (-)
	134	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1})
	135
	136	!! 0.2 Output variables
	137	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol !! Soil temperature (K)
	138	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K)
	139	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Tsol_rad (W m^{-2})
[3406]	140	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: evapot !! Soil Potential Evaporation (mm/tstep)
	141	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep)
[2581]	142	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! Surface specific humidity (kg kg^{-1})
	143	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2})
	144	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2})
	145	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1})
	146
	147
	148	!! 0.4 Local variables
	149	INTEGER(i_std) :: ier
	150
	151	!_ ================================================================================================================================
[947]	152
[2581]	153	IF (printlev>=3) WRITE (numout,*) 'enerbil_initialize : Start initillization'
	154
	155	!! 1. Allocate module variables
	156	ALLOCATE (psold(kjpindex),stat=ier)
	157	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','')
	158
	159	ALLOCATE (qsol_sat(kjpindex),stat=ier)
	160	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','')
	161
	162	ALLOCATE (pdqsold(kjpindex),stat=ier)
	163	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','')
	164
	165	ALLOCATE (psnew(kjpindex),stat=ier)
	166	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','')
	167
	168	ALLOCATE (qsol_sat_new(kjpindex),stat=ier)
	169	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable qsol_sat_new','','')
	170
	171	ALLOCATE (netrad(kjpindex),stat=ier)
	172	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable netrad','','')
	173
	174	ALLOCATE (lwabs(kjpindex),stat=ier)
	175	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable lwabs','','')
	176
	177	ALLOCATE (lwup(kjpindex),stat=ier)
	178	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable lwup','','')
	179
	180	ALLOCATE (lwnet(kjpindex),stat=ier)
	181	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable lwnet','','')
	182
	183	ALLOCATE (fluxsubli(kjpindex),stat=ier)
	184	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable fluxsubli','','')
	185
	186	ALLOCATE (qsat_air(kjpindex),stat=ier)
	187	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable qsat_air','','')
	188
	189	ALLOCATE (tair(kjpindex),stat=ier)
	190	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable tair','','')
	191
	192	ALLOCATE (q_sol_pot(kjpindex),stat=ier)
	193	IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable q_sol_pot','','')
	194
	195
	196	!! 2. Initialize variables from restart file or by default values
	197	!! The variables read are: temp_sol (surface temperature), qsurf (near surface specific humidity),
	198	!! evapot (soil potential evaporation), evapot_corr (corrected soil potential evaporation), tsolrad
	199	!! (radiative surface temperature), evapora (evaporation), fluxlat (latent heat flux), fluxsens
	200	!! (sensible heat flux) and temp_sol_new (new surface temperature).
	201	IF (printlev>=3) WRITE (numout,*) 'Read a restart file for ENERBIL variables'
	202
	203	CALL ioconf_setatt_p('UNITS', 'K')
	204	CALL ioconf_setatt_p('LONG_NAME','Surface temperature')
	205	CALL restget_p (rest_id, 'temp_sol', nbp_glo, 1, 1, kjit, .TRUE., temp_sol, "gather", nbp_glo, index_g)
	206
	207	!Config Key = ENERBIL_TSURF
	208	!Config Desc = Initial temperature if not found in restart
	209	!Config If = OK_SECHIBA
	210	!Config Def = 280.
	211	!Config Help = The initial value of surface temperature if its value is not found
	212	!Config in the restart file. This should only be used if the model is
	213	!Config started without a restart file.
	214	!Config Units = Kelvin [K]
	215	CALL setvar_p (temp_sol, val_exp,'ENERBIL_TSURF', 280._r_std)
	216
[2861]	217	! Initialize temp_sol_new with temp_sol. These variables are always equal in the beginning of a new time step.
	218	temp_sol_new(:) = temp_sol(:)
	219
[2581]	220	CALL ioconf_setatt_p('UNITS', 'g/g')
	221	CALL ioconf_setatt_p('LONG_NAME','near surface specific humidity')
	222	CALL restget_p (rest_id, 'qsurf', nbp_glo, 1, 1, kjit, .TRUE., qsurf, "gather", nbp_glo, index_g)
	223	IF ( ALL( qsurf(:) .EQ. val_exp ) ) THEN
	224	qsurf(:) = qair(:)
	225	ENDIF
	226
	227	CALL ioconf_setatt_p('UNITS', 'mm day^{-1}')
	228	CALL ioconf_setatt_p('LONG_NAME','Soil Potential Evaporation')
	229	CALL restget_p (rest_id, 'evapot', nbp_glo, 1, 1, kjit, .TRUE., evapot, "gather", nbp_glo, index_g)
	230	CALL setvar_p (evapot, val_exp, 'ENERBIL_EVAPOT', zero)
	231
	232	CALL ioconf_setatt_p('UNITS', 'mm day^{-1}')
	233	CALL ioconf_setatt_p('LONG_NAME','Corrected Soil Potential Evaporation')
	234	CALL restget_p (rest_id, 'evapot_corr', nbp_glo, 1, 1, kjit, .TRUE., evapot_corr, "gather", nbp_glo, index_g)
	235	!Config Key = ENERBIL_EVAPOT
	236	!Config Desc = Initial Soil Potential Evaporation
	237	!Config If = OK_SECHIBA
	238	!Config Def = 0.0
	239	!Config Help = The initial value of soil potential evaporation if its value
	240	!Config is not found in the restart file. This should only be used if
	241	!Config the model is started without a restart file.
	242	!Config Units =
	243	CALL setvar_p (evapot_corr, val_exp, 'ENERBIL_EVAPOT', zero)
	244
	245	CALL ioconf_setatt_p('UNITS', 'K')
	246	CALL ioconf_setatt_p('LONG_NAME','Radiative surface temperature')
	247	CALL restget_p (rest_id, 'tsolrad', nbp_glo, 1, 1, kjit, .TRUE., tsol_rad, "gather", nbp_glo, index_g)
	248	IF ( ALL( tsol_rad(:) .EQ. val_exp ) ) THEN
	249	tsol_rad(:) = temp_sol(:)
	250	ENDIF
	251
	252	!! 1.3 Set the fluxes so that we have something reasonable and not NaN on some machines
	253	CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt')
	254	CALL ioconf_setatt_p('LONG_NAME','Evaporation')
	255	CALL restget_p (rest_id, 'evapora', nbp_glo, 1, 1, kjit, .TRUE., vevapp, "gather", nbp_glo, index_g)
	256	IF ( ALL( vevapp(:) .EQ. val_exp ) ) THEN
	257	vevapp(:) = zero
	258	ENDIF
	259
	260	CALL ioconf_setatt_p('UNITS', 'W/m^2')
	261	CALL ioconf_setatt_p('LONG_NAME','Latent heat flux')
	262	CALL restget_p (rest_id, 'fluxlat', nbp_glo, 1, 1, kjit, .TRUE., fluxlat, "gather", nbp_glo, index_g)
	263	IF ( ALL( fluxlat(:) .EQ. val_exp ) ) THEN
	264	fluxlat(:) = zero
	265	ENDIF
	266
	267	CALL ioconf_setatt_p('UNITS', 'W/m^2')
	268	CALL ioconf_setatt_p('LONG_NAME','Sensible heat flux')
	269	CALL restget_p (rest_id, 'fluxsens', nbp_glo, 1, 1, kjit, .TRUE., fluxsens, "gather", nbp_glo, index_g)
	270	IF ( ALL( fluxsens(:) .EQ. val_exp ) ) THEN
	271	fluxsens(:) = zero
	272	ENDIF
	273
	274
	275	CALL ioconf_setatt_p('UNITS', 'kg/m^2')
	276	CALL ioconf_setatt_p('LONG_NAME','Potential saturated surface humidity')
	277	CALL restget_p (rest_id, 'qsolpot', nbp_glo, 1, 1, kjit, .TRUE., q_sol_pot, "gather", nbp_glo, index_g)
	278	IF ( ALL( q_sol_pot(:) .EQ. val_exp ) ) THEN
	279	q_sol_pot = qsurf
	280	ENDIF
	281
	282	END SUBROUTINE enerbil_initialize
[947]	283
[2581]	284
	285
[947]	286	!! ===========================================================================================================================
	287	!! SUBROUTINE : enerbil_main
	288	!!
	289	!!
	290	!>\BRIEF Calls each part of the energy budget calculation in sequence
	291	!!
	292	!! DESCRIPTION :
	293	!! This is the main routine for the 'enerbil' module. It is called
	294	!! once during the initialisation of ORCHIDEE, and then once for each time step. It is called a final
	295	!! time at the culmination of the last time step, for the writing of a restart file.\n
	296	!!
	297	!! The algorithm first calls 'enerbil_begin' for initialisation, followed by 'enerbil_surftemp' to
	298	!! calculate the surface static energy and the saturated surface humidity for the new time step.
	299	!! Next is the module 'enerbil_flux' which calculates the new surface temperature, the net radiation in
	300	!! the surface layer, the total evaporation and the latent and sensible heat fluxes. Finally comes
	301	!! 'enerbil_evapveg', which calculates the evaporation and transpiration from the vegetation.\n
[8]	302
[947]	303	!! \n
[8]	304	!!
[947]	305	!! RECENT CHANGE(S): None
[8]	306	!!
[947]	307	!! MAIN OUTPUT VARIABLE(S) : evapot, evapot_corr, temp_sol, qsurf, fluxsens, fluxlat, tsol_rad,
[4384]	308	!! vevapp, temp_sol_new, vevapnu, vevapsno, transpir, vevapwet
[8]	309	!!
[947]	310	!! REFERENCE(S) :
	311	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
	312	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
	313	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
	314	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
	315	!! Circulation Model. Journal of Climate, 6, pp.248-273
	316	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
	317	!! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11):
	318	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
	319	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
	320	!! sur le cycle de l'eau, PhD Thesis, available (22/12/11):
	321	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
	322	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
	323	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
	324	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
	325	!! general circulation models. Global and Planetary Change, 19, pp.261-276
	326	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
	327	!! Interscience Publishers\n
	328	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
	329	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
	330	!!
	331	!! FLOWCHART :
	332	!! \latexonly
	333	!! \includegraphics[scale=0.5]{enerbil_main_flowchart.png}
	334	!! \endlatexonly
	335	!! \n
	336	!_ ==============================================================================================================================
	337
[2591]	338	SUBROUTINE enerbil_main (kjit, kjpindex, &
[947]	339	& index, indexveg, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, &
[4197]	340	& qair, peqAcoef, peqBcoef, pb, rau, vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
[2031]	341	& emis, soilflx, soilcap, q_cdrag, humrel, fluxsens, fluxlat, &
[4384]	342	& vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, temp_sol, tsol_rad, &
[2222]	343	& temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id, &
[3959]	344	& precip_rain, pgflux, snowdz, temp_sol_add)
[2222]	345
[8]	346
[947]	347	!! 0 Variable and parameter description
	348
	349	!! 0.1 Input variables
	350
	351	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
	352	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
	353	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier (-)
	354	INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier (-)
	355	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
	356	INTEGER(i_std),DIMENSION(kjpindex*nvm), INTENT(in) :: indexveg !! Indeces of the points on the 3D map
	357	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
	358	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2})
	359	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2})
	360	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (?? J)
	361	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
	362	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
	363	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
	364	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef (see note)
	365	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef (see note)
	366	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1})
	367	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef (see note)
	368	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef (see note)
[2031]	369	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
[947]	370	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3})
	371	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-)
	372	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-)
	373	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta4 !! Bare soil resistance (-)
	374	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance
	375	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-)
[3020]	376	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Effective ground heat flux including both snow and soil (W m^{-2})
	377	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity including both snow and soil (J K^{-1])
[4146]	378	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
[2031]	379
[947]	380	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1)
	381	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Interception resistance (-)
	382	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Vegetation resistance (-)
	383	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3pot !! Vegetation resistance for potential transpiration
[2222]	384	REAL(r_std),DIMENSION (kjpindex),INTENT (in) :: precip_rain !! Rainfall
[3883]	385	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer
[2031]	386
[947]	387	!! 0.2 Output variables
	388
	389	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation (mm day^{-1})
	390	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation (mm day^{-1})
	391	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapflo !! Floodplains evaporation
	392	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration (mm day^{-1})
	393	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpot !! Potential transpiration
	394	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception (mm day^{-1})
[2222]	395	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K)
[2650]	396	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
[3883]	397
[947]	398	!! 0.3 Modified variables
	399
[3406]	400	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation (mm/tstep)
	401	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep)
[947]	402	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol !! Soil temperature (K)
	403	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: qsurf !! Surface specific humidity (kg kg^{-1})
[2868]	404	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2})
	405	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2})
[947]	406	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tsol_rad !! Tsol_rad (W m^{-2})
[2868]	407	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1})
[2222]	408	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: pgflux !! Net energy into snowpack(W/m^2)
[2031]	409
[947]	410	!! 0.4 Local variables
[3883]	411
[947]	412	REAL(r_std),DIMENSION (kjpindex) :: epot_air_new, qair_new
[3904]	413	REAL(r_std),DIMENSION (kjpindex) :: diffevap !! Difference betwence vevapp and composing fluxes (Kg/m^2/s)
[2222]	414	INTEGER(i_std) :: ji,ii,iv
[947]	415
[2222]	416
[947]	417	!_ ================================================================================================================================
	418
[2581]	419	!! 1. Computes some initialisation variables
[947]	420
	421	! Computes some initialisation variables psold (the old surface static energy), qsol_sat
	422	! (the saturated surface humidity) and pdqsold (the derivative of the saturated surface humidity,
	423	! calculated with respect to the surface temperature at the 'old' timestep)
[2650]	424	CALL enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis)
[8]	425
[2222]	426
[2581]	427	!! 2. Computes psnew (the surface static energy at the 'new' timestep)
[947]	428
	429	! Computes psnew (the surface static energy at the 'new' timestep), qsol_sat_new (the surface
	430	! saturated humidity at the 'new' timestep), temp_sol_new (the surface temperature at the 'new'
	431	! timestep), qair_new (the lowest atmospheric humidity at the 'new' timestep) and epot_air_new
	432	! (the lowest atmospheric evaporation potential at the 'new' timestep)
[2591]	433	CALL enerbil_surftemp (kjpindex, zlev, emis, &
[8]	434	& epot_air, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,&
[4197]	435	& vbeta1, vbeta5, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, &
[3959]	436	& qair_new, epot_air_new)
[8]	437
[2222]	438
[4212]	439	!! 3. Diagnose components of the energy budget
[947]	440
	441	! Diagnoses lwup (upwards long wave radiation), lwnet (net long wave radiation), tsol_rad (radiative
	442	! ground temperature), netrad (net radiation), qsurf (surface humidity), vevapp (total evaporation),
	443	! evapot (evaporation potential), evapot_corr (evaporation potential correction factor),
	444	! fluxlat (latent heat flux), fluxsubli (sublimination heat flux) and fluxsens (sensible heat flux).
	445
[4197]	446	CALL enerbil_flux (kjpindex, emis, temp_sol, rau, u, v, q_cdrag, vbeta, vbeta1, vbeta5, &
[8]	447	& qair_new, epot_air_new, psnew, qsurf, &
	448	& fluxsens , fluxlat , fluxsubli, vevapp, temp_sol_new, lwdown, swnet, &
[2222]	449	& lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr, &
[3020]	450	& precip_rain,snowdz,temp_air,pgflux, soilcap, temp_sol_add)
[8]	451
[4212]	452
	453	!! 4. Diagnoses the values for evaporation and transpiration
[947]	454
	455	! Diagnoses the values for evaporation and transpiration: vevapsno (snow evaporation), vevapnu
[2031]	456	! (bare soil evaporation), transpir (transpiration) and vevapwet (interception)
[2591]	457	CALL enerbil_evapveg (kjpindex, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
[2031]	458	& rau, u, v, q_cdrag, qair_new, humrel, vevapsno, vevapnu , vevapflo, vevapwet, &
[3959]	459	& transpir, transpot, evapot)
[8]	460
[2222]	461
[947]	462
[4384]	463	!! 5. Write diagnosics
[8]	464
[1788]	465	CALL xios_orchidee_send_field("netrad",netrad)
[3406]	466	CALL xios_orchidee_send_field("evapot",evapot/dt_sechiba)
	467	CALL xios_orchidee_send_field("evapot_corr",evapot_corr/dt_sechiba)
[1788]	468	CALL xios_orchidee_send_field("lwdown",lwabs)
	469	CALL xios_orchidee_send_field("lwnet",lwnet)
	470	CALL xios_orchidee_send_field("Qv",fluxsubli)
	471
[3904]	472	DO ji=1,kjpindex
	473	diffevap(ji) = vevapp(ji) - ( SUM(vevapwet(ji,:)) + &
	474	SUM(transpir(ji,:)) + vevapnu(ji) + vevapsno(ji) + vevapflo(ji) )
	475	ENDDO
	476	CALL xios_orchidee_send_field("diffevap",diffevap/dt_sechiba) ! mm/s
	477
[8]	478	IF ( .NOT. almaoutput ) THEN
[1078]	479	CALL histwrite_p(hist_id, 'netrad', kjit, netrad, kjpindex, index)
	480	CALL histwrite_p(hist_id, 'evapot', kjit, evapot, kjpindex, index)
	481	CALL histwrite_p(hist_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index)
	482	CALL histwrite_p(hist_id, 'lwdown', kjit, lwabs, kjpindex, index)
	483	CALL histwrite_p(hist_id, 'lwnet', kjit, lwnet, kjpindex, index)
[8]	484	IF ( hist2_id > 0 ) THEN
[1078]	485	CALL histwrite_p(hist2_id, 'netrad', kjit, netrad, kjpindex, index)
	486	CALL histwrite_p(hist2_id, 'evapot', kjit, evapot, kjpindex, index)
	487	CALL histwrite_p(hist2_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index)
	488	CALL histwrite_p(hist2_id, 'lwdown', kjit, lwabs, kjpindex, index)
	489	CALL histwrite_p(hist2_id, 'lwnet', kjit, lwnet, kjpindex, index)
[8]	490	ENDIF
	491	ELSE
[1078]	492	CALL histwrite_p(hist_id, 'LWnet', kjit, lwnet, kjpindex, index)
	493	CALL histwrite_p(hist_id, 'Qv', kjit, fluxsubli, kjpindex, index)
	494	CALL histwrite_p(hist_id, 'PotEvap', kjit, evapot_corr, kjpindex, index)
	495	CALL histwrite_p(hist_id, 'PotEvapOld', kjit, evapot, kjpindex, index)
[8]	496	IF ( hist2_id > 0 ) THEN
[1078]	497	CALL histwrite_p(hist2_id, 'LWnet', kjit, lwnet, kjpindex, index)
	498	CALL histwrite_p(hist2_id, 'Qv', kjit, fluxsubli, kjpindex, index)
	499	CALL histwrite_p(hist2_id, 'PotEvap', kjit, evapot_corr, kjpindex, index)
[8]	500	ENDIF
	501	ENDIF
	502
[2348]	503	IF (printlev>=3) WRITE (numout,*) ' enerbil_main Done '
[8]	504
	505	END SUBROUTINE enerbil_main
	506
[947]	507
[2581]	508	!! =============================================================================================================================
	509	!! SUBROUTINE: enerbil_finalize
	510	!!
	511	!>\BRIEF Write to restart file
	512	!!
	513	!! DESCRIPTION: This subroutine writes the module variables and variables calculated in enerbil
	514	!! to restart file
	515	!!
	516	!! RECENT CHANGE(S): None
	517	!!
	518	!! REFERENCE(S): None
	519	!!
	520	!! FLOWCHART: None
	521	!! \n
	522	!_ ==============================================================================================================================
	523	SUBROUTINE enerbil_finalize (kjit, kjpindex, rest_id, &
	524	evapot, evapot_corr, temp_sol, tsol_rad, &
	525	qsurf, fluxsens, fluxlat, vevapp )
	526
	527	!! 0 Variable and parameter description
[947]	528	!! 0.1 Input variables
[2581]	529	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
	530	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
	531	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (-)
[3406]	532	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm/tstep)
	533	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep)
[2581]	534	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature (K)
	535	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tsol_rad !! Tsol_rad (W m^{-2})
	536	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Surface specific humidity (kg kg^{-1})
	537	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxsens !! Sensible heat flux (W m^{-2})
	538	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxlat !! Latent heat flux (W m^{-2})
	539	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapp !! Total of evaporation (mm day^{-1})
[2031]	540
[947]	541	!_ ================================================================================================================================
[2581]	542
	543	!! 1. Write variables to restart file to be used for the next simulation
	544	IF (printlev>=3) WRITE (numout,*) 'Write restart file with ENERBIL variables'
[947]	545
[2581]	546	CALL restput_p(rest_id, 'temp_sol', nbp_glo, 1, 1, kjit, temp_sol, 'scatter', nbp_glo, index_g)
	547	CALL restput_p(rest_id, 'qsurf', nbp_glo, 1, 1, kjit, qsurf, 'scatter', nbp_glo, index_g)
	548	CALL restput_p(rest_id, 'evapot', nbp_glo, 1, 1, kjit, evapot, 'scatter', nbp_glo, index_g)
	549	CALL restput_p(rest_id, 'evapot_corr', nbp_glo, 1, 1, kjit, evapot_corr, 'scatter', nbp_glo, index_g)
	550	CALL restput_p(rest_id, 'tsolrad', nbp_glo, 1, 1, kjit, tsol_rad, 'scatter', nbp_glo, index_g)
	551	CALL restput_p(rest_id, 'evapora', nbp_glo, 1, 1, kjit, vevapp, 'scatter', nbp_glo, index_g)
	552	CALL restput_p(rest_id, 'fluxlat', nbp_glo, 1, 1, kjit, fluxlat, 'scatter', nbp_glo, index_g)
	553	CALL restput_p(rest_id, 'fluxsens', nbp_glo, 1, 1, kjit, fluxsens, 'scatter', nbp_glo, index_g)
	554	CALL restput_p(rest_id, 'qsolpot', nbp_glo, 1, 1, kjit, q_sol_pot, 'scatter', nbp_glo, index_g)
	555
	556	END SUBROUTINE enerbil_finalize
[8]	557
	558
	559
	560
[947]	561	!! =============================================================================================================================
	562	!! SUBROUTINE : enerbil_clear
	563	!!
	564	!>\BRIEF Routine deallocates clear output variables if already allocated.
	565	!!
	566	!! DESCRIPTION : This is a 'housekeeping' routine that deallocates the key output
	567	!! variables, if they have already been allocated. The variables that are deallocated are psold,
	568	!! qsol_sat, pdqsold, psnew, qsol_sat_new, netrad, lwabs, lwup, lwnet, fluxsubli, qsat_air, tair
	569	!!
	570	!! RECENT CHANGE(S) : None
	571	!!
	572	!! MAIN OUTPUT VARIABLE(S) : None
	573	!!
	574	!! REFERENCES : None
	575	!!
	576	!! FLOWCHART : None
	577	!! \n
	578	!_ ==============================================================================================================================
	579
[8]	580	SUBROUTINE enerbil_clear ()
	581	IF ( ALLOCATED (psold)) DEALLOCATE (psold)
	582	IF ( ALLOCATED (qsol_sat)) DEALLOCATE (qsol_sat)
	583	IF ( ALLOCATED (pdqsold)) DEALLOCATE (pdqsold)
	584	IF ( ALLOCATED (psnew)) DEALLOCATE (psnew)
	585	IF ( ALLOCATED (qsol_sat_new)) DEALLOCATE (qsol_sat_new)
	586	IF ( ALLOCATED (netrad)) DEALLOCATE (netrad)
	587	IF ( ALLOCATED (lwabs)) DEALLOCATE (lwabs)
	588	IF ( ALLOCATED (lwup)) DEALLOCATE (lwup)
	589	IF ( ALLOCATED (lwnet)) DEALLOCATE (lwnet)
	590	IF ( ALLOCATED (fluxsubli)) DEALLOCATE (fluxsubli)
	591	IF ( ALLOCATED (qsat_air)) DEALLOCATE (qsat_air)
	592	IF ( ALLOCATED (tair)) DEALLOCATE (tair)
[947]	593	IF ( ALLOCATED (q_sol_pot)) DEALLOCATE (q_sol_pot)
	594
	595	END SUBROUTINE enerbil_clear
[8]	596
[947]	597
	598	!! =============================================================================================================================
	599	!! SUBROUTINE : enerbil_begin
[8]	600	!!
[947]	601	!>\BRIEF Preliminary variables required for the calculation of
	602	!! the energy budget are derived.
	603	!!
	604	!! DESCRIPTION : This routines computes preliminary variables required
	605	!! for the calculation of the energy budget: the old surface static energy (psold), the surface saturation
	606	!! humidity (qsol_sat), the derivative of satured specific humidity at the old temperature (pdqsold)
	607	!! and the net radiation (netrad).
	608	!!
	609	!! MAIN OUTPUT VARIABLE(S) : psold, qsol_sat, pdqsold, netrad
	610	!!
	611	!! REFERENCE(S) :
	612	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
	613	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
	614	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
	615	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
	616	!! Circulation Model. Journal of Climate, 6, pp.248-273
	617	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
	618	!! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11):
	619	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
	620	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
	621	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
	622	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
	623	!! general circulation models. Global and Planetary Change, 19, pp.261-276
	624	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
	625	!! Interscience Publishers\n
	626	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
	627	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
	628	!!
	629	!! FLOWCHART : None
	630	!! \n
	631	!_ ==============================================================================================================================
	632
[2650]	633	SUBROUTINE enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis)
[8]	634
[947]	635	!! 0. Variable and parameter declaration
[8]	636
[947]	637	!! 0.1 Input variables
[8]	638
[947]	639	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
	640	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature (K)
	641	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2})
	642	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2})
[2031]	643	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
[947]	644	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-)
[2222]	645
[947]	646	!! 0.2 Output variables
	647
	648	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: psold !! Old surface dry static energy (J kg^{-1})
	649	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat !! Saturated specific humudity for old temperature
	650	!! (kg kg^{-1})
	651	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: pdqsold !! Derivative of satured specific humidity at the old
	652	!! temperature (kg (kg s)^{-1})
	653	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: netrad !! Net radiation (W m^{-2})
	654
	655	!! 0.3 Modified variables
	656
	657	!! 0.4 Local variables
	658
	659	INTEGER(i_std) :: ji
	660	REAL(r_std), DIMENSION(kjpindex) :: dev_qsol
	661	REAL(r_std), PARAMETER :: missing = 999998.
[2222]	662
[947]	663	!_ ================================================================================================================================
	664
	665	!! 1. Computes psold (the surface static energy for the old timestep)
	666
	667	!! We here define the surface static energy for the 'old' timestep, in terms of the surface
	668	!! temperature and heat capacity.
	669	!! \latexonly
	670	!! \input{enerbilbegin1.tex}
	671	!! \endlatexonly
[8]	672	psold(:) = temp_sol(:)*cp_air
[947]	673	!! 2. Computes qsol_sat (the surface saturated humidity).
	674
	675	!! We call the routine 'qsatcalc' from within the module 'src_parameters/constantes_veg'.
[8]	676	CALL qsatcalc (kjpindex, temp_sol, pb, qsol_sat)
	677	IF ( diag_qsat ) THEN
	678	IF ( ANY(ABS(qsol_sat(:)) .GT. missing) ) THEN
	679	DO ji = 1, kjpindex
	680	IF ( ABS(qsol_sat(ji)) .GT. missing) THEN
	681	WRITE(numout,*) 'ERROR on ji = ', ji
	682	WRITE(numout,*) 'temp_sol(ji), pb(ji) :', temp_sol(ji), pb(ji)
[1078]	683	CALL ipslerr_p (3,'enerbil_begin', &
[8]	684	& 'qsol too high ','','')
	685	ENDIF
	686	ENDDO
	687	ENDIF
	688	ENDIF
[947]	689
	690	!! 3. Computes pdqsold
	691
	692	!! Computes pdqsold (the derivative of the saturated humidity with respect to temperature
	693	!! analysed at the surface temperature at the 'old' timestep.
[2093]	694	!! We call the routine 'dev_qsatcalc' from qsat_moisture module.
[8]	695	CALL dev_qsatcalc (kjpindex, temp_sol, pb, dev_qsol)
[947]	696
	697	!! \latexonly
	698	!! \input{enerbilbegin2.tex}
	699	!! \endlatexonly
[2093]	700	pdqsold(:) = dev_qsol(:)
	701
[8]	702
	703	IF ( diag_qsat ) THEN
	704	IF ( ANY(ABS( pdqsold(:)) .GT. missing) ) THEN
	705	DO ji = 1, kjpindex
	706	IF ( ABS( pdqsold(ji)) .GT. missing ) THEN
	707	WRITE(numout,*) 'ERROR on ji = ', ji
	708	WRITE(numout,*) 'temp_sol(ji), pb(ji) :', temp_sol(ji), pb(ji)
[1078]	709	CALL ipslerr_p (3,'enerbil_begin', &
[8]	710	& 'pdqsold too high ','','')
	711	ENDIF
	712	ENDDO
	713	ENDIF
	714	ENDIF
	715
[947]	716	!! 4. Computes the net radiation and the absorbed LW radiation absorbed at the surface.
	717
	718	!! Long wave radiation absorbed by the surface is the product of emissivity and downwelling LW radiation
	719	!! \latexonly
	720	!! \input{enerbilbegin3.tex}
	721	!! \endlatexonly
[8]	722	lwabs(:) = emis(:) * lwdown(:)
[947]	723
	724	!! Net radiation is calculated as:
	725	!! \latexonly
	726	!! \input{enerbilbegin4.tex}
	727	!! \endlatexonly
[8]	728	netrad(:) = lwdown(:) + swnet (:) - (emis(:) * c_stefan * temp_sol(:)*4 + (un - emis(:)) lwdown(:))
[2222]	729
[2348]	730	IF (printlev>=3) WRITE (numout,*) ' enerbil_begin done '
[8]	731
	732	END SUBROUTINE enerbil_begin
	733
[947]	734
	735	!! =============================================================================================================================
	736	!! SUBROUTINE : enerbil_surftemp
[8]	737	!!
[947]	738	!>\BRIEF This routine computes the new surface static energy
	739	!! (psnew) and the saturated humidity at the surface (qsol_sat_new).
[8]	740	!!
[947]	741	!! DESCRIPTION : This is the key part of the enerbil module, for which
[2591]	742	!! the energy budget in the surface layer is solved and changes over the model timestep 'dt_sechiba' are
[947]	743	!! quantified for the surface static energy, surface temperature, surface humidity, the 'air' (or lowest
	744	!! level atmospheric model) temperature, the 'air' (or lowest level atmospheric model) humidity,
	745	!! according to the method that is laid out by Dufresne \& Ghattas (2009) and Shultz et al. (2001).
	746	!!
	747	!! It computes the energy balance at the surface with an implicit scheme, that is connected to the
	748	!! Richtmyer and Morton algorithm of the PBL. By computing the surface temperature and surface humidity
	749	!! the routine also implicitly estimates the various fluxes which balance in order to give the new
	750	!! temperature. Thus once the surface temperature has been obtained all the fluxes need to be diagnosed.
	751	!!
[5470]	752	!! Since the explicit snow scheme is adapted, the notion of skin layer in ORCHIDEE is abandoned and
[2222]	753	!! the first thin snow layer is adopted to solve the surface energy budget.
	754	!! Implicit method is still used for coupling to the atmosphere.
	755	!! In this new scheme, the snow temperature profile is simulatenously solved
	756	!! along with the surface snow temperature, and the details are referenced to Boone et al. (2010)
	757	!!
[947]	758	!! MAIN OUTPUT VARIABLE(S) : psnew, qsol_sat_new, temp_sol_new, qair_new, epot_air_new
	759	!!
	760	!! REFERENCE(S) :
	761	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
	762	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
	763	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
	764	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
	765	!! Circulation Model. Journal of Climate, 6, pp.248-273
	766	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
	767	!! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11):
	768	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
	769	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
	770	!! sur le cycle de l'eau, PhD Thesis, available (22/12/11):
	771	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pd
	772	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
	773	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
	774	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
	775	!! general circulation models. Global and Planetary Change, 19, pp.261-276
	776	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
	777	!! Interscience Publishers
	778	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
	779	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
	780	!!
	781	!! FLOWCHART : None
	782	!!
	783	!_ ==============================================================================================================================
	784
	785
[2591]	786	SUBROUTINE enerbil_surftemp (kjpindex, zlev, emis, epot_air, &
[8]	787	& petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,&
[4197]	788	& vbeta1, vbeta5, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, &
[3959]	789	& qair_new, epot_air_new)
[8]	790
	791
[2650]	792
[947]	793	!! 0. Variable and parameter declaration
[8]	794
[947]	795	!! 0.1 Input variables
	796
	797	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
	798	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
	799	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-)
	800	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (?? J)
	801	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef (see note)
	802	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef (see note)
	803	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1})
	804	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef (see note)
	805	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef (see note)
	806	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Soil flux (W m^{-2})
	807	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3})
	808	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind velocity by directional components
	809	!! u (Eastwards) and v (Northwards) (m s^{-1})
[4146]	810	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
[947]	811	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-)
	812	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-)
	813	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance
	814	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity (J K^{-1])
	815	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2})
	816	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2})
[2222]	817
[947]	818	!! 0.2 Output variables
	819
	820	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: psnew !! New surface static energy (J kg^{-1})
	821	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat_new !! New saturated surface air moisture (kg kg^{-1})
	822	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K)
	823	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qair_new !! New air moisture (kg kg^{-1})
	824	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: epot_air_new !! New air temperature (K)
	825
	826
	827	!! 0.4 Local variables
	828
	829	INTEGER(i_std) :: ji
[8]	830	REAL(r_std),DIMENSION (kjpindex) :: zicp
	831	REAL(r_std) :: fevap
	832	REAL(r_std) :: sensfl_old, larsub_old, lareva_old, dtheta, sum_old, sum_sns
	833	REAL(r_std) :: zikt, zikq, netrad_sns, sensfl_sns, larsub_sns, lareva_sns
	834	REAL(r_std) :: speed
[4870]	835	REAL(r_std),DIMENSION (kjpindex) :: dtes !! Diagnostic output variable for the change of the thermal energy content of surface for which the energy balance is calculated (J m-2)
[2222]	836
[947]	837	!_ ================================================================================================================================
[2222]	838
[8]	839	zicp = un / cp_air
[947]	840
[8]	841	DO ji=1,kjpindex
[947]	842	!! 1. Derivation of auxiliary variables
	843
	844	!! \latexonly
	845	!! \input{enerbilsurftemp1.tex}
	846	!! \endlatexonly
[8]	847	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
	848	!
[947]	849	!! \latexonly
	850	!! \input{enerbilsurftemp2.tex}
	851	!! \endlatexonly
[8]	852	zikt = 1/(rau(ji) * speed * q_cdrag(ji))
	853	zikq = 1/(rau(ji) * speed * q_cdrag(ji))
[947]	854
	855	!! 2. Computation of fluxes for the old surface conditions
	856
	857	!! As a reference, we first calculate the sensible and latent heat for the 'old' timestep.
	858
	859	!! 2.1 Sensible heat for the 'old' timestep
	860	!! This is equation (64) of (Dufresne & Ghattas, 2009), rewritten in terms of H_old. We make the
	861	!! approximation that (P_r/P)^{kappa}=1 and convert the A and B coefficients to the ORCHIDEE
	862	!! format (see introductory note, above). Also, the equation here is in terms of surface static
	863	!! energy, which is defined as per section 3.1, above.
	864	!! \latexonly
	865	!! \input{enerbilsurftemp3.tex}
	866	!! \endlatexonly
[8]	867	sensfl_old = (petBcoef(ji) - psold(ji)) / (zikt - petAcoef(ji))
[947]	868	!! 2.2 Latent heat for the 'old' timestep
	869	!! This is equation (70) of (Dufresne & Ghattas, 2009), rewritten in terms of {\lambda E}_{old}.
	870	!! Again we convert the A and B coefficients from the LMDZ format to the ORCHIDEE format.\n
	871	!! There are two forms of the equation - first for the latent heat as a result of sublimination
	872	!! processes and then as a result of evaporative processes:
	873	!! \latexonly
	874	!! \input{enerbilsurftemp4.tex}
	875	!! \endlatexonly
[2425]	876	larsub_old = chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) * &
	877	(peqBcoef(ji) - qsol_sat(ji)) / (zikq - vbeta1(ji) * (un - vbeta5(ji))* peqAcoef(ji))
[947]	878	!! \latexonly
	879	!! \input{enerbilsurftemp5.tex}
	880	!! \endlatexonly
	881	lareva_old = chalev0 * (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * &
[4197]	882	(peqBcoef(ji) - qsol_sat(ji)) / &
[2425]	883	(zikq - (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * peqAcoef(ji)) &
	884	+ chalev0 * vbeta5(ji) * (peqBcoef(ji) - qsol_sat(ji)) / (zikq - vbeta5(ji) * peqAcoef(ji))
[947]	885
	886	!! 3. Calculation of sensitivity terms
	887
	888	!! We here calculate the sensitivity for the different fluxes to changes in dtheta, which is the
[2591]	889	!! change in the surface static energy over the model time step (dt_sechiba).
[947]	890
	891	!! 3.1 Net radiation sensitivity
	892	!! This is an explicit-implicit representation of the Stefan-Boltzmann law - the explicit terms
	893	!! are ${ps}_{old}^3$, and it is completed when multiplied by dtheta (defined above).
	894	!! \latexonly
	895	!! \input{enerbilsurftemp6.tex}
	896	!! \endlatexonly
[8]	897	netrad_sns = zicp(ji) * quatre * emis(ji) * c_stefan * ((zicp(ji) * psold(ji))**3)
[947]	898
	899	!! 3.2 Sensible heat flux sensitivity
	900	!! This is the temperature sensitivity term N_1^h derived in equation (66) of (Dufresne & Ghattas, 2009),
	901	!! where we again assume that (P_r/P)^{kappa}=1, convert the A and B coefficients to the ORCHIDEE format,
	902	!! and rewrite in terms of surface static energy, rather than temperature (see section 3.1, above).
	903	!! \latexonly
	904	!! \input{enerbilsurftemp7.tex}
	905	!! \endlatexonly
[8]	906	sensfl_sns = un / (zikt - petAcoef(ji))
[947]	907	!! 3.3 Latent heat flux sensitivity
	908	!! This is the humidity sensitivity term N_1^q derived in equation (72) of (Dufresne & Ghattas, 2009), where
	909	!! the A and B coefficients are written in the ORCHIDEE format.
	910	!! larsub_sns is the latent heat sensitivity for sublimination. The coefficient vbeta1 is the evaporation
	911	!! coefficient. It represents the relationship between the real and potential evaporation. It is derived
	912	!! in the module 'src_sechiba/diffuco_snow'.
	913	!! \latexonly
	914	!! \input{enerbilsurftemp8.tex}
	915	!! \endlatexonly
[2425]	916	larsub_sns = chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) * zicp(ji) * &
	917	pdqsold(ji) / (zikq - vbeta1(ji) * (un - vbeta5(ji)) * peqAcoef(ji))
[947]	918
	919	!! lareva_sns is the latent heat sensitivity for evaporation. vbeta1 (src_sechiba/diffuco_snow),
[4197]	920	!! and vbeta (src_sechiba/diffuco_comb) are the evaporation
[947]	921	!! coefficients.
	922	!! \latexonly
	923	!! \input{enerbilsurftemp9.tex}
	924	!! \endlatexonly
[4197]	925	lareva_sns = chalev0 * ((un - vbeta1(ji))(un - vbeta5(ji)) vbeta(ji) + vbeta5(ji)) * &
[2425]	926	& zicp(ji) * pdqsold(ji) / &
[4197]	927	(zikq - ((un - vbeta1(ji))(un - vbeta5(ji)) vbeta(ji) + vbeta5(ji))* peqAcoef(ji))
[2222]	928
[947]	929	!! 4. Solution of the energy balance
	930
	931	!! 4.1 We calculate the total flux for the 'old' timestep.
	932	!! \latexonly
	933	!! \input{enerbilsurftemp10.tex}
	934	!! \endlatexonly
[3020]	935	sum_old = netrad(ji) + sensfl_old + larsub_old + lareva_old + soilflx(ji)
[3059]	936
[947]	937	!! 4.2 We calculate the total sensitivity dtheta (the change in the
	938	!! surface static energy over the timestep.
	939	!! \latexonly
	940	!! \input{enerbilsurftemp11.tex}
	941	!! \endlatexonly
[8]	942	sum_sns = netrad_sns + sensfl_sns + larsub_sns + lareva_sns
[947]	943	!! 4.3 Calculation of dtheta (change in surface static energy over the
	944	!! timestep.
	945	!! \latexonly
	946	!! \input{enerbilsurftemp12.tex}
	947	!! \endlatexonly
[3059]	948
[3020]	949	dtheta = dt_sechiba * sum_old / (zicp(ji) * soilcap(ji) + dt_sechiba * sum_sns)
[3059]	950
[947]	951	!! 4.4 Determination of state variables for the 'new' timestep
	952	!! No we have dtheta, we can determine the surface static energy that
	953	!! corresponds to the 'new' timestep.
	954	!! \latexonly
	955	!! \input{enerbilsurftemp13.tex}
	956	!! \endlatexonly
[8]	957	psnew(ji) = psold(ji) + dtheta
[947]	958
	959	!! The new surface saturated humidity can be calculated by equation (69)
	960	!! of (Dufresne & Ghattas, 2009), in which we substitute dtheta for the
	961	!! change between old and new temperature using the relationship from 3.1,
	962	!! above.
	963	!! \latexonly
	964	!! \input{enerbilsurftemp14.tex}
	965	!! \endlatexonly
[8]	966	qsol_sat_new(ji) = qsol_sat(ji) + zicp(ji) * pdqsold(ji) * dtheta
[947]	967	!! The new surface temperature is determined from the new surface static temperature,
	968	!! again by using the relationship in section 3.1.
	969	!! \latexonly
	970	!! \input{enerbilsurftemp15.tex}
	971	!! \endlatexonly
[8]	972	temp_sol_new(ji) = psnew(ji) / cp_air
[947]	973
	974	!! 4.5 Calculation of new evaporation potential and new evaporation latent heat
	975	!! flux (???)
	976	!! \latexonly
	977	!! \input{enerbilsurftemp16.tex}
	978	!! \endlatexonly
[8]	979	epot_air_new(ji) = zikt * (sensfl_old - sensfl_sns * dtheta) + psnew(ji)
[947]	980
	981	!! \latexonly
	982	!! \input{enerbilsurftemp17.tex}
	983	!! \endlatexonly
[8]	984	fevap = (lareva_old - lareva_sns * dtheta) + (larsub_old - larsub_sns * dtheta)
[947]	985
[8]	986	IF ( ABS(fevap) < EPSILON(un) ) THEN
[947]	987
	988	!! \latexonly
	989	!! \input{enerbilsurftemp18.tex}
	990	!! \endlatexonly
[8]	991	qair_new(ji) = qair(ji)
	992	ELSE
[947]	993	!! \latexonly
	994	!! \input{enerbilsurftemp19.tex}
	995	!! \endlatexonly
	996	qair_new(ji) = zikq * un / ( chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) + &
[4197]	997	& chalev0 * ((un - vbeta1(ji))(un - vbeta5(ji)) vbeta(ji) + vbeta5(ji)) ) &
[947]	998	& * fevap + qsol_sat_new(ji)
[8]	999	ENDIF
[2650]	1000
[4870]	1001	! Calculate dtes only for diagnostic output and send it to xios
	1002	dtes(ji) = (soilcap(ji)/cp_air)*dtheta/dt_sechiba
	1003
[8]	1004	ENDDO
	1005
[4870]	1006	! Send diagnostic variable to XIOS
	1007	CALL xios_orchidee_send_field("dtes", dtes)
	1008
[2348]	1009	IF (printlev>=3) WRITE (numout,*) ' enerbil_surftemp done '
[8]	1010
	1011	END SUBROUTINE enerbil_surftemp
	1012
	1013
[947]	1014	!! =============================================================================================================================
	1015	!! SUBROUTINE : enerbil_pottemp
	1016	!!
	1017	!>\BRIEF This subroutine computes the surface temperature and humidity should the surface been unstressed
	1018	!!
	1019	!! DESCRIPTION :
	1020	!!
	1021	!! MAIN OUTPUT VARIABLE(S) :
	1022	!!
	1023	!! REFERENCE(S) :
	1024	!!
	1025	!! FLOWCHART : None
	1026	!!
	1027	!_ ==============================================================================================================================
	1028
[2591]	1029	SUBROUTINE enerbil_pottemp (kjpindex, zlev, emis, epot_air, &
[947]	1030	& petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,&
[4197]	1031	& vbeta1, vbeta5, soilcap, lwdown, swnet, q_sol_pot, temp_sol_pot)
[947]	1032
	1033	! interface
[8]	1034	! input scalar
	1035	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
[947]	1036	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
[8]	1037	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity
[947]	1038	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy
	1039	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef
	1040	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef
	1041	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
	1042	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef
	1043	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef
	1044	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Soil flux
[8]	1045	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
[947]	1046	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind
[4146]	1047	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
[8]	1048	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient
[947]	1049	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance
	1050	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance
	1051	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity
	1052	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux
	1053	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
[8]	1054	! output fields
[947]	1055	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: q_sol_pot !! Potential surface air moisture
	1056	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_pot !! Potential soil temperature
[8]	1057
	1058
	1059	! local declaration
[947]	1060	INTEGER(i_std) :: ji
	1061	REAL(r_std),DIMENSION (kjpindex) :: zicp
	1062	REAL(r_std) :: fevap
	1063	REAL(r_std) :: sensfl_old, larsub_old, lareva_old, dtheta, sum_old, sum_sns
	1064	REAL(r_std) :: zikt, zikq, netrad_sns, sensfl_sns, larsub_sns, lareva_sns
	1065	REAL(r_std) :: speed
[8]	1066
[947]	1067	!_ ==============================================================================================================================
	1068
	1069	zicp = un / cp_air
[8]	1070	!
[947]	1071	DO ji=1,kjpindex
[8]	1072
[947]	1073	dtheta = zero
	1074	fevap = zero
	1075
	1076	temp_sol_pot(ji) = temp_sol_pot(ji) + dtheta
	1077
	1078	q_sol_pot(ji) = q_sol_pot(ji) + fevap
	1079
	1080	ENDDO
	1081
	1082	END SUBROUTINE enerbil_pottemp
	1083
	1084
	1085	!! =============================================================================================================================
	1086	!! SUBROUTINE : enerbil_flux
	1087	!!
	1088	!>\BRIEF Computes the new soil temperature, net radiation and
	1089	!! latent and sensible heat flux for the new time step.
	1090	!!
	1091	!! DESCRIPTION : This routine diagnoses, based on the new soil temperature,
	1092	!! the net radiation, the total evaporation, the latent heat flux, the sensible heat flux and the
	1093	!! sublimination flux. It also diagnoses the potential evaporation used for the fluxes (evapot) and the potential
	1094	!! as defined by Penman & Monteith (Monteith, 1965) based on the correction term developed by Chris
	1095	!! Milly (1992). This Penman-Monteith formulation is required for the estimation of bare soil evaporation
	1096	!! when the 11 layer CWRR moisture scheme is used for the soil hydrology.
	1097	!!
	1098	!! MAIN OUTPUT VARIABLE(S) : qsurf, fluxsens, fluxlat, fluxsubli, vevapp, lwup, lwnet,
	1099	!! tsol_rad, netrad, evapot, evapot_corr
	1100	!!
	1101	!! REFERENCE(S) :
	1102	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
	1103	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
	1104	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
	1105	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
	1106	!! Circulation Model. Journal of Climate, 6, pp.248-273
	1107	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
	1108	!! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11):
	1109	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
	1110	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
	1111	!! sur le cycle de l'eau, PhD Thesis, available (22/12/11):
	1112	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
	1113	!! - Monteith, JL, 1965. Evaporation and Environment, paper presented at Symposium of the Society
	1114	!! for Experimental Biology
	1115	!! - Monteith & Unsworth, 2008. Principles of Environmental Physics (third edition), published Elsevier
	1116	!! ISBN 978-0-12-505103-3
	1117	!! - Milly, P. C. D., 1992: Potential Evaporation and Soil Moisture in General Circulation Models.
	1118	!! Journal of Climate, 5, pp. 209â226.
	1119	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
	1120	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
	1121	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
	1122	!! general circulation models. Global and Planetary Change, 19, pp.261-276
	1123	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
	1124	!! Interscience Publishers
	1125	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
	1126	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
	1127	!!
	1128	!! FLOWCHART : None
	1129	!! \n
	1130	!_ ==============================================================================================================================
	1131
[4197]	1132	SUBROUTINE enerbil_flux (kjpindex, emis, temp_sol, rau, u, v, q_cdrag, vbeta, vbeta1, vbeta5, &
[947]	1133	& qair, epot_air, psnew, qsurf, fluxsens, fluxlat, fluxsubli, vevapp, temp_sol_new, &
[2222]	1134	& lwdown, swnet, lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr,&
[3020]	1135	& precip_rain,snowdz,temp_air,pgflux, soilcap, temp_sol_add)
[947]	1136
	1137	!! 0. Variable and parameter declaration
	1138
	1139	!! 0.1 Input variables
	1140
	1141	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
	1142	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-)
	1143	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Surface temperature (K)
	1144	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density (kg m^{-3})
	1145	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u,v !! Wind velocity in components u and v (m s^{-1})
[4146]	1146	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
[947]	1147	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-)
	1148	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-)
	1149	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Flood resistance
	1150	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1})
	1151	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (J)
	1152	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: psnew !! New surface static energy (J kg^{-1})
	1153	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature (K)
[2031]	1154	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
[2222]	1155	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth
	1156	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: precip_rain !! Rainfall
	1157	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: temp_air !! Air temperature
[2650]	1158	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Downward long wave radiation (W m^{-2})
	1159	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net short wave radiation (W m^{-2})
[3020]	1160	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil calorific capacity including snow and soil (J K^{-1})
[2222]	1161
[947]	1162	!! 0.2 Output variables
	1163
	1164	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! Surface specific humidity (kg kg^{-1})
	1165	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2})
	1166	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2})
	1167	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsubli !! Energy of sublimation (W m^{-2})
	1168	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1})
	1169	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwup !! Long-wave up radiation (W m^{-2})
	1170	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwnet !! Long-wave net radiation (W m^{-2})
	1171	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature (W m^{-2})
[2650]	1172	REAL(r_std), DIMENSION (kjpindex),INTENT(out) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
[2222]	1173
	1174
	1175	!! 0.3 Modified variables
[947]	1176	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: netrad !! Net radiation (W m^{-2})
[3406]	1177	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation (mm/tstep)
	1178	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep)
[2222]	1179	REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: pgflux !! Net energy into snowpack(W/m^2)
[947]	1180
	1181	!! 0.4 Local variables
	1182
	1183	INTEGER(i_std) :: ji
	1184	REAL(r_std),DIMENSION (kjpindex) :: grad_qsat
	1185	REAL(r_std) :: correction
	1186	REAL(r_std) :: speed !! Speed (m s^{-1}),
[2650]	1187	REAL(r_std) :: qc !! Surface drag coefficient (??)
[947]	1188	LOGICAL,DIMENSION (kjpindex) :: warning_correction
[2650]	1189	REAL(r_std), DIMENSION (kjpindex) :: zgflux
	1190	REAL(r_std), DIMENSION (kjpindex) :: qsol_sat_tmp
	1191	REAL(r_std), DIMENSION (kjpindex) :: zerocelcius
	1192	REAL(r_std), DIMENSION (kjpindex) :: PHPSNOW
	1193	REAL(r_std), DIMENSION (kjpindex) :: lwup_tmp
	1194	REAL(r_std), DIMENSION (kjpindex) :: netrad_tmp
	1195	REAL(r_std), DIMENSION (kjpindex) :: fluxsens_tmp
	1196	REAL(r_std), DIMENSION (kjpindex) :: fluxlat_tmp
[7084]	1197	REAL(r_std), DIMENSION (kjpindex) :: betatot !! beta used in the computation of qsurf (-)
[947]	1198	!_ ================================================================================================================================
	1199
[2650]	1200	zerocelcius(:) = tp_00
	1201	CALL qsatcalc (kjpindex, zerocelcius, pb, qsol_sat_tmp)
[2222]	1202
[8]	1203	DO ji=1,kjpindex
[2650]	1204	!! 1. Determination of 'housekeeping' variables
[947]	1205
	1206	!! The horizontal wind speed is calculated from the velocities along each axis using Pythagorus' theorem.
	1207	!! \latexonly
	1208	!! \input{enerbilflux1.tex}
	1209	!! \endlatexonly
[8]	1210	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
[947]	1211
	1212	!! From definition, the surface drag coefficient is defined:
	1213	!! \latexonly
	1214	!! \input{enerbilflux2.tex}
	1215	!! \endlatexonly
[8]	1216	qc = speed * q_cdrag(ji)
[947]	1217
	1218	!! 2. Calculation of the net upward flux of longwave radiation
	1219
	1220	!! We first of all calculate the radiation as a result of the Stefan-Boltzmann equation,
	1221	!! which is the sum of the calculated values at the surface temperature at the 'old'
	1222	!! temperature and the value that corresponds to the difference between the old temperature
	1223	!! and the temperature at the 'new' timestep.
	1224	!! \latexonly
	1225	!! \input{enerbilflux3.tex}
	1226	!! \endlatexonly
[8]	1227	lwup(ji) = emis(ji) * c_stefan * temp_sol(ji)**4 + &
	1228	& quatre * emis(ji) * c_stefan * temp_sol(ji)*3 &
[947]	1229	& (temp_sol_new(ji) - temp_sol(ji))
	1230
	1231	!! We then add to this value that from the reflected longwave radiation:
	1232	!! \latexonly
	1233	!! \input{enerbilflux4.tex}
	1234	!! \endlatexonly
[8]	1235	lwup(ji) = lwup(ji) + (un - emis(ji)) * lwdown(ji)
[947]	1236
[2183]	1237	!! The radiative surface temperature is calculated by inverting the Stefan-Boltzmann relation:
[947]	1238	!! \latexonly
	1239	!! \input{enerbilflux5.tex}
[2183]	1240	!! \endlatexonly
	1241	!! The implicit solution computes an emitted long-wave flux which is a limited Taylor expansion
	1242	!! of the future surface temperature around the current values. Thus the long-wave flux does not
	1243	!! correspond to the new surface temperature but some intermediate value. So we need to deduce the
	1244	!! radiative surface temperature to which this flux corresponds.
	1245	!!
	1246	tsol_rad(ji) = (lwup(ji)/ (emis(ji) * c_stefan)) **(1./quatre)
[947]	1247
	1248	!! qsurf (the surface specific humidity) is a simple diagnostic which will be used by the
	1249	!! GCM to compute the dependence of of the surface layer stability on moisture.
	1250	!! \latexonly
	1251	!! \input{enerbilflux6.tex}
	1252	!! \endlatexonly
[7084]	1253
	1254	betatot(ji) = vbeta1(ji) * (un - vbeta5(ji)) + vbeta5(ji) + &
	1255	& (un - vbeta1(ji))(un - vbeta5(ji)) vbeta(ji)
	1256
	1257	qsurf(ji) = betatot(ji) * qsol_sat_new(ji) + (un - betatot(ji)) * qair(ji)
	1258
[6019]	1259	! When no turbulence and no wind (qc small) the near surface air moisture is equal to qair.
	1260	! This is related to problem with interpolation in stdlevvar in LMDZ.
	1261	IF ( qc .LE. min_qc ) qsurf(ji) = qair(ji)
	1262
	1263
[947]	1264	!! \latexonly
	1265	!! \input{enerbilflux7.tex}
	1266	!! \endlatexonly
[8]	1267	qsurf(ji) = MAX(qsurf(ji), qair(ji))
	1268
[947]	1269	!! Net downward radiation is the sum of the down-welling less the up-welling long wave flux, plus
	1270	!! the short wave radiation.
	1271	!! \latexonly
	1272	!! \emissivity absorbed lw radiationinput{enerbilflux8.tex}
	1273	!! \endlatexonly
[8]	1274	netrad(ji) = lwdown(ji) + swnet(ji) - lwup(ji)
[947]	1275
	1276	!! 'vevapp' is the sum of the total evaporative processes (snow plus non-snow processes).
	1277	!! \latexonly
	1278	!! \input{enerbilflux9.tex}
	1279	!! \endlatexonly
[2591]	1280	vevapp(ji) = dt_sechiba * rau(ji) * qc * (vbeta1(ji) * (un - vbeta5(ji)) + vbeta5(ji)) * &
[947]	1281	& (qsol_sat_new(ji) - qair(ji)) + &
[2591]	1282	& dt_sechiba * rau(ji) * qc * (un - vbeta1(ji))(un-vbeta5(ji)) vbeta(ji) * &
[4197]	1283	& (qsol_sat_new(ji) - qair(ji))
[947]	1284
	1285	!! The total latent heat flux is the sum of the snow plus non-snow processes.
	1286	!! \latexonly
	1287	!! \input{enerbilflux10.tex}
	1288	!! \endlatexonly
	1289	fluxlat(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) * (un - vbeta5(ji)) * &
[8]	1290	& (qsol_sat_new(ji) - qair(ji)) + &
[947]	1291	& chalev0 * rau(ji) * qc * vbeta5(ji) *&
	1292	& (qsol_sat_new(ji) - qair(ji)) + &
	1293	& chalev0 * rau(ji) * qc * (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * &
[4197]	1294	& (qsol_sat_new(ji) - qair(ji))
[947]	1295
	1296	!! The sublimination flux concerns is calculated using vbeta1, the snow resistance.
	1297	!! \latexonly
	1298	!! \input{enerbilflux11.tex}
	1299	!! \endlatexonly
	1300	fluxsubli(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) * (un - vbeta5(ji)) * &
[8]	1301	& (qsol_sat_new(ji) - qair(ji))
[947]	1302
	1303	!! The sensible heat flux is a factor of the difference between the new surface static energy
	1304	!! and the potential energy of air.
	1305	!! \latexonly
	1306	!! \input{enerbilflux12.tex}
	1307	!! \endlatexonly
[8]	1308	fluxsens(ji) = rau(ji) * qc * (psnew(ji) - epot_air(ji))
[947]	1309
	1310	!! This is the net longwave downwards radiation.
	1311	!! \latexonly
	1312	!! \input{enerbilflux13.tex}
	1313	!! \endlatexonly
[8]	1314	lwnet(ji) = lwdown(ji) - lwup(ji)
[947]	1315
	1316	!! Diagnoses the potential evaporation used for the fluxes (evapot)
	1317	!! \latexonly
	1318	!! \input{enerbilflux14.tex}
	1319	!! \endlatexonly
[2591]	1320	evapot(ji) = MAX(zero, dt_sechiba * rau(ji) * qc * (qsol_sat_new(ji) - qair(ji)))
[947]	1321
	1322	!! From definition we can say:
	1323	!! \latexonly
	1324	!! \input{enerbilflux15.tex}
	1325	!! \endlatexonly
[8]	1326	tair(ji) = epot_air(ji) / cp_air
[2650]	1327
	1328
	1329	!! To calculate net energy flux into the snowpack
[5470]	1330	PHPSNOW(ji) = precip_rain(ji)(4.218E+3)(MAX(tp_00,temp_air(ji))-tp_00)/dt_sechiba ! (w/m2)
	1331	pgflux(ji) = netrad(ji) - fluxsens(ji) - fluxlat(ji) + PHPSNOW(ji)
	1332
[2650]	1333
	1334	!! To get the extra energy used to melt the snowpack
[5431]	1335	IF (temp_sol_new (ji) > tp_00 .AND. &
[4624]	1336	SUM(snowdz(ji,:)) .GT. zero .AND. soilcap(ji) .GT. min_sechiba) THEN
[2650]	1337
	1338	lwup_tmp(ji) = emis(ji) * c_stefan * temp_sol(ji)**4 + &
	1339	quatre * emis(ji) * c_stefan * temp_sol(ji)*3 (tp_00 - temp_sol(ji))
	1340
	1341	lwup_tmp(ji) = lwup_tmp(ji) + (un - emis(ji)) * lwdown(ji)
	1342	netrad_tmp(ji) = lwdown(ji) + swnet(ji) - lwup_tmp(ji)
	1343	fluxsens_tmp(ji) = rau(ji) * qc * cp_air * (tp_00 - epot_air(ji)/cp_air)
	1344	fluxlat_tmp(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) * (un-vbeta5(ji)) * &
	1345	& (qsol_sat_tmp(ji) - qair(ji)) + &
	1346	& chalev0 * rau(ji) * qc * vbeta5(ji) *&
	1347	& (qsol_sat_tmp(ji) - qair(ji)) + &
	1348	& chalev0 * rau(ji) * qc * (un - vbeta1(ji)) * (un-vbeta5(ji))* vbeta(ji) * &
[4197]	1349	& (qsol_sat_tmp(ji) - qair(ji))
[2650]	1350
	1351	zgflux(ji) = netrad_tmp(ji) - fluxsens_tmp(ji) - fluxlat_tmp(ji)+PHPSNOW(ji)
	1352
[3020]	1353	temp_sol_add(ji) = -(pgflux(ji) - zgflux(ji))*dt_sechiba/soilcap(ji)
[2650]	1354
	1355	pgflux(ji) = zgflux(ji)
	1356
	1357	ELSE
	1358
	1359	temp_sol_add(ji) = zero
	1360
	1361	ENDIF
	1362
[5470]	1363	ENDDO
[5431]	1364
[947]	1365	!! 3. Define qsat_air with the subroutine src_parameter:
[8]	1366
	1367	CALL qsatcalc(kjpindex, tair, pb, qsat_air)
	1368
	1369	CALL dev_qsatcalc(kjpindex, tair, pb, grad_qsat)
[2222]	1370
[2591]	1371	! grad_qsat(:)= (qsol_sat_new(:)- qsat_air(:)) / ((psnew(:) - epot_air(:)) / cp_air) ! * dt_sechiba
[947]	1372
[47]	1373	warning_correction(:)=.FALSE.
[8]	1374	DO ji=1,kjpindex
[947]	1375
	1376	!! \latexonly
	1377	!! \input{enerbilflux16.tex}
	1378	!! \endlatexonly
[8]	1379	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
[947]	1380
	1381	!! \latexonly
	1382	!! \input{enerbilflux17.tex}
	1383	!! \endlatexonly
[8]	1384	qc = speed * q_cdrag(ji)
[947]	1385
	1386	!! Derive the potential as defined by Penman & Monteith (Monteith, 1965) based on the correction
	1387	!! term developed by Chris Milly (1992).
[4624]	1388	IF ((evapot(ji) .GT. min_sechiba) .AND. ((psnew(ji) - epot_air(ji)) .NE. zero )) THEN
[947]	1389	!
	1390	!! \latexonly
	1391	!! \input{enerbilflux18.tex}
	1392	!! \endlatexonly
[8]	1393	correction = (quatre * emis(ji) * c_stefan * tair(ji)*3 + rau(ji) qc * cp_air + &
	1394	& chalev0 * rau(ji) * qc * grad_qsat(ji) * vevapp(ji) / evapot(ji) )
[947]	1395
	1396	!! \latexonly
	1397	!! \input{enerbilflux19.tex}
	1398	!! \endlatexonly
[8]	1399	IF (ABS(correction) .GT. min_sechiba) THEN
	1400	correction = chalev0 * rau(ji) * qc * grad_qsat(ji) * (un - vevapp(ji)/evapot(ji)) / correction
	1401	ELSE
[47]	1402	warning_correction(ji)=.TRUE.
[8]	1403	ENDIF
	1404	ELSE
	1405	correction = zero
	1406	ENDIF
	1407	correction = MAX (zero, correction)
[947]	1408
	1409	!! \latexonly
	1410	!! \input{enerbilflux20.tex}
	1411	!! \endlatexonly
[8]	1412	evapot_corr(ji) = evapot(ji) / (un + correction)
	1413
	1414	ENDDO
[47]	1415	IF ( ANY(warning_correction) ) THEN
	1416	DO ji=1,kjpindex
	1417	IF ( warning_correction(ji) ) THEN
	1418	WRITE(numout,*) ji,"Denominateur de la correction de milly nul! Aucune correction appliquee"
	1419	ENDIF
	1420	ENDDO
	1421	ENDIF
[2348]	1422	IF (printlev>=3) WRITE (numout,*) ' enerbil_flux done '
[8]	1423
	1424	END SUBROUTINE enerbil_flux
	1425
[947]	1426
	1427	!! ===========================================================================================================================
	1428	!! SUBROUTINE : enerbil_evapveg
[8]	1429	!!
[947]	1430	!>\BRIEF Splits total evaporation loss into individual process
	1431	!! components and calculates the assimilation.
	1432	!!
	1433	!! DESCRIPTION : Based on the estimation of the fluxes in enerbil_flux,
	1434	!! this routine splits the total evaporation into transpiration interception loss from vegetation,
	1435	!! bare soil evaporation and snow sublimation. It then calculates the assimilation.
	1436	!!
[2031]	1437	!! MAIN OUTPUT VARIABLE(S) : vevapsno, vevapnu, transpir, vevapwet
[947]	1438	!!
	1439	!! REFERENCE(S) :
	1440	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
	1441	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
	1442	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
	1443	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
	1444	!! Circulation Model. Journal of Climate, 6, pp.248-273
	1445	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
	1446	!! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11):
	1447	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
	1448	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
	1449	!! sur le cycle de l'eau, PhD Thesis, available (22/12/11):
	1450	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
	1451	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
	1452	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
	1453	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
	1454	!! general circulation models. Global and Planetary Change, 19, pp.261-276
	1455	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
	1456	!! Interscience Publishers
	1457	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
	1458	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
	1459	!!
	1460	!! FLOWCHART : None
	1461	!! \n
	1462	!_ ==============================================================================================================================
[8]	1463
[2591]	1464	SUBROUTINE enerbil_evapveg (kjpindex, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
[2031]	1465	& rau, u, v, q_cdrag, qair, humrel, vevapsno, vevapnu , vevapflo, vevapwet, &
[3959]	1466	& transpir, transpot, evapot)
[8]	1467
[947]	1468	!! 0. Variable and parameter declaration
	1469
	1470	!! 0.1 Input variables
[8]	1471
[947]	1472	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
	1473	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-)
	1474	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta4 !! Bare soil resistance (-)
	1475	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance
	1476	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density (kg m^{-3})
	1477	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind velocity in directions u and v (m s^{-1})
[4146]	1478	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
[947]	1479	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1})
[3406]	1480	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm/tstep)
[947]	1481	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1)
	1482	!!$ DS 15022011 humrel was used in a previous version of Orchidee, developped by Nathalie. Need to be discussed if it should be introduce again
	1483	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Interception resistance (-)
	1484	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Vegetation resistance (-)
	1485	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3pot !! Vegetation resistance for potential transpiration
	1486
	1487	!! 0.2 Output variables
[8]	1488
[947]	1489	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation (mm day^{-1})
	1490	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation (mm day^{-1})
[2222]	1491	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapflo !! Floodplains evaporation
[947]	1492	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration (mm day^{-1})
	1493	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpot !! Potential Transpiration
	1494	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception (mm day^{-1})
[8]	1495
[947]	1496	!! 0.3 Modified variables
[3883]	1497
[947]	1498	!! 0.4 Local variables
	1499
	1500	INTEGER(i_std) :: ji, jv
	1501	REAL(r_std), DIMENSION(kjpindex) :: xx
	1502	REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum, vbeta3sum
	1503	REAL(r_std) :: speed !! Speed (m s^{-1})
[2222]	1504
[947]	1505	!_ ==============================================================================================================================
	1506	! initialisation: utile pour calculer l'evaporation des floodplains dans lesquelles il y a de la vegetation
	1507	vbeta2sum(:) = 0.
	1508	vbeta3sum(:) = 0.
	1509	DO jv = 1, nvm
	1510	vbeta2sum(:) = vbeta2sum(:) + vbeta2(:,jv)
	1511	vbeta3sum(:) = vbeta3sum(:) + vbeta3(:,jv)
	1512	ENDDO
	1513
	1514	!! 1. Compute vevapsno (evaporation from snow) and vevapnu (bare soil evaporation)
[2222]	1515
[8]	1516	DO ji=1,kjpindex
	1517
[947]	1518	!! \latexonly
	1519	!! \input{enerbilevapveg1.tex}
	1520	!! \endlatexonly
[8]	1521	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
	1522
[947]	1523
	1524	!! 1.1 Snow sublimation
	1525	!! \latexonly
	1526	!! \input{enerbilevapveg2.tex}
	1527	!! \endlatexonly
[2591]	1528	vevapsno(ji) = (un - vbeta5(ji)) * vbeta1(ji) * dt_sechiba * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji))
[8]	1529
[947]	1530	!! 1.2 Bare soil evaporation
	1531	!! \latexonly
	1532	!! \input{enerbilevapveg3.tex}
	1533	!! \endlatexonly
[2591]	1534	vevapnu(ji) = (un - vbeta1(ji)) * (un-vbeta5(ji)) * vbeta4(ji) * dt_sechiba * rau(ji) * speed * q_cdrag(ji) &
[947]	1535	& * (qsol_sat_new(ji) - qair(ji))
[8]	1536	!
[947]	1537	! 1.3 floodplains evaporation - transpiration et interception prioritaires dans les floodplains
[8]	1538	!
[4067]	1539	vevapflo(ji) = vbeta5(ji) &
[2591]	1540	& * dt_sechiba * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji))
[8]	1541
	1542	END DO
	1543
[947]	1544	!! 2. Compute transpir (transpiration) and vevapwet (interception)
	1545
	1546	!! Preliminaries
[8]	1547	DO ji = 1, kjpindex
[947]	1548
	1549	!! \latexonly
	1550	!! \input{enerbilevapveg4.tex}
	1551	!! \endlatexonly
[8]	1552	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
[947]	1553
	1554	!! \latexonly
	1555	!! \input{enerbilevapveg5.tex}
	1556	!! \endlatexonly
[4067]	1557	xx(ji) = dt_sechiba * (un-vbeta1(ji)) * (un-vbeta5(ji)) * (qsol_sat_new(ji)-qair(ji)) * rau(ji) * speed * q_cdrag(ji)
[947]	1558
[8]	1559	ENDDO
[947]	1560
[8]	1561	DO jv=1,nvm
	1562	DO ji=1,kjpindex
[947]	1563
	1564	!! 2.1 Calculate interception loss
	1565	!! \latexonly
	1566	!! \input{enerbilevapveg6.tex}
	1567	!! \endlatexonly
[8]	1568	vevapwet(ji,jv) = xx(ji) * vbeta2(ji,jv)
	1569	!
[947]	1570	!! 2.2 Calculate transpiration
	1571	!! \latexonly
	1572	!! \input{enerbilevapveg7.tex}
	1573	!! \endlatexonly
[8]	1574	transpir (ji,jv) = xx(ji) * vbeta3(ji,jv)
[947]	1575
	1576	transpot(ji,jv) = xx(ji) * vbeta3pot(ji,jv)
	1577
[8]	1578	END DO
	1579	END DO
	1580
[947]	1581
[8]	1582
[2348]	1583	IF (printlev>=3) WRITE (numout,*) ' enerbil_evapveg done '
[8]	1584
	1585	END SUBROUTINE enerbil_evapveg
	1586
	1587	END MODULE enerbil

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

Context Navigation

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

Download in other formats: