Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_sechiba/thermosoil.f90 @ 2041

Last change on this file since 2041 was 570, checked in by didier.solyga, 13 years ago
Add labels in thermosoil. Delete old comments
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
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$, $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 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 If = OK_SECHIBA
474	!Config Help = The initial value of the temperature profile in the soil if
475	!Config its value is not found in the restart file. This should only
476	!Config be used if the model is started without a restart file. Here
477	!Config we only require one value as we will assume a constant
478	!Config throughout the column.
479	!Config Units = Kelvin [K]
480	!
481	CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std)
482
483	ENDIF
484
485	IF (long_print) WRITE (numout,*) ' thermosoil_init done '
486
487	END SUBROUTINE thermosoil_init
488
489	!! Function for distributing the levels
490	!!
491	SUBROUTINE thermosoil_clear()
492
493	l_first_thermosoil=.TRUE.
494
495	IF ( ALLOCATED (ptn)) DEALLOCATE (ptn)
496	IF ( ALLOCATED (z1)) DEALLOCATE (z1)
497	IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd)
498	IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd)
499	IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa)
500	IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa)
501	IF ( ALLOCATED (zdz1)) DEALLOCATE (zdz1)
502	IF ( ALLOCATED (zdz2)) DEALLOCATE (zdz2)
503	IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en)
504	IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg)
505	IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg)
506	IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr)
507	IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr)
508	IF ( ALLOCATED (wetdiag)) DEALLOCATE (wetdiag)
509
510	END SUBROUTINE thermosoil_clear
511	!!
512	!!
513	FUNCTION fz(rk) RESULT (fz_result)
514
515	! interface
516	! input value
517	REAL(r_std), INTENT(in) :: rk
518	! output value
519	REAL(r_std) :: fz_result
520
521	fz_result = fz1 * (zalph ** rk - un) / (zalph - un)
522
523	END FUNCTION fz
524
525	!! Thermosoil's variables initialisation
526	!!
527	SUBROUTINE thermosoil_var_init(kjpindex, zz, zz_coef, dz1, dz2, pkappa, pcapa, pcapa_en, shumdiag, stempdiag)
528
529	! interface description
530	! input scalar
531	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
532	! input fields
533	REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostic humidity profile
534	! output fields
535	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz !!
536	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz_coef
537	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz1 !!
538	REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz2 !! tailles des couches
539	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa !!
540	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa_en
541	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pkappa !!
542	REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostic temp profile
543
544	! local declaration
545	INTEGER(i_std) :: ier, ji, jg
546	REAL(r_std) :: sum
547
548	!
549	! 0. initialisation
550	!
551	cstgrnd=SQRT(one_day / pi)
552	lskin = SQRT(so_cond / so_capa * one_day / pi)
553	fz1 = 0.3_r_std * cstgrnd
554	zalph = deux
555	!
556	! 1. Computing the depth of the Temperature level, using a
557	! non dimentional variable x = z/lskin, lskin beeing
558	! the skin depth
559	!
560
561	DO jg=1,ngrnd
562	!!! This needs to be solved soon. Either we allow CPP options in SECHIBA or the VPP
563	!!! fixes its compiler !
564	!!!#ifdef VPP5000
565	dz2(jg) = fz(REAL(jg,r_std)-undemi+undemi) - fz(REAL(jg-1,r_std)-undemi+undemi)
566	!!!#else
567	!!! dz2(jg) = fz(REAL(jg,r_std)) - fz(REAL(jg-1,r_std))
568	!!!#endif
569	ENDDO
570	!
571	! 1.2 The undimentional depth is computed.
572	! ------------------------------------
573	DO jg=1,ngrnd
574	zz(jg) = fz(REAL(jg,r_std) - undemi)
575	zz_coef(jg) = fz(REAL(jg,r_std)-undemi+undemi)
576	ENDDO
577	!
578	! 1.3 Converting to meters.
579	! --------------------
580	DO jg=1,ngrnd
581	zz(jg) = zz(jg) / cstgrnd * lskin
582	zz_coef(jg) = zz_coef(jg) / cstgrnd * lskin
583	dz2(jg) = dz2(jg) / cstgrnd * lskin
584	ENDDO
585	!
586	! 1.4 Computing some usefull constants.
587	! --------------------------------
588	DO jg=1,ngrnd-1
589	dz1(jg) = un / (zz(jg+1) - zz(jg))
590	ENDDO
591	lambda = zz(1) * dz1(1)
592	!
593	! 1.5 Get the wetness profice on this grid
594	! ------------------------------------
595	!
596	CALL thermosoil_humlev(kjpindex, shumdiag)
597	!
598	! 1.6 Thermal conductivity at all levels.
599	! ----------------------------------
600	DO jg = 1,ngrnd
601	DO ji = 1,kjpindex
602	pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry)
603	pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
604	pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
605	ENDDO
606	ENDDO
607	!
608	! 2. Diagnostics.
609	! -----------
610
611	WRITE (numout,*) 'diagnostic des niveaux dans le sol'
612	WRITE (numout,*) 'niveaux intermediaires et pleins'
613	sum = zero
614	DO jg=1,ngrnd
615	sum = sum + dz2(jg)
616	WRITE (numout,*) zz(jg),sum
617	ENDDO
618
619	!
620	! 3. Compute the diagnostic temperature profile
621	!
622
623	CALL thermosoil_diaglev(kjpindex, stempdiag)
624	!
625	!
626
627	IF (long_print) WRITE (numout,*) ' thermosoil_var_init done '
628
629	END SUBROUTINE thermosoil_var_init
630
631	!! Computation of cgrnd and dgrnd coefficient at the t time-step.
632	!!
633	!! Needs previous time step values.
634	!!
635	SUBROUTINE thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,&
636	& zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa)
637
638	! interface description
639	! input scalar
640	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
641	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
642	! input fields
643	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !!
644	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !!
645	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: zz !!
646	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz1 !!
647	REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz2 !!
648	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (in) :: ptn
649	! output fields
650	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !!
651	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !!
652	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: z1 !!
653	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa !!
654	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa_en !!
655	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pkappa !!
656	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !!
657	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !!
658	REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: zdz1 !!
659	REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: zdz2 !!
660
661	! local declaration
662	INTEGER(i_std) :: ji, jg
663	REAL(r_std), DIMENSION(kjpindex) :: snow_h, zx1, zx2
664
665	!
666	! objet: computation of cgrnd and dgrnd coefficient at the t time-step.
667	! ------
668	!
669	! ---------------------------------------------------------------
670	! Computation of the Cgrd and Dgrd coefficient for the next step:
671	! ---------------------------------------------------------------
672	!
673	DO ji = 1,kjpindex
674	snow_h(ji) = snow(ji) / sn_dens
675	!
676	! Traitement special pour le premiere couche
677	!
678	IF ( snow_h(ji) .GT. zz_coef(1) ) THEN
679	pcapa(ji,1) = sn_capa
680	pcapa_en(ji,1) = sn_capa
681	pkappa(ji,1) = sn_cond
682	ELSE IF ( snow_h(ji) .GT. zero ) THEN
683	pcapa_en(ji,1) = sn_capa
684	zx1(ji) = snow_h(ji) / zz_coef(1)
685	zx2(ji) = ( zz_coef(1) - snow_h(ji)) / zz_coef(1)
686	pcapa(ji,1) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet
687	pkappa(ji,1) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet )
688	ELSE
689	pcapa(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry)
690	pkappa(ji,1) = so_cond_dry + wetdiag(ji,1)*(so_cond_wet - so_cond_dry)
691	pcapa_en(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry)
692	ENDIF
693	!
694	DO jg = 2, ngrnd - 2
695	IF ( snow_h(ji) .GT. zz_coef(jg) ) THEN
696	pcapa(ji,jg) = sn_capa
697	pkappa(ji,jg) = sn_cond
698	pcapa_en(ji,jg) = sn_capa
699	ELSE IF ( snow_h(ji) .GT. zz_coef(jg-1) ) THEN
700	zx1(ji) = (snow_h(ji) - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1))
701	zx2(ji) = ( zz_coef(jg) - snow_h(ji)) / (zz_coef(jg) - zz_coef(jg-1))
702	pcapa(ji,jg) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet
703	pkappa(ji,jg) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet )
704	pcapa_en(ji,jg) = sn_capa
705	ELSE
706	pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
707	pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry)
708	pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry)
709	ENDIF
710	ENDDO
711	!
712	!
713	ENDDO
714	!
715	DO jg=1,ngrnd
716	DO ji=1,kjpindex
717	zdz2(ji,jg)=pcapa(ji,jg) * dz2(jg)/dtradia
718	ENDDO
719	ENDDO
720	!
721	DO jg=1,ngrnd-1
722	DO ji=1,kjpindex
723	zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg)
724	ENDDO
725	ENDDO
726	!
727	DO ji = 1,kjpindex
728	z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1)
729	cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji)
730	dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji)
731	ENDDO
732
733	DO jg = ngrnd-1,2,-1
734	DO ji = 1,kjpindex
735	z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg)))
736	cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji)
737	dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji)
738	ENDDO
739	ENDDO
740
741	! ---------------------------------------------------------
742	! computation of the surface diffusive flux from ground and
743	! calorific capacity of the ground:
744	! ---------------------------------------------------------
745
746	DO ji = 1,kjpindex
747	soilflx(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1))
748	soilcap(ji) = (zdz2(ji,1) * dtradia + dtradia * (un - dgrnd(ji,1)) * zdz1(ji,1))
749	z1(ji) = lambda * (un - dgrnd(ji,1)) + un
750	soilcap(ji) = soilcap(ji) / z1(ji)
751	soilflx(ji) = soilflx(ji) + &
752	& soilcap(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dtradia
753	ENDDO
754
755	IF (long_print) WRITE (numout,*) ' thermosoil_coef done '
756
757	END SUBROUTINE thermosoil_coef
758
759	!! Computation of : the ground temperature evolution
760	!!
761	SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, ptn, stempdiag)
762
763	! interface description
764	! input scalar
765	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
766	! input fields
767	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
768	! modified fields
769	REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (inout) :: ptn !! Different levels soil temperature
770	! output fields
771	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostoc profile
772
773	! local declaration
774	INTEGER(i_std) :: ji, jg
775	!
776	! objet: computation of : the ground temperature evolution
777	! ------
778	!
779	! Method: implicit time integration
780	! -------
781	! Consecutives ground temperatures are related by:
782	! T(k+1) = C(k) + D(k)*T(k) (1)
783	! the coefficients C and D are computed at the t-dt time-step.
784	! Routine structure:
785	! -new temperatures are computed using (1)
786	!
787	!
788	! surface temperature
789	DO ji = 1,kjpindex
790	ptn(ji,1) = (lambda * cgrnd(ji,1) + temp_sol_new(ji)) / (lambda * (un - dgrnd(ji,1)) + un)
791	ENDDO
792
793	! other temperatures
794	DO jg = 1,ngrnd-1
795	DO ji = 1,kjpindex
796	ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg)
797	ENDDO
798	ENDDO
799
800	CALL thermosoil_diaglev(kjpindex, stempdiag)
801
802	IF (long_print) WRITE (numout,*) ' thermosoil_profile done '
803
804	END SUBROUTINE thermosoil_profile
805	!!
806	!!
807	!! Diagnostic soil temperature profile on new levels
808	!!
809	!!
810	SUBROUTINE thermosoil_diaglev(kjpindex, stempdiag)
811	! interface description
812	! input scalar
813	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
814	! input fields
815	!
816	! modified fields
817	!
818	! output fields
819	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostoc profile
820	!
821	! local variable
822	!
823	INTEGER(i_std) :: ji, jd, jg
824	REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog
825	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfact
826	!
827	LOGICAL, PARAMETER :: check=.FALSE.
828	!
829	!
830	IF ( .NOT. ALLOCATED(intfact)) THEN
831	!
832	ALLOCATE(intfact(nbdl, ngrnd))
833	!
834	prev_diag = zero
835	DO jd = 1, nbdl
836	lev_diag = diaglev(jd)
837	prev_prog = zero
838	DO jg = 1, ngrnd
839	IF ( jg == ngrnd .AND. (prev_prog + dz2(jg)) < lev_diag ) THEN
840	!! Just make sure we cover the deepest layers
841	lev_prog = lev_diag
842	ELSE
843	lev_prog = prev_prog + dz2(jg)
844	ENDIF
845	intfact(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
846	prev_prog = lev_prog
847	ENDDO
848	prev_diag = lev_diag
849	ENDDO
850	!
851	IF ( check ) THEN
852	WRITE(numout,*) 'thermosoil_diagev -- thermosoil_diaglev -- thermosoil_diaglev --'
853	DO jd = 1, nbdl
854	WRITE(numout,*) jd, '-', intfact(jd,1:ngrnd)
855	ENDDO
856	WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM"
857	DO jd = 1, nbdl
858	WRITE(numout,*) jd, '-', SUM(intfact(jd,1:ngrnd))
859	ENDDO
860	WRITE(numout,*) 'thermosoil_diaglev -- thermosoil_diaglev -- thermosoil_diaglev --'
861	ENDIF
862	!
863	ENDIF
864
865	stempdiag(:,:) = zero
866	DO jg = 1, ngrnd
867	DO jd = 1, nbdl
868	DO ji = 1, kjpindex
869	stempdiag(ji,jd) = stempdiag(ji,jd) + ptn(ji,jg)*intfact(jd,jg)
870	ENDDO
871	ENDDO
872	ENDDO
873
874	END SUBROUTINE thermosoil_diaglev
875	!!
876	!!
877	!! Put soil wetness on the temperature levels
878	!!
879	!!
880	SUBROUTINE thermosoil_humlev(kjpindex, shumdiag)
881	! interface description
882	! input scalar
883	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
884	! input fields
885	!
886	! modified fields
887	!
888	! output fields
889	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostoc profile
890	!
891	! local variable
892	!
893	INTEGER(i_std) :: ji, jd, jg
894	REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog
895	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfactw
896	!
897	LOGICAL, PARAMETER :: check=.FALSE.
898	!
899	!
900	IF ( .NOT. ALLOCATED(intfactw)) THEN
901	!
902	ALLOCATE(intfactw(ngrnd, nbdl))
903	!
904	prev_diag = zero
905	DO jd = 1, ngrnd
906	lev_diag = prev_diag + dz2(jd)
907	prev_prog = zero
908	DO jg = 1, nbdl
909	IF ( jg == nbdl .AND. diaglev(jg) < lev_diag ) THEN
910	!! Just make sure we cover the deepest layers
911	lev_prog = lev_diag
912	ELSE
913	lev_prog = diaglev(jg)
914	ENDIF
915	intfactw(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
916	prev_prog = lev_prog
917	ENDDO
918	prev_diag = lev_diag
919	ENDDO
920	!
921	IF ( check ) THEN
922	WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --'
923	DO jd = 1, ngrnd
924	WRITE(numout,*) jd, '-', intfactw(jd,1:nbdl)
925	ENDDO
926	WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM"
927	DO jd = 1, ngrnd
928	WRITE(numout,*) jd, '-', SUM(intfactw(jd,1:nbdl))
929	ENDDO
930	WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --'
931	ENDIF
932	!
933	ENDIF
934
935	wetdiag(:,:) = zero
936	DO jg = 1, nbdl
937	DO jd = 1, ngrnd
938	DO ji = 1, kjpindex
939	wetdiag(ji,jd) = wetdiag(ji,jd) + shumdiag(ji,jg)*intfactw(jd,jg)
940	ENDDO
941	ENDDO
942	ENDDO
943
944	END SUBROUTINE thermosoil_humlev
945
946	SUBROUTINE thermosoil_energy(kjpindex, temp_sol_new, soilcap, first_call)
947	! interface description
948	! input scalar
949	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
950	LOGICAL, INTENT (in) :: first_call !!
951	! input fields
952	!
953	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
954	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil capacity
955	!
956	! modified fields
957	!
958	! output fields
959	!
960	! local variable
961	!
962	INTEGER(i_std) :: ji, jg
963	!
964	!
965	IF (first_call) THEN
966
967	DO ji = 1, kjpindex
968	surfheat_incr(ji) = zero
969	coldcont_incr(ji) = zero
970	temp_sol_beg(ji) = temp_sol_new(ji)
971	!
972	DO jg = 1, ngrnd
973	ptn_beg(ji,jg) = ptn(ji,jg)
974	ENDDO
975	!
976	ENDDO
977
978	RETURN
979
980	ENDIF
981
982	DO ji = 1, kjpindex
983	surfheat_incr(ji) = zero
984	coldcont_incr(ji) = zero
985	ENDDO
986	!
987	! Sum up the energy content of all layers in the soil.
988	!
989	DO ji = 1, kjpindex
990	!
991	IF (pcapa_en(ji,1) .LE. sn_capa) THEN
992	!
993	! Verify the energy conservation in the surface layer
994	!
995	coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
996	surfheat_incr(ji) = zero
997	ELSE
998	!
999	! Verify the energy conservation in the surface layer
1000	!
1001	surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
1002	coldcont_incr(ji) = zero
1003	ENDIF
1004	ENDDO
1005
1006	ptn_beg(:,:) = ptn(:,:)
1007	temp_sol_beg(:) = temp_sol_new(:)
1008
1009	END SUBROUTINE thermosoil_energy
1010
1011	END MODULE thermosoil

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

Download in other formats: