Context Navigation

soil.F @ 634

Last change on this file since 634 was 634, checked in by Laurent Fairhead, 19 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 8.0 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 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	REAL lambda
84	cym SAVE dz1,dz2,zc,zd,lambda
85	SAVE dz1,dz2,lambda
86	LOGICAL firstcall, firstsurf(nbsrf)
87	SAVE firstcall, firstsurf
88	REAL isol,isno,iice
89	SAVE isol,isno,iice
90
91	DATA firstcall/.true./
92	DATA firstsurf/.TRUE.,.TRUE.,.TRUE.,.TRUE./
93
94	DATA isol,isno,iice/2000.,2000.,2000./
95	LOGICAL,SAVE :: First=.true.
96	c-----------------------------------------------------------------------
97	c Depthts:
98	c --------
99
100	REAL fz,rk,fz1,rk1,rk2
101	fz(rk)=fz1(dalph_soil*rk-1.)/(dalph_soil-1.)
102	pfluxgrd(:) = 0.
103	c calcul de l'inertie thermique a partir de la variable rnat.
104	c on initialise a iice meme au-dessus d'un point de mer au cas
105	c ou le point de mer devienne point de glace au pas suivant
106	c on corrige si on a un point de terre avec ou sans glace
107	c
108	IF (first) THEN
109	allocate(zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf))
110	first=.false.
111	ENDIF
112
113	IF (indice.EQ.is_sic) THEN
114	DO ig = 1, knon
115	ztherm_i(ig) = iice
116	IF (snow(ig).GT.0.0) ztherm_i(ig) = isno
117	ENDDO
118	ELSE IF (indice.EQ.is_lic) THEN
119	DO ig = 1, knon
120	ztherm_i(ig) = iice
121	IF (snow(ig).GT.0.0) ztherm_i(ig) = isno
122	ENDDO
123	ELSE IF (indice.EQ.is_ter) THEN
124	DO ig = 1, knon
125	ztherm_i(ig) = isol
126	IF (snow(ig).GT.0.0) ztherm_i(ig) = isno
127	ENDDO
128	ELSE IF (indice.EQ.is_oce) THEN
129	DO ig = 1, knon
130	ztherm_i(ig) = iice
131	ENDDO
132	ELSE
133	PRINT*, "valeur d indice non prevue", indice
134	CALL abort
135	ENDIF
136
137
138	c$$$ IF (firstcall) THEN
139	IF (firstsurf(indice)) THEN
140
141	c-----------------------------------------------------------------------
142	c ground levels
143	c grnd=z/l where l is the skin depth of the diurnal cycle:
144	c --------------------------------------------------------
145
146	min_period=1800. ! en secondes
147	dalph_soil=2. ! rapport entre les epaisseurs de 2 couches succ.
148
149	OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
150	READ(99,*) min_period
151	READ(99,*) dalph_soil
152	PRINT*,'Discretization for the soil model'
153	PRINT*,'First level e-folding depth',min_period,
154	s ' dalph',dalph_soil
155	CLOSE(99)
156	9999 CONTINUE
157
158	c la premiere couche represente un dixieme de cycle diurne
159	fz1=sqrt(min_period/3.14)
160
161	DO jk=1,nsoilmx
162	rk1=jk
163	rk2=jk-1
164	dz2(jk)=fz(rk1)-fz(rk2)
165	ENDDO
166	DO jk=1,nsoilmx-1
167	rk1=jk+.5
168	rk2=jk-.5
169	dz1(jk)=1./(fz(rk1)-fz(rk2))
170	ENDDO
171	lambda=fz(.5)*dz1(1)
172	PRINT*,'full layers, intermediate layers (seconds)'
173	DO jk=1,nsoilmx
174	rk=jk
175	rk1=jk+.5
176	rk2=jk-.5
177	PRINT *,'fz=',
178	. fz(rk1)fz(rk2)3.14,fz(rk)fz(rk)3.14
179	ENDDO
180	C PB
181	firstsurf(indice) = .FALSE.
182	c$$$ firstcall =.false.
183
184	c Initialisations:
185	c ----------------
186
187	ELSE !--not firstcall
188	c-----------------------------------------------------------------------
189	c Computation of the soil temperatures using the Cgrd and Dgrd
190	c coefficient computed at the previous time-step:
191	c -----------------------------------------------
192
193	c surface temperature
194	DO ig=1,knon
195	ptsoil(ig,1)=(lambda*zc(ig,1,indice)+ptsrf(ig))/
196	s (lambda*(1.-zd(ig,1,indice))+1.)
197	ENDDO
198
199	c other temperatures
200	DO jk=1,nsoilmx-1
201	DO ig=1,knon
202	ptsoil(ig,jk+1)=zc(ig,jk,indice)+zd(ig,jk,indice)
203	$ *ptsoil(ig,jk)
204	ENDDO
205	ENDDO
206
207	ENDIF !--not firstcall
208	c-----------------------------------------------------------------------
209	c Computation of the Cgrd and Dgrd coefficient for the next step:
210	c ---------------------------------------------------------------
211
212	c$$$ PB ajout pour cas glace de mer
213	IF (indice .EQ. is_sic) THEN
214	DO ig = 1 , knon
215	ptsoil(ig,nsoilmx) = RTT - 1.8
216	END DO
217	ENDIF
218
219	DO jk=1,nsoilmx
220	zdz2(jk)=dz2(jk)/ptimestep
221	ENDDO
222
223	DO ig=1,knon
224	z1(ig,indice)=zdz2(nsoilmx)+dz1(nsoilmx-1)
225	zc(ig,nsoilmx-1,indice)=
226	$ zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig,indice)
227	zd(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1(ig,indice)
228	ENDDO
229
230	DO jk=nsoilmx-1,2,-1
231	DO ig=1,knon
232	z1(ig,indice)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)
233	$ *(1.-zd(ig,jk,indice)))
234	zc(ig,jk-1,indice)=
235	s (ptsoil(ig,jk)zdz2(jk)+dz1(jk)zc(ig,jk,indice))
236	$ *z1(ig,indice)
237	zd(ig,jk-1,indice)=dz1(jk-1)*z1(ig,indice)
238	ENDDO
239	ENDDO
240
241	c-----------------------------------------------------------------------
242	c computation of the surface diffusive flux from ground and
243	c calorific capacity of the ground:
244	c ---------------------------------
245
246	DO ig=1,knon
247	pfluxgrd(ig)=ztherm_i(ig)dz1(1)
248	s (zc(ig,1,indice)+(zd(ig,1,indice)-1.)*ptsoil(ig,1))
249	pcapcal(ig)=ztherm_i(ig)*
250	s (dz2(1)+ptimestep(1.-zd(ig,1,indice))dz1(1))
251	z1(ig,indice)=lambda*(1.-zd(ig,1,indice))+1.
252	pcapcal(ig)=pcapcal(ig)/z1(ig,indice)
253	pfluxgrd(ig) = pfluxgrd(ig)
254	s + pcapcal(ig) * (ptsoil(ig,1) * z1(ig,indice)
255	$ - lambda * zc(ig,1,indice)
256	$ - ptsrf(ig))
257	s /ptimestep
258	ENDDO
259
260	RETURN
261	END

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

Download in other formats:

Original Format