Context Navigation

enerbil.f90 @ 5925

Last change on this file since 5925 was 3313, checked in by josefine.ghattas, 8 years ago

Move the part for ok_exlicitsnow from enerbil_fusion to the the end of thermosoil_main and thermosoilc_main. This does not make any change in the order of calculations.
temp_sol_new is now modified in thermosoil_main/thermosoilc_main (intent(inout)) but only for the case ok_explicitsnow
call enerbil_fusion after hydrol_main as done before changeset [3269]
protected variable fusion from writing for case ok_explicitsnow. This variable is calculated by enerbil_fusion

See also ticket #229

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

File size: 93.9 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 day^{-1})
145	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: evapot_corr !! Soil Potential Evaporation Correction (mm day^{-1})
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, snowrho, 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
402	!! 0.2 Output variables
403
404	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation (mm day^{-1})
405	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation (mm day^{-1})
406	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapflo !! Floodplains evaporation
407	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration (mm day^{-1})
408	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpot !! Potential transpiration
409	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception (mm day^{-1})
410	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: t2mdiag !! Diagnostic temperature at 2m above ground (K)
411	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K)
412	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
413	!! 0.3 Modified variables
414
415	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation (mm day^{-1})
416	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction (mm day^{-1})
417	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol !! Soil temperature (K)
418	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: qsurf !! Surface specific humidity (kg kg^{-1})
419	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2})
420	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2})
421	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tsol_rad !! Tsol_rad (W m^{-2})
422	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1})
423	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: pgflux !! Net energy into snowpack(W/m^2)
424	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density at each snow layer
425	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth at each snow layer
426
427	!! 0.4 Local variables
428	REAL(r_std),DIMENSION (kjpindex) :: epot_air_new, qair_new
429	INTEGER(i_std) :: ji,ii,iv
430
431
432	!_ ================================================================================================================================
433
434	!! 1. Computes some initialisation variables
435
436	! Computes some initialisation variables psold (the old surface static energy), qsol_sat
437	! (the saturated surface humidity) and pdqsold (the derivative of the saturated surface humidity,
438	! calculated with respect to the surface temperature at the 'old' timestep)
439	CALL enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis)
440
441
442	!! 2. Computes psnew (the surface static energy at the 'new' timestep)
443
444	! Computes psnew (the surface static energy at the 'new' timestep), qsol_sat_new (the surface
445	! saturated humidity at the 'new' timestep), temp_sol_new (the surface temperature at the 'new'
446	! timestep), qair_new (the lowest atmospheric humidity at the 'new' timestep) and epot_air_new
447	! (the lowest atmospheric evaporation potential at the 'new' timestep)
448	CALL enerbil_surftemp (kjpindex, zlev, emis, &
449	& epot_air, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,&
450	& valpha, vbeta1, vbeta5, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, &
451	& qair_new, epot_air_new, &
452	& snowdz)
453
454
455	! 3. computes surface temperature and humidity for a saturated surface
456	!
457	CALL enerbil_pottemp (kjpindex, zlev, emis, &
458	& epot_air, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,&
459	& valpha, vbeta1, vbeta5, soilcap, lwdown, swnet, q_sol_pot, temp_sol_pot)
460
461
462	!! 4. Diagnoses components of the energy budget
463
464	! Diagnoses lwup (upwards long wave radiation), lwnet (net long wave radiation), tsol_rad (radiative
465	! ground temperature), netrad (net radiation), qsurf (surface humidity), vevapp (total evaporation),
466	! evapot (evaporation potential), evapot_corr (evaporation potential correction factor),
467	! fluxlat (latent heat flux), fluxsubli (sublimination heat flux) and fluxsens (sensible heat flux).
468
469	CALL enerbil_flux (kjpindex, emis, temp_sol, rau, u, v, q_cdrag, vbeta, valpha, vbeta1, vbeta5, &
470	& qair_new, epot_air_new, psnew, qsurf, &
471	& fluxsens , fluxlat , fluxsubli, vevapp, temp_sol_new, lwdown, swnet, &
472	& lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr, &
473	& precip_rain,snowdz,temp_air,pgflux, soilcap, temp_sol_add)
474
475	!! 5. Diagnoses the values for evaporation and transpiration
476
477	! Diagnoses the values for evaporation and transpiration: vevapsno (snow evaporation), vevapnu
478	! (bare soil evaporation), transpir (transpiration) and vevapwet (interception)
479	CALL enerbil_evapveg (kjpindex, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
480	& rau, u, v, q_cdrag, qair_new, humrel, vevapsno, vevapnu , vevapflo, vevapwet, &
481	& transpir, transpot, evapot, &
482	& snowdz,snowrho,evapot_corr)
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*one_day/dt_sechiba)
494	CALL xios_orchidee_send_field("evapot_alma",evapot/dt_sechiba)
495	CALL xios_orchidee_send_field("evapot_corr",evapot_corr*one_day/dt_sechiba)
496	CALL xios_orchidee_send_field("evapot_corr_alma",evapot_corr/dt_sechiba)
497	CALL xios_orchidee_send_field("lwdown",lwabs)
498	CALL xios_orchidee_send_field("lwnet",lwnet)
499	CALL xios_orchidee_send_field("Qv",fluxsubli)
500	CALL xios_orchidee_send_field("PotSurfT",temp_sol_pot)
501
502
503	IF ( .NOT. almaoutput ) THEN
504	CALL histwrite_p(hist_id, 'netrad', kjit, netrad, kjpindex, index)
505	CALL histwrite_p(hist_id, 'evapot', kjit, evapot, kjpindex, index)
506	CALL histwrite_p(hist_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index)
507	CALL histwrite_p(hist_id, 'lwdown', kjit, lwabs, kjpindex, index)
508	CALL histwrite_p(hist_id, 'lwnet', kjit, lwnet, kjpindex, index)
509	IF ( hist2_id > 0 ) THEN
510	CALL histwrite_p(hist2_id, 'netrad', kjit, netrad, kjpindex, index)
511	CALL histwrite_p(hist2_id, 'evapot', kjit, evapot, kjpindex, index)
512	CALL histwrite_p(hist2_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index)
513	CALL histwrite_p(hist2_id, 'lwdown', kjit, lwabs, kjpindex, index)
514	CALL histwrite_p(hist2_id, 'lwnet', kjit, lwnet, kjpindex, index)
515	ENDIF
516	ELSE
517	CALL histwrite_p(hist_id, 'LWnet', kjit, lwnet, kjpindex, index)
518	CALL histwrite_p(hist_id, 'Qv', kjit, fluxsubli, kjpindex, index)
519	CALL histwrite_p(hist_id, 'PotEvap', kjit, evapot_corr, kjpindex, index)
520	CALL histwrite_p(hist_id, 'PotEvapOld', kjit, evapot, kjpindex, index)
521	CALL histwrite_p(hist_id, 'PotSurfT', kjit, temp_sol_pot, kjpindex, index)
522	IF ( hist2_id > 0 ) THEN
523	CALL histwrite_p(hist2_id, 'LWnet', kjit, lwnet, kjpindex, index)
524	CALL histwrite_p(hist2_id, 'Qv', kjit, fluxsubli, kjpindex, index)
525	CALL histwrite_p(hist2_id, 'PotEvap', kjit, evapot_corr, kjpindex, index)
526	ENDIF
527	ENDIF
528
529	IF (printlev>=3) WRITE (numout,*) ' enerbil_main Done '
530
531	END SUBROUTINE enerbil_main
532
533
534	!! =============================================================================================================================
535	!! SUBROUTINE: enerbil_finalize
536	!!
537	!>\BRIEF Write to restart file
538	!!
539	!! DESCRIPTION: This subroutine writes the module variables and variables calculated in enerbil
540	!! to restart file
541	!!
542	!! RECENT CHANGE(S): None
543	!!
544	!! REFERENCE(S): None
545	!!
546	!! FLOWCHART: None
547	!! \n
548	!_ ==============================================================================================================================
549	SUBROUTINE enerbil_finalize (kjit, kjpindex, rest_id, &
550	evapot, evapot_corr, temp_sol, tsol_rad, &
551	qsurf, fluxsens, fluxlat, vevapp )
552
553	!! 0 Variable and parameter description
554	!! 0.1 Input variables
555	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
556	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
557	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (-)
558	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm day^{-1})
559	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation Correction (mm day^{-1})
560	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature (K)
561	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tsol_rad !! Tsol_rad (W m^{-2})
562	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Surface specific humidity (kg kg^{-1})
563	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxsens !! Sensible heat flux (W m^{-2})
564	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxlat !! Latent heat flux (W m^{-2})
565	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapp !! Total of evaporation (mm day^{-1})
566
567	!_ ================================================================================================================================
568
569	!! 1. Write variables to restart file to be used for the next simulation
570	IF (printlev>=3) WRITE (numout,*) 'Write restart file with ENERBIL variables'
571
572	CALL restput_p(rest_id, 'temp_sol', nbp_glo, 1, 1, kjit, temp_sol, 'scatter', nbp_glo, index_g)
573	CALL restput_p(rest_id, 'qsurf', nbp_glo, 1, 1, kjit, qsurf, 'scatter', nbp_glo, index_g)
574	CALL restput_p(rest_id, 'evapot', nbp_glo, 1, 1, kjit, evapot, 'scatter', nbp_glo, index_g)
575	CALL restput_p(rest_id, 'evapot_corr', nbp_glo, 1, 1, kjit, evapot_corr, 'scatter', nbp_glo, index_g)
576	CALL restput_p(rest_id, 'tsolrad', nbp_glo, 1, 1, kjit, tsol_rad, 'scatter', nbp_glo, index_g)
577	CALL restput_p(rest_id, 'evapora', nbp_glo, 1, 1, kjit, vevapp, 'scatter', nbp_glo, index_g)
578	CALL restput_p(rest_id, 'fluxlat', nbp_glo, 1, 1, kjit, fluxlat, 'scatter', nbp_glo, index_g)
579	CALL restput_p(rest_id, 'fluxsens', nbp_glo, 1, 1, kjit, fluxsens, 'scatter', nbp_glo, index_g)
580	CALL restput_p(rest_id, 'tempsolpot', nbp_glo, 1, 1, kjit, temp_sol_pot, 'scatter', nbp_glo, index_g)
581	CALL restput_p(rest_id, 'qsolpot', nbp_glo, 1, 1, kjit, q_sol_pot, 'scatter', nbp_glo, index_g)
582
583	END SUBROUTINE enerbil_finalize
584
585
586
587
588	!! =============================================================================================================================
589	!! SUBROUTINE : enerbil_clear
590	!!
591	!>\BRIEF Routine deallocates clear output variables if already allocated.
592	!!
593	!! DESCRIPTION : This is a 'housekeeping' routine that deallocates the key output
594	!! variables, if they have already been allocated. The variables that are deallocated are psold,
595	!! qsol_sat, pdqsold, psnew, qsol_sat_new, netrad, lwabs, lwup, lwnet, fluxsubli, qsat_air, tair
596	!!
597	!! RECENT CHANGE(S) : None
598	!!
599	!! MAIN OUTPUT VARIABLE(S) : None
600	!!
601	!! REFERENCES : None
602	!!
603	!! FLOWCHART : None
604	!! \n
605	!_ ==============================================================================================================================
606
607	SUBROUTINE enerbil_clear ()
608	IF ( ALLOCATED (psold)) DEALLOCATE (psold)
609	IF ( ALLOCATED (qsol_sat)) DEALLOCATE (qsol_sat)
610	IF ( ALLOCATED (pdqsold)) DEALLOCATE (pdqsold)
611	IF ( ALLOCATED (psnew)) DEALLOCATE (psnew)
612	IF ( ALLOCATED (qsol_sat_new)) DEALLOCATE (qsol_sat_new)
613	IF ( ALLOCATED (netrad)) DEALLOCATE (netrad)
614	IF ( ALLOCATED (lwabs)) DEALLOCATE (lwabs)
615	IF ( ALLOCATED (lwup)) DEALLOCATE (lwup)
616	IF ( ALLOCATED (lwnet)) DEALLOCATE (lwnet)
617	IF ( ALLOCATED (fluxsubli)) DEALLOCATE (fluxsubli)
618	IF ( ALLOCATED (qsat_air)) DEALLOCATE (qsat_air)
619	IF ( ALLOCATED (tair)) DEALLOCATE (tair)
620	IF ( ALLOCATED (q_sol_pot)) DEALLOCATE (q_sol_pot)
621	IF ( ALLOCATED (temp_sol_pot)) DEALLOCATE (temp_sol_pot)
622
623	END SUBROUTINE enerbil_clear
624
625
626	!! =============================================================================================================================
627	!! SUBROUTINE : enerbil_begin
628	!!
629	!>\BRIEF Preliminary variables required for the calculation of
630	!! the energy budget are derived.
631	!!
632	!! DESCRIPTION : This routines computes preliminary variables required
633	!! for the calculation of the energy budget: the old surface static energy (psold), the surface saturation
634	!! humidity (qsol_sat), the derivative of satured specific humidity at the old temperature (pdqsold)
635	!! and the net radiation (netrad).
636	!!
637	!! MAIN OUTPUT VARIABLE(S) : psold, qsol_sat, pdqsold, netrad
638	!!
639	!! REFERENCE(S) :
640	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
641	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
642	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
643	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
644	!! Circulation Model. Journal of Climate, 6, pp.248-273
645	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
646	!! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11):
647	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
648	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
649	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
650	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
651	!! general circulation models. Global and Planetary Change, 19, pp.261-276
652	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
653	!! Interscience Publishers\n
654	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
655	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
656	!!
657	!! FLOWCHART : None
658	!! \n
659	!_ ==============================================================================================================================
660
661	SUBROUTINE enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis)
662
663	!! 0. Variable and parameter declaration
664
665	!! 0.1 Input variables
666
667	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
668	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature (K)
669	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2})
670	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2})
671	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
672	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-)
673
674	!! 0.2 Output variables
675
676	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: psold !! Old surface dry static energy (J kg^{-1})
677	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat !! Saturated specific humudity for old temperature
678	!! (kg kg^{-1})
679	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: pdqsold !! Derivative of satured specific humidity at the old
680	!! temperature (kg (kg s)^{-1})
681	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: netrad !! Net radiation (W m^{-2})
682
683	!! 0.3 Modified variables
684
685	!! 0.4 Local variables
686
687	INTEGER(i_std) :: ji
688	REAL(r_std), DIMENSION(kjpindex) :: dev_qsol
689	REAL(r_std), PARAMETER :: missing = 999998.
690
691	!_ ================================================================================================================================
692
693	!! 1. Computes psold (the surface static energy for the old timestep)
694
695	!! We here define the surface static energy for the 'old' timestep, in terms of the surface
696	!! temperature and heat capacity.
697	!! \latexonly
698	!! \input{enerbilbegin1.tex}
699	!! \endlatexonly
700	psold(:) = temp_sol(:)*cp_air
701	!! 2. Computes qsol_sat (the surface saturated humidity).
702
703	!! We call the routine 'qsatcalc' from within the module 'src_parameters/constantes_veg'.
704	CALL qsatcalc (kjpindex, temp_sol, pb, qsol_sat)
705	IF ( diag_qsat ) THEN
706	IF ( ANY(ABS(qsol_sat(:)) .GT. missing) ) THEN
707	DO ji = 1, kjpindex
708	IF ( ABS(qsol_sat(ji)) .GT. missing) THEN
709	WRITE(numout,*) 'ERROR on ji = ', ji
710	WRITE(numout,*) 'temp_sol(ji), pb(ji) :', temp_sol(ji), pb(ji)
711	CALL ipslerr_p (3,'enerbil_begin', &
712	& 'qsol too high ','','')
713	ENDIF
714	ENDDO
715	ENDIF
716	ENDIF
717
718	!! 3. Computes pdqsold
719
720	!! Computes pdqsold (the derivative of the saturated humidity with respect to temperature
721	!! analysed at the surface temperature at the 'old' timestep.
722	!! We call the routine 'dev_qsatcalc' from qsat_moisture module.
723	CALL dev_qsatcalc (kjpindex, temp_sol, pb, dev_qsol)
724
725	!! \latexonly
726	!! \input{enerbilbegin2.tex}
727	!! \endlatexonly
728	pdqsold(:) = dev_qsol(:)
729
730
731	IF ( diag_qsat ) THEN
732	IF ( ANY(ABS( pdqsold(:)) .GT. missing) ) THEN
733	DO ji = 1, kjpindex
734	IF ( ABS( pdqsold(ji)) .GT. missing ) THEN
735	WRITE(numout,*) 'ERROR on ji = ', ji
736	WRITE(numout,*) 'temp_sol(ji), pb(ji) :', temp_sol(ji), pb(ji)
737	CALL ipslerr_p (3,'enerbil_begin', &
738	& 'pdqsold too high ','','')
739	ENDIF
740	ENDDO
741	ENDIF
742	ENDIF
743
744	!! 4. Computes the net radiation and the absorbed LW radiation absorbed at the surface.
745
746	!! Long wave radiation absorbed by the surface is the product of emissivity and downwelling LW radiation
747	!! \latexonly
748	!! \input{enerbilbegin3.tex}
749	!! \endlatexonly
750	lwabs(:) = emis(:) * lwdown(:)
751
752	!! Net radiation is calculated as:
753	!! \latexonly
754	!! \input{enerbilbegin4.tex}
755	!! \endlatexonly
756	netrad(:) = lwdown(:) + swnet (:) - (emis(:) * c_stefan * temp_sol(:)*4 + (un - emis(:)) lwdown(:))
757
758	IF (printlev>=3) WRITE (numout,*) ' enerbil_begin done '
759
760	END SUBROUTINE enerbil_begin
761
762
763	!! =============================================================================================================================
764	!! SUBROUTINE : enerbil_surftemp
765	!!
766	!>\BRIEF This routine computes the new surface static energy
767	!! (psnew) and the saturated humidity at the surface (qsol_sat_new).
768	!!
769	!! DESCRIPTION : This is the key part of the enerbil module, for which
770	!! the energy budget in the surface layer is solved and changes over the model timestep 'dt_sechiba' are
771	!! quantified for the surface static energy, surface temperature, surface humidity, the 'air' (or lowest
772	!! level atmospheric model) temperature, the 'air' (or lowest level atmospheric model) humidity,
773	!! according to the method that is laid out by Dufresne \& Ghattas (2009) and Shultz et al. (2001).
774	!!
775	!! It computes the energy balance at the surface with an implicit scheme, that is connected to the
776	!! Richtmyer and Morton algorithm of the PBL. By computing the surface temperature and surface humidity
777	!! the routine also implicitly estimates the various fluxes which balance in order to give the new
778	!! temperature. Thus once the surface temperature has been obtained all the fluxes need to be diagnosed.
779	!!
780	!! If ok_explicitsnow is used, the notion of skin layer in ORCHIDEE is abandoned and
781	!! the first thin snow layer is adopted to solve the surface energy budget.
782	!! Implicit method is still used for coupling to the atmosphere.
783	!! In this new scheme, the snow temperature profile is simulatenously solved
784	!! along with the surface snow temperature, and the details are referenced to Boone et al. (2010)
785	!!
786	!! MAIN OUTPUT VARIABLE(S) : psnew, qsol_sat_new, temp_sol_new, qair_new, epot_air_new
787	!!
788	!! REFERENCE(S) :
789	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
790	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
791	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
792	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
793	!! Circulation Model. Journal of Climate, 6, pp.248-273
794	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
795	!! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11):
796	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
797	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
798	!! sur le cycle de l'eau, PhD Thesis, available (22/12/11):
799	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pd
800	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
801	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
802	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
803	!! general circulation models. Global and Planetary Change, 19, pp.261-276
804	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
805	!! Interscience Publishers
806	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
807	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
808	!!
809	!! FLOWCHART : None
810	!!
811	!_ ==============================================================================================================================
812
813
814	SUBROUTINE enerbil_surftemp (kjpindex, zlev, emis, epot_air, &
815	& petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,&
816	& valpha, vbeta1, vbeta5, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, &
817	& qair_new, epot_air_new, &
818	& snowdz)
819
820
821
822	!! 0. Variable and parameter declaration
823
824	!! 0.1 Input variables
825
826	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
827	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
828	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-)
829	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (?? J)
830	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef (see note)
831	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef (see note)
832	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1})
833	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef (see note)
834	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef (see note)
835	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Soil flux (W m^{-2})
836	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3})
837	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind velocity by directional components
838	!! u (Eastwards) and v (Northwards) (m s^{-1})
839	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag (m s^{-1})
840	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-)
841	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: valpha !! Resistance coefficient (-)
842	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-)
843	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance
844	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity (J K^{-1])
845	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2})
846	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2})
847	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth
848
849	!! 0.2 Output variables
850
851	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: psnew !! New surface static energy (J kg^{-1})
852	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat_new !! New saturated surface air moisture (kg kg^{-1})
853	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K)
854	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qair_new !! New air moisture (kg kg^{-1})
855	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: epot_air_new !! New air temperature (K)
856
857
858	!! 0.4 Local variables
859
860	INTEGER(i_std) :: ji
861	REAL(r_std),DIMENSION (kjpindex) :: zicp
862	REAL(r_std) :: fevap
863	REAL(r_std) :: sensfl_old, larsub_old, lareva_old, dtheta, sum_old, sum_sns
864	REAL(r_std) :: zikt, zikq, netrad_sns, sensfl_sns, larsub_sns, lareva_sns
865	REAL(r_std) :: speed
866
867	!_ ================================================================================================================================
868
869	zicp = un / cp_air
870
871	DO ji=1,kjpindex
872	!! 1. Derivation of auxiliary variables
873
874	!! \latexonly
875	!! \input{enerbilsurftemp1.tex}
876	!! \endlatexonly
877	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
878	!
879	!! \latexonly
880	!! \input{enerbilsurftemp2.tex}
881	!! \endlatexonly
882	zikt = 1/(rau(ji) * speed * q_cdrag(ji))
883	zikq = 1/(rau(ji) * speed * q_cdrag(ji))
884
885	!! 2. Computation of fluxes for the old surface conditions
886
887	!! As a reference, we first calculate the sensible and latent heat for the 'old' timestep.
888
889	!! 2.1 Sensible heat for the 'old' timestep
890	!! This is equation (64) of (Dufresne & Ghattas, 2009), rewritten in terms of H_old. We make the
891	!! approximation that (P_r/P)^{kappa}=1 and convert the A and B coefficients to the ORCHIDEE
892	!! format (see introductory note, above). Also, the equation here is in terms of surface static
893	!! energy, which is defined as per section 3.1, above.
894	!! \latexonly
895	!! \input{enerbilsurftemp3.tex}
896	!! \endlatexonly
897	sensfl_old = (petBcoef(ji) - psold(ji)) / (zikt - petAcoef(ji))
898	!! 2.2 Latent heat for the 'old' timestep
899	!! This is equation (70) of (Dufresne & Ghattas, 2009), rewritten in terms of {\lambda E}_{old}.
900	!! Again we convert the A and B coefficients from the LMDZ format to the ORCHIDEE format.\n
901	!! There are two forms of the equation - first for the latent heat as a result of sublimination
902	!! processes and then as a result of evaporative processes:
903	!! \latexonly
904	!! \input{enerbilsurftemp4.tex}
905	!! \endlatexonly
906	larsub_old = chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) * &
907	(peqBcoef(ji) - qsol_sat(ji)) / (zikq - vbeta1(ji) * (un - vbeta5(ji))* peqAcoef(ji))
908	!! \latexonly
909	!! \input{enerbilsurftemp5.tex}
910	!! \endlatexonly
911	lareva_old = chalev0 * (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * &
912	(peqBcoef(ji) - valpha(ji) * qsol_sat(ji)) / &
913	(zikq - (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * peqAcoef(ji)) &
914	+ chalev0 * vbeta5(ji) * (peqBcoef(ji) - qsol_sat(ji)) / (zikq - vbeta5(ji) * peqAcoef(ji))
915
916	!! 3. Calculation of sensitivity terms
917
918	!! We here calculate the sensitivity for the different fluxes to changes in dtheta, which is the
919	!! change in the surface static energy over the model time step (dt_sechiba).
920
921	!! 3.1 Net radiation sensitivity
922	!! This is an explicit-implicit representation of the Stefan-Boltzmann law - the explicit terms
923	!! are ${ps}_{old}^3$, and it is completed when multiplied by dtheta (defined above).
924	!! \latexonly
925	!! \input{enerbilsurftemp6.tex}
926	!! \endlatexonly
927	netrad_sns = zicp(ji) * quatre * emis(ji) * c_stefan * ((zicp(ji) * psold(ji))**3)
928
929	!! 3.2 Sensible heat flux sensitivity
930	!! This is the temperature sensitivity term N_1^h derived in equation (66) of (Dufresne & Ghattas, 2009),
931	!! where we again assume that (P_r/P)^{kappa}=1, convert the A and B coefficients to the ORCHIDEE format,
932	!! and rewrite in terms of surface static energy, rather than temperature (see section 3.1, above).
933	!! \latexonly
934	!! \input{enerbilsurftemp7.tex}
935	!! \endlatexonly
936	sensfl_sns = un / (zikt - petAcoef(ji))
937	!! 3.3 Latent heat flux sensitivity
938	!! This is the humidity sensitivity term N_1^q derived in equation (72) of (Dufresne & Ghattas, 2009), where
939	!! the A and B coefficients are written in the ORCHIDEE format.
940	!! larsub_sns is the latent heat sensitivity for sublimination. The coefficient vbeta1 is the evaporation
941	!! coefficient. It represents the relationship between the real and potential evaporation. It is derived
942	!! in the module 'src_sechiba/diffuco_snow'.
943	!! \latexonly
944	!! \input{enerbilsurftemp8.tex}
945	!! \endlatexonly
946	larsub_sns = chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) * zicp(ji) * &
947	pdqsold(ji) / (zikq - vbeta1(ji) * (un - vbeta5(ji)) * peqAcoef(ji))
948
949	!! lareva_sns is the latent heat sensitivity for evaporation. vbeta1 (src_sechiba/diffuco_snow),
950	!! vbeta (src_sechiba/diffuco_comb) and valpha (src_sechiba/diffuco_comb) are the evaporation
951	!! coefficients.
952	!! \latexonly
953	!! \input{enerbilsurftemp9.tex}
954	!! \endlatexonly
955	lareva_sns = chalev0 * ((un - vbeta1(ji))(un - vbeta5(ji)) vbeta(ji) * valpha(ji) + vbeta5(ji)) * &
956	& zicp(ji) * pdqsold(ji) / &
957	(zikq - ((un - vbeta1(ji))(un - vbeta5(ji)) vbeta(ji) * valpha(ji) + vbeta5(ji))* peqAcoef(ji))
958
959	!! 4. Solution of the energy balance
960
961	!! 4.1 We calculate the total flux for the 'old' timestep.
962	!! \latexonly
963	!! \input{enerbilsurftemp10.tex}
964	!! \endlatexonly
965	sum_old = netrad(ji) + sensfl_old + larsub_old + lareva_old + soilflx(ji)
966
967	!! 4.2 We calculate the total sensitivity dtheta (the change in the
968	!! surface static energy over the timestep.
969	!! \latexonly
970	!! \input{enerbilsurftemp11.tex}
971	!! \endlatexonly
972	sum_sns = netrad_sns + sensfl_sns + larsub_sns + lareva_sns
973	!! 4.3 Calculation of dtheta (change in surface static energy over the
974	!! timestep.
975	!! \latexonly
976	!! \input{enerbilsurftemp12.tex}
977	!! \endlatexonly
978
979	dtheta = dt_sechiba * sum_old / (zicp(ji) * soilcap(ji) + dt_sechiba * sum_sns)
980
981	!! 4.4 Determination of state variables for the 'new' timestep
982	!! No we have dtheta, we can determine the surface static energy that
983	!! corresponds to the 'new' timestep.
984	!! \latexonly
985	!! \input{enerbilsurftemp13.tex}
986	!! \endlatexonly
987	psnew(ji) = psold(ji) + dtheta
988
989	!! The new surface saturated humidity can be calculated by equation (69)
990	!! of (Dufresne & Ghattas, 2009), in which we substitute dtheta for the
991	!! change between old and new temperature using the relationship from 3.1,
992	!! above.
993	!! \latexonly
994	!! \input{enerbilsurftemp14.tex}
995	!! \endlatexonly
996	qsol_sat_new(ji) = qsol_sat(ji) + zicp(ji) * pdqsold(ji) * dtheta
997	!! The new surface temperature is determined from the new surface static temperature,
998	!! again by using the relationship in section 3.1.
999	!! \latexonly
1000	!! \input{enerbilsurftemp15.tex}
1001	!! \endlatexonly
1002	temp_sol_new(ji) = psnew(ji) / cp_air
1003
1004	!! 4.5 Calculation of new evaporation potential and new evaporation latent heat
1005	!! flux (???)
1006	!! \latexonly
1007	!! \input{enerbilsurftemp16.tex}
1008	!! \endlatexonly
1009	epot_air_new(ji) = zikt * (sensfl_old - sensfl_sns * dtheta) + psnew(ji)
1010
1011	!! \latexonly
1012	!! \input{enerbilsurftemp17.tex}
1013	!! \endlatexonly
1014	fevap = (lareva_old - lareva_sns * dtheta) + (larsub_old - larsub_sns * dtheta)
1015
1016	IF ( ABS(fevap) < EPSILON(un) ) THEN
1017
1018	!! \latexonly
1019	!! \input{enerbilsurftemp18.tex}
1020	!! \endlatexonly
1021	qair_new(ji) = qair(ji)
1022	ELSE
1023	!! \latexonly
1024	!! \input{enerbilsurftemp19.tex}
1025	!! \endlatexonly
1026	qair_new(ji) = zikq * un / ( chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) + &
1027	& chalev0 * ((un - vbeta1(ji))(un - vbeta5(ji)) vbeta(ji) * valpha(ji) + vbeta5(ji)) ) &
1028	& * fevap + qsol_sat_new(ji)
1029	ENDIF
1030
1031	ENDDO
1032
1033	IF (printlev>=3) WRITE (numout,*) ' enerbil_surftemp done '
1034
1035	END SUBROUTINE enerbil_surftemp
1036
1037
1038	!! =============================================================================================================================
1039	!! SUBROUTINE : enerbil_pottemp
1040	!!
1041	!>\BRIEF This subroutine computes the surface temperature and humidity should the surface been unstressed
1042	!!
1043	!! DESCRIPTION :
1044	!!
1045	!! MAIN OUTPUT VARIABLE(S) :
1046	!!
1047	!! REFERENCE(S) :
1048	!!
1049	!! FLOWCHART : None
1050	!!
1051	!_ ==============================================================================================================================
1052
1053	SUBROUTINE enerbil_pottemp (kjpindex, zlev, emis, epot_air, &
1054	& petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,&
1055	& valpha, vbeta1, vbeta5, soilcap, lwdown, swnet, q_sol_pot, temp_sol_pot)
1056
1057	! interface
1058	! input scalar
1059	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1060	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
1061	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity
1062	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy
1063	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef
1064	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef
1065	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
1066	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef
1067	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef
1068	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Soil flux
1069	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
1070	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind
1071	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !!
1072	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient
1073	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: valpha !! Resistance coefficient
1074	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance
1075	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance
1076	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity
1077	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux
1078	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
1079	! output fields
1080	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: q_sol_pot !! Potential surface air moisture
1081	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_pot !! Potential soil temperature
1082
1083
1084	! local declaration
1085	INTEGER(i_std) :: ji
1086	REAL(r_std),DIMENSION (kjpindex) :: zicp
1087	REAL(r_std) :: fevap
1088	REAL(r_std) :: sensfl_old, larsub_old, lareva_old, dtheta, sum_old, sum_sns
1089	REAL(r_std) :: zikt, zikq, netrad_sns, sensfl_sns, larsub_sns, lareva_sns
1090	REAL(r_std) :: speed
1091
1092	!_ ==============================================================================================================================
1093
1094	zicp = un / cp_air
1095	!
1096	DO ji=1,kjpindex
1097
1098	dtheta = zero
1099	fevap = zero
1100
1101	temp_sol_pot(ji) = temp_sol_pot(ji) + dtheta
1102
1103	q_sol_pot(ji) = q_sol_pot(ji) + fevap
1104
1105	ENDDO
1106
1107	END SUBROUTINE enerbil_pottemp
1108
1109
1110	!! =============================================================================================================================
1111	!! SUBROUTINE : enerbil_flux
1112	!!
1113	!>\BRIEF Computes the new soil temperature, net radiation and
1114	!! latent and sensible heat flux for the new time step.
1115	!!
1116	!! DESCRIPTION : This routine diagnoses, based on the new soil temperature,
1117	!! the net radiation, the total evaporation, the latent heat flux, the sensible heat flux and the
1118	!! sublimination flux. It also diagnoses the potential evaporation used for the fluxes (evapot) and the potential
1119	!! as defined by Penman & Monteith (Monteith, 1965) based on the correction term developed by Chris
1120	!! Milly (1992). This Penman-Monteith formulation is required for the estimation of bare soil evaporation
1121	!! when the 11 layer CWRR moisture scheme is used for the soil hydrology.
1122	!!
1123	!! MAIN OUTPUT VARIABLE(S) : qsurf, fluxsens, fluxlat, fluxsubli, vevapp, lwup, lwnet,
1124	!! tsol_rad, netrad, evapot, evapot_corr
1125	!!
1126	!! REFERENCE(S) :
1127	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
1128	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
1129	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1130	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1131	!! Circulation Model. Journal of Climate, 6, pp.248-273
1132	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
1133	!! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11):
1134	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
1135	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1136	!! sur le cycle de l'eau, PhD Thesis, available (22/12/11):
1137	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1138	!! - Monteith, JL, 1965. Evaporation and Environment, paper presented at Symposium of the Society
1139	!! for Experimental Biology
1140	!! - Monteith & Unsworth, 2008. Principles of Environmental Physics (third edition), published Elsevier
1141	!! ISBN 978-0-12-505103-3
1142	!! - Milly, P. C. D., 1992: Potential Evaporation and Soil Moisture in General Circulation Models.
1143	!! Journal of Climate, 5, pp. 209â226.
1144	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
1145	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
1146	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
1147	!! general circulation models. Global and Planetary Change, 19, pp.261-276
1148	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
1149	!! Interscience Publishers
1150	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
1151	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
1152	!!
1153	!! FLOWCHART : None
1154	!! \n
1155	!_ ==============================================================================================================================
1156
1157	SUBROUTINE enerbil_flux (kjpindex, emis, temp_sol, rau, u, v, q_cdrag, vbeta, valpha, vbeta1, vbeta5, &
1158	& qair, epot_air, psnew, qsurf, fluxsens, fluxlat, fluxsubli, vevapp, temp_sol_new, &
1159	& lwdown, swnet, lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr,&
1160	& precip_rain,snowdz,temp_air,pgflux, soilcap, temp_sol_add)
1161
1162	!! 0. Variable and parameter declaration
1163
1164	!! 0.1 Input variables
1165
1166	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1167	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-)
1168	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Surface temperature (K)
1169	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density (kg m^{-3})
1170	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u,v !! Wind velocity in components u and v (m s^{-1})
1171	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag (m s^{-1})
1172	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-)
1173	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: valpha !! Resistance coefficient (-)
1174	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-)
1175	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Flood resistance
1176	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1})
1177	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (J)
1178	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: psnew !! New surface static energy (J kg^{-1})
1179	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature (K)
1180	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
1181	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth
1182	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: precip_rain !! Rainfall
1183	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: temp_air !! Air temperature
1184	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Downward long wave radiation (W m^{-2})
1185	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net short wave radiation (W m^{-2})
1186	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil calorific capacity including snow and soil (J K^{-1})
1187
1188	!! 0.2 Output variables
1189
1190	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! Surface specific humidity (kg kg^{-1})
1191	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2})
1192	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2})
1193	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsubli !! Energy of sublimation (W m^{-2})
1194	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1})
1195	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwup !! Long-wave up radiation (W m^{-2})
1196	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwnet !! Long-wave net radiation (W m^{-2})
1197	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature (W m^{-2})
1198	REAL(r_std), DIMENSION (kjpindex),INTENT(out) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
1199
1200
1201	!! 0.3 Modified variables
1202	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: netrad !! Net radiation (W m^{-2})
1203	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation (mm day^{-1})
1204	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction
1205	!! (?? mm day^{-1})
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	!- Penser a sortir evapot en meme temps qu'evapot_corr tdo.
1394
1395	warning_correction(:)=.FALSE.
1396	DO ji=1,kjpindex
1397
1398	!! \latexonly
1399	!! \input{enerbilflux16.tex}
1400	!! \endlatexonly
1401	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
1402
1403	!! \latexonly
1404	!! \input{enerbilflux17.tex}
1405	!! \endlatexonly
1406	qc = speed * q_cdrag(ji)
1407
1408	!! Derive the potential as defined by Penman & Monteith (Monteith, 1965) based on the correction
1409	!! term developed by Chris Milly (1992).
1410	IF ((evapot(ji) .GT. zero) .AND. ((psnew(ji) - epot_air(ji)) .NE. zero )) THEN
1411	!
1412	!! \latexonly
1413	!! \input{enerbilflux18.tex}
1414	!! \endlatexonly
1415	correction = (quatre * emis(ji) * c_stefan * tair(ji)*3 + rau(ji) qc * cp_air + &
1416	& chalev0 * rau(ji) * qc * grad_qsat(ji) * vevapp(ji) / evapot(ji) )
1417
1418	!! \latexonly
1419	!! \input{enerbilflux19.tex}
1420	!! \endlatexonly
1421	IF (ABS(correction) .GT. min_sechiba) THEN
1422	correction = chalev0 * rau(ji) * qc * grad_qsat(ji) * (un - vevapp(ji)/evapot(ji)) / correction
1423	ELSE
1424	warning_correction(ji)=.TRUE.
1425	ENDIF
1426	ELSE
1427	correction = zero
1428	ENDIF
1429	correction = MAX (zero, correction)
1430
1431	!! \latexonly
1432	!! \input{enerbilflux20.tex}
1433	!! \endlatexonly
1434	evapot_corr(ji) = evapot(ji) / (un + correction)
1435
1436	ENDDO
1437	IF ( ANY(warning_correction) ) THEN
1438	DO ji=1,kjpindex
1439	IF ( warning_correction(ji) ) THEN
1440	WRITE(numout,*) ji,"Denominateur de la correction de milly nul! Aucune correction appliquee"
1441	ENDIF
1442	ENDDO
1443	ENDIF
1444	IF (printlev>=3) WRITE (numout,*) ' enerbil_flux done '
1445
1446	END SUBROUTINE enerbil_flux
1447
1448
1449	!! ===========================================================================================================================
1450	!! SUBROUTINE : enerbil_evapveg
1451	!!
1452	!>\BRIEF Splits total evaporation loss into individual process
1453	!! components and calculates the assimilation.
1454	!!
1455	!! DESCRIPTION : Based on the estimation of the fluxes in enerbil_flux,
1456	!! this routine splits the total evaporation into transpiration interception loss from vegetation,
1457	!! bare soil evaporation and snow sublimation. It then calculates the assimilation.
1458	!!
1459	!! MAIN OUTPUT VARIABLE(S) : vevapsno, vevapnu, transpir, vevapwet
1460	!!
1461	!! REFERENCE(S) :
1462	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
1463	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
1464	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1465	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1466	!! Circulation Model. Journal of Climate, 6, pp.248-273
1467	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
1468	!! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11):
1469	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
1470	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1471	!! sur le cycle de l'eau, PhD Thesis, available (22/12/11):
1472	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1473	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
1474	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
1475	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
1476	!! general circulation models. Global and Planetary Change, 19, pp.261-276
1477	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
1478	!! Interscience Publishers
1479	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
1480	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
1481	!!
1482	!! FLOWCHART : None
1483	!! \n
1484	!_ ==============================================================================================================================
1485
1486	SUBROUTINE enerbil_evapveg (kjpindex, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
1487	& rau, u, v, q_cdrag, qair, humrel, vevapsno, vevapnu , vevapflo, vevapwet, &
1488	& transpir, transpot, evapot, &
1489	& snowdz,snowrho,evapot_corr)
1490
1491	!! 0. Variable and parameter declaration
1492
1493	!! 0.1 Input variables
1494
1495	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1496	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-)
1497	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta4 !! Bare soil resistance (-)
1498	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance
1499	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density (kg m^{-3})
1500	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind velocity in directions u and v (m s^{-1})
1501	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag (m s^{-1})
1502	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1})
1503	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm day^{-1})
1504	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1)
1505	!!$ DS 15022011 humrel was used in a previous version of Orchidee, developped by Nathalie. Need to be discussed if it should be introduce again
1506	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Interception resistance (-)
1507	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Vegetation resistance (-)
1508	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3pot !! Vegetation resistance for potential transpiration
1509	REAL(r_std), DIMENSION(kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density
1510
1511	!! 0.2 Output variables
1512
1513	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation (mm day^{-1})
1514	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation (mm day^{-1})
1515	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapflo !! Floodplains evaporation
1516	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration (mm day^{-1})
1517	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpot !! Potential Transpiration
1518	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception (mm day^{-1})
1519
1520	!! 0.3 Modified variables
1521	REAL(r_std), DIMENSION(kjpindex,nsnow),INTENT(inout):: snowdz !! Snow depth
1522	REAL(r_std), DIMENSION(kjpindex),INTENT(inout) :: evapot_corr !! Correction term in evapotranspiration
1523	!! 0.4 Local variables
1524
1525	INTEGER(i_std) :: ji, jv
1526	REAL(r_std), DIMENSION(kjpindex) :: xx
1527	REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum, vbeta3sum
1528	REAL(r_std) :: speed !! Speed (m s^{-1})
1529	REAL(r_std), DIMENSION(kjpindex) :: ZSNOWEVAPS,ZEVAPCOR,ZSNOWDZ
1530
1531	!_ ==============================================================================================================================
1532	! initialisation: utile pour calculer l'evaporation des floodplains dans lesquelles il y a de la vegetation
1533	vbeta2sum(:) = 0.
1534	vbeta3sum(:) = 0.
1535	DO jv = 1, nvm
1536	vbeta2sum(:) = vbeta2sum(:) + vbeta2(:,jv)
1537	vbeta3sum(:) = vbeta3sum(:) + vbeta3(:,jv)
1538	ENDDO
1539
1540	!! 1. Compute vevapsno (evaporation from snow) and vevapnu (bare soil evaporation)
1541
1542	DO ji=1,kjpindex
1543
1544	!! \latexonly
1545	!! \input{enerbilevapveg1.tex}
1546	!! \endlatexonly
1547	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
1548
1549
1550	!! 1.1 Snow sublimation
1551	!! \latexonly
1552	!! \input{enerbilevapveg2.tex}
1553	!! \endlatexonly
1554	vevapsno(ji) = (un - vbeta5(ji)) * vbeta1(ji) * dt_sechiba * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji))
1555
1556
1557	! For the case of ok_explicitsnow activated, remove mass from uppermost snow layer in
1558	! response to sublimation. Do not take into account evaporation. The evaporation on the
1559	! snow surface could be possible in this case since the new snow model allows liquid water to be existed.
1560	IF ( ok_explicitsnow ) THEN
1561	! Sublimation of snow ice:
1562	! Reduce snow layer thickness and total snow depth if sublimation.
1563	! Add to correction term if potential sublimation exceeds available snow cover.
1564
1565	ZSNOWEVAPS(ji) = vevapsno(ji)/(chalsu0*snowrho(ji,1))
1566	ZSNOWDZ(ji) = snowdz(ji,1) - ZSNOWEVAPS(ji)
1567
1568	ZEVAPCOR(ji) = MAX(0.0,-ZSNOWDZ(ji))*snowrho(ji,1)/dt_sechiba
1569
1570	! Update evaporation correction term:
1571	evapot_corr(ji) = evapot_corr(ji) + ZEVAPCOR(ji)
1572	ENDIF
1573
1574	!! 1.2 Bare soil evaporation
1575	!! \latexonly
1576	!! \input{enerbilevapveg3.tex}
1577	!! \endlatexonly
1578	vevapnu(ji) = (un - vbeta1(ji)) * (un-vbeta5(ji)) * vbeta4(ji) * dt_sechiba * rau(ji) * speed * q_cdrag(ji) &
1579	& * (qsol_sat_new(ji) - qair(ji))
1580	!
1581	! 1.3 floodplains evaporation - transpiration et interception prioritaires dans les floodplains
1582	!
1583	vevapflo(ji) = vbeta5(ji) * (1 - vbeta2sum(ji) - vbeta3sum(ji)) &
1584	& * dt_sechiba * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji))
1585
1586	END DO
1587
1588	!! 2. Compute transpir (transpiration) and vevapwet (interception)
1589
1590	!! Preliminaries
1591	DO ji = 1, kjpindex
1592
1593	!! \latexonly
1594	!! \input{enerbilevapveg4.tex}
1595	!! \endlatexonly
1596	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
1597
1598	!! \latexonly
1599	!! \input{enerbilevapveg5.tex}
1600	!! \endlatexonly
1601	xx(ji) = dt_sechiba * (un-vbeta1(ji)) * (qsol_sat_new(ji)-qair(ji)) * rau(ji) * speed * q_cdrag(ji)
1602
1603	ENDDO
1604
1605	DO jv=1,nvm
1606	DO ji=1,kjpindex
1607
1608	!! 2.1 Calculate interception loss
1609	!! \latexonly
1610	!! \input{enerbilevapveg6.tex}
1611	!! \endlatexonly
1612	vevapwet(ji,jv) = xx(ji) * vbeta2(ji,jv)
1613	!
1614	!! 2.2 Calculate transpiration
1615	!! \latexonly
1616	!! \input{enerbilevapveg7.tex}
1617	!! \endlatexonly
1618	transpir (ji,jv) = xx(ji) * vbeta3(ji,jv)
1619
1620	transpot(ji,jv) = xx(ji) * vbeta3pot(ji,jv)
1621
1622	END DO
1623	END DO
1624
1625
1626
1627	IF (printlev>=3) WRITE (numout,*) ' enerbil_evapveg done '
1628
1629	END SUBROUTINE enerbil_evapveg
1630
1631
1632	!! =============================================================================================================================
1633	!! SUBROUTINE : enerbil_fusion
1634	!!
1635	!>\BRIEF Computes new soil temperature due to
1636	!! ice and snow melt.
1637	!!
1638	!! DESCRIPTION : This routine computes new soil temperature due to
1639	!! ice and snow melt. It is the second part of main routine for enerbil module, and is called every
1640	!! time step.
1641	!!
1642	!! MAIN OUTPUT VARIABLE(S) : temp_sol_new, fusion
1643	!!
1644	!! REFERENCE(S) :
1645	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
1646	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
1647	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1648	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1649	!! Circulation Model. Journal of Climate, 6, pp.248-273
1650	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
1651	!! interface avec la surface planetaire dans LMDZ. Technical note, available (22/12/11):
1652	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
1653	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1654	!! sur le cycle de l'eau. PhD Thesis, available (22/12/11):
1655	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1656	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
1657	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
1658	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
1659	!! general circulation models. Global and Planetary Change, 19, pp.261-276
1660	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
1661	!! Interscience Publishers
1662	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
1663	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
1664	!!
1665	!! FLOWCHART : None
1666	!! \n
1667	!_ ==============================================================================================================================
1668
1669	SUBROUTINE enerbil_fusion (kjpindex, tot_melt, soilcap, snowdz, &
1670	temp_sol_new, fusion)
1671
1672
1673	!! 0. Variable and parameter declaration
1674
1675	!! 0.1 Input variables
1676
1677	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1678	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tot_melt !! Total melt (??)
1679	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil heat capacity (J K^{-1])
1680	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow layer thickness (m)
1681
1682	!! 0.2 Output variables
1683
1684	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fusion !! Fusion (??)
1685
1686	!! 0.3 Modified variables
1687
1688	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! New soil temperature (K)
1689
1690	!! 0.4 Local variables
1691	INTEGER(i_std) :: ji !! (-)
1692
1693	!_ ==============================================================================================================================
1694
1695	!! Initialisation
1696	IF (printlev>=3) WRITE (numout,*) ' enerbil_fusion start ', MINVAL(soilcap), MINLOC(soilcap),&
1697	& MAXVAL(soilcap), MAXLOC(soilcap)
1698
1699	!! 1. Calculate new soil temperature due to ice and snow melt
1700
1701	IF (.NOT. ok_explicitsnow) THEN
1702	!! Surface temperature is reduced if there is snow melt
1703
1704	DO ji=1,kjpindex
1705	!! \latexonly
1706	!! \input{enerbilfusion1.tex}
1707	!! \endlatexonly
1708	fusion(ji) = tot_melt(ji) * chalfu0 / dt_sechiba
1709
1710	!! \latexonly
1711	!! \input{enerbilfusion2.tex}
1712	!! \endlatexonly
1713	temp_sol_new(ji) = temp_sol_new(ji) - ((tot_melt(ji) * chalfu0) / soilcap(ji))
1714	END DO
1715	ENDIF
1716
1717	IF (printlev>=3) WRITE (numout,*) ' enerbil_fusion done '
1718
1719	END SUBROUTINE enerbil_fusion
1720
1721
1722	!! =============================================================================================================================
1723	!! SUBROUTINE : enerbil_t2mdiag
1724	!!
1725	!> BRIEF This routine diagnoses the air temperature at 2
1726	!! metres above ground.
1727	!!
1728	!! DESCRIPTION : The 2 metre above ground air temperature is currently
1729	!! set to be the same as the general air temperaturem, as the surface model currently simulates the
1730	!! above ground atmosphere as a uniform box (and the lowest layer of the LMDZ atmospheric model).
1731	!!
1732	!! MAIN OUTPUT VARIABLE(S) : t2mdiag
1733	!!
1734	!! REFERENCE(S) :
1735	!! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled
1736	!! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278
1737	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1738	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1739	!! Circulation Model. Journal of Climate, 6, pp.248-273
1740	!! - Dufresne, J-L & Ghattas, J, 2009. Description du schÃ©ma de la couche limite turbulente et la
1741	!! interface avec la surface planetaire dans LMDZ, Technical note, available from:
1742	!! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf
1743	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1744	!! sur le cycle de l'eau, PhD Thesis, available from:
1745	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1746	!! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J
1747	!! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D &
1748	!! Xue, Y, 1998. A proposal for a general interface between land surface schemes and
1749	!! general circulation models. Global and Planetary Change, 19, pp.261-276
1750	!! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems.
1751	!! Interscience Publishers\n
1752	!! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere
1753	!! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663.
1754	!!
1755	!! FLOWCHART : None
1756	!! \n
1757	!_ ==============================================================================================================================
1758
1759	SUBROUTINE enerbil_t2mdiag (kjpindex, temp_air, t2mdiag)
1760
1761	!! 0. Variable and parameter declaration
1762
1763	!! 0.1 Input variables
1764
1765	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1766	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
1767
1768	!! 0.2 Output variables
1769	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: t2mdiag !! 2-metre temperature (K)
1770
1771	!! 0.3 Modified variables
1772
1773	!! 0.4 Local variables
1774
1775	!_ ==============================================================================================================================
1776
1777	!! The 2 metre air temperature is assumed to be the same of the value of temp_air - the atmospheric air
1778	!! temperature, which also taken as the temperature at the lowest level of the LMDZ. 't2mdiag' is just
1779	!! a value used for diagnostic purposes, admd it is no applied amywhere else.
1780	!! \latexonly
1781	!! \input{enerbilt2mdiag1.tex}
1782	!! \endlatexonly
1783	t2mdiag(:) = temp_air(:)
1784
1785	IF (printlev>=3) WRITE (numout,*) ' enerbil_t2mdiag done '
1786
1787	END SUBROUTINE enerbil_t2mdiag
1788
1789
1790
1791	END MODULE enerbil

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

Context Navigation

source: branches/publications/ORCHIDEE-LEAK-r5919/src_sechiba/enerbil.f90 @ 5925

Download in other formats: