SURFEX v8.1
General documentation of Surfex
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,PSOC,K,NP,PFRACSOC,PWSAT,PWFC,PWWILT,KPATCH)
7 ! ########################################################################
8 !
9 !!**** *ISBA_SOC_PARAMETERS*
10 !!
11 !! PURPOSE
12 !! -------
13 !
14 ! ISBA parameterizations for soil thermal and hydraulic properties
15 ! are modified to accommodate both mineral and organic carbon soils
16 ! according to observations from Boelter (1969).
17 ! Disctinction is made for Fibric soil (76.8038 % of fiber content)
18 ! and Sapric soil (21.7815 % of fiber content)
19 !
20 !!** METHOD
21 !! ------
22 !
23 ! Direct calculation
24 !
25 !! EXTERNAL
26 !! --------
27 !
28 ! None
29 !!
30 !! IMPLICIT ARGUMENTS
31 !! ------------------
32 !!
33 !!
34 !! REFERENCE
35 !! ---------
36 !!
37 !! AUTHOR
38 !! ------
39 !! B. Decharme
40 !!
41 !! MODIFICATIONS
42 !! -------------
43 !! Original 10/12/11
44 !! (B. Decharme) 04/2013 ksat anisotropy factor
45 !-------------------------------------------------------------------------------
46 !
47 USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_np_t
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) :: PSOC
64 !
65 TYPE(isba_k_t), INTENT(INOUT) :: K
66 TYPE(isba_np_t), INTENT(INOUT) :: NP
67 !
68 REAL, DIMENSION(:,:), INTENT(OUT) :: PFRACSOC
69 INTEGER, INTENT(IN) :: KPATCH
70 !
71 REAL, DIMENSION(:,:), INTENT(INOUT) :: PWSAT
72 REAL, DIMENSION(:,:), INTENT(INOUT) :: PWFC
73 REAL, DIMENSION(:,:), INTENT(INOUT) :: PWWILT
74 !
75 !* 0.2 declarations of local parameter
76 !
77 TYPE(isba_p_t), POINTER :: PK
78 !
79 REAL, DIMENSION(2), PARAMETER :: ZCONDSAT = (/24.192,0.00864/) !Peatland hydraulic conductivity (m/day)
80  !from Letts et al. (2000)
81 !
82 REAL, DIMENSION(2), PARAMETER :: ZBCOEF = (/2.7,12.0/) !Peatland b coef (-)
83  !from Letts et al. (2000)
84 !
85 REAL, DIMENSION(2), PARAMETER :: ZMPOTSAT = (/-0.0103,-0.0101/) !Peatland matric potential (m)
86  !from Letts et al. (2000)
87 !
88 REAL, DIMENSION(2), PARAMETER :: ZWSAT = (/0.930,0.845/) !Peatland porosity (-)
89  !from Boelter (1969) PTF for
90  !Fibric soil = 76.8038 % of fiber content
91  !Sapric soil = 21.7815 % of fiber content
92 
93 REAL, DIMENSION(2), PARAMETER :: ZWFC = (/0.369,0.719/) !Peatland field capacity (-)
94  !Water potential at -0.1 bar given by
95  !Boelter (1969) PTF
96 !
97 REAL, DIMENSION(2), PARAMETER :: ZWWILT = (/0.073,0.222/) !Peatland wilting point (-)
98  !Water potential at -15 bar given by
99  !Boelter (1969) PTF
100 !
101 REAL, DIMENSION(2), PARAMETER :: ZWD0 = (/0.212,0.716/) !Peatland Topmodel D0 water equivalent (-)
102  !using hydro cond at 0.1mm/days
103 !
104 REAL, DIMENSION(2), PARAMETER :: ZANISO = (/2.0,48.0/) !Peatland ksat anisotropy factor (-)
105 !
106 !
107 !HWSD data profile
108 REAL, PARAMETER :: ZDGHWSD_TOP = 0.3
109 REAL, PARAMETER :: ZDGHWSD_SUB = 1.0
110 REAL, PARAMETER :: ZDGHWSD_INF = 1000.
111 !
112 !* 0.3 declarations of local variables
113 !
114 REAL, DIMENSION(SIZE(PWSAT,1)) :: ZPEAT_PROFILE, ZMOSS_DEPTH
115 !
116 REAL, DIMENSION(SIZE(PWSAT,1)) :: ZMASK, ZRHO_TOP, ZRHO_SUB, ZRHO_INF
117 !
118 REAL, DIMENSION(SIZE(PWSAT,1),SIZE(NP%AL(1)%XDG,2)) :: ZDG_SOIL, ZDZG_SOIL, ZRHO_SOC, ZMID_SOIL
119 !
120 REAL, DIMENSION(SIZE(PWSAT,1),SIZE(NP%AL(1)%XDG,2)) :: ZPEAT_BCOEF,ZPEAT_MPOTSAT,&
121  ZPEAT_WSAT,ZPEAT_WFC, &
122  ZPEAT_WWILT,ZPEAT_WD0, &
123  ZPEAT_ANISO, ZPEAT_RHO
124 !
125 REAL, DIMENSION(SIZE(PWSAT,1),SIZE(NP%AL(1)%XDG,2),KPATCH) :: ZPEAT_CONDSAT, ZMID_CONDSAT
126 !
127 REAL, DIMENSION(SIZE(PWSAT,1)) :: ZREFDEPTH,ZF_BCOEF,ZF_MPOTSAT, &
128  ZLOG_MOSS,ZLOG_PEAT_DEPTH , &
129  ZF_WSAT,ZF_CONDSAT,ZF_WFC, &
130  ZF_WWILT, ZF_WD0, ZF_ANISO
131 
132 REAL :: ZA, ZB, ZLOG1, ZLOG2, ZMOSS_DENSITY, &
133  ZTOP, ZSUB, ZFTOP, ZFSUB
134 !
135 REAL, DIMENSION(2) :: ZLOG_CONDSAT,ZLOG_BCOEF,ZLOG_MPOTSAT, &
136  ZLOG_WSAT,ZLOG_WFC,ZLOG_WWILT,ZLOG_WD0,&
137  ZLOG_ANISO
138 !
139 INTEGER :: INI, INL, INP, JI, JL, JP, IMASK
140 !
141 REAL(KIND=JPRB) :: ZHOOK_HANDLE
142 !
143 !-------------------------------------------------------------------------------
144 !
145 IF (lhook) CALL dr_hook('ISBA_SOC_PARAMETERS',0,zhook_handle)
146 !
147 ini = SIZE(pwsat,1)
148 inl = SIZE(np%AL(1)%XDG,2)
149 inp = kpatch
150 !
151 zmask(:) = 0.0
152 zrho_top(:) = 0.0
153 zrho_sub(:) = 0.0
154 zrho_inf(:) = 0.0
155 !
156 zrho_soc(:,:) = 0.0
157 !
158 zpeat_rho(:,: )=0.0
159 zpeat_bcoef(:,: )=0.0
160 zpeat_mpotsat(:,: )=0.0
161 zpeat_wsat(:,: )=0.0
162 zpeat_wfc(:,: )=0.0
163 zpeat_wwilt(:,: )=0.0
164 zpeat_wd0(:,: )=0.0
165 zpeat_aniso(:,: )=0.0
166 zpeat_condsat(:,:,:)=0.0
167 !
168 pfracsoc(:,:)=xundef
169 !
170 !-------------------------------------------------------------------------------
171 !
172 zdg_soil(:,:) = 0.0
173 DO jp = 1,inp
174  pk => np%AL(jp)
175  DO ji=1,pk%NSIZE_P
176  imask = pk%NR_P(ji)
177  IF(pk%XPATCH(ji)>0.0)THEN
178  zmask(imask)=zmask(imask)+pk%XPATCH(ji)
179  DO jl=1,inl
180  zdg_soil(imask,jl)= zdg_soil(imask,jl) + pk%XDG(ji,jl)*pk%XPATCH(ji)
181  ENDDO
182  ENDIF
183  ENDDO
184 ENDDO
185 !
186 DO jl = 1,inl
187  WHERE (zmask(:)/=0.) zdg_soil(:,jl) = zdg_soil(:,jl)/zmask(:)
188 ENDDO
189 !
190 zdzg_soil(:,1)=zdg_soil(:,1)
191 DO jl=2,inl
192  DO ji=1,ini
193  zdzg_soil(ji,jl)=zdg_soil(ji,jl)-zdg_soil(ji,jl-1)
194  ENDDO
195 ENDDO
196 !
197 zmid_soil(:,1)=0.5*zdg_soil(:,1)
198 DO jl=2,inl
199  DO ji=1,ini
200  zmid_soil(ji,jl)=0.5*(zdg_soil(ji,jl)+zdg_soil(ji,jl-1))
201  ENDDO
202 ENDDO
203 !
204 DO jp=1,inp
205  pk => np%AL(jp)
206  DO jl=1,inl
207  DO ji=1,pk%NSIZE_P
208  imask = pk%NR_P(ji)
209  IF(pk%XPATCH(ji)/=xundef)THEN
210  zmid_condsat(imask,jl,jp)=zmid_soil(imask,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 (ji>np%AL(jp)%NSIZE_P) cycle
324  imask = np%AL(jp)%NR_P(ji)
325  IF(np%AL(jp)%XPATCH(ji)/=xundef)THEN
326  zrefdepth(imask)=min(zpeat_profile(imask),max(zmoss_depth(imask),zmid_condsat(imask,jl,jp)))
327  zrefdepth(imask)=log(zrefdepth(imask))-zlog_moss(imask)
328  zpeat_condsat(imask,jl,jp)=zcondsat(1)*exp(zf_condsat(imask)*zrefdepth(imask))/xday
329  ENDIF
330  ENDDO
331 !
332  ENDIF
333  ENDDO
334 ENDDO
335 !
336 !-------------------------------------------------------------------------------
337 !
338 DO jl=1,inl
339  DO ji=1,ini
340  IF(zmask(ji)>0.0)THEN
341 ! Soil organic carbon fraction
342  pfracsoc(ji,jl) = min(1.0,zrho_soc(ji,jl)/zpeat_rho(ji,jl))
343 ! New soil thermal properties
344  k%XHCAPSOIL(ji,jl) = (1.0-pfracsoc(ji,jl))*k%XHCAPSOIL(ji,jl) + pfracsoc(ji,jl)*xomrho*xomsph
345  k%XCONDDRY (ji,jl) = (k%XCONDDRY(ji,jl)**(1.0-pfracsoc(ji,jl))) * (xomconddry**pfracsoc(ji,jl))
346  k%XCONDSLD (ji,jl) = (k%XCONDSLD(ji,jl)**(1.0-pfracsoc(ji,jl))) * (xomcondsld**pfracsoc(ji,jl))
347 ! New soil hydraulic properties
348  k%XBCOEF (ji,jl) = (1.0-pfracsoc(ji,jl))*k%XBCOEF (ji,jl) + pfracsoc(ji,jl)*zpeat_bcoef(ji,jl)
349  k%XMPOTSAT(ji,jl) = (1.0-pfracsoc(ji,jl))*k%XMPOTSAT(ji,jl) + pfracsoc(ji,jl)*zpeat_mpotsat(ji,jl)
350  pwsat(ji,jl) = (1.0-pfracsoc(ji,jl))*pwsat(ji,jl) + pfracsoc(ji,jl)*zpeat_wsat(ji,jl)
351  pwfc(ji,jl) = (1.0-pfracsoc(ji,jl))*pwfc(ji,jl) + pfracsoc(ji,jl)*zpeat_wfc(ji,jl)
352  pwwilt(ji,jl) = (1.0-pfracsoc(ji,jl))*pwwilt(ji,jl) + pfracsoc(ji,jl)*zpeat_wwilt(ji,jl)
353  ENDIF
354  ENDDO
355 ENDDO
356 !
357 DO jp=1,inp
358  pk => np%AL(jp)
359  DO jl=1,inl
360  DO ji=1,pk%NSIZE_P
361  imask = pk%NR_P(ji)
362  IF(pk%XPATCH(ji)/=xundef .AND. zmask(imask)>0.0)THEN
363  pk%XCONDSAT (ji,jl) = pk%XCONDSAT(ji,jl)**(1.0-pfracsoc(imask,jl)) * &
364  zpeat_condsat(imask,jl,jp)**pfracsoc(imask,jl)
365  ENDIF
366  ENDDO
367  ENDDO
368 ENDDO
369 !
370 IF(hrunoff=='SGH')THEN
371  DO jl=1,inl
372  DO ji=1,ini
373  IF(zmask(ji)>0.0)THEN
374  k%XWD0 (ji,jl) = (1.0-pfracsoc(ji,jl))*k%XWD0 (ji,jl) + pfracsoc(ji,jl)*zpeat_wd0(ji,jl)
375  k%XKANISO(ji,jl) = (1.0-pfracsoc(ji,jl))*k%XKANISO(ji,jl) + pfracsoc(ji,jl)*zpeat_aniso(ji,jl)
376  ENDIF
377  ENDDO
378  ENDDO
379 ENDIF
380 !
381 !-------------------------------------------------------------------------------
382 !
383 IF (lhook) CALL dr_hook('ISBA_SOC_PARAMETERS',1,zhook_handle)
384 !
385 END SUBROUTINE isba_soc_parameters
386 
387 
388 
389 
390 
modd_surf_par
Definition: modd_surf_par.F90:6
modd_isba_n
Definition: modd_isban.F90:6
modd_isba_n::isba_k_t
Definition: modd_isban.F90:139
yomhook::dr_hook
Definition: yomhook.F90:20
isba_soc_parameters
subroutine isba_soc_parameters(HRUNOFF, PSOC, K, NP, PFRACSOC, PWSAT, PWFC, PWWILT, KPATCH)
Definition: isba_soc_parameters.F90:7
modd_surf_par::xundef
real, parameter xundef
Definition: modd_surf_par.F90:43
parkind1::jprb
integer, parameter jprb
Definition: parkind1.F90:32
modd_csts::xday
real, save xday
Definition: modd_csts.F90:45
yomhook::lhook
logical lhook
Definition: yomhook.F90:15
modd_csts
Definition: modd_csts.F90:6
modd_isba_n::isba_np_t
Definition: modd_isban.F90:479
modd_isba_n::isba_p_t
Definition: modd_isban.F90:219
Generated on Tue Jan 16 2018 16:23:29 for SURFEX v8.1 by   doxygen 1.8.13