Context Navigation

surface.F90 @ 4245

Last change on this file since 4245 was 4245, checked in by dubos, 4 years ago
simple_physics : turn zc, zd, capcal, fluxgrd into local temporaries
File size: 11.3 KB

Line
1	MODULE surface
2
3	#include "use_logging.h"
4
5	IMPLICIT NONE
6	PRIVATE
7	SAVE
8
9	REAL, PARAMETER :: pi=2.*ASIN(1.)
10
11	! common variables
12	REAL, PUBLIC :: I_mer,I_ter,Cd_mer,Cd_ter, &
13	& alb_mer,alb_ter,emi_mer,emi_ter
14
15	! precomputed variables
16	REAL :: lambda
17	REAL, ALLOCATABLE :: dz1(:),dz2(:)
18	!$OMP THREADPRIVATE(dz1,dz2)
19	REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:)
20	!$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)
21
22	! internal state, written to / read from disk at checkpoint / restart
23	REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:)
24	!$OMP THREADPRIVATE(tsurf, tsoil)
25	! variables below should be temporary arrays, not persistent
26	REAL, ALLOCATABLE :: zc(:,:),zd(:,:), capcal(:), fluxgrd(:)
27	!$OMP THREADPRIVATE(zc,zd, capcal, fluxgrd)
28
29	PUBLIC :: init_soil, &
30	soil, soil_new, soil_forward, soil_backward, &
31	zc, zd, &
32	rnatur, albedo, emissiv, z0, inertie, &
33	tsurf, tsoil, capcal, fluxgrd
34
35	CONTAINS
36
37	SUBROUTINE init_soil(nsoil)
38	INTEGER, INTENT(IN) :: nsoil
39	REAL :: min_period,dalph_soil, rk,fz1,rk1,rk2
40	INTEGER :: jk
41
42	!-----------------------------------------------------------------------
43	! ground levels
44	! grnd=z/l where l is the skin depth of the diurnal cycle:
45	! --------------------------------------------------------
46
47	WRITELOG(,) 'nsoil,firstcall=',nsoil, .TRUE.
48
49	ALLOCATE(dz1(nsoil),dz2(nsoil))
50
51	min_period=20000.
52	dalph_soil=2.
53
54	! la premiere couche represente un dixieme de cycle diurne
55	fz1=sqrt(min_period/pi)
56
57	DO jk=1,nsoil
58	rk1=jk
59	rk2=jk-1
60	dz2(jk)=fz(rk1)-fz(rk2)
61	ENDDO
62	DO jk=1,nsoil-1
63	rk1=jk+.5
64	rk2=jk-.5
65	dz1(jk)=1./(fz(rk1)-fz(rk2))
66	ENDDO
67	lambda=fz(.5)*dz1(1)
68
69	WRITELOG(,) 'full layers, intermediate layers (secoonds)'
70	DO jk=1,nsoil
71	rk=jk
72	rk1=jk+.5
73	rk2=jk-.5
74	WRITELOG(,) fz(rk1)fz(rk2)pi, &
75	& fz(rk)fz(rk)pi
76	ENDDO
77	LOG_INFO('init_soil')
78
79	CONTAINS
80
81	FUNCTION fz(rk) RESULT(val)
82	REAL :: val, rk
83	val = fz1(dalph_soil*rk-1.)/(dalph_soil-1.)
84	END FUNCTION fz
85
86	END SUBROUTINE init_soil
87
88	PURE SUBROUTINE soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
89	INTEGER, INTENT(IN) :: ngrid, nsoil ! number of columns, of soil layers
90	REAL, INTENT(IN) :: zc(ngrid, nsoil), zd(ngrid, nsoil) ! LU factorization
91	REAL, INTENT(IN) :: ptsrf(ngrid) ! new surface temperature
92	REAL, INTENT(INOUT) :: ptsoil(ngrid,nsoil) ! soil temperature
93	INTEGER :: ig, jk
94
95	!-----------------------------------------------------------------------
96	! Computation of the soil temperatures using the Cgrd and Dgrd
97	! coefficient computed during the forward sweep
98	! -----------------------------------------------
99
100	! surface temperature => temperature in first soil layer
101	DO ig=1,ngrid
102	ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ &
103	& (lambda*(1.-zd(ig,1))+1.)
104	ENDDO
105
106	! other temperatures
107	DO jk=1,nsoil-1
108	DO ig=1,ngrid
109	ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
110	ENDDO
111	ENDDO
112	END SUBROUTINE Soil_backward
113
114	PURE SUBROUTINE soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
115	& zc, zd, pcapcal, pfluxgrd)
116	INTEGER, INTENT(IN) :: ngrid, nsoil ! number of columns, of soil layers
117	REAL, INTENT(IN) :: ptimestep, & ! time step
118	& ptherm_i(ngrid), & ! thermal inertia ??
119	& ptsrf(ngrid), & ! surface temperature before heat conduction
120	& ptsoil(ngrid, nsoil) ! soil temperature before heat conduction
121	REAL, INTENT(OUT) :: zc(ngrid,nsoil), &
122	& zd(ngrid, nsoil), & ! LU factorization for backward sweep
123	& pcapcal(ngrid), & ! effective calorific capacity
124	& pfluxgrd(ngrid) ! conductive heat flux at the ground
125	REAL :: z1, zdz2(ngrid)
126	INTEGER :: jk, ig
127
128	!-----------------------------------------------------------------------
129	! Computation of the Cgrd and Dgrd coefficients the backward sweep :
130	! ---------------------------------------------------------------
131
132	DO jk=1,nsoil
133	zdz2(jk)=dz2(jk)/ptimestep
134	ENDDO
135
136	DO ig=1,ngrid
137	z1=zdz2(nsoil)+dz1(nsoil-1)
138	zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1
139	zd(ig,nsoil-1)=dz1(nsoil-1)/z1
140	ENDDO
141
142	DO jk=nsoil-1,2,-1
143	DO ig=1,ngrid
144	z1=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
145	zc(ig,jk-1)= &
146	& (ptsoil(ig,jk)zdz2(jk)+dz1(jk)zc(ig,jk))*z1
147	zd(ig,jk-1)=dz1(jk-1)*z1
148	ENDDO
149	ENDDO
150
151	!-----------------------------------------------------------------------
152	! computation of the surface diffusive flux from ground and
153	! calorific capacity of the ground:
154	! ---------------------------------
155
156	DO ig=1,ngrid
157	pfluxgrd(ig)=ptherm_i(ig)dz1(1) &
158	& (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
159	z1=lambda*(1.-zd(ig,1))+1.
160	pcapcal(ig)=ptherm_i(ig)* &
161	& ptimestep(zdz2(1)+(1.-zd(ig,1))dz1(1))/z1
162	pfluxgrd(ig)=pfluxgrd(ig) &
163	& +pcapcal(ig)(ptsoil(ig,1)z1-lambda*zc(ig,1)-ptsrf(ig)) &
164	& /ptimestep
165	ENDDO
166	END SUBROUTINE soil_forward
167
168	SUBROUTINE soil_new(ngrid,nsoil,ptimestep,ptherm_i, ptsrf,ptsoil, pcapcal,pfluxgrd)
169	INTEGER, INTENT(IN) :: ngrid, nsoil ! number of columns, of soil layers
170	REAL, INTENT(IN) :: ptimestep, & ! time step
171	& ptherm_i(ngrid) ! thermal inertia ??
172	REAL, INTENT(INOUT) :: ptsrf(ngrid), & ! surface temperature
173	& ptsoil(ngrid,nsoil) ! soil temperature
174	REAL, INTENT(OUT) :: pcapcal(ngrid), & ! effective calorific capacity
175	& pfluxgrd(ngrid) ! conductive heat flux at the ground
176	CALL soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
177	CALL soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
178	& zc, zd, pcapcal, pfluxgrd)
179	END SUBROUTINE soil_new
180
181	SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i, &
182	& ptimestep,ptsrf,ptsoil, &
183	& pcapcal,pfluxgrd)
184
185	!=======================================================================
186	!
187	! Auteur: Frederic Hourdin 30/01/92
188	! -------
189	!
190	! objet: computation of : the soil temperature evolution
191	! ------ the surfacic heat capacity "Capcal"
192	! the surface conduction flux pcapcal
193	!
194	!
195	! Method: implicit time integration
196	! -------
197	! Consecutive ground temperatures are related by:
198	! T(k+1) = C(k) + D(k)*T(k) (1)
199	! the coefficients C and D are computed at the t-dt time-step.
200	! Routine structure:
201	! 1)new temperatures are computed using (1)
202	! 2)C and D coefficients are computed from the new temperature
203	! profile for the t+dt time-step
204	! 3)the coefficients A and B are computed where the diffusive
205	! fluxes at the t+dt time-step is given by
206	! Fdiff = A + B Ts(t+dt)
207	! or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
208	! with F0 = A + B (Ts(t))
209	! Capcal = B*dt
210	!
211	! Interface:
212	! ----------
213	!
214	! Arguments:
215	! ----------
216	! ngrid number of grid-points
217	! ptimestep physical timestep (s)
218	! pto(ngrid,nsoil) temperature at time-step t (K)
219	! ptn(ngrid,nsoil) temperature at time step t+dt (K)
220	! pcapcal(ngrid) specific heat (Wm-2s*K-1)
221	! pfluxgrd(ngrid) surface diffusive flux from ground (Wm-2)
222	!
223	!=======================================================================
224	! declarations:
225	! -------------
226
227
228	!-----------------------------------------------------------------------
229	! arguments
230	! ---------
231
232	INTEGER ngrid,nsoil
233	REAL ptimestep
234	REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid)
235	REAL pcapcal(ngrid),pfluxgrd(ngrid)
236	LOGICAL firstcall
237
238
239	!-----------------------------------------------------------------------
240	! local arrays
241	! ------------
242
243	INTEGER ig,jk
244	REAL zdz2(nsoil),z1(ngrid)
245
246	IF (firstcall) THEN
247	! init_soil is now called by iniphyparam
248	! CALL init_soil(ngrid, nsoil)
249	ELSE
250	IF(.FALSE.) THEN
251	!-----------------------------------------------------------------------
252	! Computation of the soil temperatures using the Cgrd and Dgrd
253	! coefficient computed at the previous time-step:
254	! -----------------------------------------------
255
256	! surface temperature
257	DO ig=1,ngrid
258	ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ &
259	& (lambda*(1.-zd(ig,1))+1.)
260	ENDDO
261
262	! other temperatures
263	DO jk=1,nsoil-1
264	DO ig=1,ngrid
265	ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
266	ENDDO
267	ENDDO
268	ELSE
269	CALL soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
270	END IF
271
272	ENDIF
273
274	IF(.FALSE.) THEN
275	!-----------------------------------------------------------------------
276	! Computation of the Cgrd and Dgrd coefficient for the next step:
277	! ---------------------------------------------------------------
278
279	DO jk=1,nsoil
280	zdz2(jk)=dz2(jk)/ptimestep
281	ENDDO
282
283	DO ig=1,ngrid
284	z1(ig)=zdz2(nsoil)+dz1(nsoil-1)
285	zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)
286	zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)
287	ENDDO
288
289	DO jk=nsoil-1,2,-1
290	DO ig=1,ngrid
291	z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
292	zc(ig,jk-1)= &
293	& (ptsoil(ig,jk)zdz2(jk)+dz1(jk)zc(ig,jk))*z1(ig)
294	zd(ig,jk-1)=dz1(jk-1)*z1(ig)
295	ENDDO
296	ENDDO
297
298	!-----------------------------------------------------------------------
299	! computation of the surface diffusive flux from ground and
300	! calorific capacity of the ground:
301	! ---------------------------------
302
303	DO ig=1,ngrid
304	pfluxgrd(ig)=ptherm_i(ig)dz1(1) &
305	& (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
306	z1(ig)=lambda*(1.-zd(ig,1))+1.
307	pcapcal(ig)=ptherm_i(ig)* &
308	& ptimestep(zdz2(1)+(1.-zd(ig,1))dz1(1))/z1(ig)
309	pfluxgrd(ig)=pfluxgrd(ig) &
310	& +pcapcal(ig)(ptsoil(ig,1)z1(ig)-lambda*zc(ig,1)-ptsrf(ig)) &
311	& /ptimestep
312	ENDDO
313	ELSE
314	CALL soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
315	& zc, zd, pcapcal, pfluxgrd)
316	END IF
317	END SUBROUTINE soil
318
319	END MODULE surface
320

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

Download in other formats:

Original Format