source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/soil.F90 @ 3814

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

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