Context Navigation

enerbil.f90 @ 7760

Last change on this file since 7760 was 4493, checked in by jan.polcher, 8 years ago

Modifications :

routing_reg has now a minimal initialisation process which allows to build/read the routing network so we know on which processors the river outflow (estuaries) are. This allows to define an OASIS partition for the rivers.
The routing now also stores the outflow of the largest rivers to the ocean internaly.
A function exists in order to transfer this information into the OASIS interface routine (orchideeoasis).
There a oasis_put is made to send off the information to the ocean model. The flow values as well as the lon/lat of the estuaries are sent.
An attempt is implemented in driver2oasis to put some data into XIOS. It does not work yet !
Correction of some variables in file_def_orchidee.xml

Remaining issues :

Synchronise of driver2oasis and orchidee before the xios_orchidee_finalize or oasis_terminate.

Property svn:keywords set to Revision Date HeadURL Date Author Revision

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_sechiba/enerbil.f90 @ 7760

Download in other formats: