source: LMDZ4/branches/LMDZ4_par_0/libf/phylmd/soil.F @ 5179

Last change on this file since 5179 was 634, checked in by Laurent Fairhead, 20 years ago

Modifications faites à la physique pour la rendre parallele YM
Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester
les modifs pleinement avant leurs inclusions dans le tronc principal
LF

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