SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
isba_soc_parameters.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 ! #########################
6 SUBROUTINE isba_soc_parameters (HRUNOFF,PPATCH,PDG,PSOC,PBCOEF,PMPOTSAT, &
7  pcondsat,pwsat,phcapsoil,pconddry,pcondsld,&
8  pwfc,pwwilt,pwd0,paniso,pfracsoc )
9 ! ########################################################################
10 !
11 !!**** *ISBA_SOC_PARAMETERS*
12 !!
13 !! PURPOSE
14 !! -------
15 !
16 ! ISBA parameterizations for soil thermal and hydraulic properties
17 ! are modified to accommodate both mineral and organic carbon soils
18 ! according to observations from Boelter (1969).
19 ! Disctinction is made for Fibric soil (76.8038 % of fiber content)
20 ! and Sapric soil (21.7815 % of fiber content)
21 !
22 !!** METHOD
23 !! ------
24 !
25 ! Direct calculation
26 !
27 !! EXTERNAL
28 !! --------
29 !
30 ! None
31 !!
32 !! IMPLICIT ARGUMENTS
33 !! ------------------
34 !!
35 !!
36 !! REFERENCE
37 !! ---------
38 !!
39 !! AUTHOR
40 !! ------
41 !! B. Decharme
42 !!
43 !! MODIFICATIONS
44 !! -------------
45 !! Original 10/12/11
46 !! (B. Decharme) 04/2013 ksat anisotropy factor
47 !-------------------------------------------------------------------------------
48 !
49 USE modd_surf_par, ONLY : xundef
50 USE modd_csts, ONLY : xday
51 USE modd_isba_par, ONLY : xomrho, xomsph, xomconddry, xomcondsld
52 !
53 USE yomhook ,ONLY : lhook, dr_hook
54 USE parkind1 ,ONLY : jprb
55 
56 !* 0.1 declarations of arguments
57 !
58 !
59 IMPLICIT NONE
60 !
61  CHARACTER(LEN=4), INTENT(IN) :: hrunoff
62 !
63 REAL, DIMENSION(:,:,:),INTENT(IN) :: pdg
64 !
65 REAL, DIMENSION(:,:), INTENT(IN) :: ppatch
66 !
67 REAL, DIMENSION(:,:), INTENT(IN) :: psoc
68 !
69 REAL, DIMENSION(:,:,:),INTENT(INOUT) :: pcondsat
70 !
71 REAL, DIMENSION(:,:), INTENT(INOUT) :: pbcoef,pmpotsat, &
72  phcapsoil,pconddry, &
73  pcondsld
74 !
75 REAL, DIMENSION(:,:), INTENT(INOUT) :: pwsat,pwfc,pwwilt,pwd0
76 !
77 REAL, DIMENSION(:,:), INTENT(INOUT) :: paniso
78 !
79 REAL, DIMENSION(:,:), INTENT(OUT) :: pfracsoc
80 !
81 !* 0.2 declarations of local parameter
82 !
83 REAL, DIMENSION(2), PARAMETER :: zcondsat = (/24.192,0.00864/) !Peatland hydraulic conductivity (m/day)
84  !from Letts et al. (2000)
85 !
86 REAL, DIMENSION(2), PARAMETER :: zbcoef = (/2.7,12.0/) !Peatland b coef (-)
87  !from Letts et al. (2000)
88 !
89 REAL, DIMENSION(2), PARAMETER :: zmpotsat = (/-0.0103,-0.0101/) !Peatland matric potential (m)
90  !from Letts et al. (2000)
91 !
92 REAL, DIMENSION(2), PARAMETER :: zwsat = (/0.930,0.845/) !Peatland porosity (-)
93  !from Boelter (1969) PTF for
94  !Fibric soil = 76.8038 % of fiber content
95  !Sapric soil = 21.7815 % of fiber content
96 
97 REAL, DIMENSION(2), PARAMETER :: zwfc = (/0.369,0.719/) !Peatland field capacity (-)
98  !Water potential at -0.1 bar given by
99  !Boelter (1969) PTF
100 !
101 REAL, DIMENSION(2), PARAMETER :: zwwilt = (/0.073,0.222/) !Peatland wilting point (-)
102  !Water potential at -15 bar given by
103  !Boelter (1969) PTF
104 !
105 REAL, DIMENSION(2), PARAMETER :: zwd0 = (/0.212,0.716/) !Peatland Topmodel D0 water equivalent (-)
106  !using hydro cond at 0.1mm/days
107 !
108 REAL, DIMENSION(2), PARAMETER :: zaniso = (/2.0,48.0/) !Peatland ksat anisotropy factor (-)
109 !
110 !
111 !HWSD data profile
112 REAL, PARAMETER :: zdghwsd_top = 0.3
113 REAL, PARAMETER :: zdghwsd_sub = 1.0
114 REAL, PARAMETER :: zdghwsd_inf = 1000.
115 !
116 !* 0.3 declarations of local variables
117 !
118 REAL, DIMENSION(SIZE(PDG,1)) :: zpeat_profile, zmoss_depth
119 !
120 REAL, DIMENSION(SIZE(PDG,1)) :: zmask, zrho_top, zrho_sub, zrho_inf
121 !
122 REAL, DIMENSION(SIZE(PDG,1),SIZE(PDG,2)) :: zdg_soil, zdzg_soil, zrho_soc, zmid_soil
123 !
124 REAL, DIMENSION(SIZE(PDG,1),SIZE(PDG,2)) :: zpeat_bcoef,zpeat_mpotsat,&
125  zpeat_wsat,zpeat_wfc, &
126  zpeat_wwilt,zpeat_wd0, &
127  zpeat_aniso, zpeat_rho
128 !
129 REAL, DIMENSION(SIZE(PDG,1),SIZE(PDG,2),SIZE(PDG,3)) :: zpeat_condsat, zmid_condsat
130 !
131 REAL, DIMENSION(SIZE(PDG,1)) :: zrefdepth,zf_bcoef,zf_mpotsat, &
132  zlog_moss,zlog_peat_depth , &
133  zf_wsat,zf_condsat,zf_wfc, &
134  zf_wwilt, zf_wd0, zf_aniso
135 
136 REAL :: za, zb, zlog1, zlog2, zmoss_density, &
137  ztop, zsub, zftop, zfsub
138 !
139 REAL, DIMENSION(2) :: zlog_condsat,zlog_bcoef,zlog_mpotsat, &
140  zlog_wsat,zlog_wfc,zlog_wwilt,zlog_wd0,&
141  zlog_aniso
142 !
143 INTEGER :: ini, inl, inp, ji, jl, jp
144 !
145 REAL(KIND=JPRB) :: zhook_handle
146 !
147 !-------------------------------------------------------------------------------
148 !
149 IF (lhook) CALL dr_hook('ISBA_SOC_PARAMETERS',0,zhook_handle)
150 !
151 ini=SIZE(pdg,1)
152 inl=SIZE(pdg,2)
153 inp=SIZE(pdg,3)
154 !
155 zmask(:) = 0.0
156 zrho_top(:) = 0.0
157 zrho_sub(:) = 0.0
158 zrho_inf(:) = 0.0
159 !
160 zdg_soil(:,:) = 0.0
161 zrho_soc(:,:) = 0.0
162 !
163 zpeat_rho(:,: )=0.0
164 zpeat_bcoef(:,: )=0.0
165 zpeat_mpotsat(:,: )=0.0
166 zpeat_wsat(:,: )=0.0
167 zpeat_wfc(:,: )=0.0
168 zpeat_wwilt(:,: )=0.0
169 zpeat_wd0(:,: )=0.0
170 zpeat_aniso(:,: )=0.0
171 zpeat_condsat(:,:,:)=0.0
172 !
173 pfracsoc(:,:)=xundef
174 !
175 !-------------------------------------------------------------------------------
176 !
177 DO jp=1,inp
178  DO ji=1,ini
179  zmask(ji)=zmask(ji)+ppatch(ji,jp)
180  ENDDO
181 ENDDO
182 !
183 DO jl=1,inl
184  DO ji=1,ini
185  IF(zmask(ji)>0.0)THEN
186  zdg_soil(ji,jl)=sum(pdg(ji,jl,:)*ppatch(ji,:),ppatch(ji,:)>0.0) &
187  /sum(ppatch(ji,:),ppatch(ji,:)>0.0)
188  ENDIF
189  ENDDO
190 ENDDO
191 !
192 zdzg_soil(:,1)=zdg_soil(:,1)
193 DO jl=2,inl
194  DO ji=1,ini
195  zdzg_soil(ji,jl)=zdg_soil(ji,jl)-zdg_soil(ji,jl-1)
196  ENDDO
197 ENDDO
198 !
199 zmid_soil(:,1)=0.5*zdg_soil(:,1)
200 DO jl=2,inl
201  DO ji=1,ini
202  zmid_soil(ji,jl)=0.5*(zdg_soil(ji,jl)+zdg_soil(ji,jl-1))
203  ENDDO
204 ENDDO
205 !
206 DO jp=1,inp
207  DO jl=1,inl
208  DO ji=1,ini
209  IF(ppatch(ji,jp)/=xundef)THEN
210  zmid_condsat(ji,jl,jp)=zmid_soil(ji,jl)
211  ENDIF
212  ENDDO
213  ENDDO
214 ENDDO
215 !
216 !-------------------------------------------------------------------------------
217 !
218 ! Compute the SOC density distribution (kg.m-3)
219 !
220 zlog1=log(zdghwsd_top/zdghwsd_sub)
221 zlog2=log(zdghwsd_inf/zdghwsd_sub)
222 DO ji=1,ini
223  IF(zmask(ji)>0.0)THEN
224  zrho_top(ji) = psoc(ji,1)/zdghwsd_top
225  zrho_sub(ji) = psoc(ji,2)/(zdghwsd_sub-zdghwsd_top)
226  IF(zrho_top(ji)>zrho_sub(ji))THEN
227  zb = log(psoc(ji,1)/(psoc(ji,1)+psoc(ji,2)))/zlog1
228  za = (psoc(ji,1)+psoc(ji,2))/(zdghwsd_inf-zdghwsd_sub)
229  zrho_inf(ji) = za*(exp(zb*zlog2)-1.0)
230  ELSE
231  zrho_inf(ji) = zrho_sub(ji)
232  ENDIF
233  ENDIF
234 ENDDO
235 !
236 ! Compute the SOC density distribution (kg.m-3)
237 !
238 !
239 DO ji=1,ini
240  ztop=0.0
241  zsub=0.0
242  DO jl=1,inl
243  ztop=zsub
244  zsub=zsub+zdzg_soil(ji,jl)
245  IF(zsub<=zdghwsd_top)THEN
246  zrho_soc(ji,jl)=zrho_top(ji)
247  ELSEIF(ztop>=zdghwsd_top.AND.zsub<=zdghwsd_sub)THEN
248  zrho_soc(ji,jl)=zrho_sub(ji)
249  ELSEIF(ztop>=zdghwsd_sub)THEN
250  zrho_soc(ji,jl)=zrho_inf(ji)
251  ELSEIF(ztop<zdghwsd_top.AND.zsub>zdghwsd_top)THEN
252  zftop=min(1.0,max(0.0,zdghwsd_top-ztop))/(zsub-ztop)
253  zfsub=min(1.0,max(0.0,zsub-zdghwsd_top))/(zsub-ztop)
254  zrho_soc(ji,jl)=zftop*zrho_top(ji)+zfsub*zrho_sub(ji)
255  ELSEIF(ztop<zdghwsd_sub.AND.zsub>zdghwsd_sub)THEN
256  zftop=min(1.0,max(0.0,zdghwsd_sub-ztop))/(zsub-ztop)
257  zfsub=min(1.0,max(0.0,zsub-zdghwsd_sub))/(zsub-ztop)
258  zrho_soc(ji,jl)=zftop*zrho_sub(ji)+zfsub*zrho_inf(ji)
259  ENDIF
260  ENDDO
261 ENDDO
262 !
263 !-------------------------------------------------------------------------------
264 !
265 ! Define the Peatland soil properties
266 !
267 zlog_condsat(:) = log(zcondsat(:))
268 zlog_bcoef(:) = log(zbcoef(:))
269 zlog_mpotsat(:) = log(-zmpotsat(:))
270 zlog_wsat(:) = log(zwsat(:))
271 zlog_wfc(:) = log(zwfc(:))
272 zlog_wwilt(:) = log(zwwilt(:))
273 zlog_wd0(:) = log(zwd0(:))
274 zlog_aniso(:) = log(zaniso(:))
275 !
276 zpeat_profile(:) = 1.0
277 !
278 zmoss_density=(1.0-zwsat(1))*xomrho
279 !
280 WHERE(zrho_top(:)<zmoss_density)
281  zmoss_depth(:) = 2.5e-3 ! => Small fibric soil at surface (<=> moss=2.5mm)
282 ELSEWHERE
283  zmoss_depth(:) = 0.01 ! => Fibric soil at surface (<=> moss=1cm)
284 ENDWHERE
285 !
286 WHERE(zmask(:)>0.0)
287 !
288  zlog_moss(:) = log(zmoss_depth(:))
289  zlog_peat_depth(:) = log(zpeat_profile(:))
290 !
291  zf_condsat(:) =(zlog_condsat(2)-zlog_condsat(1))/(zlog_peat_depth(:)-zlog_moss(:))
292  zf_bcoef(:) =(zlog_bcoef(2)-zlog_bcoef(1))/(zlog_peat_depth(:)-zlog_moss(:))
293  zf_mpotsat(:) =(zlog_mpotsat(2)-zlog_mpotsat(1))/(zlog_peat_depth(:)-zlog_moss(:))
294  zf_wsat(:) =(zlog_wsat(2)-zlog_wsat(1))/(zlog_peat_depth(:)-zlog_moss(:))
295  zf_wfc(:) =(zlog_wfc(2)-zlog_wfc(1))/(zlog_peat_depth(:)-zlog_moss(:))
296  zf_wwilt(:) =(zlog_wwilt(2)-zlog_wwilt(1))/(zlog_peat_depth(:)-zlog_moss(:))
297  zf_wd0(:) =(zlog_wd0(2)-zlog_wd0(1))/(zlog_peat_depth(:)-zlog_moss(:))
298  zf_aniso(:) =(zlog_aniso(2)-zlog_aniso(1))/(zlog_peat_depth(:)-zlog_moss(:))
299 !
300 ENDWHERE
301 !
302 !-------------------------------------------------------------------------------
303 !
304 ! Compute the Peatland soil properties profile
305 !
306 DO jl=1,inl
307  DO ji=1,ini
308  IF(zmask(ji)>0.0)THEN
309 !
310  zrefdepth(ji)=min(zpeat_profile(ji),max(zmoss_depth(ji),zmid_soil(ji,jl)))
311  zrefdepth(ji)=log(zrefdepth(ji))-zlog_moss(ji)
312  zpeat_mpotsat(ji,jl)=zmpotsat(1)*exp(zf_mpotsat(ji)*zrefdepth(ji))
313  zpeat_wsat(ji,jl)=zwsat(1)*exp(zf_wsat(ji)*zrefdepth(ji))
314  zpeat_bcoef(ji,jl)=zbcoef(1)*exp(zf_bcoef(ji)*zrefdepth(ji))
315  zpeat_wwilt(ji,jl)=zwwilt(1)*exp(zf_wwilt(ji)*zrefdepth(ji))
316  zpeat_wd0(ji,jl)=zwd0(1)*exp(zf_wd0(ji)*zrefdepth(ji))
317  zpeat_aniso(ji,jl)=zaniso(1)*exp(zf_aniso(ji)*zrefdepth(ji))
318  zpeat_wfc(ji,jl)=zwfc(1)*exp(zf_wfc(ji)*zrefdepth(ji))
319 !
320  zpeat_rho(ji,jl)=(1.0-zpeat_wsat(ji,jl))*xomrho
321 !
322  DO jp=1,inp
323  IF(ppatch(ji,jp)/=xundef)THEN
324  zrefdepth(ji)=min(zpeat_profile(ji),max(zmoss_depth(ji),zmid_condsat(ji,jl,jp)))
325  zrefdepth(ji)=log(zrefdepth(ji))-zlog_moss(ji)
326  zpeat_condsat(ji,jl,jp)=zcondsat(1)*exp(zf_condsat(ji)*zrefdepth(ji))/xday
327  ENDIF
328  ENDDO
329 !
330  ENDIF
331  ENDDO
332 ENDDO
333 !
334 !-------------------------------------------------------------------------------
335 !
336 DO jl=1,inl
337  DO ji=1,ini
338  IF(zmask(ji)>0.0)THEN
339 ! Soil organic carbon fraction
340  pfracsoc(ji,jl ) = min(1.0,zrho_soc(ji,jl)/zpeat_rho(ji,jl))
341 ! New soil thermal properties
342  phcapsoil(ji,jl ) = (1.0-pfracsoc(ji,jl))*phcapsoil(ji,jl) + pfracsoc(ji,jl)*xomrho*xomsph
343  pconddry(ji,jl ) = (pconddry(ji,jl)**(1.0-pfracsoc(ji,jl))) * (xomconddry**pfracsoc(ji,jl))
344  pcondsld(ji,jl ) = (pcondsld(ji,jl)**(1.0-pfracsoc(ji,jl))) * (xomcondsld**pfracsoc(ji,jl))
345 ! New soil hydraulic properties
346  pbcoef(ji,jl ) = (1.0-pfracsoc(ji,jl))*pbcoef(ji,jl) + pfracsoc(ji,jl)*zpeat_bcoef(ji,jl)
347  pmpotsat(ji,jl ) = (1.0-pfracsoc(ji,jl))*pmpotsat(ji,jl) + pfracsoc(ji,jl)*zpeat_mpotsat(ji,jl)
348  pwsat(ji,jl ) = (1.0-pfracsoc(ji,jl))*pwsat(ji,jl) + pfracsoc(ji,jl)*zpeat_wsat(ji,jl)
349  pwfc(ji,jl ) = (1.0-pfracsoc(ji,jl))*pwfc(ji,jl) + pfracsoc(ji,jl)*zpeat_wfc(ji,jl)
350  pwwilt(ji,jl ) = (1.0-pfracsoc(ji,jl))*pwwilt(ji,jl) + pfracsoc(ji,jl)*zpeat_wwilt(ji,jl)
351  DO jp=1,inp
352  IF(ppatch(ji,jp)/=xundef)THEN
353  pcondsat(ji,jl,jp) = pcondsat(ji,jl,jp)**(1.0-pfracsoc(ji,jl))*zpeat_condsat(ji,jl,jp)**pfracsoc(ji,jl)
354  ENDIF
355  ENDDO
356  ENDIF
357  ENDDO
358 ENDDO
359 !
360 IF(hrunoff=='SGH')THEN
361  DO jl=1,inl
362  DO ji=1,ini
363  IF(zmask(ji)>0.0)THEN
364  pwd0(ji,jl) = (1.0-pfracsoc(ji,jl))*pwd0(ji,jl) + pfracsoc(ji,jl)*zpeat_wd0(ji,jl)
365  paniso(ji,jl) = (1.0-pfracsoc(ji,jl))*paniso(ji,jl) + pfracsoc(ji,jl)*zpeat_aniso(ji,jl)
366  ENDIF
367  ENDDO
368  ENDDO
369 ENDIF
370 !
371 !-------------------------------------------------------------------------------
372 !
373 IF (lhook) CALL dr_hook('ISBA_SOC_PARAMETERS',1,zhook_handle)
374 !
375 END SUBROUTINE isba_soc_parameters
376 
377 
378 
379 
380 
subroutine isba_soc_parameters(HRUNOFF, PPATCH, PDG, PSOC, PBCOEF, PMPOTSAT, PCONDSAT, PWSAT, PHCAPSOIL, PCONDDRY, PCONDSLD, PWFC, PWWILT, PWD0, PANISO, PFRACSOC)