Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/thermosoil.f90 @ 326

Last change on this file since 326 was 257, checked in by didier.solyga, 13 years ago
Externalized version merged with the trunk
File size: 37.6 KB

Line
1	!!
2	!! This module computes soil thermodynamic
3	!!
4	!! @author Marie-Alice Foujols and Jan Polcher
5	!! @Version : $Revision: 45 $, $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $
6	!!
7	!< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/thermosoil.f90 $
8	!< $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $
9	!< $Author: mmaipsl $
10	!< $Revision: 45 $
11	!! IPSL (2006)
12	!! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
13	!!
14	MODULE thermosoil
15
16	! routines called : restput, restget
17	!
18	USE ioipsl
19	!
20	! modules used :
21	USE constantes
22	USE sechiba_io
23	USE grid
24	USE parallel
25
26	IMPLICIT NONE
27
28	! public routines :
29	! thermosoil_main
30	PRIVATE
31	PUBLIC :: thermosoil_main,thermosoil_clear
32
33	!
34	! variables used inside thermosoil module : declaration and initialisation
35	!
36	LOGICAL, SAVE :: l_first_thermosoil=.TRUE. !! Initialisation has to be done one time
37
38	CHARACTER(LEN=80) , SAVE :: var_name !! To store variables names for I/O
39	REAL(r_std), SAVE :: lambda, cstgrnd, lskin, fz1, zalph
40
41	! two dimensions array allocated, computed, saved and got in thermosoil module
42
43	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn !! Different levels soil temperature
44
45	! one dimension array allocated, computed and used in thermosoil module exclusively
46
47	REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz !!
48	REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz_coef
49	REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz1 !!
50	REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz2 !!
51
52	! one dimension array allocated, computed and used in thermosoil module exclusively
53
54	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z1 !!
55
56	! two dimensions arrays allocated, computed and used in thermosoil module exclusively
57
58	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !!
59	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !!
60	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !!
61	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !!
62	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz1 !!
63	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz2 !!
64	!
65	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! Capacity used for energy_incr calculation
66	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg !! Temperature at the beginning of the time step
67	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at the beginning of the timestep
68	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat
69	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow cold content
70	!!
71	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: wetdiag !! Soil weetness on the thermodynamical levels
72	!!
73	CONTAINS
74
75	!! Main routine for thermosoil module.
76	!!
77	!! Algorithm:
78	!! - call thermosoil_var_init to initialise variables
79	!! - call thermosoil_coef for coefficient
80	!! - call thermosoil_profile for soil profiling
81	!!
82	!! @call thermosoil_var_init
83	!! @call thermosoil_coef
84	!! @call thermosoil_profile
85	!!
86	SUBROUTINE thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexgrnd, &
87	& temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id)
88
89	! interface description
90	! input scalar
91	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
92	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
93	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and history file identifier
94	INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier
95	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
96	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
97	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
98	! input fields
99	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
100	INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in) :: indexgrnd !! Indeces of the points on the 3D map
101	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
102	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow quantity
103	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostic of relative humidity
104	! output fields
105	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! Soil capacity
106	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx
107	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (inout):: stempdiag !! diagnostic temp profile
108
109	REAL(r_std),DIMENSION (kjpindex,ngrnd) :: temp
110	REAL(r_std),DIMENSION (kjpindex,ngrnd-1) :: temp1
111	REAL(r_std),DIMENSION (kjpindex) :: temp2
112	!
113	! do initialisation
114	!
115	IF (l_first_thermosoil) THEN
116
117	IF (long_print) WRITE (numout,*) ' l_first_thermosoil : call thermosoil_init '
118
119	!
120	! do needed allocation
121	!
122	CALL thermosoil_init (kjit, ldrestart_read, kjpindex, index, rest_id)
123
124	!
125	! computes some physical constantes and array depending soil level depth
126	!
127	CALL thermosoil_var_init (kjpindex, zz, zz_coef, dz1, dz2, pkappa, pcapa, pcapa_en, shumdiag, stempdiag)
128
129	!
130	! computes cgrd and dgrd coefficient from previous time step (restart)
131	!
132	CALL thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,&
133	& zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa)
134
135	CALL thermosoil_energy (kjpindex, temp_sol_new, soilcap, .TRUE.)
136
137	IF (ldrestart_read) THEN
138	IF (long_print) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables'
139
140	var_name= 'cgrnd'
141	CALL ioconf_setatt('UNITS', '-')
142	CALL ioconf_setatt('LONG_NAME','Cgrnd coefficient.')
143	IF ( ok_var(var_name) ) THEN
144	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
145	IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
146	cgrnd(:,:)=temp1(:,:)
147	ENDIF
148	ENDIF
149
150	var_name= 'dgrnd'
151	CALL ioconf_setatt('UNITS', '-')
152	CALL ioconf_setatt('LONG_NAME','Dgrnd coefficient.')
153	IF ( ok_var(var_name) ) THEN
154	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
155	IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
156	dgrnd(:,:)=temp1(:,:)
157	ENDIF
158	ENDIF
159
160	var_name= 'z1'
161	CALL ioconf_setatt('UNITS', '-')
162	CALL ioconf_setatt('LONG_NAME','?.')
163	IF ( ok_var(var_name) ) THEN
164	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g)
165	IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN
166	z1(:)=temp2(:)
167	ENDIF
168	ENDIF
169
170	var_name= 'pcapa'
171	CALL ioconf_setatt('UNITS', '-')
172	CALL ioconf_setatt('LONG_NAME','?.')
173	IF ( ok_var(var_name) ) THEN
174	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
175	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
176	pcapa(:,:)=temp(:,:)
177	ENDIF
178	ENDIF
179
180	var_name= 'pcapa_en'
181	CALL ioconf_setatt('UNITS', '-')
182	CALL ioconf_setatt('LONG_NAME','?.')
183	IF ( ok_var(var_name) ) THEN
184	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
185	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
186	pcapa_en(:,:)=temp(:,:)
187	ENDIF
188	ENDIF
189
190	var_name= 'pkappa'
191	CALL ioconf_setatt('UNITS', '-')
192	CALL ioconf_setatt('LONG_NAME','?.')
193	IF ( ok_var(var_name) ) THEN
194	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
195	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
196	pkappa(:,:)=temp(:,:)
197	ENDIF
198	ENDIF
199
200	var_name= 'zdz1'
201	CALL ioconf_setatt('UNITS', '-')
202	CALL ioconf_setatt('LONG_NAME','?.')
203	IF ( ok_var(var_name) ) THEN
204	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
205	IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
206	zdz1(:,:)=temp1(:,:)
207	ENDIF
208	ENDIF
209
210	var_name= 'zdz2'
211	CALL ioconf_setatt('UNITS', '-')
212	CALL ioconf_setatt('LONG_NAME','?.')
213	IF ( ok_var(var_name) ) THEN
214	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
215	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
216	zdz2(:,:)=temp(:,:)
217	ENDIF
218	ENDIF
219
220	var_name='temp_sol_beg'
221	CALL ioconf_setatt('UNITS', 'K')
222	CALL ioconf_setatt('LONG_NAME','Old Surface temperature')
223	IF ( ok_var(var_name) ) THEN
224	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g)
225	IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN
226	temp_sol_beg(:) = temp2(:)
227	ENDIF
228	ENDIF
229
230	ENDIF
231
232	RETURN
233
234	ENDIF
235
236	!
237	! prepares restart file for the next simulation
238	!
239	IF (ldrestart_write) THEN
240
241	IF (long_print) WRITE (numout,*) ' we have to complete restart file with THERMOSOIL variables'
242
243	var_name= 'ptn'
244	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, ptn, 'scatter', nbp_glo, index_g)
245
246	var_name= 'cgrnd'
247	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g)
248	var_name= 'dgrnd'
249	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g)
250
251	var_name= 'z1'
252	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, z1, 'scatter', nbp_glo, index_g)
253
254	var_name= 'pcapa'
255	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pcapa, 'scatter', nbp_glo, index_g)
256
257	var_name= 'pcapa_en'
258	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pcapa_en, 'scatter', nbp_glo, index_g)
259
260	var_name= 'pkappa'
261	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pkappa, 'scatter', nbp_glo, index_g)
262
263	var_name= 'zdz1'
264	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, zdz1, 'scatter', nbp_glo, index_g)
265
266	var_name= 'zdz2'
267	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, zdz2, 'scatter', nbp_glo, index_g)
268
269	var_name= 'temp_sol_beg'
270	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol_beg, 'scatter', nbp_glo, index_g)
271
272	var_name= 'soilcap'
273	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g)
274
275	var_name= 'soilflx'
276	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g)
277
278	! read in enerbil
279	var_name= 'temp_sol_new'
280	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol_new, 'scatter', nbp_glo, index_g)
281
282	RETURN
283
284	END IF
285
286	!
287	! Put the soil wetnesss diagnostic on the levels of the soil temprature
288	!
289	CALL thermosoil_humlev(kjpindex, shumdiag)
290
291	!
292	! computes profile with previous cgrd and dgrd coefficient
293	!
294	CALL thermosoil_profile (kjpindex, temp_sol_new, ptn, stempdiag)
295	CALL thermosoil_energy (kjpindex, temp_sol_new, soilcap, .FALSE.)
296	!
297	IF ( .NOT. almaoutput ) THEN
298	CALL histwrite(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
299	ELSE
300	CALL histwrite(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
301	CALL histwrite(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
302	CALL histwrite(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
303	CALL histwrite(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
304	ENDIF
305	IF ( hist2_id > 0 ) THEN
306	IF ( .NOT. almaoutput ) THEN
307	CALL histwrite(hist2_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
308	ELSE
309	CALL histwrite(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
310	CALL histwrite(hist2_id, 'Qg', kjit, soilflx, kjpindex, index)
311	CALL histwrite(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
312	CALL histwrite(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
313	ENDIF
314	ENDIF
315	!
316	! computes cgrd and dgrd coefficient
317	!
318	CALL thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,&
319	& zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa)
320
321	IF (long_print) WRITE (numout,*) ' thermosoil_main done '
322
323	END SUBROUTINE thermosoil_main
324
325	!! Initialisation for thermosoil module.
326	!! - does dynamic allocation for local arrays
327	!! - reads _restart_ file or set initial values to a raisonable value
328	!!
329	SUBROUTINE thermosoil_init(kjit, ldrestart_read, kjpindex, index, rest_id)
330
331	! interface description
332	! input scalar
333	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
334	LOGICAL,INTENT (in) :: ldrestart_read !! Logical for restart file to read
335	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
336	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
337
338	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier
339	! input fields
340	! output fields
341
342	! local declaration
343	INTEGER(i_std) :: ier
344
345	! initialisation
346	IF (l_first_thermosoil) THEN
347	l_first_thermosoil=.FALSE.
348	ELSE
349	WRITE (numout,*) ' l_first_thermosoil false . we stop '
350	STOP 'thermosoil_init'
351	ENDIF
352
353	! two dimensions array allocation
354
355	ALLOCATE (ptn(kjpindex,ngrnd),stat=ier)
356	IF (ier.NE.0) THEN
357	WRITE (numout,*) ' error in ptn allocation. We stop. We need ',kjpindex,' fois ',ngrnd,' words = '&
358	& , kjpindex*ngrnd
359	STOP 'thermosoil_init'
360	END IF
361
362	! one dimension array allocation
363
364	ALLOCATE (z1(kjpindex),stat=ier)
365	IF (ier.NE.0) THEN
366	WRITE (numout,*) ' error in z1 allocation. We STOP. We need ',kjpindex,' words '
367	STOP 'thermosoil_init'
368	END IF
369
370	! two dimension array allocation
371
372	ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier)
373	IF (ier.NE.0) THEN
374	WRITE (numout,*) ' error in cgrnd allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '&
375	& , kjpindex*(ngrnd-1)
376	STOP 'thermosoil_init'
377	END IF
378
379	ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier)
380	IF (ier.NE.0) THEN
381	WRITE (numout,*) ' error in dgrnd allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '&
382	& , kjpindex*(ngrnd-1)
383	STOP 'thermosoil_init'
384	END IF
385
386	ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier)
387	IF (ier.NE.0) THEN
388	WRITE (numout,*) ' error in pcapa allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
389	& , kjpindex*ngrnd
390	STOP 'thermosoil_init'
391	END IF
392
393	ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier)
394	IF (ier.NE.0) THEN
395	WRITE (numout,*) ' error in pkappa allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
396	& , kjpindex*ngrnd
397	STOP 'thermosoil_init'
398	END IF
399
400	ALLOCATE (zdz1(kjpindex,ngrnd-1),stat=ier)
401	IF (ier.NE.0) THEN
402	WRITE (numout,*) ' error in zdz1 allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '&
403	& , kjpindex*(ngrnd-1)
404	STOP 'thermosoil_init'
405	END IF
406
407	ALLOCATE (zdz2(kjpindex,ngrnd),stat=ier)
408	IF (ier.NE.0) THEN
409	WRITE (numout,*) ' error in zdz2 allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
410	& , kjpindex*ngrnd
411	STOP 'thermosoil_init'
412	END IF
413
414	ALLOCATE (surfheat_incr(kjpindex),stat=ier)
415	IF (ier.NE.0) THEN
416	WRITE (numout,*) ' error in surfheat_incr allocation. We STOP. We need ',kjpindex,' words = '&
417	& , kjpindex
418	STOP 'thermosoil_init'
419	END IF
420
421	ALLOCATE (coldcont_incr(kjpindex),stat=ier)
422	IF (ier.NE.0) THEN
423	WRITE (numout,*) ' error in coldcont_incr allocation. We STOP. We need ',kjpindex,' words = '&
424	& , kjpindex
425	STOP 'thermosoil_init'
426	END IF
427
428	ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier)
429	IF (ier.NE.0) THEN
430	WRITE (numout,*) ' error in pcapa_en allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
431	& , kjpindex*ngrnd
432	STOP 'thermosoil_init'
433	END IF
434
435	ALLOCATE (ptn_beg(kjpindex,ngrnd),stat=ier)
436	IF (ier.NE.0) THEN
437	WRITE (numout,*) ' error in ptn_beg allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
438	& , kjpindex*ngrnd
439	STOP 'thermosoil_init'
440	END IF
441
442	ALLOCATE (temp_sol_beg(kjpindex),stat=ier)
443	IF (ier.NE.0) THEN
444	WRITE (numout,*) ' error in temp_sol_beg allocation. We STOP. We need ',kjpindex,' words = '&
445	& , kjpindex
446	STOP 'thermosoil_init'
447	END IF
448
449	ALLOCATE (wetdiag(kjpindex,ngrnd),stat=ier)
450	IF (ier.NE.0) THEN
451	WRITE (numout,*) ' error in wetdiag allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
452	& , kjpindex*ngrnd
453	STOP 'thermosoil_init'
454	END IF
455	!
456	! open restart input file done by sechiba_init
457	! and read data from restart input file for THERMOSOIL process
458
459	IF (ldrestart_read) THEN
460	IF (long_print) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables'
461
462	var_name= 'ptn'
463	CALL ioconf_setatt('UNITS', 'K')
464	CALL ioconf_setatt('LONG_NAME','Soil Temperature profile')
465	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., ptn, "gather", nbp_glo, index_g)
466	!
467	! change restart If they were not found in the restart file
468	!
469	!Config Key = THERMOSOIL_TPRO
470	!Config Desc = Initial soil temperature profile if not found in restart
471	!Config Def = 280.
472	!Config Help = The initial value of the temperature profile in the soil if
473	!Config its value is not found in the restart file. This should only
474	!Config be used if the model is started without a restart file. Here
475	!Config we only require one value as we will assume a constant
476	!Config throughout the column.
477	!
478	CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std)
479
480	ENDIF
481
482	IF (long_print) WRITE (numout,*) ' thermosoil_init done '
483
484	END SUBROUTINE thermosoil_init
485
486	!! Function for distributing the levels
487	!!
488	SUBROUTINE thermosoil_clear()
489
490	l_first_thermosoil=.TRUE.
491
492	IF ( ALLOCATED (ptn)) DEALLOCATE (ptn)
493	IF ( ALLOCATED (z1)) DEALLOCATE (z1)
494	IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd)
495	IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd)
496	IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa)
497	IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa)
498	IF ( ALLOCATED (zdz1)) DEALLOCATE (zdz1)
499	IF ( ALLOCATED (zdz2)) DEALLOCATE (zdz2)
500	IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en)
501	IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg)
502	IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg)
503	IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr)
504	IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr)
505	IF ( ALLOCATED (wetdiag)) DEALLOCATE (wetdiag)
506
507	END SUBROUTINE thermosoil_clear
508	!!
509	!!
510	FUNCTION fz(rk) RESULT (fz_result)
511
512	! interface
513	! input value
514	REAL(r_std), INTENT(in) :: rk
515	! output value
516	REAL(r_std) :: fz_result
517
518	fz_result = fz1 * (zalph ** rk - un) / (zalph - un)
519
520	END FUNCTION fz
521
522	!! Thermosoil's variables initialisation
523	!!
524	SUBROUTINE thermosoil_var_init(kjpindex, zz, zz_coef, dz1, dz2, pkappa, pcapa, pcapa_en, shumdiag, stempdiag)
525
526	! interface description
527	! input scalar
528	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
529	! input fields
530	REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostic humidity profile
531	! output fields
532	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz !!
533	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz_coef
534	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz1 !!
535	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz2 !! tailles des couches
536	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa !!
537	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa_en
538	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pkappa !!
539	REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostic temp profile
540
541	! local declaration
542	INTEGER(i_std) :: ier, ji, jg
543	REAL(r_std) :: sum
544
545	!
546	! 0. initialisation
547	!
548	cstgrnd=SQRT(one_day / pi)
549	lskin = SQRT(so_cond / so_capa * one_day / pi)
550	fz1 = 0.3_r_std * cstgrnd
551	zalph = deux
552	!
553	! 1. Computing the depth of the Temperature level, using a
554	! non dimentional variable x = z/lskin, lskin beeing
555	! the skin depth
556	!
557
558	DO jg=1,ngrnd
559	!!! This needs to be solved soon. Either we allow CPP options in SECHIBA or the VPP
560	!!! fixes its compiler !
561	!!!#ifdef VPP5000
562	dz2(jg) = fz(REAL(jg,r_std)-undemi+undemi) - fz(REAL(jg-1,r_std)-undemi+undemi)
563	!!!#else
564	!!! dz2(jg) = fz(REAL(jg,r_std)) - fz(REAL(jg-1,r_std))
565	!!!#endif
566	ENDDO
567	!
568	! 1.2 The undimentional depth is computed.
569	! ------------------------------------
570	DO jg=1,ngrnd
571	zz(jg) = fz(REAL(jg,r_std) - undemi)
572	zz_coef(jg) = fz(REAL(jg,r_std)-undemi+undemi)
573	ENDDO
574	!
575	! 1.3 Converting to meters.
576	! --------------------
577	DO jg=1,ngrnd
578	zz(jg) = zz(jg) / cstgrnd * lskin
579	zz_coef(jg) = zz_coef(jg) / cstgrnd * lskin
580	dz2(jg) = dz2(jg) / cstgrnd * lskin
581	ENDDO
582	!
583	! 1.4 Computing some usefull constants.
584	! --------------------------------
585	DO jg=1,ngrnd-1
586	dz1(jg) = un / (zz(jg+1) - zz(jg))
587	ENDDO
588	lambda = zz(1) * dz1(1)
589	!
590	! 1.5 Get the wetness profice on this grid
591	! ------------------------------------
592	!
593	CALL thermosoil_humlev(kjpindex, shumdiag)
594	!
595	! 1.6 Thermal conductivity at all levels.
596	! ----------------------------------
597	DO jg = 1,ngrnd
598	DO ji = 1,kjpindex
599	pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry)
600	pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
601	pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
602	ENDDO
603	ENDDO
604	!
605	! 2. Diagnostics.
606	! -----------
607
608	WRITE (numout,*) 'diagnostic des niveaux dans le sol'
609	WRITE (numout,*) 'niveaux intermediaires et pleins'
610	sum = zero
611	DO jg=1,ngrnd
612	sum = sum + dz2(jg)
613	WRITE (numout,*) zz(jg),sum
614	ENDDO
615
616	!
617	! 3. Compute the diagnostic temperature profile
618	!
619
620	CALL thermosoil_diaglev(kjpindex, stempdiag)
621	!
622	!
623
624	IF (long_print) WRITE (numout,*) ' thermosoil_var_init done '
625
626	END SUBROUTINE thermosoil_var_init
627
628	!! Computation of cgrnd and dgrnd coefficient at the t time-step.
629	!!
630	!! Needs previous time step values.
631	!!
632	SUBROUTINE thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,&
633	& zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa)
634
635	! interface description
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) :: temp_sol_new !!
641	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !!
642	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: zz !!
643	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz1 !!
644	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz2 !!
645	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (in) :: ptn
646	! output fields
647	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !!
648	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !!
649	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: z1 !!
650	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa !!
651	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa_en !!
652	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pkappa !!
653	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !!
654	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !!
655	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: zdz1 !!
656	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: zdz2 !!
657
658	! local declaration
659	INTEGER(i_std) :: ji, jg
660	REAL(r_std), DIMENSION(kjpindex) :: snow_h, zx1, zx2
661
662	!
663	! objet: computation of cgrnd and dgrnd coefficient at the t time-step.
664	! ------
665	!
666	! ---------------------------------------------------------------
667	! Computation of the Cgrd and Dgrd coefficient for the next step:
668	! ---------------------------------------------------------------
669	!
670	DO ji = 1,kjpindex
671	snow_h(ji) = snow(ji) / sn_dens
672	!
673	! Traitement special pour le premiere couche
674	!
675	IF ( snow_h(ji) .GT. zz_coef(1) ) THEN
676	pcapa(ji,1) = sn_capa
677	pcapa_en(ji,1) = sn_capa
678	pkappa(ji,1) = sn_cond
679	ELSE IF ( snow_h(ji) .GT. zero ) THEN
680	pcapa_en(ji,1) = sn_capa
681	zx1(ji) = snow_h(ji) / zz_coef(1)
682	zx2(ji) = ( zz_coef(1) - snow_h(ji)) / zz_coef(1)
683	pcapa(ji,1) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet
684	pkappa(ji,1) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet )
685	ELSE
686	pcapa(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry)
687	pkappa(ji,1) = so_cond_dry + wetdiag(ji,1)*(so_cond_wet - so_cond_dry)
688	pcapa_en(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry)
689	ENDIF
690	!
691	DO jg = 2, ngrnd - 2
692	IF ( snow_h(ji) .GT. zz_coef(jg) ) THEN
693	pcapa(ji,jg) = sn_capa
694	pkappa(ji,jg) = sn_cond
695	pcapa_en(ji,jg) = sn_capa
696	ELSE IF ( snow_h(ji) .GT. zz_coef(jg-1) ) THEN
697	zx1(ji) = (snow_h(ji) - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1))
698	zx2(ji) = ( zz_coef(jg) - snow_h(ji)) / (zz_coef(jg) - zz_coef(jg-1))
699	pcapa(ji,jg) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet
700	pkappa(ji,jg) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet )
701	pcapa_en(ji,jg) = sn_capa
702	ELSE
703	pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
704	pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry)
705	pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
706	ENDIF
707	ENDDO
708	!
709	!
710	ENDDO
711	!
712	DO jg=1,ngrnd
713	DO ji=1,kjpindex
714	zdz2(ji,jg)=pcapa(ji,jg) * dz2(jg)/dtradia
715	ENDDO
716	ENDDO
717	!
718	DO jg=1,ngrnd-1
719	DO ji=1,kjpindex
720	zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg)
721	ENDDO
722	ENDDO
723	!
724	DO ji = 1,kjpindex
725	z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1)
726	cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji)
727	dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji)
728	ENDDO
729
730	DO jg = ngrnd-1,2,-1
731	DO ji = 1,kjpindex
732	z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg)))
733	cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji)
734	dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji)
735	ENDDO
736	ENDDO
737
738	! ---------------------------------------------------------
739	! computation of the surface diffusive flux from ground and
740	! calorific capacity of the ground:
741	! ---------------------------------------------------------
742
743	DO ji = 1,kjpindex
744	soilflx(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1))
745	soilcap(ji) = (zdz2(ji,1) * dtradia + dtradia * (un - dgrnd(ji,1)) * zdz1(ji,1))
746	z1(ji) = lambda * (un - dgrnd(ji,1)) + un
747	soilcap(ji) = soilcap(ji) / z1(ji)
748	soilflx(ji) = soilflx(ji) + &
749	& soilcap(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dtradia
750	ENDDO
751
752	IF (long_print) WRITE (numout,*) ' thermosoil_coef done '
753
754	END SUBROUTINE thermosoil_coef
755
756	!! Computation of : the ground temperature evolution
757	!!
758	SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, ptn, stempdiag)
759
760	! interface description
761	! input scalar
762	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
763	! input fields
764	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
765	! modified fields
766	REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (inout) :: ptn !! Different levels soil temperature
767	! output fields
768	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostoc profile
769
770	! local declaration
771	INTEGER(i_std) :: ji, jg
772	!
773	! objet: computation of : the ground temperature evolution
774	! ------
775	!
776	! Method: implicit time integration
777	! -------
778	! Consecutives ground temperatures are related by:
779	! T(k+1) = C(k) + D(k)*T(k) (1)
780	! the coefficients C and D are computed at the t-dt time-step.
781	! Routine structure:
782	! -new temperatures are computed using (1)
783	!
784	!
785	! surface temperature
786	DO ji = 1,kjpindex
787	ptn(ji,1) = (lambda * cgrnd(ji,1) + temp_sol_new(ji)) / (lambda * (un - dgrnd(ji,1)) + un)
788	ENDDO
789
790	! other temperatures
791	DO jg = 1,ngrnd-1
792	DO ji = 1,kjpindex
793	ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg)
794	ENDDO
795	ENDDO
796
797	CALL thermosoil_diaglev(kjpindex, stempdiag)
798
799	IF (long_print) WRITE (numout,*) ' thermosoil_profile done '
800
801	END SUBROUTINE thermosoil_profile
802	!!
803	!!
804	!! Diagnostic soil temperature profile on new levels
805	!!
806	!!
807	SUBROUTINE thermosoil_diaglev(kjpindex, stempdiag)
808	! interface description
809	! input scalar
810	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
811	! input fields
812	!
813	! modified fields
814	!
815	! output fields
816	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostoc profile
817	!
818	! local variable
819	!
820	INTEGER(i_std) :: ji, jd, jg
821	REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog
822	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfact
823	!
824	LOGICAL, PARAMETER :: check=.FALSE.
825	!
826	!
827	IF ( .NOT. ALLOCATED(intfact)) THEN
828	!
829	ALLOCATE(intfact(nbdl, ngrnd))
830	!
831	prev_diag = zero
832	DO jd = 1, nbdl
833	lev_diag = diaglev(jd)
834	prev_prog = zero
835	DO jg = 1, ngrnd
836	IF ( jg == ngrnd .AND. (prev_prog + dz2(jg)) < lev_diag ) THEN
837	!! Just make sure we cover the deepest layers
838	lev_prog = lev_diag
839	ELSE
840	lev_prog = prev_prog + dz2(jg)
841	ENDIF
842	intfact(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
843	prev_prog = lev_prog
844	ENDDO
845	prev_diag = lev_diag
846	ENDDO
847	!
848	IF ( check ) THEN
849	WRITE(numout,*) 'thermosoil_diagev -- thermosoil_diaglev -- thermosoil_diaglev --'
850	DO jd = 1, nbdl
851	WRITE(numout,*) jd, '-', intfact(jd,1:ngrnd)
852	ENDDO
853	WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM"
854	DO jd = 1, nbdl
855	WRITE(numout,*) jd, '-', SUM(intfact(jd,1:ngrnd))
856	ENDDO
857	WRITE(numout,*) 'thermosoil_diaglev -- thermosoil_diaglev -- thermosoil_diaglev --'
858	ENDIF
859	!
860	ENDIF
861
862	stempdiag(:,:) = zero
863	DO jg = 1, ngrnd
864	DO jd = 1, nbdl
865	DO ji = 1, kjpindex
866	stempdiag(ji,jd) = stempdiag(ji,jd) + ptn(ji,jg)*intfact(jd,jg)
867	ENDDO
868	ENDDO
869	ENDDO
870
871	END SUBROUTINE thermosoil_diaglev
872	!!
873	!!
874	!! Put soil wetness on the temperature levels
875	!!
876	!!
877	SUBROUTINE thermosoil_humlev(kjpindex, shumdiag)
878	! interface description
879	! input scalar
880	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
881	! input fields
882	!
883	! modified fields
884	!
885	! output fields
886	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostoc profile
887	!
888	! local variable
889	!
890	INTEGER(i_std) :: ji, jd, jg
891	REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog
892	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfactw
893	!
894	LOGICAL, PARAMETER :: check=.FALSE.
895	!
896	!
897	IF ( .NOT. ALLOCATED(intfactw)) THEN
898	!
899	ALLOCATE(intfactw(ngrnd, nbdl))
900	!
901	prev_diag = zero
902	DO jd = 1, ngrnd
903	lev_diag = prev_diag + dz2(jd)
904	prev_prog = zero
905	DO jg = 1, nbdl
906	IF ( jg == nbdl .AND. diaglev(jg) < lev_diag ) THEN
907	!! Just make sure we cover the deepest layers
908	lev_prog = lev_diag
909	ELSE
910	lev_prog = diaglev(jg)
911	ENDIF
912	intfactw(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
913	prev_prog = lev_prog
914	ENDDO
915	prev_diag = lev_diag
916	ENDDO
917	!
918	IF ( check ) THEN
919	WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --'
920	DO jd = 1, ngrnd
921	WRITE(numout,*) jd, '-', intfactw(jd,1:nbdl)
922	ENDDO
923	WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM"
924	DO jd = 1, ngrnd
925	WRITE(numout,*) jd, '-', SUM(intfactw(jd,1:nbdl))
926	ENDDO
927	WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --'
928	ENDIF
929	!
930	ENDIF
931
932	wetdiag(:,:) = zero
933	DO jg = 1, nbdl
934	DO jd = 1, ngrnd
935	DO ji = 1, kjpindex
936	wetdiag(ji,jd) = wetdiag(ji,jd) + shumdiag(ji,jg)*intfactw(jd,jg)
937	ENDDO
938	ENDDO
939	ENDDO
940
941	END SUBROUTINE thermosoil_humlev
942
943	SUBROUTINE thermosoil_energy(kjpindex, temp_sol_new, soilcap, first_call)
944	! interface description
945	! input scalar
946	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
947	LOGICAL, INTENT (in) :: first_call !!
948	! input fields
949	!
950	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
951	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil capacity
952	!
953	! modified fields
954	!
955	! output fields
956	!
957	! local variable
958	!
959	INTEGER(i_std) :: ji, jg
960	!
961	!
962	IF (first_call) THEN
963
964	DO ji = 1, kjpindex
965	surfheat_incr(ji) = zero
966	coldcont_incr(ji) = zero
967	temp_sol_beg(ji) = temp_sol_new(ji)
968	!
969	DO jg = 1, ngrnd
970	ptn_beg(ji,jg) = ptn(ji,jg)
971	ENDDO
972	!
973	ENDDO
974
975	RETURN
976
977	ENDIF
978
979	DO ji = 1, kjpindex
980	surfheat_incr(ji) = zero
981	coldcont_incr(ji) = zero
982	ENDDO
983	!
984	! Sum up the energy content of all layers in the soil.
985	!
986	DO ji = 1, kjpindex
987	!
988	IF (pcapa_en(ji,1) .LE. sn_capa) THEN
989	!
990	! Verify the energy conservation in the surface layer
991	!
992	coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
993	surfheat_incr(ji) = zero
994	ELSE
995	!
996	! Verify the energy conservation in the surface layer
997	!
998	surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
999	coldcont_incr(ji) = zero
1000	ENDIF
1001	ENDDO
1002
1003	ptn_beg(:,:) = ptn(:,:)
1004	temp_sol_beg(:) = temp_sol_new(:)
1005
1006	END SUBROUTINE thermosoil_energy
1007
1008	END MODULE thermosoil

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

Download in other formats: