Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_sechiba/enerbil.f90 @ 6268

Last change on this file since 6268 was 566, checked in by didier.solyga, 13 years ago
Update Parameters labels : add Config Units and Config If for all parameters
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 46.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: