source: LMDZ4/trunk/libf/phytherm/soil.F @ 815

Last change on this file since 815 was 814, checked in by Laurent Fairhead, 17 years ago

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