source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_sechiba/enerbil.f90 @ 7541

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