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

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

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