@@ -2511,7 +2511,19 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso ,
25112511 real (kind=kind_phys), dimension(-nsnow+1: 0) :: tksno !snow thermal conductivity (j/m3/k)
25122512 real (kind=kind_phys), dimension( 1:nsoil) :: sice !soil ice content
25132513 real (kind=kind_phys), parameter :: sbeta = -2.0
2514+ real (kind=kind_phys), dimension(4,20) :: soil_carbon ! soil carbon content [kg/m3]
2515+ real (kind=kind_phys), parameter :: soil_carbon_df = 0.25 ! soil carbon therm cond (Lawrence and Slater)
2516+ real (kind=kind_phys), parameter :: soil_carbon_hcpct = 2.5e6 ! soil carbon heat capacity (Lawrence and Slater)
25142517! --------------------------------------------------------------------------------------------------
2518+ ! soil carbon [kg/m3] by vegetation type estimated from global PNNL soil carbon dataset
2519+ ! and VIIRS surface type
2520+
2521+ soil_carbon(1,:) = (/90,65,90,65,90,40,50,50,40,50,90,60,60,60,0,20,0,90,90,60/)
2522+ soil_carbon(2,:) = (/40,30,40,30,40,25,30,30,25,30,40,30,30,30,0,15,0,60,60,40/)
2523+ soil_carbon(3,:) = (/20,15,20,15,20,15,20,15,15,15,25,20,20,20,0,10,0,40,40,30/)
2524+ soil_carbon(4,:) = (/15,10,15,10,15,10,15,10,10,10,20,10,10,10,0,10,0,40,30,20/)
2525+
2526+ soil_carbon = soil_carbon / 130.0 ! convert to soil carbon relative to peat
25152527
25162528! compute snow thermal conductivity and heat capacity
25172529
@@ -2530,6 +2542,11 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso ,
25302542 hcpct(iz) = sh2o(iz)*cwat + (1.0-parameters%smcmax(iz))*parameters%csoil &
25312543 + (parameters%smcmax(iz)-smc(iz))*cpair + sice(iz)*cice
25322544 call tdfcnd (parameters,iz,df(iz), smc(iz), sh2o(iz))
2545+
2546+ ! adjust for soil carbon organic content
2547+
2548+ ! hcpct(iz) = (1.0 - soil_carbon(iz,vegtyp)) * hcpct(iz) + soil_carbon(iz,vegtyp) * soil_carbon_hcpct
2549+ df(iz) = (1.0 - soil_carbon(iz,vegtyp)) * df(iz) + soil_carbon(iz,vegtyp) * soil_carbon_df
25332550 end do
25342551
25352552 if ( parameters%urban_flag ) then
@@ -3003,7 +3020,11 @@ subroutine albedo (parameters,vegtyp ,ist ,ice ,nsoil , & !in
30033020 if (ib.eq.1) fsun = 0.
30043021 end do
30053022
3006- if(cosz <= 0) goto 100
3023+ ! snow age
3024+
3025+ call snow_age (parameters,dt,tg,sneqvo,sneqv,tauss,fage)
3026+
3027+ if(cosz > 0) then
30073028
30083029! weight reflectance/transmittance by lai and sai
30093030
@@ -3015,10 +3036,6 @@ subroutine albedo (parameters,vegtyp ,ist ,ice ,nsoil , & !in
30153036 tau(ib) = max(parameters%taul(ib)*wl+parameters%taus(ib)*ws, mpe)
30163037 end do
30173038
3018- ! snow age
3019-
3020- call snow_age (parameters,dt,tg,sneqvo,sneqv,tauss,fage)
3021-
30223039! snow albedos: only if cosz > 0 and fsno > 0
30233040
30243041 if(opt_alb == 1) &
@@ -3067,8 +3084,7 @@ subroutine albedo (parameters,vegtyp ,ist ,ice ,nsoil , & !in
30673084 wl = ext
30683085 end if
30693086 fsun = wl
3070-
3071- 100 continue
3087+ end if
30723088
30733089 end subroutine albedo
30743090
@@ -9126,9 +9142,7 @@ subroutine groundwater(parameters,nsnow ,nsoil ,dt ,sice ,zsoil , & !in
91269142
91279143! recharge rate qin to groundwater
91289144
9129- ! ka = hk(iwt)
9130- ! harmonic average, c.he changed based on gy niu's update
9131- ka = 2.0*(hk(iwt)*parameters%dksat(iwt)*1.0e3) / (hk(iwt)+parameters%dksat(iwt)*1.0e3)
9145+ ka = 0.5*(hk(iwt)+parameters%dksat(iwt)*1.0e3)
91329146
91339147 wh_zwt = - zwt * 1.e3 !(mm)
91349148 wh = smpfz - znode(iwt)*1.e3 !(mm)
0 commit comments