source: LMDZ4/trunk/libf/phylmd/soil.F @ 540

Last change on this file since 540 was 524, checked in by lmdzadmin, 20 years ago

Initial revision

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