Skip to content

Commit 46ea27e

Browse files
jjokellajjokella
committed
Soil Hydraulic Parameter Perturbations for eCLM
1 parent 860a68c commit 46ea27e

1 file changed

Lines changed: 20 additions & 24 deletions

File tree

‎src/clm5/biogeophys/SoilStateInitTimeConstMod.F90‎

Lines changed: 20 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -150,7 +150,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
150150
real(r8) ,pointer :: gti (:) ! read in - fmax
151151
real(r8) ,pointer :: sand3d (:,:) ! read in - soil texture: percent sand (needs to be a pointer for use in ncdio)
152152
real(r8) ,pointer :: clay3d (:,:) ! read in - soil texture: percent clay (needs to be a pointer for use in ncdio)
153-
#ifdef USE_PDAF
153+
! SHP start
154154
real(r8) ,pointer :: psis_sat (:,:) ! read in - soil parameter: sucsat (needs to be a pointer for use in ncdio)
155155
real(r8) ,pointer :: shape_param (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio)
156156
real(r8) ,pointer :: thetas (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio)
@@ -160,16 +160,16 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
160160
real(r8) ,pointer :: shape_param_adj (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio)
161161
real(r8) ,pointer :: thetas_adj (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio)
162162
real(r8) ,pointer :: ks_adj (:,:) ! read in - soil parameter: xksat (needs to be a pointer for use in ncdio)
163-
#endif
163+
! SHP end
164164
real(r8) ,pointer :: organic3d (:,:) ! read in - organic matter: kg/m3 (needs to be a pointer for use in ncdio)
165165
character(len=256) :: locfn ! local filename
166166
integer :: ipedof
167167
integer :: begp, endp
168168
integer :: begc, endc
169169
integer :: begg, endg
170-
#ifdef USE_PDAF
170+
! SHP start
171171
logical :: parameters_in_file, parameters_in_file_adj
172-
#endif
172+
! SHP end
173173
!-----------------------------------------------------------------------
174174

175175
begp = bounds%begp; endp= bounds%endp
@@ -239,7 +239,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
239239

240240
allocate(sand3d(begg:endg,nlevsoifl))
241241
allocate(clay3d(begg:endg,nlevsoifl))
242-
#ifdef USE_PDAF
242+
! SHP start
243243
allocate(thetas(begg:endg,nlevsoifl))
244244
allocate(shape_param(begg:endg,nlevsoifl))
245245
allocate(psis_sat(begg:endg,nlevsoifl))
@@ -249,7 +249,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
249249
allocate(shape_param_adj(begg:endg,nlevgrnd))
250250
allocate(psis_sat_adj(begg:endg,nlevgrnd))
251251
allocate(ks_adj(begg:endg,nlevgrnd))
252-
#endif
252+
! SHP end
253253

254254
! Determine organic_max from parameter file
255255

@@ -287,7 +287,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
287287
call endrun(msg=' ERROR: PCT_CLAY NOT on surfdata file'//errMsg(sourcefile, __LINE__))
288288
end if
289289

290-
#ifdef USE_PDAF
290+
! SHP start
291291
! include option to also read hydraulic parameters from file. Keep it variable so that the code also works for surface files that were
292292
! generated without parameter perturbation and parameter as input variables
293293

@@ -340,7 +340,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
340340
parameters_in_file_adj = .False.
341341
end if
342342
end if
343-
#endif
343+
! SHP end
344344
do p = begp,endp
345345
g = patch%gridcell(p)
346346
if ( sand3d(g,1)+clay3d(g,1) == 0.0_r8 )then
@@ -533,7 +533,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
533533
call pedotransf(ipedof, sand, clay, &
534534
soilstate_inst%watsat_col(c,lev), soilstate_inst%bsw_col(c,lev), soilstate_inst%sucsat_col(c,lev), xksat)
535535

536-
#ifdef USE_PDAF
536+
! SHP start
537537
! if parameters are included in the file, watsat,... are overwritten with the values from there. If not, the pedotransfer
538538
! function is used
539539

@@ -552,7 +552,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
552552
xksat = ks(col%gridcell(c), nlevsoifl)
553553
end if
554554
end if
555-
#endif
555+
! SHP end
556556
om_watsat = max(0.93_r8 - 0.1_r8 *(zsoi(lev)/zsapric), 0.83_r8)
557557
om_b = min(2.7_r8 + 9.3_r8 *(zsoi(lev)/zsapric), 12.0_r8)
558558
om_sucsat = min(10.3_r8 - 0.2_r8 *(zsoi(lev)/zsapric), 10.1_r8)
@@ -561,15 +561,13 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
561561
soilstate_inst%bd_col(c,lev) = (1._r8 - soilstate_inst%watsat_col(c,lev))*2.7e3_r8
562562
soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac) * soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac
563563
tkm = (1._r8-om_frac) * (8.80_r8*sand+2.92_r8*clay)/(sand+clay)+om_tkm*om_frac ! W/(m K)
564-
#ifndef USE_PDAF
565-
soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b
566-
#else
564+
! SHP adapt start
567565
if (parameters_in_file) then
568566
soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b
569567
else
570-
soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b
568+
soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b
571569
end if
572-
#endif
570+
! SHP adapt end
573571
soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac) * soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac
574572
soilstate_inst%hksat_min_col(c,lev) = xksat
575573

@@ -593,15 +591,15 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
593591
end if
594592
soilstate_inst%hksat_col(c,lev) = uncon_frac*uncon_hksat + (perc_frac*om_frac)*om_hksat
595593

596-
#ifdef USE_PDAF
594+
! SHP start
597595
if (parameters_in_file_adj) then
598596
! Use values from the file for the soil layers
599597
soilstate_inst%watsat_col(c,lev) = thetas_adj(col%gridcell(c), lev)
600598
soilstate_inst%bsw_col(c,lev) = shape_param_adj(col%gridcell(c), lev)
601599
soilstate_inst%sucsat_col(c,lev) = psis_sat_adj(col%gridcell(c), lev)
602600
soilstate_inst%hksat_col(c,lev) = ks_adj(col%gridcell(c), lev) ! mm/s
603601
end if
604-
#endif
602+
! SHP end
605603
soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev))
606604

607605
soilstate_inst%tksatu_col(c,lev) = soilstate_inst%tkmg_col(c,lev)*0.57_r8**soilstate_inst%watsat_col(c,lev)
@@ -675,15 +673,13 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
675673
soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac)*soilstate_inst%watsat_col(c,lev) + om_watsat_lake * om_frac
676674

677675
tkm = (1._r8-om_frac)*(8.80_r8*sand+2.92_r8*clay)/(sand+clay) + om_tkm * om_frac ! W/(m K)
678-
#ifndef USE_PDAF
679-
soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake
680-
#else
676+
! SHP adapt start
681677
if (parameters_in_file .or. parameters_in_file_adj) then
682678
soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake
683679
else
684-
soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake
680+
soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake
685681
end if
686-
#endif
682+
! SHP adapt end
687683

688684
soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac)*soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac
689685

@@ -747,10 +743,10 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
747743

748744
deallocate(sand3d, clay3d, organic3d)
749745
deallocate(zisoifl, zsoifl, dzsoifl)
750-
#ifdef USE_PDAF
746+
! SHP start
751747
deallocate(thetas, shape_param, psis_sat, ks)
752748
deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj)
753-
#endif
749+
! SHP end
754750

755751
end subroutine SoilStateInitTimeConst
756752

0 commit comments

Comments
 (0)