Context Navigation

soil.F @ 3855

Last change on this file since 3855 was 704, checked in by Laurent Fairhead, 18 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 8.2 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil,
5	s pcapcal, pfluxgrd)
6	use dimphy
7	IMPLICIT NONE
8
9	c=======================================================================
10	c
11	c Auteur: Frederic Hourdin 30/01/92
12	c -------
13	c
14	c objet: computation of : the soil temperature evolution
15	c ------ the surfacic heat capacity "Capcal"
16	c the surface conduction flux pcapcal
17	c
18	c
19	c Method: implicit time integration
20	c -------
21	c Consecutive ground temperatures are related by:
22	c T(k+1) = C(k) + D(k)*T(k) (1)
23	c the coefficients C and D are computed at the t-dt time-step.
24	c Routine structure:
25	c 1)new temperatures are computed using (1)
26	c 2)C and D coefficients are computed from the new temperature
27	c profile for the t+dt time-step
28	c 3)the coefficients A and B are computed where the diffusive
29	c fluxes at the t+dt time-step is given by
30	c Fdiff = A + B Ts(t+dt)
31	c or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
32	c with F0 = A + B (Ts(t))
33	c Capcal = B*dt
34	c
35	c Interface:
36	c ----------
37	c
38	c Arguments:
39	c ----------
40	c ptimestep physical timestep (s)
41	c indice sub-surface index
42	c snow(klon,nbsrf) snow
43	c ptsrf(klon) surface temperature at time-step t (K)
44	c ptsoil(klon,nsoilmx) temperature inside the ground (K)
45	c pcapcal(klon) surfacic specific heat (Wm-2s*K-1)
46	c pfluxgrd(klon) surface diffusive flux from ground (Wm-2)
47	c
48	c=======================================================================
49	c declarations:
50	c -------------
51
52	cym#include "dimensions.h"
53	#include "YOMCST.h"
54	cym#include "dimphy.h"
55	#include "dimsoil.h"
56	#include "indicesol.h"
57
58	c-----------------------------------------------------------------------
59	c arguments
60	c ---------
61
62	REAL ptimestep
63	INTEGER indice, knon
64	REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon)
65	REAL pcapcal(klon),pfluxgrd(klon)
66
67	c-----------------------------------------------------------------------
68	c local arrays
69	c ------------
70
71	INTEGER ig,jk
72	c@$$ REAL zdz2(nsoilmx),z1(klon)
73	REAL zdz2(nsoilmx),z1(klon,nbsrf)
74	REAL,SAVE :: min_period,dalph_soil
75	REAL ztherm_i(klon)
76
77	c local saved variables:
78	c ----------------------
79	REAL dz1(nsoilmx),dz2(nsoilmx)
80	c@$$ REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
81	cym REAL zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf)
82	REAL,ALLOCATABLE,SAVE :: zc(:,:,:),zd(:,:,:)
83	c$OMP THREADPRIVATE(zc,zd)
84	REAL lambda
85	cym SAVE dz1,dz2,zc,zd,lambda
86	SAVE dz1,dz2,lambda
87	c$OMP THREADPRIVATE(dz1,dz2,lambda)
88	LOGICAL firstcall, firstsurf(nbsrf)
89	SAVE firstcall, firstsurf
90	c$OMP THREADPRIVATE(firstcall, firstsurf)
91	REAL isol,isno,iice
92	SAVE isol,isno,iice
93	c$OMP THREADPRIVATE(isol,isno,iice)
94	DATA firstcall/.true./
95	DATA firstsurf/.TRUE.,.TRUE.,.TRUE.,.TRUE./
96
97	DATA isol,isno,iice/2000.,2000.,2000./
98	LOGICAL,SAVE :: First=.true.
99	c$OMP THREADPRIVATE(First)
100	c-----------------------------------------------------------------------
101	c Depthts:
102	c --------
103
104	REAL fz,rk,fz1,rk1,rk2
105	fz(rk)=fz1(dalph_soil*rk-1.)/(dalph_soil-1.)
106	pfluxgrd(:) = 0.
107	c calcul de l'inertie thermique a partir de la variable rnat.
108	c on initialise a iice meme au-dessus d'un point de mer au cas
109	c ou le point de mer devienne point de glace au pas suivant
110	c on corrige si on a un point de terre avec ou sans glace
111	c
112	IF (first) THEN
113	allocate(zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf))
114	first=.false.
115	ENDIF
116
117	IF (indice.EQ.is_sic) THEN
118	DO ig = 1, knon
119	ztherm_i(ig) = iice
120	IF (snow(ig).GT.0.0) ztherm_i(ig) = isno
121	ENDDO
122	ELSE IF (indice.EQ.is_lic) THEN
123	DO ig = 1, knon
124	ztherm_i(ig) = iice
125	IF (snow(ig).GT.0.0) ztherm_i(ig) = isno
126	ENDDO
127	ELSE IF (indice.EQ.is_ter) THEN
128	DO ig = 1, knon
129	ztherm_i(ig) = isol
130	IF (snow(ig).GT.0.0) ztherm_i(ig) = isno
131	ENDDO
132	ELSE IF (indice.EQ.is_oce) THEN
133	DO ig = 1, knon
134	ztherm_i(ig) = iice
135	ENDDO
136	ELSE
137	PRINT*, "valeur d indice non prevue", indice
138	CALL abort
139	ENDIF
140
141
142	c@$$ IF (firstcall) THEN
143	IF (firstsurf(indice)) THEN
144
145	c-----------------------------------------------------------------------
146	c ground levels
147	c grnd=z/l where l is the skin depth of the diurnal cycle:
148	c --------------------------------------------------------
149
150	min_period=1800. ! en secondes
151	dalph_soil=2. ! rapport entre les epaisseurs de 2 couches succ.
152	c$OMP MASTER
153	OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
154	READ(99,*) min_period
155	READ(99,*) dalph_soil
156	PRINT*,'Discretization for the soil model'
157	PRINT*,'First level e-folding depth',min_period,
158	s ' dalph',dalph_soil
159	CLOSE(99)
160	9999 CONTINUE
161	c$OMP END MASTER
162	c$OMP BARRIER
163
164	c la premiere couche represente un dixieme de cycle diurne
165	fz1=sqrt(min_period/3.14)
166
167	DO jk=1,nsoilmx
168	rk1=jk
169	rk2=jk-1
170	dz2(jk)=fz(rk1)-fz(rk2)
171	ENDDO
172	DO jk=1,nsoilmx-1
173	rk1=jk+.5
174	rk2=jk-.5
175	dz1(jk)=1./(fz(rk1)-fz(rk2))
176	ENDDO
177	lambda=fz(.5)*dz1(1)
178	PRINT*,'full layers, intermediate layers (seconds)'
179	DO jk=1,nsoilmx
180	rk=jk
181	rk1=jk+.5
182	rk2=jk-.5
183	PRINT *,'fz=',
184	. fz(rk1)fz(rk2)3.14,fz(rk)fz(rk)3.14
185	ENDDO
186	C PB
187	firstsurf(indice) = .FALSE.
188	c@$$ firstcall =.false.
189
190	c Initialisations:
191	c ----------------
192
193	ELSE !--not firstcall
194	c-----------------------------------------------------------------------
195	c Computation of the soil temperatures using the Cgrd and Dgrd
196	c coefficient computed at the previous time-step:
197	c -----------------------------------------------
198
199	c surface temperature
200	DO ig=1,knon
201	ptsoil(ig,1)=(lambda*zc(ig,1,indice)+ptsrf(ig))/
202	s (lambda*(1.-zd(ig,1,indice))+1.)
203	ENDDO
204
205	c other temperatures
206	DO jk=1,nsoilmx-1
207	DO ig=1,knon
208	ptsoil(ig,jk+1)=zc(ig,jk,indice)+zd(ig,jk,indice)
209	$ *ptsoil(ig,jk)
210	ENDDO
211	ENDDO
212
213	ENDIF !--not firstcall
214	c-----------------------------------------------------------------------
215	c Computation of the Cgrd and Dgrd coefficient for the next step:
216	c ---------------------------------------------------------------
217
218	c@$$ PB ajout pour cas glace de mer
219	IF (indice .EQ. is_sic) THEN
220	DO ig = 1 , knon
221	ptsoil(ig,nsoilmx) = RTT - 1.8
222	END DO
223	ENDIF
224
225	DO jk=1,nsoilmx
226	zdz2(jk)=dz2(jk)/ptimestep
227	ENDDO
228
229	DO ig=1,knon
230	z1(ig,indice)=zdz2(nsoilmx)+dz1(nsoilmx-1)
231	zc(ig,nsoilmx-1,indice)=
232	$ zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig,indice)
233	zd(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1(ig,indice)
234	ENDDO
235
236	DO jk=nsoilmx-1,2,-1
237	DO ig=1,knon
238	z1(ig,indice)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)
239	$ *(1.-zd(ig,jk,indice)))
240	zc(ig,jk-1,indice)=
241	s (ptsoil(ig,jk)zdz2(jk)+dz1(jk)zc(ig,jk,indice))
242	$ *z1(ig,indice)
243	zd(ig,jk-1,indice)=dz1(jk-1)*z1(ig,indice)
244	ENDDO
245	ENDDO
246
247	c-----------------------------------------------------------------------
248	c computation of the surface diffusive flux from ground and
249	c calorific capacity of the ground:
250	c ---------------------------------
251
252	DO ig=1,knon
253	pfluxgrd(ig)=ztherm_i(ig)dz1(1)
254	s (zc(ig,1,indice)+(zd(ig,1,indice)-1.)*ptsoil(ig,1))
255	pcapcal(ig)=ztherm_i(ig)*
256	s (dz2(1)+ptimestep(1.-zd(ig,1,indice))dz1(1))
257	z1(ig,indice)=lambda*(1.-zd(ig,1,indice))+1.
258	pcapcal(ig)=pcapcal(ig)/z1(ig,indice)
259	pfluxgrd(ig) = pfluxgrd(ig)
260	s + pcapcal(ig) * (ptsoil(ig,1) * z1(ig,indice)
261	$ - lambda * zc(ig,1,indice)
262	$ - ptsrf(ig))
263	s /ptimestep
264	ENDDO
265
266	RETURN
267	END

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

Download in other formats:

Original Format