Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

enerbil.f90

Last change on this file was 5323, checked in by chunjing.qiu, 7 years ago
Add oldpeat
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 107.3 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_sechiba/enerbil.f90

Download in other formats: