source: LMDZ5/trunk/libf/phylmd/soil.F90 @ 1748

Last change on this file since 1748 was 1575, checked in by jghattas, 13 years ago
  • Added suffix _mpi_rank to lmdz.out text file. Each processus now write into seperate file but only if lunout/=6 (lunout is set in run.def).
  • Change some print* into write(lunout,*)
  • phytrac.F90 : always include ini_histrac and write_histrac. The file histrac.nc is written if ecrit_tra> 0 (set in physiq.def). Change default value of ecrit_tra into 0.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.9 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  INCLUDE "iniprint.h"
57!-----------------------------------------------------------------------
58! Arguments
59! ---------
60  REAL, INTENT(IN)                     :: ptimestep
61  INTEGER, INTENT(IN)                  :: indice, knon
62  REAL, DIMENSION(klon), INTENT(IN)    :: snow
63  REAL, DIMENSION(klon), INTENT(IN)    :: ptsrf
64 
65  REAL, DIMENSION(klon,nsoilmx), INTENT(INOUT) :: ptsoil
66  REAL, DIMENSION(klon), INTENT(OUT)           :: pcapcal
67  REAL, DIMENSION(klon), INTENT(OUT)           :: pfluxgrd
68
69!-----------------------------------------------------------------------
70! Local variables
71! ---------------
72  INTEGER                             :: ig, jk, ierr
73  REAL                                :: min_period,dalph_soil
74  REAL, DIMENSION(nsoilmx)            :: zdz2
75  REAL                                :: z1s
76  REAL, DIMENSION(klon)               :: ztherm_i
77  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: C_coef, D_coef
78
79! Local saved variables
80! ---------------------
81  REAL, SAVE                     :: lambda
82!$OMP THREADPRIVATE(lambda)
83  REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2
84!$OMP THREADPRIVATE(dz1,dz2)
85  LOGICAL, SAVE                  :: firstcall=.TRUE.
86!$OMP THREADPRIVATE(firstcall)
87   
88!-----------------------------------------------------------------------
89!   Depthts:
90!   --------
91  REAL fz,rk,fz1,rk1,rk2
92  fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
93
94
95!-----------------------------------------------------------------------
96! Calculation of some constants
97! NB! These constants do not depend on the sub-surfaces
98!-----------------------------------------------------------------------
99
100  IF (firstcall) THEN
101!-----------------------------------------------------------------------
102!   ground levels
103!   grnd=z/l where l is the skin depth of the diurnal cycle:
104!-----------------------------------------------------------------------
105
106     min_period=1800. ! en secondes
107     dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
108!$OMP MASTER
109     IF (is_mpi_root) THEN
110        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
111        IF (ierr == 0) THEN ! Read file only if it exists
112           READ(99,*) min_period
113           READ(99,*) dalph_soil
114           WRITE(lunout,*)'Discretization for the soil model'
115           WRITE(lunout,*)'First level e-folding depth',min_period, &
116                '   dalph',dalph_soil
117           CLOSE(99)
118        END IF
119     ENDIF
120!$OMP END MASTER
121     CALL bcast(min_period)
122     CALL bcast(dalph_soil)
123
124!   la premiere couche represente un dixieme de cycle diurne
125     fz1=SQRT(min_period/3.14)
126     
127     DO jk=1,nsoilmx
128        rk1=jk
129        rk2=jk-1
130        dz2(jk)=fz(rk1)-fz(rk2)
131     ENDDO
132     DO jk=1,nsoilmx-1
133        rk1=jk+.5
134        rk2=jk-.5
135        dz1(jk)=1./(fz(rk1)-fz(rk2))
136     ENDDO
137     lambda=fz(.5)*dz1(1)
138     WRITE(lunout,*)'full layers, intermediate layers (seconds)'
139     DO jk=1,nsoilmx
140        rk=jk
141        rk1=jk+.5
142        rk2=jk-.5
143        WRITE(lunout,*)'fz=', &
144             fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
145     ENDDO
146
147     firstcall =.FALSE.
148  END IF
149
150
151!-----------------------------------------------------------------------
152!   Calcul de l'inertie thermique a partir de la variable rnat.
153!   on initialise a inertie_ice meme au-dessus d'un point de mer au cas
154!   ou le point de mer devienne point de glace au pas suivant
155!   on corrige si on a un point de terre avec ou sans glace
156!
157!-----------------------------------------------------------------------
158  IF (indice == is_sic) THEN
159     DO ig = 1, knon
160        ztherm_i(ig)   = inertie_ice
161        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
162     ENDDO
163  ELSE IF (indice == is_lic) THEN
164     DO ig = 1, knon
165        ztherm_i(ig)   = inertie_ice
166        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
167     ENDDO
168  ELSE IF (indice == is_ter) THEN
169     DO ig = 1, knon
170        ztherm_i(ig)   = inertie_sol
171        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
172     ENDDO
173  ELSE IF (indice == is_oce) THEN
174     DO ig = 1, knon
175        ztherm_i(ig)   = inertie_ice
176     ENDDO
177  ELSE
178     WRITE(lunout,*) "valeur d indice non prevue", indice
179     CALL abort
180  ENDIF
181
182
183!-----------------------------------------------------------------------
184! 1)
185! Calculation of Cgrf and Dgrd coefficients using soil temperature from
186! previous time step.
187!
188! These variables are recalculated on the local compressed grid instead
189! of saved in restart file.
190!-----------------------------------------------------------------------
191  DO jk=1,nsoilmx
192     zdz2(jk)=dz2(jk)/ptimestep
193  ENDDO
194 
195  DO ig=1,knon
196     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
197     C_coef(ig,nsoilmx-1,indice)= &
198          zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
199     D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s
200  ENDDO
201 
202  DO jk=nsoilmx-1,2,-1
203     DO ig=1,knon
204        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
205             *(1.-D_coef(ig,jk,indice)))
206        C_coef(ig,jk-1,indice)= &
207             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
208        D_coef(ig,jk-1,indice)=dz1(jk-1)*z1s
209     ENDDO
210  ENDDO
211
212!-----------------------------------------------------------------------
213! 2)
214! Computation of the soil temperatures using the Cgrd and Dgrd
215! coefficient computed above
216!
217!-----------------------------------------------------------------------
218
219!    Surface temperature
220  DO ig=1,knon
221     ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/  &
222          (lambda*(1.-D_coef(ig,1,indice))+1.)
223  ENDDO
224 
225!   Other temperatures
226  DO jk=1,nsoilmx-1
227     DO ig=1,knon
228        ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) &
229             *ptsoil(ig,jk)
230     ENDDO
231  ENDDO
232
233  IF (indice == is_sic) THEN
234     DO ig = 1 , knon
235        ptsoil(ig,nsoilmx) = RTT - 1.8
236     END DO
237  ENDIF
238
239!-----------------------------------------------------------------------
240! 3)
241! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil
242! temperature
243!-----------------------------------------------------------------------
244  DO ig=1,knon
245     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
246     C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
247     D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s
248  ENDDO
249 
250  DO jk=nsoilmx-1,2,-1
251     DO ig=1,knon
252        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
253             *(1.-D_coef(ig,jk,indice)))
254        C_coef(ig,jk-1,indice) = &
255             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
256        D_coef(ig,jk-1,indice) = dz1(jk-1)*z1s
257     ENDDO
258  ENDDO
259
260!-----------------------------------------------------------------------
261! 4)
262! Computation of the surface diffusive flux from ground and
263! calorific capacity of the ground
264!-----------------------------------------------------------------------
265  DO ig=1,knon
266     pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* &
267          (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1))
268     pcapcal(ig)  = ztherm_i(ig)* &
269          (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1))
270     z1s = lambda*(1.-D_coef(ig,1,indice))+1.
271     pcapcal(ig)  = pcapcal(ig)/z1s
272     pfluxgrd(ig) = pfluxgrd(ig) &
273          + pcapcal(ig) * (ptsoil(ig,1) * z1s &
274          - lambda * C_coef(ig,1,indice) &
275          - ptsrf(ig)) &
276          /ptimestep
277  ENDDO
278   
279END SUBROUTINE soil
Note: See TracBrowser for help on using the repository browser.