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

Last change on this file since 4011 was 2915, checked in by jbmadeleine, 7 years ago

Improved the way thermal inertia works in simulations using imposed SST and sea ice.
Removed inertie_ice which was confusing; Added inertie_sic for sea-ice and inertie_lic for
land-ice. When snow is present, inertie_sno is used instead of inertie_sol, inertie_sic or
inertie_lic. For sea-ice, added a flag called iflag_sic:
If iflag_sic=0, thermal inertia is changed over sea-ice when snow is present
If iflag_sic=1, thermal inertia is kept constant and equal to inertie_sic

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.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  USE indice_sol_mod
10  USE print_control_mod, ONLY: lunout
11
12  IMPLICIT NONE
13
14!=======================================================================
15!
16!   Auteur:  Frederic Hourdin     30/01/92
17!   -------
18!
19!   Object:  Computation of : the soil temperature evolution
20!   -------                   the surfacic heat capacity "Capcal"
21!                            the surface conduction flux pcapcal
22!
23!
24!   Method: Implicit time integration
25!   -------
26!   Consecutive ground temperatures are related by:
27!           T(k+1) = C(k) + D(k)*T(k)  (*)
28!   The coefficients C and D are computed at the t-dt time-step.
29!   Routine structure:
30!   1) C and D coefficients are computed from the old temperature
31!   2) new temperatures are computed using (*)
32!   3) C and D coefficients are computed from the new temperature
33!      profile for the t+dt time-step
34!   4) the coefficients A and B are computed where the diffusive
35!      fluxes at the t+dt time-step is given by
36!             Fdiff = A + B Ts(t+dt)
37!      or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
38!             with F0 = A + B (Ts(t))
39!                 Capcal = B*dt
40!           
41!   Interface:
42!   ----------
43!
44!   Arguments:
45!   ----------
46!   ptimestep            physical timestep (s)
47!   indice               sub-surface index
48!   snow(klon)           snow
49!   ptsrf(klon)          surface temperature at time-step t (K)
50!   ptsoil(klon,nsoilmx) temperature inside the ground (K)
51!   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
52!   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
53!   
54!=======================================================================
55  INCLUDE "YOMCST.h"
56  INCLUDE "dimsoil.h"
57  INCLUDE "comsoil.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_sic 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!   iophys can be used to write the ztherm_i variable in a phys.nc file
159!   and check the results; to do so, add "CALL iophys_ini" in physiq_mod
160!   and add knindex to the list of inputs in all the calls to soil.F90
161!   (and to soil.F90 itself !)
162!-----------------------------------------------------------------------
163
164  IF (indice == is_sic) THEN
165     DO ig = 1, knon
166        ztherm_i(ig)   = inertie_sic
167     ENDDO
168     IF (iflag_sic == 0) THEN
169       DO ig = 1, knon
170         IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
171       ENDDO
172!      Otherwise sea-ice keeps the same inertia, even when covered by snow
173     ENDIF
174!    CALL iophys_ecrit_index('ztherm_sic', 1, 'ztherm_sic', 'USI', &
175!      knon, knindex, ztherm_i)
176  ELSE IF (indice == is_lic) THEN
177     DO ig = 1, knon
178        ztherm_i(ig)   = inertie_lic
179        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
180     ENDDO
181!    CALL iophys_ecrit_index('ztherm_lic', 1, 'ztherm_lic', 'USI', &
182!      knon, knindex, ztherm_i)
183  ELSE IF (indice == is_ter) THEN
184     DO ig = 1, knon
185        ztherm_i(ig)   = inertie_sol
186        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
187     ENDDO
188!    CALL iophys_ecrit_index('ztherm_ter', 1, 'ztherm_ter', 'USI', &
189!      knon, knindex, ztherm_i)
190  ELSE IF (indice == is_oce) THEN
191     DO ig = 1, knon
192!       This is just in case, but SST should be used by the model anyway
193        ztherm_i(ig)   = inertie_sic
194     ENDDO
195!    CALL iophys_ecrit_index('ztherm_oce', 1, 'ztherm_oce', 'USI', &
196!      knon, knindex, ztherm_i)
197  ELSE
198     WRITE(lunout,*) "valeur d indice non prevue", indice
199     call abort_physic("soil", "", 1)
200  ENDIF
201
202
203!-----------------------------------------------------------------------
204! 1)
205! Calculation of Cgrf and Dgrd coefficients using soil temperature from
206! previous time step.
207!
208! These variables are recalculated on the local compressed grid instead
209! of saved in restart file.
210!-----------------------------------------------------------------------
211  DO jk=1,nsoilmx
212     zdz2(jk)=dz2(jk)/ptimestep
213  ENDDO
214 
215  DO ig=1,knon
216     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
217     C_coef(ig,nsoilmx-1,indice)= &
218          zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
219     D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s
220  ENDDO
221 
222  DO jk=nsoilmx-1,2,-1
223     DO ig=1,knon
224        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
225             *(1.-D_coef(ig,jk,indice)))
226        C_coef(ig,jk-1,indice)= &
227             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
228        D_coef(ig,jk-1,indice)=dz1(jk-1)*z1s
229     ENDDO
230  ENDDO
231
232!-----------------------------------------------------------------------
233! 2)
234! Computation of the soil temperatures using the Cgrd and Dgrd
235! coefficient computed above
236!
237!-----------------------------------------------------------------------
238
239!    Surface temperature
240  DO ig=1,knon
241     ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/  &
242          (lambda*(1.-D_coef(ig,1,indice))+1.)
243  ENDDO
244 
245!   Other temperatures
246  DO jk=1,nsoilmx-1
247     DO ig=1,knon
248        ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) &
249             *ptsoil(ig,jk)
250     ENDDO
251  ENDDO
252
253  IF (indice == is_sic) THEN
254     DO ig = 1 , knon
255        ptsoil(ig,nsoilmx) = RTT - 1.8
256     END DO
257  ENDIF
258
259!-----------------------------------------------------------------------
260! 3)
261! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil
262! temperature
263!-----------------------------------------------------------------------
264  DO ig=1,knon
265     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
266     C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
267     D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s
268  ENDDO
269 
270  DO jk=nsoilmx-1,2,-1
271     DO ig=1,knon
272        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
273             *(1.-D_coef(ig,jk,indice)))
274        C_coef(ig,jk-1,indice) = &
275             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
276        D_coef(ig,jk-1,indice) = dz1(jk-1)*z1s
277     ENDDO
278  ENDDO
279
280!-----------------------------------------------------------------------
281! 4)
282! Computation of the surface diffusive flux from ground and
283! calorific capacity of the ground
284!-----------------------------------------------------------------------
285  DO ig=1,knon
286     pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* &
287          (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1))
288     pcapcal(ig)  = ztherm_i(ig)* &
289          (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1))
290     z1s = lambda*(1.-D_coef(ig,1,indice))+1.
291     pcapcal(ig)  = pcapcal(ig)/z1s
292     pfluxgrd(ig) = pfluxgrd(ig) &
293          + pcapcal(ig) * (ptsoil(ig,1) * z1s &
294          - lambda * C_coef(ig,1,indice) &
295          - ptsrf(ig)) &
296          /ptimestep
297  ENDDO
298   
299END SUBROUTINE soil
Note: See TracBrowser for help on using the repository browser.