source: tags/ORCHIDEE_1_9_5/ORCHIDEE/src_sechiba/enerbil.f90 @ 8

Last change on this file since 8 was 8, checked in by orchidee, 14 years ago

import first tag equivalent to CVS orchidee_1_9_5 + OOL_1_9_5
File size: 46.3 KB
Line 
1!!
2!! This module computes energy bilan on continental surface
3!!
4!! @author Marie-Alice Foujols and Jan Polcher
5!! @Version : $Revision: 1.24 $, $Date: 2009/01/07 13:39:45 $
6!!
7!! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/enerbil.f90,v 1.24 2009/01/07 13:39:45 ssipsl Exp $
8!! IPSL (2006)
9!!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
10!!
11MODULE enerbil
12
13  ! routines called : restput, restget
14  !
15  USE ioipsl
16  !
17  ! modules used :
18  USE constantes
19  USE constantes_veg
20  USE sechiba_io
21  USE parallel
22!  USE write_field_p, only : WriteFieldI_p 
23  IMPLICIT NONE
24
25  ! public routines :
26  ! enerbil_main
27  ! enerbil_fusion
28  PRIVATE
29  PUBLIC :: enerbil_main, enerbil_fusion,enerbil_clear
30
31  !
32  ! variables used inside enerbil module : declaration and initialisation
33  !
34  LOGICAL, SAVE                              :: l_first_enerbil=.TRUE.  !! Initialisation has to be done one time
35
36  CHARACTER(LEN=80), SAVE                    :: var_name                !! To store variables names for I/O
37
38  ! one dimension array allocated, computed and used in enerbil module exclusively
39
40  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: psold                   !! Old surface dry static energy
41  !! saturated specific humudity for old temperature
42  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: qsol_sat
43  !! derivative of satured specific humidity at the old temperature
44  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: pdqsold
45  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: psnew                   !! New surface static energy
46  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: qsol_sat_new            !! New saturated surface air moisture
47  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: netrad                  !! Net radiation
48  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: lwabs                   !! LW radiation absorbed by the surface
49  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: lwup                    !! Long-wave up radiation
50  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: lwnet                   !! Net Long-wave radiation
51  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: fluxsubli               !! Energy of sublimation
52  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: qsat_air                !!
53  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tair                    !!
54CONTAINS
55
56  !!
57  !! Main routine for *enerbil* module
58  !! - called only one time for initialisation
59  !! - called every time step
60  !! - called one more time at last time step for writing _restart_ file
61  !!
62  !! Algorithm:
63  !! - call enerbil_begin for initialisation
64  !! - call enerbil_surftemp for psnew and qsol_sat_new
65  !! - call enerbil_flux for tsol_new, netrad, vevapp, fluxlat and fluxsens
66  !! - call enerbil_evapveg for evaporation and transpiration
67  !!
68  !! @call enerbil_begin
69  !! @call enerbil_surftemp
70  !! @call enerbil_flux
71  !! @call enerbil_evapveg
72  !!
73  SUBROUTINE enerbil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, &
74     & index, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, &
75     & qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, &
76     & cimean, ccanopy, emis, soilflx, soilcap, q_cdrag, humrel, fluxsens, fluxlat, &
77     & vevapp, transpir, gpp, vevapnu, vevapwet, vevapsno, t2mdiag, temp_sol, tsol_rad, &
78     & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id ) 
79
80    ! interface description
81    ! input scalar
82    INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number
83    INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size
84    INTEGER(i_std),INTENT (in)                         :: rest_id,hist_id  !! _Restart_ file and _history_ file identifier
85    INTEGER(i_std),INTENT (in)                         :: hist2_id         !! _history_ file 2 identifier
86    REAL(r_std), INTENT (in)                           :: dtradia          !! Time step in seconds
87    LOGICAL, INTENT(in)                               :: ldrestart_read   !! Logical for _restart_ file to read
88    LOGICAL, INTENT(in)                               :: ldrestart_write  !! Logical for _restart_ file to write
89    ! input fields
90    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)   :: index            !! Indeces of the points on the map
91    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: zlev             !! Height of first layer
92    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: lwdown           !! Down-welling long-wave flux
93    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: swnet            !! Net surface short-wave flux
94    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: epot_air         !! Air potential energy
95    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: temp_air         !! Air temperature in Kelvin
96    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: u                !! zonal wind (m/s)
97    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: v                !! north-south wind (m/s)
98    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: petAcoef         !! PetAcoef
99    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: petBcoef         !! PetBcoef
100    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: qair             !! Lowest level specific humidity
101    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: peqAcoef         !! PeqAcoef
102    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: peqBcoef         !! PeqBcoef
103    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: pb               !! Lowest level pressure
104    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: rau              !! Density
105    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: vbeta            !! Resistance coefficient
106    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: valpha           !! Resistance coefficient
107    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: vbeta1           !! Snow resistance
108    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: vbeta4           !! Bare soil resistance
109    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: emis             !! Emissivity
110    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: soilflx          !! Soil flux
111    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: soilcap          !! Soil calorific capacity
112    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: q_cdrag          !! This is the cdrag without the wind multiplied
113    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ccanopy          !! CO2 concentration in the canopy
114    REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel           !! Relative humidity
115    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: vbeta2           !! Interception resistance
116    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: vbeta3           !! Vegetation resistance
117    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: vbetaco2         !! Vegetation resistance to CO2
118    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: cimean           !! mean Ci
119    ! modified fields
120    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: evapot           !! Soil Potential Evaporation
121    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: evapot_corr !! Soil Potential Evaporation Correction
122    ! output fields
123    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: fluxsens         !! Sensible chaleur flux
124    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: fluxlat          !! Latent chaleur flux
125    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: vevapp           !! Total of evaporation
126    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: vevapnu          !! Bare soil evaporation
127    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: vevapsno         !! Snow evaporation
128    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: tsol_rad         !! Tsol_rad
129    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: temp_sol_new     !! New soil temperature
130    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: temp_sol         !! Soil temperature
131    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: qsurf            !! Surface specific humidity
132    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir         !! Transpiration
133    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp              !! Assimilation, gC/m**2 total area.
134    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet         !! Interception
135    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: t2mdiag          !! 2-meter temperature
136    !
137    ! LOCAL
138    !
139    REAL(r_std),DIMENSION (kjpindex) :: epot_air_new, qair_new
140    !
141    ! do initialisation
142    !
143    IF (l_first_enerbil) THEN
144
145        IF (long_print) WRITE (numout,*) ' l_first_enerbil : call enerbil_init '
146        CALL enerbil_init (kjit, ldrestart_read, kjpindex, index, rest_id, qair, temp_sol, temp_sol_new,  &
147             &             qsurf, tsol_rad, vevapp, fluxsens, fluxlat, gpp, evapot, evapot_corr)
148
149        CALL enerbil_var_init (kjpindex, temp_air, t2mdiag)
150
151        RETURN
152
153    ENDIF
154
155    !
156    ! prepares restart file for the next simulation
157    !
158    IF (ldrestart_write) THEN
159
160        IF (long_print) WRITE (numout,*) ' we have to complete restart file with ENERBIL variables '
161
162        var_name= 'temp_sol' 
163        CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit,  temp_sol, 'scatter',  nbp_glo, index_g)
164
165        var_name= 'qsurf' 
166        CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit,  qsurf, 'scatter',  nbp_glo, index_g)
167
168        var_name= 'evapot' 
169        CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit,  evapot, 'scatter',  nbp_glo, index_g)
170
171        var_name= 'evapot_corr' 
172        CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit,  evapot_corr, 'scatter',  nbp_glo, index_g)
173
174        var_name= 'tsolrad'
175        CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit,  tsol_rad, 'scatter',  nbp_glo, index_g)
176
177        var_name= 'evapora'
178        CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit,  vevapp, 'scatter',  nbp_glo, index_g)
179
180        var_name= 'fluxlat'
181        CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit,  fluxlat, 'scatter',  nbp_glo, index_g)
182
183        var_name= 'fluxsens'
184        CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit,  fluxsens, 'scatter',  nbp_glo, index_g)
185
186        IF ( control%stomate_watchout .OR. control%ok_co2 ) THEN
187
188          ! The gpp could in principle be recalculated at the beginning of the run.
189          ! However, we would need several variables that are not stored in the restart files.
190          !
191          var_name= 'gpp'
192          CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, gpp, 'scatter',  nbp_glo, index_g)
193
194        ENDIF
195
196        RETURN
197
198    END IF
199
200    !
201    !
202    ! 1. computes some initialisation: psold, qsol_sat and pdqsold
203    !
204    CALL enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis)
205
206    !
207    ! 2. computes psnew, qsol_sat_new, temp_sol_new, qair_new and epot_air_new
208    !
209    CALL enerbil_surftemp (kjpindex, dtradia, zlev, emis, &
210       & epot_air, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,&
211       & valpha, vbeta1, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, &
212       & qair_new, epot_air_new)
213
214    !
215    ! 3. computes lwup, lwnet, tsol_rad, netrad, qsurf, vevapp, evapot, evapot_corr, fluxlat, fluxsubli and fluxsens
216    !
217
218    CALL enerbil_flux (kjpindex, dtradia, emis, temp_sol, rau, u, v, q_cdrag, vbeta, valpha, vbeta1, &
219       & qair_new, epot_air_new, psnew, qsurf, &
220       & fluxsens , fluxlat , fluxsubli, vevapp, temp_sol_new, lwdown, swnet, &
221       & lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr)
222
223    !
224    ! 4. computes in details evaporation and transpiration
225    !
226
227    CALL enerbil_evapveg (kjpindex, dtradia, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, cimean, &
228       & ccanopy, rau, u, v, q_cdrag, qair_new, humrel, vevapsno, vevapnu , vevapwet, transpir, gpp, evapot)
229
230    !
231    ! 5. diagnose 2-meter temperatures
232    !
233
234    CALL enerbil_t2mdiag (kjpindex, temp_air, t2mdiag)
235
236    IF ( .NOT. almaoutput ) THEN
237       CALL histwrite(hist_id, 'netrad', kjit, netrad, kjpindex, index)
238       CALL histwrite(hist_id, 'evapot', kjit, evapot, kjpindex, index)
239       CALL histwrite(hist_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index)
240       CALL histwrite(hist_id, 'lwdown', kjit, lwabs,  kjpindex, index)
241       CALL histwrite(hist_id, 'lwnet',  kjit, lwnet,  kjpindex, index)
242       IF ( hist2_id > 0 ) THEN
243          CALL histwrite(hist2_id, 'netrad', kjit, netrad, kjpindex, index)
244          CALL histwrite(hist2_id, 'evapot', kjit, evapot, kjpindex, index)
245          CALL histwrite(hist2_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index)
246          CALL histwrite(hist2_id, 'lwdown', kjit, lwabs,  kjpindex, index)
247          CALL histwrite(hist2_id, 'lwnet',  kjit, lwnet,  kjpindex, index)
248       ENDIF
249    ELSE
250       CALL histwrite(hist_id, 'LWnet', kjit, lwnet, kjpindex, index)
251       CALL histwrite(hist_id, 'Qv', kjit, fluxsubli, kjpindex, index)
252       CALL histwrite(hist_id, 'PotEvap', kjit, evapot_corr, kjpindex, index)
253       IF ( hist2_id > 0 ) THEN
254          CALL histwrite(hist2_id, 'LWnet', kjit, lwnet, kjpindex, index)
255          CALL histwrite(hist2_id, 'Qv', kjit, fluxsubli, kjpindex, index)
256          CALL histwrite(hist2_id, 'PotEvap', kjit, evapot_corr, kjpindex, index)
257       ENDIF
258    ENDIF
259
260    IF (long_print) WRITE (numout,*) ' enerbil_main Done '
261
262  END SUBROUTINE enerbil_main
263
264  !! Algorithm:
265  !! - dynamic allocation for local array
266  !!
267  SUBROUTINE enerbil_init (kjit, ldrestart_read, kjpindex, index, rest_id, qair, temp_sol, temp_sol_new, &
268       &                   qsurf, tsol_rad, vevapp, fluxsens, fluxlat, gpp, evapot, evapot_corr)
269
270    ! interface description
271    ! input scalar
272    INTEGER(i_std), INTENT (in)                              :: kjit               !! Time step number
273    LOGICAL ,INTENT (in)                                    :: ldrestart_read     !! Logical for _restart_ file to read
274    INTEGER(i_std), INTENT (in)                              :: kjpindex           !! Domain size
275    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index              !! Indeces of the points on the map
276    INTEGER(i_std), INTENT (in)                              :: rest_id            !! _Restart_ file identifier
277    ! input fields
278    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair               !! Lowest level specific humidity
279    ! output scalar
280    ! output fields, they need to initialized somehow for the model forcing ORCHIDEE.
281    !
282    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: temp_sol           !! Soil temperature
283    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: temp_sol_new       !! New soil temperature
284    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsurf              !! near surface specific humidity
285    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: tsol_rad           !! Tsol_rad
286    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapp           !! Total of evaporation
287    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsens         !! Sensible chaleur flux
288    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat          !! Latent chaleur flux
289    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: gpp                !! Assimilation, gC/m**2 total area.
290    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: evapot             !! Soil Potential Evaporation
291    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: evapot_corr        !! Soil Potential Evaporation Correction
292
293    ! local declaration
294    INTEGER(i_std)                                          :: ier
295
296    ! initialisation
297    IF (l_first_enerbil) THEN
298        l_first_enerbil=.FALSE.
299    ELSE
300        WRITE (numout,*) ' l_first_enerbil false . we stop '
301        STOP 'enerbil_init'
302    ENDIF
303
304    ALLOCATE (psold(kjpindex),stat=ier)
305    IF (ier.NE.0) THEN
306        WRITE (numout,*) ' error in psold allocation. We stop.We need kjpindex words = ',kjpindex
307        STOP 'enerbil_init'
308    END IF
309
310    ALLOCATE (qsol_sat(kjpindex),stat=ier)
311    IF (ier.NE.0) THEN
312        WRITE (numout,*) ' error in qsol_sat allocation. We stop. We need kjpindex words = ',kjpindex
313        STOP 'enerbil_init'
314    END IF
315
316    ALLOCATE (pdqsold(kjpindex),stat=ier)
317    IF (ier.NE.0) THEN
318        WRITE (numout,*) ' error in pdqsold allocation. We stop. We need kjpindex words = ',kjpindex
319        STOP 'enerbil_init'
320    END IF
321
322    ALLOCATE (psnew(kjpindex),stat=ier)
323    IF (ier.NE.0) THEN
324        WRITE (numout,*) ' error in psnew allocation. We stop. We need kjpindex words = ',kjpindex
325        STOP 'enerbil_init'
326    END IF
327
328    ALLOCATE (qsol_sat_new(kjpindex),stat=ier)
329    IF (ier.NE.0) THEN
330        WRITE (numout,*) ' error in qsol_sat_new allocation. We stop. We need kjpindex words = ',kjpindex
331        STOP 'enerbil_init'
332    END IF
333
334    ALLOCATE (netrad(kjpindex),stat=ier)
335    IF (ier.NE.0) THEN
336        WRITE (numout,*) ' error in netrad allocation. We stop. We need kjpindex words = ',kjpindex
337        STOP 'enerbil_init'
338    END IF
339
340    ALLOCATE (lwabs(kjpindex),stat=ier)
341    IF (ier.NE.0) THEN
342        WRITE (numout,*) ' error in lwabs allocation. We stop. We need kjpindex words = ',kjpindex
343        STOP 'enerbil_init'
344    END IF
345
346    ALLOCATE (lwup(kjpindex),stat=ier)
347    IF (ier.NE.0) THEN
348        WRITE (numout,*) ' error in lwup allocation. We stop. We need kjpindex words = ',kjpindex
349        STOP 'enerbil_init'
350    END IF
351
352    ALLOCATE (lwnet(kjpindex),stat=ier)
353    IF (ier.NE.0) THEN
354        WRITE (numout,*) ' error in lwnet allocation. We stop. We need kjpindex words = ',kjpindex
355        STOP 'enerbil_init'
356    END IF
357
358    ALLOCATE (fluxsubli(kjpindex),stat=ier)
359    IF (ier.NE.0) THEN
360        WRITE (numout,*) ' error in fluxsubli allocation. We stop. We need kjpindex words = ',kjpindex
361        STOP 'enerbil_init'
362    END IF
363
364    ALLOCATE (qsat_air(kjpindex),stat=ier)
365    IF (ier.NE.0) THEN
366        WRITE (numout,*) ' error in qsat_air allocation. We stop. We need kjpindex words = ',kjpindex
367        STOP 'enerbil_init'
368    END IF
369
370    ALLOCATE (tair(kjpindex),stat=ier)
371    IF (ier.NE.0) THEN
372        WRITE (numout,*) ' error in tair allocation. We stop. We need kjpindex words = ',kjpindex
373        STOP 'enerbil_init'
374    END IF
375
376    ! open restart input file done by enerbil_init
377    ! and read data from restart input file for ENERBIL process
378
379    IF (ldrestart_read) THEN
380
381        IF (long_print) WRITE (numout,*) ' we have to read a restart file for HYDROLOGIC variables'
382
383        var_name='temp_sol'
384        CALL ioconf_setatt('UNITS', 'K')
385        CALL ioconf_setatt('LONG_NAME','Surface temperature')
386        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp_sol, "gather", nbp_glo, index_g)
387        !
388        !Config Key  = ENERBIL_TSURF
389        !Config Desc = Initial temperature if not found in restart
390        !Config Def  = 280.
391        !Config Help = The initial value of surface temperature if its value is not found
392        !Config        in the restart file. This should only be used if the model is
393        !Config        started without a restart file.
394        !
395        CALL setvar_p (temp_sol, val_exp,'ENERBIL_TSURF', 280._r_std)
396        !
397        var_name= 'qsurf'
398        CALL ioconf_setatt('UNITS', 'g/g')
399        CALL ioconf_setatt('LONG_NAME','near surface specific humidity')
400        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., qsurf, "gather", nbp_glo, index_g)
401        IF ( ALL( qsurf(:) .EQ. val_exp ) ) THEN
402           qsurf(:) = qair(:)
403        ENDIF
404        !
405        var_name= 'evapot'
406        CALL ioconf_setatt('UNITS', 'mm/d')
407        CALL ioconf_setatt('LONG_NAME','Soil Potential Evaporation')
408        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., evapot, "gather", nbp_glo, index_g)
409        !
410        var_name= 'evapot_corr'
411        CALL ioconf_setatt('UNITS', 'mm/d')
412        CALL ioconf_setatt('LONG_NAME','Corrected Soil Potential Evaporation')
413        IF ( ok_var(var_name) ) THEN
414           CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., evapot_corr, "gather", nbp_glo, index_g)
415        ENDIF
416        !
417        !Config Key  = ENERBIL_EVAPOT
418        !Config Desc = Initial Soil Potential Evaporation
419        !Config Def  = 0.0
420        !Config Help = The initial value of soil potential evaporation if its value
421        !Config        is not found in the restart file. This should only be used if
422        !Config        the model is started without a restart file.
423        !
424        CALL setvar_p (evapot, val_exp, 'ENERBIL_EVAPOT', 0.0_r_std)
425        IF ( ok_var("evapot_corr") ) THEN
426           CALL setvar_p (evapot_corr, val_exp, 'ENERBIL_EVAPOT', 0.0_r_std)
427        ENDIF
428        !
429        var_name= 'tsolrad'
430        CALL ioconf_setatt('UNITS', 'K')
431        CALL ioconf_setatt('LONG_NAME','Radiative surface temperature')
432        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., tsol_rad, "gather", nbp_glo, index_g)
433        IF ( ALL( tsol_rad(:) .EQ. val_exp ) ) THEN
434           tsol_rad(:) = temp_sol(:)
435        ENDIF
436        !
437        ! Set the fluxes so that we have something reasonable and not NaN on some machines
438        !
439        var_name= 'evapora'
440        CALL ioconf_setatt('UNITS', 'Kg/m^2/dt')
441        CALL ioconf_setatt('LONG_NAME','Evaporation')
442        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vevapp, "gather", nbp_glo, index_g)
443        IF ( ALL( vevapp(:) .EQ. val_exp ) ) THEN
444           vevapp(:) = 0.0_r_std
445        ENDIF
446        !
447        var_name= 'fluxlat'
448        CALL ioconf_setatt('UNITS', 'W/m^2')
449        CALL ioconf_setatt('LONG_NAME','Latent heat flux')
450        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., fluxlat, "gather", nbp_glo, index_g)
451        IF ( ALL( fluxlat(:) .EQ. val_exp ) ) THEN
452           fluxlat(:) = 0.0_r_std
453        ENDIF
454        !
455        var_name= 'fluxsens'
456        CALL ioconf_setatt('UNITS', 'W/m^2')
457        CALL ioconf_setatt('LONG_NAME','Sensible heat flux')
458        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., fluxsens, "gather", nbp_glo, index_g)
459        IF ( ALL( fluxsens(:) .EQ. val_exp ) ) THEN
460           fluxsens(:) = 0.0_r_std
461        ENDIF
462        !
463        ! If we are with STOMATE
464        !
465        IF ( control%stomate_watchout .OR. control%ok_co2 ) THEN
466
467          ! The gpp could in principle be recalculated at the beginning of the run.
468          ! However, we would need several variables that are not stored in the restart files.
469          var_name= 'gpp'
470          CALL ioconf_setatt('UNITS', 'gC/m**2/time step')
471          CALL ioconf_setatt('LONG_NAME','Gross primary productivity')
472          CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., gpp, "gather", nbp_glo, index_g)
473
474          IF ( ALL( gpp(:,:) .EQ. val_exp ) ) THEN
475             gpp(:,:) = zero
476          ENDIF
477
478        ENDIF
479
480        ! initialises temp_sol_new
481        ! saved in thermosoil
482        var_name='temp_sol_new'
483        CALL ioconf_setatt('UNITS', 'K')
484        CALL ioconf_setatt('LONG_NAME','New Surface temperature')
485        IF ( ok_var(var_name) ) THEN
486           CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp_sol_new, "gather", nbp_glo, index_g)
487           IF ( ALL( temp_sol_new(:) .EQ. val_exp ) ) THEN
488              temp_sol_new(:) = temp_sol(:)
489           ENDIF
490        ELSE
491           ! initialises temp_sol_new
492           temp_sol_new(:) = temp_sol(:)
493        ENDIF
494
495    ELSE
496       ! initialises temp_sol_new
497       temp_sol_new(:) = temp_sol(:)
498    ENDIF
499
500    IF (long_print) WRITE (numout,*) ' enerbil_init done '
501
502  END SUBROUTINE enerbil_init
503
504  SUBROUTINE enerbil_clear ()
505    l_first_enerbil=.TRUE.
506    IF ( ALLOCATED (psold)) DEALLOCATE (psold)
507    IF ( ALLOCATED (qsol_sat)) DEALLOCATE (qsol_sat)
508    IF ( ALLOCATED (pdqsold)) DEALLOCATE (pdqsold)
509    IF ( ALLOCATED (psnew)) DEALLOCATE (psnew)
510    IF ( ALLOCATED (qsol_sat_new)) DEALLOCATE (qsol_sat_new)
511    IF ( ALLOCATED (netrad)) DEALLOCATE (netrad)
512    IF ( ALLOCATED (lwabs)) DEALLOCATE (lwabs)
513    IF ( ALLOCATED (lwup)) DEALLOCATE (lwup)
514    IF ( ALLOCATED (lwnet)) DEALLOCATE (lwnet)
515    IF ( ALLOCATED (fluxsubli)) DEALLOCATE (fluxsubli)
516    IF ( ALLOCATED (qsat_air)) DEALLOCATE (qsat_air)
517    IF ( ALLOCATED (tair)) DEALLOCATE (tair)
518   !
519   ! open restart input file done by enerbil_init
520    ! and read data from restart input file for ENERBIL process
521
522 
523  END SUBROUTINE enerbil_clear
524 
525  SUBROUTINE enerbil_var_init (kjpindex, temp_air, t2mdiag)
526
527    ! interface description
528    ! input scalar
529    INTEGER(i_std), INTENT(in)                         :: kjpindex       !! Domain size
530    ! input fields
531    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: temp_air       !! Air temperature in Kelvin
532    ! modified fields
533    ! output fields
534    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: t2mdiag        !! 2-meter temperature
535
536    CALL enerbil_t2mdiag (kjpindex, temp_air, t2mdiag)
537
538    IF (long_print) WRITE (numout,*) ' enerbil_var_init done '
539
540  END SUBROUTINE enerbil_var_init
541
542  !! This routines computes psold, qsol_sat, pdqsold and netrad
543  !!
544  SUBROUTINE enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis)
545
546    ! interface description
547    INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size
548    ! input fields
549    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: temp_sol         !! Soil temperature in Kelvin
550    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: lwdown           !! Down-welling long-wave flux
551    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: swnet            !! Net surface short-wave flux
552    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: pb               !! Lowest level pressure
553    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: emis             !! Emissivity
554    ! output fields
555    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: psold            !! Old surface dry static energy
556    !! Saturated specific humudity for old temperature
557    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: qsol_sat
558    !! Derivative of satured specific humidity at the old temperature
559    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: pdqsold
560    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: netrad           !! Net radiation
561
562    ! local declaration
563    INTEGER(i_std)                                :: ji
564    REAL(r_std), DIMENSION(kjpindex)               :: dev_qsol
565    REAL(r_std), PARAMETER                         :: missing = 999998.
566    ! initialisation
567
568    !
569    ! 1. computes psold
570    !
571    psold(:) = temp_sol(:)*cp_air
572    !
573    ! 2. computes qsol_sat
574    !
575    CALL qsatcalc (kjpindex, temp_sol, pb, qsol_sat)
576
577    IF ( diag_qsat ) THEN
578      IF ( ANY(ABS(qsol_sat(:)) .GT. missing) ) THEN
579        DO ji = 1, kjpindex
580          IF ( ABS(qsol_sat(ji)) .GT. missing) THEN
581            WRITE(numout,*) 'ERROR on ji = ', ji
582            WRITE(numout,*) 'temp_sol(ji),  pb(ji) :', temp_sol(ji),  pb(ji)
583            CALL ipslerr (3,'enerbil_begin', &
584 &           'qsol too high ','','')
585          ENDIF
586        ENDDO
587      ENDIF
588    ENDIF
589
590    !
591    ! 3. computes pdqsold
592    !
593    CALL dev_qsatcalc (kjpindex, temp_sol, pb, dev_qsol)
594    DO ji = 1, kjpindex
595      pdqsold(ji) = dev_qsol(ji) * ( pb(ji)**kappa ) / cp_air
596    ENDDO
597
598    IF ( diag_qsat ) THEN
599      IF ( ANY(ABS( pdqsold(:)) .GT. missing) ) THEN
600        DO ji = 1, kjpindex
601          IF ( ABS( pdqsold(ji)) .GT. missing ) THEN
602            WRITE(numout,*) 'ERROR on ji = ', ji
603            WRITE(numout,*) 'temp_sol(ji),  pb(ji) :', temp_sol(ji),  pb(ji)
604            CALL ipslerr (3,'enerbil_begin', &
605 &           'pdqsold too high ','','')
606          ENDIF
607        ENDDO
608      ENDIF
609    ENDIF
610
611    !
612    ! 4. computes netrad and absorbed LW radiation absorbed at the surface
613    !
614    lwabs(:) = emis(:) * lwdown(:)
615    netrad(:) = lwdown(:) + swnet (:) - (emis(:) * c_stefan * temp_sol(:)**4 + (un - emis(:)) * lwdown(:)) 
616    !
617    IF (long_print) WRITE (numout,*) ' enerbil_begin done '
618
619  END SUBROUTINE enerbil_begin
620
621  !! This routine computes psnew and qsol_sat_new
622  !!
623  !! Computes the energy balance at the surface with an implicit scheme
624  !! that is connected to the Richtmyer and Morton algorithm of the PBL.
625  !!
626  SUBROUTINE enerbil_surftemp (kjpindex, dtradia, zlev, emis, epot_air, &
627     & petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,&
628     & valpha, vbeta1, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, &
629     & qair_new, epot_air_new) 
630
631    ! interface
632    ! input scalar
633    INTEGER(i_std), INTENT(in)                               :: kjpindex      !! Domain size
634    REAL(r_std), INTENT(in)                                  :: dtradia       !! Time step in seconds
635    ! input fields
636    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: zlev          !! Height of first layer
637    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: emis          !! Emissivity
638    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: epot_air      !! Air potential energy
639    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petAcoef      !! PetAcoef
640    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: petBcoef      !! PetBcoef
641    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair          !! Lowest level specific humidity
642    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqAcoef      !! PeqAcoef
643    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: peqBcoef      !! PeqBcoef
644    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: soilflx       !! Soil flux
645    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: rau           !! Density
646    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u, v          !! Wind
647    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q_cdrag       !!
648    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: vbeta         !! Resistance coefficient
649    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: valpha        !! Resistance coefficient
650    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: vbeta1        !! Snow resistance
651    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: soilcap       !! Soil calorific capacity
652    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: lwdown        !! Down-welling long-wave flux
653    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swnet         !! Net surface short-wave flux
654    ! output fields
655    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: psnew         !! New surface static energy
656    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsol_sat_new  !! New saturated surface air moisture
657    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: temp_sol_new  !! New soil temperature
658    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qair_new      !! New air moisture
659    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: epot_air_new  !! New air temperature
660
661
662    ! local declaration
663    INTEGER(i_std)                  :: ji
664    REAL(r_std),DIMENSION (kjpindex) :: zicp
665    REAL(r_std)                      :: fevap
666    REAL(r_std)                      :: sensfl_old, larsub_old, lareva_old, dtheta, sum_old, sum_sns
667    REAL(r_std)                      :: zikt, zikq, netrad_sns, sensfl_sns, larsub_sns, lareva_sns
668    REAL(r_std)                      :: speed
669
670    zicp = un / cp_air
671    !
672    DO ji=1,kjpindex
673      !
674      !
675      ! Help variables
676      !
677      !
678      speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))
679      !
680      zikt = 1/(rau(ji) * speed * q_cdrag(ji))
681      zikq = 1/(rau(ji) * speed * q_cdrag(ji))
682      !
683      !
684      !  The first step is to compute the fluxes for the old surface conditions
685      !
686      !
687      sensfl_old = (petBcoef(ji) -  psold(ji)) / (zikt -  petAcoef(ji))
688      larsub_old = chalsu0 * vbeta1(ji) * (peqBcoef(ji) -  qsol_sat(ji)) / (zikq - peqAcoef(ji))
689      lareva_old = chalev0 * (un - vbeta1(ji)) * vbeta(ji) * &
690           & (peqBcoef(ji) -  valpha(ji) * qsol_sat(ji)) / (zikq - peqAcoef(ji))
691      !
692      !
693      !  Next the sensitivity terms are computed
694      !
695      !
696      netrad_sns = zicp(ji) * quatre * emis(ji) * c_stefan * ((zicp(ji) * psold(ji))**3)
697      sensfl_sns = un / (zikt -  petAcoef(ji))
698      larsub_sns = chalsu0 * vbeta1(ji) * zicp(ji) * pdqsold(ji) / (zikq - peqAcoef(ji))
699      lareva_sns = chalev0 * (un - vbeta1(ji)) * vbeta(ji) * valpha(ji) * &
700           & zicp(ji) * pdqsold(ji) / (zikq - peqAcoef(ji))
701      !
702      !
703      !  Now we are solving the energy balance
704      !
705      !
706      sum_old = netrad(ji) + sensfl_old + larsub_old + lareva_old + soilflx(ji)
707      sum_sns = netrad_sns + sensfl_sns + larsub_sns + lareva_sns
708      dtheta = dtradia * sum_old / (zicp(ji) * soilcap(ji) + dtradia * sum_sns)
709      !
710      !
711      psnew(ji) =  psold(ji) + dtheta
712      !
713      qsol_sat_new(ji) = qsol_sat(ji) + zicp(ji) * pdqsold(ji) * dtheta
714      !
715      temp_sol_new(ji) = psnew(ji) / cp_air
716      !
717      epot_air_new(ji) = zikt * (sensfl_old - sensfl_sns * dtheta) + psnew(ji)
718      !
719      fevap = (lareva_old - lareva_sns * dtheta) + (larsub_old - larsub_sns * dtheta)
720      IF ( ABS(fevap) < EPSILON(un) ) THEN
721        qair_new(ji) = qair(ji)
722      ELSE
723        qair_new(ji) = zikq * un / (chalev0 * (un - vbeta1(ji)) * vbeta(ji) * valpha(ji) + &
724           & chalsu0 *  vbeta1(ji)) * fevap + qsol_sat_new(ji)
725      ENDIF
726      !
727      !
728    ENDDO
729
730    IF (long_print) WRITE (numout,*) ' enerbil_surftemp done '
731
732  END SUBROUTINE enerbil_surftemp
733
734  !! This routine computes tsol_new, netrad, vevapp, fluxlat, fluxsubli and fluxsens
735  !!
736  SUBROUTINE enerbil_flux (kjpindex, dtradia, emis, temp_sol, rau, u, v, q_cdrag, vbeta, valpha, &
737     & vbeta1, qair, epot_air, psnew, qsurf, fluxsens, fluxlat, fluxsubli, vevapp, temp_sol_new, &
738     & lwdown, swnet, lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr)
739
740    ! interface description
741    ! input scalar
742    INTEGER(i_std), INTENT(in)                               :: kjpindex      !! Domain size
743    REAL(r_std), INTENT(in)                                  :: dtradia       !! Time step in seconds
744    ! input fields
745    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: emis          !! Emissivity
746    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_sol      !! Soil temperature
747    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: rau           !! Density
748    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u,v           !! wind
749    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q_cdrag       !!
750    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: vbeta         !! Resistance coefficient
751    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: valpha        !! Resistance coefficient
752    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: vbeta1        !! Snow resistance
753    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair          !! Lowest level specific humidity
754    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: epot_air      !! Air potential energy
755    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: psnew         !! New surface static energy
756    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_sol_new  !! New soil temperature
757    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: pb            !! Lowest level pressure
758    ! output fields
759    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsurf         !! Surface specific humidity
760    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsens      !! Sensible chaleur flux
761    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat       !! Latent chaleur flux
762    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsubli     !! Energy of sublimation
763    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapp        !! Total of evaporation
764    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: lwdown        !! Downward Long-wave radiation
765    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swnet         !! Net SW radiation
766    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: lwup          !! Long-wave up radiation
767    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: lwnet         !! Long-wave net radiation
768    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: tsol_rad      !! Radiative surface temperature
769    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: netrad        !! Net radiation
770    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: evapot        !! Soil Potential Evaporation
771    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: evapot_corr   !! Soil Potential Evaporation Correction
772
773
774    ! local declaration
775    INTEGER(i_std)                                 :: ji
776    REAL(r_std),DIMENSION (kjpindex)                :: grad_qsat
777    REAL(r_std)                                     :: correction
778    REAL(r_std)                                     :: speed, qc
779    ! initialisation
780
781    !
782    ! 1. computes temp_sol_new, netrad, vevapp, fluxlat, fluxsubli, fluxsens
783    !
784
785    DO ji=1,kjpindex
786
787      speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))
788      qc = speed * q_cdrag(ji)
789
790      lwup(ji) = emis(ji) * c_stefan * temp_sol(ji)**4 + &
791           &     quatre * emis(ji) * c_stefan * temp_sol(ji)**3 * &
792           &     (temp_sol_new(ji) - temp_sol(ji)) 
793!!
794!!    Add the reflected LW radiation
795!!
796      lwup(ji) = lwup(ji)  +  (un - emis(ji)) * lwdown(ji)
797
798!!      tsol_rad(ji) =  (lwup(ji)/ (emis(ji) * c_stefan)) **(1./quatre)
799!!  Need to check the equations
800!!
801      tsol_rad(ji) = emis(ji) * c_stefan * temp_sol(ji)**4 + lwup(ji)
802!!
803!!  This is a simple diagnostic which will be used by the GCM to compute the dependence of
804!!  of the surface layer stability on moisture.
805!!
806      qsurf(ji) = vbeta1(ji) * qsol_sat_new(ji) + &
807           & (un - vbeta1(ji)) * vbeta(ji) * valpha(ji) * qsol_sat_new(ji)
808      qsurf(ji) = MAX(qsurf(ji), qair(ji))
809     
810      netrad(ji) = lwdown(ji) + swnet(ji) - lwup(ji) 
811
812      vevapp(ji) = dtradia * rau(ji) * qc * vbeta1(ji) * (qsol_sat_new(ji) - qair(ji)) + &
813           & dtradia * rau(ji) * qc * (un - vbeta1(ji)) * vbeta(ji) * &
814           & (valpha(ji) * qsol_sat_new(ji) - qair(ji))
815
816      fluxlat(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) *&
817           & (qsol_sat_new(ji) - qair(ji)) + &
818           &  chalev0 * rau(ji) * qc * (un - vbeta1(ji)) * vbeta(ji) * &
819           & (valpha(ji) * qsol_sat_new(ji) - qair(ji))
820
821      fluxsubli(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) *&
822           & (qsol_sat_new(ji) - qair(ji)) 
823
824      fluxsens(ji) =  rau(ji) * qc * (psnew(ji) - epot_air(ji))
825
826      lwnet(ji) = lwdown(ji) - lwup(ji)
827 
828      evapot(ji) =  MAX(zero, dtradia * rau(ji) * qc * (qsol_sat_new(ji) - qair(ji)))
829
830      tair(ji)  =  epot_air(ji) / cp_air
831
832     ENDDO 
833
834! define qsat_air avec subroutine de src_parameter:
835
836    CALL qsatcalc(kjpindex, tair, pb, qsat_air)
837
838    CALL dev_qsatcalc(kjpindex, tair, pb, grad_qsat)
839!    grad_qsat(:)= (qsol_sat_new(:)- qsat_air(:)) / ((psnew(:) - epot_air(:)) / cp_air) ! * dtradia
840    !- Penser a sortir evapot en meme temps qu'evapot_corr tdo.
841    DO ji=1,kjpindex
842
843      speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))
844      qc = speed * q_cdrag(ji)
845
846       IF ((evapot(ji) .GT. zero) .AND. ((psnew(ji) - epot_air(ji)) .NE. zero )) THEN
847
848          correction =  (quatre * emis(ji) * c_stefan * tair(ji)**3 + rau(ji) * qc * cp_air + &
849               &                  chalev0 * rau(ji) * qc * grad_qsat(ji) * vevapp(ji) / evapot(ji) )
850          IF (ABS(correction) .GT. min_sechiba) THEN
851             correction = chalev0 * rau(ji) * qc * grad_qsat(ji) * (un - vevapp(ji)/evapot(ji)) / correction
852          ELSE
853             WRITE(numout,*) "Denominateur de la correction de milly nul! Aucune correction appliquee"
854          ENDIF
855       ELSE
856          correction = zero
857       ENDIF
858       correction = MAX (zero, correction)
859       
860       evapot_corr(ji) = evapot(ji) / (un + correction)
861       
862    ENDDO
863
864    IF (long_print) WRITE (numout,*) ' enerbil_flux done '
865
866  END SUBROUTINE enerbil_flux
867
868  !! This routine computes evaporation and transpiration
869  !!
870  SUBROUTINE enerbil_evapveg (kjpindex, dtradia, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, cimean, &
871     & ccanopy, rau, u, v, q_cdrag, qair, humrel, vevapsno, vevapnu , vevapwet, transpir, gpp, evapot)
872
873    ! interface description
874    ! input scalar
875    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size
876    REAL(r_std), INTENT(in)                                  :: dtradia           !! Time step in seconds
877    ! input fields
878    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: vbeta1           !! Snow resistance
879    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: vbeta4           !! Bare soil resistance
880    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: rau              !! Density
881    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u, v             !! Wind
882    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q_cdrag          !!
883    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair             !! Lowest level specific humidity
884    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ccanopy          !! CO2 concentration in the canopy
885    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: evapot           !! Soil Potential Evaporation
886    REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in)       :: humrel           !! Relative humidity
887    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: vbeta2           !! Interception resistance
888    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: vbeta3           !! Vegetation resistance
889    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: vbetaco2         !! Vegetation resistance to CO2
890    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: cimean           !! mean Ci
891    ! output fields
892    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapsno         !! Snow evaporation
893    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapnu          !! Bare soil evaporation
894    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: transpir         !! Transpiration
895    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: gpp              !! Assimilation, gC/m**2 total area.
896    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: vevapwet         !! Interception
897
898    ! local declaration
899    INTEGER(i_std)                                :: ji, jv
900    REAL(r_std), DIMENSION(kjpindex)               :: xx
901    REAL(r_std)                                    :: speed
902
903    ! initialisation
904
905    !
906    ! 1. computes vevapsno, vevapnu
907    !
908
909    DO ji=1,kjpindex 
910
911      speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))
912
913      !
914      ! 1.1 snow sublimation
915      !
916      vevapsno(ji) = vbeta1(ji) * dtradia * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji))
917
918      !
919      ! 1.2 bare soil evaporation
920      !
921      vevapnu(ji) = (un - vbeta1(ji)) * vbeta4(ji) * dtradia * rau(ji) * speed * q_cdrag(ji) &
922         & * (qsol_sat_new(ji) - qair(ji))
923
924    END DO
925
926    !
927    ! 2. computes transpir and vevapwet
928    !
929
930    DO ji = 1, kjpindex
931       !
932       speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))
933       !
934       xx(ji) = dtradia * (un-vbeta1(ji)) * (qsol_sat_new(ji)-qair(ji)) * rau(ji) * speed * q_cdrag(ji)
935       !
936    ENDDO
937    !
938    DO jv=1,nvm 
939      DO ji=1,kjpindex 
940        !
941        ! 2.1 Interception loss
942        !
943        vevapwet(ji,jv) = xx(ji) * vbeta2(ji,jv)
944        !
945        ! 2.2 Transpiration
946        !
947        transpir (ji,jv) = xx(ji) * vbeta3(ji,jv)
948        !
949      END DO
950    END DO
951
952    !
953    ! 3 STOMATE: Assimilation
954    !
955
956    IF ( control%ok_co2 ) THEN
957
958      DO jv = 1, nvm
959         DO ji = 1, kjpindex
960        speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))
961
962        gpp(ji,jv) = vbetaco2(ji,jv) * dtradia * rau(ji) * speed * q_cdrag(ji) * &
963                    (ccanopy(ji) - cimean(ji,jv)) * 12.e-6
964        ENDDO
965      ENDDO
966
967    ELSEIF ( control%stomate_watchout ) THEN
968
969      gpp(:,:) = 0.0
970
971    ENDIF
972
973    IF (long_print) WRITE (numout,*) ' enerbil_evapveg done '
974
975  END SUBROUTINE enerbil_evapveg
976
977  !! Second part of main routine for enerbil module
978  !! - called every time step
979  !!
980  !! Algorithm:
981  !! computes new soil temperature due to ice and snow melt
982  !!
983  SUBROUTINE enerbil_fusion (kjpindex, dtradia, tot_melt, soilcap, temp_sol_new, fusion )
984
985    ! interface description
986    ! input scalar
987    INTEGER(i_std), INTENT(in)                               :: kjpindex       !! Domain size
988    REAL(r_std), INTENT(in)                                  :: dtradia        !! Time step in seconds
989    ! input fields
990    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: tot_melt       !! Total melt
991    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: soilcap        !! Soil calorific capacity
992    ! modified fields
993    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: temp_sol_new   !! New soil temperature
994    ! output fields
995    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fusion         !! Fusion
996
997    ! local declaration
998    INTEGER(i_std)                                :: ji
999
1000    ! initialisation
1001    IF (long_print) WRITE (numout,*) ' enerbil_fusion start ', MINVAL(soilcap), MINLOC(soilcap),&
1002         & MAXVAL(soilcap), MAXLOC(soilcap)
1003    !
1004    ! 1. computes new soil temperature due to ice and snow melt
1005    !
1006    DO ji=1,kjpindex 
1007
1008      fusion(ji) = tot_melt(ji) * chalfu0 / dtradia
1009
1010      temp_sol_new(ji) = temp_sol_new(ji) - ((tot_melt(ji) * chalfu0) / soilcap(ji))
1011
1012    END DO
1013
1014    IF (long_print) WRITE (numout,*) ' enerbil_fusion done '
1015
1016  END SUBROUTINE enerbil_fusion
1017
1018  !! Diagnose 2 meter air temperature
1019  !!
1020  SUBROUTINE enerbil_t2mdiag (kjpindex, temp_air, t2mdiag)
1021
1022    ! interface description
1023    ! input scalar
1024    INTEGER(i_std), INTENT(in)                         :: kjpindex       !! Domain size
1025    ! input fields
1026    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: temp_air       !! Air temperature in Kelvin
1027    ! modified fields
1028    ! output fields
1029    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: t2mdiag        !! 2-meter temperature
1030
1031    t2mdiag(:) = temp_air(:)
1032
1033    IF (long_print) WRITE (numout,*) ' enerbil_t2mdiag done '
1034
1035  END SUBROUTINE enerbil_t2mdiag
1036
1037END MODULE enerbil
Note: See TracBrowser for help on using the repository browser.