source: LMDZ.3.3/branches/rel-LF/libf/phylmd/soil.F @ 456

Last change on this file since 456 was 235, checked in by lmdzadmin, 24 years ago

Integration des modifs pour fder de Pasb
LF

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