source: LMDZ.3.3/branches/rel-1-0-patch/libf/phylmd/soil.F @ 2167

Last change on this file since 2167 was 253, checked in by (none), 23 years ago

This commit was manufactured by cvs2svn to create branch
'rel-1-0-patch'.
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 6.9 KB
Line 
1      SUBROUTINE soil(ptimestep, indice, 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 "dimphy.h"
50#include "dimsoil.h"
51#include "indicesol.h"
52
53c-----------------------------------------------------------------------
54c  arguments
55c  ---------
56
57      REAL ptimestep
58      INTEGER indice
59      REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon)
60      REAL pcapcal(klon),pfluxgrd(klon)
61
62c-----------------------------------------------------------------------
63c  local arrays
64c  ------------
65
66      INTEGER ig,jk
67      REAL zdz2(nsoilmx),z1(klon)
68      REAL min_period,dalph_soil
69      REAL ztherm_i(klon)
70
71c   local saved variables:
72c   ----------------------
73      REAL dz1(nsoilmx),dz2(nsoilmx)
74      REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
75      REAL lambda
76      SAVE dz1,dz2,zc,zd,lambda
77      LOGICAL firstcall
78      SAVE firstcall
79      REAL isol,isno,iice
80      SAVE isol,isno,iice
81
82      DATA firstcall/.true./
83
84      DATA isol,isno,iice/2000.,2000.,2000./
85
86c-----------------------------------------------------------------------
87c   Depthts:
88c   --------
89
90      REAL fz,rk,fz1,rk1,rk2
91      fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
92
93c   calcul de l'inertie thermique a partir de la variable rnat.
94c   on initialise a iice meme au-dessus d'un point de mer au cas
95c   ou le point de mer devienne point de glace au pas suivant
96c   on corrige si on a un point de terre avec ou sans glace
97c
98      IF (indice.EQ.is_sic) THEN
99         DO ig = 1, klon
100            ztherm_i(ig)   = iice
101            IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno
102         ENDDO
103      ELSE IF (indice.EQ.is_lic) THEN
104         DO ig = 1, klon
105            ztherm_i(ig)   = iice
106            IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno
107         ENDDO
108      ELSE IF (indice.EQ.is_ter) THEN
109         DO ig = 1, klon
110            ztherm_i(ig)   = isol
111            IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno
112         ENDDO
113      ELSE IF (indice.EQ.is_oce) THEN
114         DO ig = 1, klon
115            ztherm_i(ig)   = iice
116         ENDDO
117      ELSE
118         PRINT*, "valeur d indice non prevue", indice
119         CALL abort
120      ENDIF
121
122
123      IF (firstcall) THEN
124
125c-----------------------------------------------------------------------
126c   ground levels
127c   grnd=z/l where l is the skin depth of the diurnal cycle:
128c   --------------------------------------------------------
129
130         min_period=1800. ! en secondes
131         dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
132
133         OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
134         READ(99,*) min_period
135         READ(99,*) dalph_soil
136         PRINT*,'Discretization for the soil model'
137         PRINT*,'First level e-folding depth',min_period,
138     s   '   dalph',dalph_soil
139         CLOSE(99)
1409999     CONTINUE
141
142c   la premiere couche represente un dixieme de cycle diurne
143         fz1=sqrt(min_period/3.14)
144
145         DO jk=1,nsoilmx
146            rk1=jk
147            rk2=jk-1
148            dz2(jk)=fz(rk1)-fz(rk2)
149         ENDDO
150         DO jk=1,nsoilmx-1
151            rk1=jk+.5
152            rk2=jk-.5
153            dz1(jk)=1./(fz(rk1)-fz(rk2))
154         ENDDO
155         lambda=fz(.5)*dz1(1)
156         PRINT*,'full layers, intermediate layers (seconds)'
157         DO jk=1,nsoilmx
158            rk=jk
159            rk1=jk+.5
160            rk2=jk-.5
161            PRINT *,'fz=',
162     .               fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
163         ENDDO
164
165         firstcall =.false.
166
167c   Initialisations:
168c   ----------------
169
170      ELSE   !--not firstcall
171c-----------------------------------------------------------------------
172c   Computation of the soil temperatures using the Cgrd and Dgrd
173c  coefficient computed at the previous time-step:
174c  -----------------------------------------------
175
176c    surface temperature
177         DO ig=1,klon
178            ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/
179     s      (lambda*(1.-zd(ig,1))+1.)
180         ENDDO
181
182c   other temperatures
183         DO jk=1,nsoilmx-1
184            DO ig=1,klon
185               ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
186            ENDDO
187         ENDDO
188
189      ENDIF !--not firstcall
190c-----------------------------------------------------------------------
191c   Computation of the Cgrd and Dgrd coefficient for the next step:
192c   ---------------------------------------------------------------
193
194      DO jk=1,nsoilmx
195         zdz2(jk)=dz2(jk)/ptimestep
196      ENDDO
197
198      DO ig=1,klon
199         z1(ig)=zdz2(nsoilmx)+dz1(nsoilmx-1)
200         zc(ig,nsoilmx-1)=zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig)
201         zd(ig,nsoilmx-1)=dz1(nsoilmx-1)/z1(ig)
202      ENDDO
203
204      DO jk=nsoilmx-1,2,-1
205         DO ig=1,klon
206            z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
207            zc(ig,jk-1)=
208     s      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)
209            zd(ig,jk-1)=dz1(jk-1)*z1(ig)
210         ENDDO
211      ENDDO
212
213c-----------------------------------------------------------------------
214c   computation of the surface diffusive flux from ground and
215c   calorific capacity of the ground:
216c   ---------------------------------
217
218      DO ig=1,klon
219         pfluxgrd(ig)=ztherm_i(ig)*dz1(1)*
220     s   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
221         pcapcal(ig)=ztherm_i(ig)*
222     s   (dz2(1)+ptimestep*(1.-zd(ig,1))*dz1(1))
223         z1(ig)=lambda*(1.-zd(ig,1))+1.
224         pcapcal(ig)=pcapcal(ig)/z1(ig)
225         pfluxgrd(ig)=pfluxgrd(ig)
226     s   +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))
227     s   /ptimestep
228      ENDDO
229
230      RETURN
231      END
Note: See TracBrowser for help on using the repository browser.