Context Navigation

source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE/src_sechiba/thermosoil.f90

Last change on this file was 45, checked in by mmaipsl, 14 years ago

MM: Tests with lf95 compiler : correct f95 strict norm problems.

There is no change in numerical result after these commits.

Property svn:keywords set to Revision Date HeadURL Date Author Revision

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$, $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 constantes_soil
23	USE sechiba_io
24	USE grid
25	USE parallel
26
27	IMPLICIT NONE
28
29	! public routines :
30	! thermosoil_main
31	PRIVATE
32	PUBLIC :: thermosoil_main,thermosoil_clear
33
34	!
35	! variables used inside thermosoil module : declaration and initialisation
36	!
37	LOGICAL, SAVE :: l_first_thermosoil=.TRUE. !! Initialisation has to be done one time
38
39	CHARACTER(LEN=80) , SAVE :: var_name !! To store variables names for I/O
40	REAL(r_std), SAVE :: lambda, cstgrnd, lskin, fz1, zalph
41
42	! two dimensions array allocated, computed, saved and got in thermosoil module
43
44	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn !! Different levels soil temperature
45
46	! one dimension array allocated, computed and used in thermosoil module exclusively
47
48	REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz !!
49	REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz_coef
50	REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz1 !!
51	REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz2 !!
52
53	! one dimension array allocated, computed and used in thermosoil module exclusively
54
55	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z1 !!
56
57	! two dimensions arrays allocated, computed and used in thermosoil module exclusively
58
59	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !!
60	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !!
61	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !!
62	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !!
63	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz1 !!
64	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz2 !!
65	!
66	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! Capacity used for energy_incr calculation
67	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg !! Temperature at the beginning of the time step
68	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at the beginning of the timestep
69	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat
70	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow cold content
71	!!
72	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: wetdiag !! Soil weetness on the thermodynamical levels
73	!!
74	CONTAINS
75
76	!! Main routine for thermosoil module.
77	!!
78	!! Algorithm:
79	!! - call thermosoil_var_init to initialise variables
80	!! - call thermosoil_coef for coefficient
81	!! - call thermosoil_profile for soil profiling
82	!!
83	!! @call thermosoil_var_init
84	!! @call thermosoil_coef
85	!! @call thermosoil_profile
86	!!
87	SUBROUTINE thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexgrnd, &
88	& temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id)
89
90	! interface description
91	! input scalar
92	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
93	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
94	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and history file identifier
95	INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier
96	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
97	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
98	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
99	! input fields
100	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
101	INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in) :: indexgrnd !! Indeces of the points on the 3D map
102	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
103	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow quantity
104	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostic of relative humidity
105	! output fields
106	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! Soil capacity
107	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx
108	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (inout):: stempdiag !! diagnostic temp profile
109
110	REAL(r_std),DIMENSION (kjpindex,ngrnd) :: temp
111	REAL(r_std),DIMENSION (kjpindex,ngrnd-1) :: temp1
112	REAL(r_std),DIMENSION (kjpindex) :: temp2
113	!
114	! do initialisation
115	!
116	IF (l_first_thermosoil) THEN
117
118	IF (long_print) WRITE (numout,*) ' l_first_thermosoil : call thermosoil_init '
119
120	!
121	! do needed allocation
122	!
123	CALL thermosoil_init (kjit, ldrestart_read, kjpindex, index, rest_id)
124
125	!
126	! computes some physical constantes and array depending soil level depth
127	!
128	CALL thermosoil_var_init (kjpindex, zz, zz_coef, dz1, dz2, pkappa, pcapa, pcapa_en, shumdiag, stempdiag)
129
130	!
131	! computes cgrd and dgrd coefficient from previous time step (restart)
132	!
133	CALL thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,&
134	& zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa)
135
136	CALL thermosoil_energy (kjpindex, temp_sol_new, soilcap, .TRUE.)
137
138	IF (ldrestart_read) THEN
139	IF (long_print) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables'
140
141	var_name= 'cgrnd'
142	CALL ioconf_setatt('UNITS', '-')
143	CALL ioconf_setatt('LONG_NAME','Cgrnd coefficient.')
144	IF ( ok_var(var_name) ) THEN
145	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
146	IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
147	cgrnd(:,:)=temp1(:,:)
148	ENDIF
149	ENDIF
150
151	var_name= 'dgrnd'
152	CALL ioconf_setatt('UNITS', '-')
153	CALL ioconf_setatt('LONG_NAME','Dgrnd coefficient.')
154	IF ( ok_var(var_name) ) THEN
155	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
156	IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
157	dgrnd(:,:)=temp1(:,:)
158	ENDIF
159	ENDIF
160
161	var_name= 'z1'
162	CALL ioconf_setatt('UNITS', '-')
163	CALL ioconf_setatt('LONG_NAME','?.')
164	IF ( ok_var(var_name) ) THEN
165	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g)
166	IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN
167	z1(:)=temp2(:)
168	ENDIF
169	ENDIF
170
171	var_name= 'pcapa'
172	CALL ioconf_setatt('UNITS', '-')
173	CALL ioconf_setatt('LONG_NAME','?.')
174	IF ( ok_var(var_name) ) THEN
175	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
176	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
177	pcapa(:,:)=temp(:,:)
178	ENDIF
179	ENDIF
180
181	var_name= 'pcapa_en'
182	CALL ioconf_setatt('UNITS', '-')
183	CALL ioconf_setatt('LONG_NAME','?.')
184	IF ( ok_var(var_name) ) THEN
185	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
186	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
187	pcapa_en(:,:)=temp(:,:)
188	ENDIF
189	ENDIF
190
191	var_name= 'pkappa'
192	CALL ioconf_setatt('UNITS', '-')
193	CALL ioconf_setatt('LONG_NAME','?.')
194	IF ( ok_var(var_name) ) THEN
195	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
196	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
197	pkappa(:,:)=temp(:,:)
198	ENDIF
199	ENDIF
200
201	var_name= 'zdz1'
202	CALL ioconf_setatt('UNITS', '-')
203	CALL ioconf_setatt('LONG_NAME','?.')
204	IF ( ok_var(var_name) ) THEN
205	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g)
206	IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN
207	zdz1(:,:)=temp1(:,:)
208	ENDIF
209	ENDIF
210
211	var_name= 'zdz2'
212	CALL ioconf_setatt('UNITS', '-')
213	CALL ioconf_setatt('LONG_NAME','?.')
214	IF ( ok_var(var_name) ) THEN
215	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
216	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
217	zdz2(:,:)=temp(:,:)
218	ENDIF
219	ENDIF
220
221	var_name='temp_sol_beg'
222	CALL ioconf_setatt('UNITS', 'K')
223	CALL ioconf_setatt('LONG_NAME','Old Surface temperature')
224	IF ( ok_var(var_name) ) THEN
225	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g)
226	IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN
227	temp_sol_beg(:) = temp2(:)
228	ENDIF
229	ENDIF
230
231	ENDIF
232
233	RETURN
234
235	ENDIF
236
237	!
238	! prepares restart file for the next simulation
239	!
240	IF (ldrestart_write) THEN
241
242	IF (long_print) WRITE (numout,*) ' we have to complete restart file with THERMOSOIL variables'
243
244	var_name= 'ptn'
245	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, ptn, 'scatter', nbp_glo, index_g)
246
247	var_name= 'cgrnd'
248	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g)
249	var_name= 'dgrnd'
250	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g)
251
252	var_name= 'z1'
253	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, z1, 'scatter', nbp_glo, index_g)
254
255	var_name= 'pcapa'
256	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pcapa, 'scatter', nbp_glo, index_g)
257
258	var_name= 'pcapa_en'
259	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pcapa_en, 'scatter', nbp_glo, index_g)
260
261	var_name= 'pkappa'
262	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pkappa, 'scatter', nbp_glo, index_g)
263
264	var_name= 'zdz1'
265	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, zdz1, 'scatter', nbp_glo, index_g)
266
267	var_name= 'zdz2'
268	CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, zdz2, 'scatter', nbp_glo, index_g)
269
270	var_name= 'temp_sol_beg'
271	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol_beg, 'scatter', nbp_glo, index_g)
272
273	var_name= 'soilcap'
274	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g)
275
276	var_name= 'soilflx'
277	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g)
278
279	! read in enerbil
280	var_name= 'temp_sol_new'
281	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol_new, 'scatter', nbp_glo, index_g)
282
283	RETURN
284
285	END IF
286
287	!
288	! Put the soil wetnesss diagnostic on the levels of the soil temprature
289	!
290	CALL thermosoil_humlev(kjpindex, shumdiag)
291
292	!
293	! computes profile with previous cgrd and dgrd coefficient
294	!
295	CALL thermosoil_profile (kjpindex, temp_sol_new, ptn, stempdiag)
296	CALL thermosoil_energy (kjpindex, temp_sol_new, soilcap, .FALSE.)
297	!
298	IF ( .NOT. almaoutput ) THEN
299	CALL histwrite(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
300	ELSE
301	CALL histwrite(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
302	CALL histwrite(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
303	CALL histwrite(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
304	CALL histwrite(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
305	ENDIF
306	IF ( hist2_id > 0 ) THEN
307	IF ( .NOT. almaoutput ) THEN
308	CALL histwrite(hist2_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
309	ELSE
310	CALL histwrite(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
311	CALL histwrite(hist2_id, 'Qg', kjit, soilflx, kjpindex, index)
312	CALL histwrite(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
313	CALL histwrite(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
314	ENDIF
315	ENDIF
316	!
317	! computes cgrd and dgrd coefficient
318	!
319	CALL thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,&
320	& zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa)
321
322	IF (long_print) WRITE (numout,*) ' thermosoil_main done '
323
324	END SUBROUTINE thermosoil_main
325
326	!! Initialisation for thermosoil module.
327	!! - does dynamic allocation for local arrays
328	!! - reads _restart_ file or set initial values to a raisonable value
329	!!
330	SUBROUTINE thermosoil_init(kjit, ldrestart_read, kjpindex, index, rest_id)
331
332	! interface description
333	! input scalar
334	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
335	LOGICAL,INTENT (in) :: ldrestart_read !! Logical for restart file to read
336	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
337	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
338
339	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier
340	! input fields
341	! output fields
342
343	! local declaration
344	INTEGER(i_std) :: ier
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: