Context Navigation

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

Last change on this file since 461 was 461, checked in by didier.solyga, 13 years ago
Correct bcast functions. In sechiba modules, move var_name from save variables to local variables. Correct a save attribute to parameter attribute in constantes_mtc.f90
File size: 37.7 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	REAL(r_std), SAVE :: lambda, cstgrnd, lskin, fz1, zalph
38
39	! two dimensions array allocated, computed, saved and got in thermosoil module
40
41	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn !! Different levels soil temperature
42
43	! one dimension array allocated, computed and used in thermosoil module exclusively
44
45	REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz !!
46	REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz_coef
47	REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz1 !!
48	REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz2 !!
49
50	! one dimension array allocated, computed and used in thermosoil module exclusively
51
52	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z1 !!
53
54	! two dimensions arrays allocated, computed and used in thermosoil module exclusively
55
56	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !!
57	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !!
58	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !!
59	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !!
60	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz1 !!
61	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz2 !!
62	!
63	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! Capacity used for energy_incr calculation
64	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg !! Temperature at the beginning of the time step
65	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at the beginning of the timestep
66	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat
67	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow cold content
68	!!
69	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: wetdiag !! Soil weetness on the thermodynamical levels
70	!!
71	CONTAINS
72
73	!! Main routine for thermosoil module.
74	!!
75	!! Algorithm:
76	!! - call thermosoil_var_init to initialise variables
77	!! - call thermosoil_coef for coefficient
78	!! - call thermosoil_profile for soil profiling
79	!!
80	!! @call thermosoil_var_init
81	!! @call thermosoil_coef
82	!! @call thermosoil_profile
83	!!
84	SUBROUTINE thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexgrnd, &
85	& temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id)
86
87	! interface description
88	! input scalar
89	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
90	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
91	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and history file identifier
92	INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier
93	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
94	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
95	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
96	! input fields
97	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
98	INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in) :: indexgrnd !! Indeces of the points on the 3D map
99	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
100	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow quantity
101	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostic of relative humidity
102	! output fields
103	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! Soil capacity
104	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx
105	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (inout):: stempdiag !! diagnostic temp profile
106	! local
107	REAL(r_std),DIMENSION (kjpindex,ngrnd) :: temp
108	REAL(r_std),DIMENSION (kjpindex,ngrnd-1) :: temp1
109	REAL(r_std),DIMENSION (kjpindex) :: temp2
110	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
111
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	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
345
346	! initialisation
347	IF (l_first_thermosoil) THEN
348	l_first_thermosoil=.FALSE.
349	ELSE
350	WRITE (numout,*) ' l_first_thermosoil false . we stop '
351	STOP 'thermosoil_init'
352	ENDIF
353
354	! two dimensions array allocation
355
356	ALLOCATE (ptn(kjpindex,ngrnd),stat=ier)
357	IF (ier.NE.0) THEN
358	WRITE (numout,*) ' error in ptn allocation. We stop. We need ',kjpindex,' fois ',ngrnd,' words = '&
359	& , kjpindex*ngrnd
360	STOP 'thermosoil_init'
361	END IF
362
363	! one dimension array allocation
364
365	ALLOCATE (z1(kjpindex),stat=ier)
366	IF (ier.NE.0) THEN
367	WRITE (numout,*) ' error in z1 allocation. We STOP. We need ',kjpindex,' words '
368	STOP 'thermosoil_init'
369	END IF
370
371	! two dimension array allocation
372
373	ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier)
374	IF (ier.NE.0) THEN
375	WRITE (numout,*) ' error in cgrnd allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '&
376	& , kjpindex*(ngrnd-1)
377	STOP 'thermosoil_init'
378	END IF
379
380	ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier)
381	IF (ier.NE.0) THEN
382	WRITE (numout,*) ' error in dgrnd allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '&
383	& , kjpindex*(ngrnd-1)
384	STOP 'thermosoil_init'
385	END IF
386
387	ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier)
388	IF (ier.NE.0) THEN
389	WRITE (numout,*) ' error in pcapa allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
390	& , kjpindex*ngrnd
391	STOP 'thermosoil_init'
392	END IF
393
394	ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier)
395	IF (ier.NE.0) THEN
396	WRITE (numout,*) ' error in pkappa allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
397	& , kjpindex*ngrnd
398	STOP 'thermosoil_init'
399	END IF
400
401	ALLOCATE (zdz1(kjpindex,ngrnd-1),stat=ier)
402	IF (ier.NE.0) THEN
403	WRITE (numout,*) ' error in zdz1 allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '&
404	& , kjpindex*(ngrnd-1)
405	STOP 'thermosoil_init'
406	END IF
407
408	ALLOCATE (zdz2(kjpindex,ngrnd),stat=ier)
409	IF (ier.NE.0) THEN
410	WRITE (numout,*) ' error in zdz2 allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
411	& , kjpindex*ngrnd
412	STOP 'thermosoil_init'
413	END IF
414
415	ALLOCATE (surfheat_incr(kjpindex),stat=ier)
416	IF (ier.NE.0) THEN
417	WRITE (numout,*) ' error in surfheat_incr allocation. We STOP. We need ',kjpindex,' words = '&
418	& , kjpindex
419	STOP 'thermosoil_init'
420	END IF
421
422	ALLOCATE (coldcont_incr(kjpindex),stat=ier)
423	IF (ier.NE.0) THEN
424	WRITE (numout,*) ' error in coldcont_incr allocation. We STOP. We need ',kjpindex,' words = '&
425	& , kjpindex
426	STOP 'thermosoil_init'
427	END IF
428
429	ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier)
430	IF (ier.NE.0) THEN
431	WRITE (numout,*) ' error in pcapa_en allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
432	& , kjpindex*ngrnd
433	STOP 'thermosoil_init'
434	END IF
435
436	ALLOCATE (ptn_beg(kjpindex,ngrnd),stat=ier)
437	IF (ier.NE.0) THEN
438	WRITE (numout,*) ' error in ptn_beg allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
439	& , kjpindex*ngrnd
440	STOP 'thermosoil_init'
441	END IF
442
443	ALLOCATE (temp_sol_beg(kjpindex),stat=ier)
444	IF (ier.NE.0) THEN
445	WRITE (numout,*) ' error in temp_sol_beg allocation. We STOP. We need ',kjpindex,' words = '&
446	& , kjpindex
447	STOP 'thermosoil_init'
448	END IF
449
450	ALLOCATE (wetdiag(kjpindex,ngrnd),stat=ier)
451	IF (ier.NE.0) THEN
452	WRITE (numout,*) ' error in wetdiag allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '&
453	& , kjpindex*ngrnd
454	STOP 'thermosoil_init'
455	END IF
456	!
457	! open restart input file done by sechiba_init
458	! and read data from restart input file for THERMOSOIL process
459
460	IF (ldrestart_read) THEN
461	IF (long_print) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables'
462
463	var_name= 'ptn'
464	CALL ioconf_setatt('UNITS', 'K')
465	CALL ioconf_setatt('LONG_NAME','Soil Temperature profile')
466	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., ptn, "gather", nbp_glo, index_g)
467	!
468	! change restart If they were not found in the restart file
469	!
470	!Config Key = THERMOSOIL_TPRO
471	!Config Desc = Initial soil temperature profile if not found in restart
472	!Config Def = 280.
473	!Config Help = The initial value of the temperature profile in the soil if
474	!Config its value is not found in the restart file. This should only
475	!Config be used if the model is started without a restart file. Here
476	!Config we only require one value as we will assume a constant
477	!Config throughout the column.
478	!
479	CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std)
480
481	ENDIF
482
483	IF (long_print) WRITE (numout,*) ' thermosoil_init done '
484
485	END SUBROUTINE thermosoil_init
486
487	!! Function for distributing the levels
488	!!
489	SUBROUTINE thermosoil_clear()
490
491	l_first_thermosoil=.TRUE.
492
493	IF ( ALLOCATED (ptn)) DEALLOCATE (ptn)
494	IF ( ALLOCATED (z1)) DEALLOCATE (z1)
495	IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd)
496	IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd)
497	IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa)
498	IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa)
499	IF ( ALLOCATED (zdz1)) DEALLOCATE (zdz1)
500	IF ( ALLOCATED (zdz2)) DEALLOCATE (zdz2)
501	IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en)
502	IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg)
503	IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg)
504	IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr)
505	IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr)
506	IF ( ALLOCATED (wetdiag)) DEALLOCATE (wetdiag)
507
508	END SUBROUTINE thermosoil_clear
509	!!
510	!!
511	FUNCTION fz(rk) RESULT (fz_result)
512
513	! interface
514	! input value
515	REAL(r_std), INTENT(in) :: rk
516	! output value
517	REAL(r_std) :: fz_result
518
519	fz_result = fz1 * (zalph ** rk - un) / (zalph - un)
520
521	END FUNCTION fz
522
523	!! Thermosoil's variables initialisation
524	!!
525	SUBROUTINE thermosoil_var_init(kjpindex, zz, zz_coef, dz1, dz2, pkappa, pcapa, pcapa_en, shumdiag, stempdiag)
526
527	! interface description
528	! input scalar
529	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
530	! input fields
531	REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostic humidity profile
532	! output fields
533	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz !!
534	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz_coef
535	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz1 !!
536	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz2 !! tailles des couches
537	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa !!
538	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa_en
539	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pkappa !!
540	REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostic temp profile
541
542	! local declaration
543	INTEGER(i_std) :: ier, ji, jg
544	REAL(r_std) :: sum
545
546	!
547	! 0. initialisation
548	!
549	cstgrnd=SQRT(one_day / pi)
550	lskin = SQRT(so_cond / so_capa * one_day / pi)
551	fz1 = 0.3_r_std * cstgrnd
552	zalph = deux
553	!
554	! 1. Computing the depth of the Temperature level, using a
555	! non dimentional variable x = z/lskin, lskin beeing
556	! the skin depth
557	!
558
559	DO jg=1,ngrnd
560	!!! This needs to be solved soon. Either we allow CPP options in SECHIBA or the VPP
561	!!! fixes its compiler !
562	!!!#ifdef VPP5000
563	dz2(jg) = fz(REAL(jg,r_std)-undemi+undemi) - fz(REAL(jg-1,r_std)-undemi+undemi)
564	!!!#else
565	!!! dz2(jg) = fz(REAL(jg,r_std)) - fz(REAL(jg-1,r_std))
566	!!!#endif
567	ENDDO
568	!
569	! 1.2 The undimentional depth is computed.
570	! ------------------------------------
571	DO jg=1,ngrnd
572	zz(jg) = fz(REAL(jg,r_std) - undemi)
573	zz_coef(jg) = fz(REAL(jg,r_std)-undemi+undemi)
574	ENDDO
575	!
576	! 1.3 Converting to meters.
577	! --------------------
578	DO jg=1,ngrnd
579	zz(jg) = zz(jg) / cstgrnd * lskin
580	zz_coef(jg) = zz_coef(jg) / cstgrnd * lskin
581	dz2(jg) = dz2(jg) / cstgrnd * lskin
582	ENDDO
583	!
584	! 1.4 Computing some usefull constants.
585	! --------------------------------
586	DO jg=1,ngrnd-1
587	dz1(jg) = un / (zz(jg+1) - zz(jg))
588	ENDDO
589	lambda = zz(1) * dz1(1)
590	!
591	! 1.5 Get the wetness profice on this grid
592	! ------------------------------------
593	!
594	CALL thermosoil_humlev(kjpindex, shumdiag)
595	!
596	! 1.6 Thermal conductivity at all levels.
597	! ----------------------------------
598	DO jg = 1,ngrnd
599	DO ji = 1,kjpindex
600	pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry)
601	pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
602	pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
603	ENDDO
604	ENDDO
605	!
606	! 2. Diagnostics.
607	! -----------
608
609	WRITE (numout,*) 'diagnostic des niveaux dans le sol'
610	WRITE (numout,*) 'niveaux intermediaires et pleins'
611	sum = zero
612	DO jg=1,ngrnd
613	sum = sum + dz2(jg)
614	WRITE (numout,*) zz(jg),sum
615	ENDDO
616
617	!
618	! 3. Compute the diagnostic temperature profile
619	!
620
621	CALL thermosoil_diaglev(kjpindex, stempdiag)
622	!
623	!
624
625	IF (long_print) WRITE (numout,*) ' thermosoil_var_init done '
626
627	END SUBROUTINE thermosoil_var_init
628
629	!! Computation of cgrnd and dgrnd coefficient at the t time-step.
630	!!
631	!! Needs previous time step values.
632	!!
633	SUBROUTINE thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,&
634	& zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa)
635
636	! interface description
637	! input scalar
638	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
639	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
640	! input fields
641	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !!
642	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !!
643	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: zz !!
644	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz1 !!
645	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz2 !!
646	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (in) :: ptn
647	! output fields
648	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !!
649	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !!
650	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: z1 !!
651	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa !!
652	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa_en !!
653	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pkappa !!
654	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !!
655	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !!
656	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: zdz1 !!
657	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: zdz2 !!
658
659	! local declaration
660	INTEGER(i_std) :: ji, jg
661	REAL(r_std), DIMENSION(kjpindex) :: snow_h, zx1, zx2
662
663	!
664	! objet: computation of cgrnd and dgrnd coefficient at the t time-step.
665	! ------
666	!
667	! ---------------------------------------------------------------
668	! Computation of the Cgrd and Dgrd coefficient for the next step:
669	! ---------------------------------------------------------------
670	!
671	DO ji = 1,kjpindex
672	snow_h(ji) = snow(ji) / sn_dens
673	!
674	! Traitement special pour le premiere couche
675	!
676	IF ( snow_h(ji) .GT. zz_coef(1) ) THEN
677	pcapa(ji,1) = sn_capa
678	pcapa_en(ji,1) = sn_capa
679	pkappa(ji,1) = sn_cond
680	ELSE IF ( snow_h(ji) .GT. zero ) THEN
681	pcapa_en(ji,1) = sn_capa
682	zx1(ji) = snow_h(ji) / zz_coef(1)
683	zx2(ji) = ( zz_coef(1) - snow_h(ji)) / zz_coef(1)
684	pcapa(ji,1) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet
685	pkappa(ji,1) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet )
686	ELSE
687	pcapa(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry)
688	pkappa(ji,1) = so_cond_dry + wetdiag(ji,1)*(so_cond_wet - so_cond_dry)
689	pcapa_en(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry)
690	ENDIF
691	!
692	DO jg = 2, ngrnd - 2
693	IF ( snow_h(ji) .GT. zz_coef(jg) ) THEN
694	pcapa(ji,jg) = sn_capa
695	pkappa(ji,jg) = sn_cond
696	pcapa_en(ji,jg) = sn_capa
697	ELSE IF ( snow_h(ji) .GT. zz_coef(jg-1) ) THEN
698	zx1(ji) = (snow_h(ji) - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1))
699	zx2(ji) = ( zz_coef(jg) - snow_h(ji)) / (zz_coef(jg) - zz_coef(jg-1))
700	pcapa(ji,jg) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet
701	pkappa(ji,jg) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet )
702	pcapa_en(ji,jg) = sn_capa
703	ELSE
704	pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
705	pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry)
706	pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
707	ENDIF
708	ENDDO
709	!
710	!
711	ENDDO
712	!
713	DO jg=1,ngrnd
714	DO ji=1,kjpindex
715	zdz2(ji,jg)=pcapa(ji,jg) * dz2(jg)/dtradia
716	ENDDO
717	ENDDO
718	!
719	DO jg=1,ngrnd-1
720	DO ji=1,kjpindex
721	zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg)
722	ENDDO
723	ENDDO
724	!
725	DO ji = 1,kjpindex
726	z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1)
727	cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji)
728	dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji)
729	ENDDO
730
731	DO jg = ngrnd-1,2,-1
732	DO ji = 1,kjpindex
733	z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg)))
734	cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji)
735	dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji)
736	ENDDO
737	ENDDO
738
739	! ---------------------------------------------------------
740	! computation of the surface diffusive flux from ground and
741	! calorific capacity of the ground:
742	! ---------------------------------------------------------
743
744	DO ji = 1,kjpindex
745	soilflx(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1))
746	soilcap(ji) = (zdz2(ji,1) * dtradia + dtradia * (un - dgrnd(ji,1)) * zdz1(ji,1))
747	z1(ji) = lambda * (un - dgrnd(ji,1)) + un
748	soilcap(ji) = soilcap(ji) / z1(ji)
749	soilflx(ji) = soilflx(ji) + &
750	& soilcap(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dtradia
751	ENDDO
752
753	IF (long_print) WRITE (numout,*) ' thermosoil_coef done '
754
755	END SUBROUTINE thermosoil_coef
756
757	!! Computation of : the ground temperature evolution
758	!!
759	SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, ptn, stempdiag)
760
761	! interface description
762	! input scalar
763	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
764	! input fields
765	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
766	! modified fields
767	REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (inout) :: ptn !! Different levels soil temperature
768	! output fields
769	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostoc profile
770
771	! local declaration
772	INTEGER(i_std) :: ji, jg
773	!
774	! objet: computation of : the ground temperature evolution
775	! ------
776	!
777	! Method: implicit time integration
778	! -------
779	! Consecutives ground temperatures are related by:
780	! T(k+1) = C(k) + D(k)*T(k) (1)
781	! the coefficients C and D are computed at the t-dt time-step.
782	! Routine structure:
783	! -new temperatures are computed using (1)
784	!
785	!
786	! surface temperature
787	DO ji = 1,kjpindex
788	ptn(ji,1) = (lambda * cgrnd(ji,1) + temp_sol_new(ji)) / (lambda * (un - dgrnd(ji,1)) + un)
789	ENDDO
790
791	! other temperatures
792	DO jg = 1,ngrnd-1
793	DO ji = 1,kjpindex
794	ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg)
795	ENDDO
796	ENDDO
797
798	CALL thermosoil_diaglev(kjpindex, stempdiag)
799
800	IF (long_print) WRITE (numout,*) ' thermosoil_profile done '
801
802	END SUBROUTINE thermosoil_profile
803	!!
804	!!
805	!! Diagnostic soil temperature profile on new levels
806	!!
807	!!
808	SUBROUTINE thermosoil_diaglev(kjpindex, stempdiag)
809	! interface description
810	! input scalar
811	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
812	! input fields
813	!
814	! modified fields
815	!
816	! output fields
817	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostoc profile
818	!
819	! local variable
820	!
821	INTEGER(i_std) :: ji, jd, jg
822	REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog
823	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfact
824	!
825	LOGICAL, PARAMETER :: check=.FALSE.
826	!
827	!
828	IF ( .NOT. ALLOCATED(intfact)) THEN
829	!
830	ALLOCATE(intfact(nbdl, ngrnd))
831	!
832	prev_diag = zero
833	DO jd = 1, nbdl
834	lev_diag = diaglev(jd)
835	prev_prog = zero
836	DO jg = 1, ngrnd
837	IF ( jg == ngrnd .AND. (prev_prog + dz2(jg)) < lev_diag ) THEN
838	!! Just make sure we cover the deepest layers
839	lev_prog = lev_diag
840	ELSE
841	lev_prog = prev_prog + dz2(jg)
842	ENDIF
843	intfact(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
844	prev_prog = lev_prog
845	ENDDO
846	prev_diag = lev_diag
847	ENDDO
848	!
849	IF ( check ) THEN
850	WRITE(numout,*) 'thermosoil_diagev -- thermosoil_diaglev -- thermosoil_diaglev --'
851	DO jd = 1, nbdl
852	WRITE(numout,*) jd, '-', intfact(jd,1:ngrnd)
853	ENDDO
854	WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM"
855	DO jd = 1, nbdl
856	WRITE(numout,*) jd, '-', SUM(intfact(jd,1:ngrnd))
857	ENDDO
858	WRITE(numout,*) 'thermosoil_diaglev -- thermosoil_diaglev -- thermosoil_diaglev --'
859	ENDIF
860	!
861	ENDIF
862
863	stempdiag(:,:) = zero
864	DO jg = 1, ngrnd
865	DO jd = 1, nbdl
866	DO ji = 1, kjpindex
867	stempdiag(ji,jd) = stempdiag(ji,jd) + ptn(ji,jg)*intfact(jd,jg)
868	ENDDO
869	ENDDO
870	ENDDO
871
872	END SUBROUTINE thermosoil_diaglev
873	!!
874	!!
875	!! Put soil wetness on the temperature levels
876	!!
877	!!
878	SUBROUTINE thermosoil_humlev(kjpindex, shumdiag)
879	! interface description
880	! input scalar
881	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
882	! input fields
883	!
884	! modified fields
885	!
886	! output fields
887	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostoc profile
888	!
889	! local variable
890	!
891	INTEGER(i_std) :: ji, jd, jg
892	REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog
893	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfactw
894	!
895	LOGICAL, PARAMETER :: check=.FALSE.
896	!
897	!
898	IF ( .NOT. ALLOCATED(intfactw)) THEN
899	!
900	ALLOCATE(intfactw(ngrnd, nbdl))
901	!
902	prev_diag = zero
903	DO jd = 1, ngrnd
904	lev_diag = prev_diag + dz2(jd)
905	prev_prog = zero
906	DO jg = 1, nbdl
907	IF ( jg == nbdl .AND. diaglev(jg) < lev_diag ) THEN
908	!! Just make sure we cover the deepest layers
909	lev_prog = lev_diag
910	ELSE
911	lev_prog = diaglev(jg)
912	ENDIF
913	intfactw(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
914	prev_prog = lev_prog
915	ENDDO
916	prev_diag = lev_diag
917	ENDDO
918	!
919	IF ( check ) THEN
920	WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --'
921	DO jd = 1, ngrnd
922	WRITE(numout,*) jd, '-', intfactw(jd,1:nbdl)
923	ENDDO
924	WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM"
925	DO jd = 1, ngrnd
926	WRITE(numout,*) jd, '-', SUM(intfactw(jd,1:nbdl))
927	ENDDO
928	WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --'
929	ENDIF
930	!
931	ENDIF
932
933	wetdiag(:,:) = zero
934	DO jg = 1, nbdl
935	DO jd = 1, ngrnd
936	DO ji = 1, kjpindex
937	wetdiag(ji,jd) = wetdiag(ji,jd) + shumdiag(ji,jg)*intfactw(jd,jg)
938	ENDDO
939	ENDDO
940	ENDDO
941
942	END SUBROUTINE thermosoil_humlev
943
944	SUBROUTINE thermosoil_energy(kjpindex, temp_sol_new, soilcap, first_call)
945	! interface description
946	! input scalar
947	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
948	LOGICAL, INTENT (in) :: first_call !!
949	! input fields
950	!
951	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
952	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil capacity
953	!
954	! modified fields
955	!
956	! output fields
957	!
958	! local variable
959	!
960	INTEGER(i_std) :: ji, jg
961	!
962	!
963	IF (first_call) THEN
964
965	DO ji = 1, kjpindex
966	surfheat_incr(ji) = zero
967	coldcont_incr(ji) = zero
968	temp_sol_beg(ji) = temp_sol_new(ji)
969	!
970	DO jg = 1, ngrnd
971	ptn_beg(ji,jg) = ptn(ji,jg)
972	ENDDO
973	!
974	ENDDO
975
976	RETURN
977
978	ENDIF
979
980	DO ji = 1, kjpindex
981	surfheat_incr(ji) = zero
982	coldcont_incr(ji) = zero
983	ENDDO
984	!
985	! Sum up the energy content of all layers in the soil.
986	!
987	DO ji = 1, kjpindex
988	!
989	IF (pcapa_en(ji,1) .LE. sn_capa) THEN
990	!
991	! Verify the energy conservation in the surface layer
992	!
993	coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
994	surfheat_incr(ji) = zero
995	ELSE
996	!
997	! Verify the energy conservation in the surface layer
998	!
999	surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
1000	coldcont_incr(ji) = zero
1001	ENDIF
1002	ENDDO
1003
1004	ptn_beg(:,:) = ptn(:,:)
1005	temp_sol_beg(:) = temp_sol_new(:)
1006
1007	END SUBROUTINE thermosoil_energy
1008
1009	END MODULE thermosoil

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

Download in other formats: