Context Navigation

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

Last change on this file since 64 was 64, checked in by didier.solyga, 13 years ago
Import first version of ORCHIDEE_EXT
File size: 37.4 KB

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

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

Download in other formats: