Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

enerbil.f90 @ 7329

Last change on this file since 7329 was 7191, checked in by haicheng.zhang, 3 years ago
Implementing latral transfers of sediment and POC from land to ocean through river in ORCHILEAK
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 89.1 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE-Clateral/src_sechiba/enerbil.f90 @ 7329

Download in other formats: