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

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

modifications pour faire de l'aquaplanète 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#include "comsoil.h"
59
60c-----------------------------------------------------------------------
61c  arguments
62c  ---------
63
64      REAL ptimestep
65      INTEGER indice, knon
66      REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon)
67      REAL pcapcal(klon),pfluxgrd(klon)
68
69c-----------------------------------------------------------------------
70c  local arrays
71c  ------------
72
73      INTEGER ig,jk
74c@$$      REAL zdz2(nsoilmx),z1(klon)
75      REAL zdz2(nsoilmx),z1(klon,nbsrf)
76      REAL,SAVE :: min_period,dalph_soil
77c$OMP THREADPRIVATE( min_period,dalph_soil)     
78      REAL ztherm_i(klon)
79
80c   local saved variables:
81c   ----------------------
82      REAL dz1(nsoilmx),dz2(nsoilmx)
83c@$$          REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
84cym      REAL zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf)
85      REAL,ALLOCATABLE,SAVE ::  zc(:,:,:),zd(:,:,:)
86c$OMP THREADPRIVATE(zc,zd)
87      REAL lambda
88cym      SAVE dz1,dz2,zc,zd,lambda
89      SAVE dz1,dz2,lambda
90c$OMP THREADPRIVATE(dz1,dz2,lambda)
91      LOGICAL firstcall, firstsurf(nbsrf)
92      SAVE firstcall, firstsurf
93c$OMP THREADPRIVATE(firstcall, firstsurf)
94      DATA firstcall/.true./
95      DATA firstsurf/.TRUE.,.TRUE.,.TRUE.,.TRUE./
96
97      LOGICAL,SAVE :: First=.true.
98c$OMP THREADPRIVATE(First)
99c-----------------------------------------------------------------------
100c   Depthts:
101c   --------
102
103      REAL fz,rk,fz1,rk1,rk2
104      fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
105      pfluxgrd(:) = 0.
106c   calcul de l'inertie thermique a partir de la variable rnat.
107c   on initialise a inertie_ice meme au-dessus d'un point de mer au cas
108c   ou le point de mer devienne point de glace au pas suivant
109c   on corrige si on a un point de terre avec ou sans glace
110c
111      IF (first) THEN
112        allocate(zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf))
113        first=.false.
114      ENDIF
115     
116      IF (indice.EQ.is_sic) THEN
117         DO ig = 1, knon
118            ztherm_i(ig)   = inertie_ice
119            IF (snow(ig).GT.0.0) ztherm_i(ig)   = inertie_sno
120         ENDDO
121      ELSE IF (indice.EQ.is_lic) THEN
122         DO ig = 1, knon
123            ztherm_i(ig)   = inertie_ice
124            IF (snow(ig).GT.0.0) ztherm_i(ig)   = inertie_sno
125         ENDDO
126      ELSE IF (indice.EQ.is_ter) THEN
127         DO ig = 1, knon
128            ztherm_i(ig)   = inertie_sol
129            IF (snow(ig).GT.0.0) ztherm_i(ig)   = inertie_sno
130         ENDDO
131      ELSE IF (indice.EQ.is_oce) THEN
132         DO ig = 1, knon
133            ztherm_i(ig)   = inertie_ice
134         ENDDO
135      ELSE
136         PRINT*, "valeur d indice non prevue", indice
137         CALL abort
138      ENDIF
139
140
141c@$$      IF (firstcall) THEN
142      IF (firstsurf(indice)) THEN
143
144c-----------------------------------------------------------------------
145c   ground levels
146c   grnd=z/l where l is the skin depth of the diurnal cycle:
147c   --------------------------------------------------------
148
149         min_period=1800. ! en secondes
150         dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
151c$OMP MASTER
152         IF (is_mpi_root) THEN
153         OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
154         READ(99,*) min_period
155         READ(99,*) dalph_soil
156         PRINT*,'Discretization for the soil model'
157         PRINT*,'First level e-folding depth',min_period,
158     s   '   dalph',dalph_soil
159         CLOSE(99)
1609999     CONTINUE
161         ENDIF  ! is_mpi_root
162c$OMP END MASTER
163         CALL bcast(min_period)
164         CALL bcast(dalph_soil)
165
166c   la premiere couche represente un dixieme de cycle diurne
167         fz1=sqrt(min_period/3.14)
168
169         DO jk=1,nsoilmx
170            rk1=jk
171            rk2=jk-1
172            dz2(jk)=fz(rk1)-fz(rk2)
173         ENDDO
174         DO jk=1,nsoilmx-1
175            rk1=jk+.5
176            rk2=jk-.5
177            dz1(jk)=1./(fz(rk1)-fz(rk2))
178         ENDDO
179         lambda=fz(.5)*dz1(1)
180         PRINT*,'full layers, intermediate layers (seconds)'
181         DO jk=1,nsoilmx
182            rk=jk
183            rk1=jk+.5
184            rk2=jk-.5
185            PRINT *,'fz=',
186     .               fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
187         ENDDO
188C PB
189         firstsurf(indice) = .FALSE.
190c@$$         firstcall =.false.
191
192c   Initialisations:
193c   ----------------
194
195      ELSE   !--not firstcall
196c-----------------------------------------------------------------------
197c   Computation of the soil temperatures using the Cgrd and Dgrd
198c  coefficient computed at the previous time-step:
199c  -----------------------------------------------
200
201c    surface temperature
202         DO ig=1,knon
203            ptsoil(ig,1)=(lambda*zc(ig,1,indice)+ptsrf(ig))/
204     s      (lambda*(1.-zd(ig,1,indice))+1.)
205         ENDDO
206
207c   other temperatures
208         DO jk=1,nsoilmx-1
209            DO ig=1,knon
210               ptsoil(ig,jk+1)=zc(ig,jk,indice)+zd(ig,jk,indice)
211     $            *ptsoil(ig,jk)
212            ENDDO
213         ENDDO
214
215      ENDIF !--not firstcall
216c-----------------------------------------------------------------------
217c   Computation of the Cgrd and Dgrd coefficient for the next step:
218c   ---------------------------------------------------------------
219
220c@$$  PB ajout pour cas glace de mer
221      IF (indice .EQ. is_sic) THEN
222          DO ig = 1 , knon
223            ptsoil(ig,nsoilmx) = RTT - 1.8
224          END DO
225      ENDIF
226
227      DO jk=1,nsoilmx
228         zdz2(jk)=dz2(jk)/ptimestep
229      ENDDO
230
231      DO ig=1,knon
232         z1(ig,indice)=zdz2(nsoilmx)+dz1(nsoilmx-1)
233         zc(ig,nsoilmx-1,indice)=
234     $       zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig,indice)
235         zd(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1(ig,indice)
236      ENDDO
237
238      DO jk=nsoilmx-1,2,-1
239         DO ig=1,knon
240            z1(ig,indice)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)
241     $         *(1.-zd(ig,jk,indice)))
242            zc(ig,jk-1,indice)=
243     s      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice))
244     $          *z1(ig,indice)
245            zd(ig,jk-1,indice)=dz1(jk-1)*z1(ig,indice)
246         ENDDO
247      ENDDO
248
249c-----------------------------------------------------------------------
250c   computation of the surface diffusive flux from ground and
251c   calorific capacity of the ground:
252c   ---------------------------------
253
254      DO ig=1,knon
255         pfluxgrd(ig)=ztherm_i(ig)*dz1(1)*
256     s   (zc(ig,1,indice)+(zd(ig,1,indice)-1.)*ptsoil(ig,1))
257         pcapcal(ig)=ztherm_i(ig)*
258     s   (dz2(1)+ptimestep*(1.-zd(ig,1,indice))*dz1(1))
259         z1(ig,indice)=lambda*(1.-zd(ig,1,indice))+1.
260         pcapcal(ig)=pcapcal(ig)/z1(ig,indice)
261         pfluxgrd(ig) = pfluxgrd(ig)
262     s   + pcapcal(ig) * (ptsoil(ig,1) * z1(ig,indice)
263     $       - lambda * zc(ig,1,indice)
264     $       - ptsrf(ig))
265     s   /ptimestep
266      ENDDO
267
268      RETURN
269      END
Note: See TracBrowser for help on using the repository browser.