source: perso/abdelouhab.djerrah/ORCHIDEE_ORIG/src_sechiba/enerbil.f90 @ 938

Last change on this file since 938 was 47, checked in by mmaipsl, 14 years ago

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