source: LMDZ6/trunk/libf/phylmd/soil_hetero.F90 @ 5627

Last change on this file since 5627 was 5627, checked in by amaison, 6 weeks ago

Representation of heterogeneous continental subsurfaces with parameter or flux aggregation in the simplified surface model (bucket) for 1D case studies.
File size: 8.5 KB
Line 
1!
2! $Header$
3!
4SUBROUTINE soil_hetero(ptimestep, indice, knon, snow, ptsrf, qsol, &
5     lon, lat, ptsoil, pcapcal, pfluxgrd, ztherm_i, conv_ratio)
6
7  USE YOMCST, ONLY: RTT, RPI
8  USE dimsoil_mod_h, ONLY: nsoilmx
9  USE comsoil_mod_h
10  USE compbl_mod_h
11  USE dimpft_mod_h
12  USE dimphy
13  USE mod_phys_lmdz_para
14  USE indice_sol_mod
15  USE print_control_mod, ONLY: lunout
16  USE phys_state_var_mod, ONLY: alpha_soil_tersrf, period_tersrf
17  USE surf_param_mod, ONLY: eff_surf_param
18
19  IMPLICIT NONE
20
21!=======================================================================
22!
23!   Auteur:  Frederic Hourdin     30/01/92
24!   -------
25!
26!   Object:  Computation of : the soil temperature evolution
27!   -------                   the surfacic heat capacity "Capcal"
28!                            the surface conduction flux pcapcal
29!
30!   Update: 2021/07 : soil thermal inertia, formerly a constant value,
31!   ------   can also be now a function of soil moisture (F Cheruy's idea)
32!            depending on iflag_inertie, read from physiq.def via conf_phys_m.F90
33!            ("Stage L3" Eve Rebouillat, with E Vignon, A Sima, F Cheruy)
34!           2025/04 : A. Maison, adapting the routine for heterogeneous continental sub-surfaces
35!
36!   Method: Implicit time integration
37!   -------
38!   Consecutive ground temperatures are related by:
39!           T(k+1) = C(k) + D(k)*T(k)  (*)
40!   The coefficients C and D are computed at the t-dt time-step.
41!   Routine structure:
42!   1) C and D coefficients are computed from the old temperature
43!   2) new temperatures are computed using (*)
44!   3) C and D coefficients are computed from the new temperature
45!      profile for the t+dt time-step
46!   4) the coefficients A and B are computed where the diffusive
47!      fluxes at the t+dt time-step is given by
48!             Fdiff = A + B Ts(t+dt)
49!      or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
50!             with F0 = A + B (Ts(t))
51!                 Capcal = B*dt
52!           
53!   Interface:
54!   ----------
55!
56!   Arguments:
57!   ----------
58!   ptimestep            physical timestep (s)
59!   indice               sub-surface index
60!   snow(klon)           snow
61!   ptsrf(klon)          surface temperature at time-step t (K)
62!   qsol(klon)           soil moisture (kg/m2 or mm)
63!   lon(klon)            longitude in radian
64!   lat(klon)            latitude in radian
65!   ptsoil(klon,nsoilmx) temperature inside the ground (K)
66!   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
67!   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
68!   ztherm_i(klon)       soil thermal inertia (J.m-2.K.s-1/2)
69!   conv_ratio(klon)     ratio to convert soil depths in meters (-)
70!   
71!=======================================================================
72!-----------------------------------------------------------------------
73! Arguments
74! ---------
75  REAL, INTENT(IN)                     :: ptimestep
76  INTEGER, INTENT(IN)                  :: indice, knon !, knindex
77  REAL, DIMENSION(klon), INTENT(IN)    :: snow
78  REAL, DIMENSION(klon), INTENT(IN)    :: ptsrf
79  REAL, DIMENSION(klon), INTENT(IN)    :: qsol
80  REAL, DIMENSION(klon), INTENT(IN)    :: lon
81  REAL, DIMENSION(klon), INTENT(IN)    :: lat
82  REAL, DIMENSION(klon), INTENT(IN)    :: ztherm_i
83  REAL, DIMENSION(klon), INTENT(IN)    :: conv_ratio
84
85  REAL, DIMENSION(klon,nsoilmx), INTENT(INOUT) :: ptsoil
86  REAL, DIMENSION(klon), INTENT(OUT)           :: pcapcal
87  REAL, DIMENSION(klon), INTENT(OUT)           :: pfluxgrd
88
89!-----------------------------------------------------------------------
90! Local variables
91! ---------------
92  INTEGER                             :: ig, jk, ierr
93  REAL, DIMENSION(nsoilmx)            :: zdz2
94  REAL                                :: z1s
95  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: C_coef, D_coef
96
97! Local saved variables
98! ---------------------
99  REAL, SAVE                     :: lambda
100!$OMP THREADPRIVATE(lambda)
101  REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2
102!$OMP THREADPRIVATE(dz1,dz2)
103  LOGICAL, SAVE                  :: firstcall=.TRUE.
104!$OMP THREADPRIVATE(firstcall)
105   
106!-----------------------------------------------------------------------
107!   Depthts:
108!   --------
109  REAL fz,rk,fz1,rk1,rk2
110  fz(rk)=fz1*(alpha_soil_tersrf**rk-1.)/(alpha_soil_tersrf-1.)
111!
112!-----------------------------------------------------------------------
113! Calculation of some constants
114! NB! These constants do not depend on the sub-surfaces
115!-----------------------------------------------------------------------
116
117  IF (firstcall) THEN
118!-----------------------------------------------------------------------
119!   ground levels
120!   grnd=z/l where l is the skin depth of the diurnal cycle:
121!-----------------------------------------------------------------------
122!   la premiere couche represente un dixieme de cycle diurne
123     fz1=SQRT(period_tersrf/RPI)
124     
125     DO jk=1,nsoilmx
126        rk1=jk
127        rk2=jk-1
128        dz2(jk)=fz(rk1)-fz(rk2)
129     ENDDO
130     DO jk=1,nsoilmx-1
131        rk1=jk+.5
132        rk2=jk-.5
133        dz1(jk)=1./(fz(rk1)-fz(rk2))
134     ENDDO
135     lambda=fz(.5)*dz1(1)
136     WRITE(lunout,*) 'surface index:', indice
137     WRITE(lunout,*)'full layers, intermediate layers (seconds)'
138     DO jk=1,nsoilmx
139        rk=jk
140        rk1=jk+.5
141        rk2=jk-.5
142        WRITE(lunout,*)'fz=', &
143             fz(rk1)*fz(rk2)*RPI,fz(rk)*fz(rk)*RPI
144     ENDDO
145     WRITE(lunout,*)'full layers, intermediate layers (meters)'
146     DO jk=1,nsoilmx
147        rk=jk
148        rk2=jk-.5
149        WRITE(lunout,*)'fz=', &
150             fz(rk2)*conv_ratio, fz(rk)*conv_ratio
151     ENDDO
152     firstcall =.FALSE.
153  END IF
154
155!-----------------------------------------------------------------------
156! 1)
157! Calculation of Cgrf and Dgrd coefficients using soil temperature from
158! previous time step.
159!
160! These variables are recalculated on the local compressed grid instead
161! of saved in restart file.
162!-----------------------------------------------------------------------
163  DO jk=1,nsoilmx
164     zdz2(jk)=dz2(jk)/ptimestep
165  ENDDO
166 
167  DO ig=1,knon
168     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
169     C_coef(ig,nsoilmx-1,indice)= &
170          zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
171     D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s
172  ENDDO
173 
174  DO jk=nsoilmx-1,2,-1
175     DO ig=1,knon
176        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
177             *(1.-D_coef(ig,jk,indice)))
178        C_coef(ig,jk-1,indice)= &
179             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
180        D_coef(ig,jk-1,indice)=dz1(jk-1)*z1s
181     ENDDO
182  ENDDO
183
184!-----------------------------------------------------------------------
185! 2)
186! Computation of the soil temperatures using the Cgrd and Dgrd
187! coefficient computed above
188!
189!-----------------------------------------------------------------------
190
191!    Surface temperature
192  DO ig=1,knon
193     ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/  &
194          (lambda*(1.-D_coef(ig,1,indice))+1.)
195  ENDDO
196 
197!   Other temperatures
198  DO jk=1,nsoilmx-1
199     DO ig=1,knon
200        ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) &
201             *ptsoil(ig,jk)
202     ENDDO
203  ENDDO
204
205  IF (indice == is_sic) THEN
206     DO ig = 1 , knon
207        ptsoil(ig,nsoilmx) = RTT - 1.8
208     END DO
209  ENDIF
210
211!-----------------------------------------------------------------------
212! 3)
213! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil
214! temperature
215!-----------------------------------------------------------------------
216  DO ig=1,knon
217     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
218     C_coef(ig,nsoilmx-1,indice) = 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! 4)
234! Computation of the surface diffusive flux from ground and
235! calorific capacity of the ground
236!-----------------------------------------------------------------------
237  DO ig=1,knon
238     pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* &
239          (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1))
240     pcapcal(ig)  = ztherm_i(ig)* &
241          (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1))
242     z1s = lambda*(1.-D_coef(ig,1,indice))+1.
243     pcapcal(ig)  = pcapcal(ig)/z1s
244     pfluxgrd(ig) = pfluxgrd(ig) &
245          + pcapcal(ig) * (ptsoil(ig,1) * z1s &
246          - lambda * C_coef(ig,1,indice) &
247          - ptsrf(ig)) &
248          /ptimestep
249  ENDDO
250   
251END SUBROUTINE soil_hetero
Note: See TracBrowser for help on using the repository browser.