Context Navigation

source: tags/ORCHIDEE/src_sechiba/thermosoil.f90 @ 6

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

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

Download in other formats: