source: LMDZ4/trunk/libf/phylmd/soil.F90 @ 5068

Last change on this file since 5068 was 996, checked in by lsce, 16 years ago
  • Modifications liées au calcul des nouveau sous-fractions
  • Nettoyage de ocean slab : il reste uniquement la version avec glace de mer forcé
  • Nouveaux variables pour distiguer la version et type d'ocean : type_ocean=force/slab/couple, version_ocean=opa8/nemo pour couplé ou version_ocean=sicOBS pour slab

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