Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

enerbil.f90 @ 7746

Last change on this file since 7746 was 7366, checked in by simon.bowring, 3 years ago
Biochar version
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 101.9 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_Biochar/src_sechiba/enerbil.f90 @ 7746

Download in other formats: