source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/enerbil.f90 @ 112

Last change on this file since 112 was 106, checked in by didier.solyga, 14 years ago

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