diff --git a/BRAMS/Template/RAMSIN b/BRAMS/Template/RAMSIN index ebf07b19c..6090e294f 100644 --- a/BRAMS/Template/RAMSIN +++ b/BRAMS/Template/RAMSIN @@ -2525,16 +2525,23 @@ $ED2_INFO ! often will never allow fires. ! ! 2. Fire will be triggered with enough biomass and the total soil ! ! water at the top 75 cm falls below a threshold. ! + ! 3. Fire will be triggered with enough biomass and accumulated ! + ! 30-day water deficit exceeds the threshold given by SM_FIRE. ! + ! This is soil independent. ! ! FIRE_PARAMETER -- If fire happens, this will control the intensity of the disturbance ! ! given the amount of fuel (currently the total above-ground ! ! biomass). ! - ! SM_FIRE -- This is used only when INCLUDE_FIRE = 2. The sign here matters. ! - ! >= 0. - Minimum relative soil moisture above dry air of the top 1m ! - ! that will prevent fires to happen. ! - ! < 0. - Minimum mean soil moisture potential in MPa of the top 1m ! - ! that will prevent fires to happen. The dry air soil ! - ! potential is defined as -3.1 MPa, so make sure SM_FIRE is ! - ! greater than this value. ! + ! SM_FIRE -- This is used only when INCLUDE_FIRE = 2 or 3, and it has different ! + ! meanings. The sign here matters. ! + ! When INCLUDE_FIRE = 2: ! + ! >= 0. - Minimum relative soil moisture above dry air of the top ! + ! 1m that will prevent fires to happen. ! + ! < 0. - Minimum mean soil moisture potential in MPa of the top ! + ! 1m that will prevent fires to happen. The dry air soil ! + ! potential is defined as -3.1 MPa, so make sure SM_FIRE ! + ! is greater than this value. ! + ! When INCLUDE_FIRE = 3, only positive values are allowed. This is ! + ! the minimum water deficit, in kg/m2/30 days, to trigger fires. ! !---------------------------------------------------------------------------------------! INCLUDE_FIRE = 2, FIRE_PARAMETER = 0.5, diff --git a/BRAMS/run/RAMSIN b/BRAMS/run/RAMSIN index 0cbc294b0..8ff1749af 100644 --- a/BRAMS/run/RAMSIN +++ b/BRAMS/run/RAMSIN @@ -2547,16 +2547,23 @@ $ED2_INFO ! often will never allow fires. ! ! 2. Fire will be triggered with enough biomass and the total soil ! ! water at the top 75 cm falls below a threshold. ! + ! 3. Fire will be triggered with enough biomass and accumulated ! + ! 30-day water deficit exceeds the threshold given by SM_FIRE. ! + ! This is soil independent. ! ! FIRE_PARAMETER -- If fire happens, this will control the intensity of the disturbance ! ! given the amount of fuel (currently the total above-ground ! ! biomass). ! - ! SM_FIRE -- This is used only when INCLUDE_FIRE = 2. The sign here matters. ! - ! >= 0. - Minimum relative soil moisture above dry air of the top 1m ! - ! that will prevent fires to happen. ! - ! < 0. - Minimum mean soil moisture potential in MPa of the top 1m ! - ! that will prevent fires to happen. The dry air soil ! - ! potential is defined as -3.1 MPa, so make sure SM_FIRE is ! - ! greater than this value. ! + ! SM_FIRE -- This is used only when INCLUDE_FIRE = 2 or 3, and it has different ! + ! meanings. The sign here matters. ! + ! When INCLUDE_FIRE = 2: ! + ! >= 0. - Minimum relative soil moisture above dry air of the top ! + ! 1m that will prevent fires to happen. ! + ! < 0. - Minimum mean soil moisture potential in MPa of the top ! + ! 1m that will prevent fires to happen. The dry air soil ! + ! potential is defined as -3.1 MPa, so make sure SM_FIRE ! + ! is greater than this value. ! + ! When INCLUDE_FIRE = 3, only positive values are allowed. This is ! + ! the minimum water deficit, in kg/m2/30 days, to trigger fires. ! !---------------------------------------------------------------------------------------! INCLUDE_FIRE = 2, FIRE_PARAMETER = 0.5, diff --git a/ED/Template/Template/ED2IN b/ED/Template/Template/ED2IN index 42940dda8..fa58f6b8b 100644 --- a/ED/Template/Template/ED2IN +++ b/ED/Template/Template/ED2IN @@ -1045,16 +1045,23 @@ $ED_NL ! often will never allow fires. ! ! 2. Fire will be triggered with enough biomass and the total soil ! ! water at the top 75 cm falls below a threshold. ! + ! 3. Fire will be triggered with enough biomass and accumulated ! + ! 30-day water deficit exceeds the threshold given by SM_FIRE. ! + ! This is soil independent. ! ! FIRE_PARAMETER -- If fire happens, this will control the intensity of the disturbance ! ! given the amount of fuel (currently the total above-ground ! ! biomass). ! - ! SM_FIRE -- This is used only when INCLUDE_FIRE = 2. The sign here matters. ! - ! >= 0. - Minimum relative soil moisture above dry air of the top 1m ! - ! that will prevent fires to happen. ! - ! < 0. - Minimum mean soil moisture potential in MPa of the top 1m ! - ! that will prevent fires to happen. The dry air soil ! - ! potential is defined as -3.1 MPa, so make sure SM_FIRE is ! - ! greater than this value. ! + ! SM_FIRE -- This is used only when INCLUDE_FIRE = 2 or 3, and it has different ! + ! meanings. The sign here matters. ! + ! When INCLUDE_FIRE = 2: ! + ! >= 0. - Minimum relative soil moisture above dry air of the top ! + ! 1m that will prevent fires to happen. ! + ! < 0. - Minimum mean soil moisture potential in MPa of the top ! + ! 1m that will prevent fires to happen. The dry air soil ! + ! potential is defined as -3.1 MPa, so make sure SM_FIRE ! + ! is greater than this value. ! + ! When INCLUDE_FIRE = 3, only positive values are allowed. This is ! + ! the minimum water deficit, in kg/m2/day, to trigger fires. ! !---------------------------------------------------------------------------------------! NL%INCLUDE_FIRE = myfire NL%FIRE_PARAMETER = myfuel diff --git a/ED/Template/Template/plot_monthly.r b/ED/Template/Template/plot_monthly.r index 06e3dcf17..c5f2e064e 100644 --- a/ED/Template/Template/plot_monthly.r +++ b/ED/Template/Template/plot_monthly.r @@ -335,6 +335,7 @@ for (place in myplaces){ mmsqu.wflxgc = NULL mmsqu.evap = NULL mmsqu.transp = NULL + water.deficit = NULL #----- Cohort level lists. -------------------------------------------------------------# lightco = list() @@ -828,6 +829,13 @@ for (place in myplaces){ #--------------------------------------------------------------------------------# + #--------------------------------------------------------------------------------# + # Get the water deficit. # + #--------------------------------------------------------------------------------# + water.deficit = c(water.deficit,sum(mymont$AVG.MONTHLY.WATERDEF * areapa)) + #--------------------------------------------------------------------------------# + + #--------------------------------------------------------------------------------# # If this is a biomass initialisation, or a run with anthropogenic # # disturbance, we must jitter the age so we can distinguish the patches. # @@ -2243,20 +2251,32 @@ for (place in myplaces){ obsnow = paste("obs.",iata,sep="") }#end if + #------------------------------------------------------------------------------------# + # Last check to see if we should plot it or not. # + #------------------------------------------------------------------------------------# plotit = plotit && obsnow %in% ls() + if (plotit){ + thisobs = get(obsnow) + obswhen = thisobs$tomonth + sel = thismonth >= min(obswhen) & thismonth <= max(obswhen) + plotit = any(sel) + }#end if + #------------------------------------------------------------------------------------# + #------------------------------------------------------------------------------------# + # Enter here only if there is any overlap of time between observations and # + # model. # + #------------------------------------------------------------------------------------# if (plotit){ #---------------------------------------------------------------------------------# # Copy the observations to a scratch variable. # #---------------------------------------------------------------------------------# - thisobs = get(obsnow) mnvar = paste("emean",vname,sep=".") obsmean = thisobs[[mnvar]] - obswhen = thisobs$tomonth #---------------------------------------------------------------------------------# @@ -2272,15 +2292,14 @@ for (place in myplaces){ #----- Define the number of layers. ----------------------------------------------# - sel = thismonth >= min(obswhen) & thismonth <= max(obswhen) - thiswhen = thismonth - thismean = get(vname) + thiswhen = thismonth [sel] + thismean = get(vname)[sel] #---------------------------------------------------------------------------------# #----- Find the plot range. ------------------------------------------------------# - ylimit = range(c(thismean[sel],obsmean),na.rm=TRUE) + ylimit = range(c(thismean,obsmean),na.rm=TRUE) #----- Expand the upper range in so the legend doesn't hide things. --------------# if (ylimit[1] == ylimit[2] & ylimit[1] == 0){ ylimit[1] = -1 @@ -2411,6 +2430,28 @@ for (place in myplaces){ + #---------------------------------------------------------------------------------# + # Some observations do not have enough measurements to make a full year. If # + # this is the case, then we must split the observations into smaller intervals so # + # the polygon works. In case no observation is available, make the vectors NULL # + # so we will not plot observations at all. # + #---------------------------------------------------------------------------------# + if (all(is.na(obsmean+obssdev))){ + obs.x = NULL + obs.ylow = NULL + obs.yhigh = NULL + }else{ + #------ Find the periods with continous data. ---------------------------------# + ok = is.finite(obsmean+obssdev) + obs.x = montmont[ok] + obs.ylow = obsmean [ok] - obssdev[ok] + obs.yhigh = obsmean [ok] + obssdev[ok] + #------------------------------------------------------------------------------# + }#end if + #---------------------------------------------------------------------------------# + + + #---------------------------------------------------------------------------------# # Check whether the time series directory exists. If not, create it. # #---------------------------------------------------------------------------------# @@ -2421,28 +2462,27 @@ for (place in myplaces){ - #----- Define the number of layers. ----------------------------------------------# - thismean = mont12mn[[vname]] - thissdev = mont12sd[[vname]] - #---------------------------------------------------------------------------------# - - - #---------------------------------------------------------------------------------# - # Some variables have no standard deviation in the model. Make them 0 if this # - # is the case. # + # Define the number of layers. Some variables have no standard deviation in # + # the model, so Make them 0 if this is the case. # #---------------------------------------------------------------------------------# - if (length(thissdev) == 0){ + thismean = mont12mn[[vname]] + thissdev = mont12sd[[vname]] + if (length(mont12sd[[vname]]) == 0){ thissdev = 0. * thismean + }else{ + thissdev = mont12sd[[vname]] }#end if + mod.x = montmont + mod.ylow = thismean - thissdev + mod.yhigh = thismean + thissdev #---------------------------------------------------------------------------------# #----- Find the plot range. ------------------------------------------------------# if (plotsd){ - ylimit = range(c(thismean + thissdev ,thismean - thissdev - ,obsmean + obssdev ,obsmean - obssdev ),na.rm=TRUE) + ylimit = range(c(mod.ylow,mod.yhigh,obs.ylow,obs.yhigh),na.rm=TRUE) }else{ ylimit = range(c(thismean,obsmean),na.rm=TRUE) }#end if @@ -2491,11 +2531,17 @@ for (place in myplaces){ abline(v=montplot$levels,h=axTicks(side=2),col="gray52",lty="solid") }#end if if (plotsd){ - err.x = c(montmont,rev(montmont),NA,montmont,rev(montmont)) - err.y = c(thismean + thissdev,rev(thismean) - rev(thissdev),NA - ,obsmean + obssdev ,rev(obsmean ) - rev(obssdev ) ) - polygon(x=err.x,y=err.y,col=errcolours,angle=angle,density=dens - ,lty="solid",lwd=shwd) + if (is.null(obs.x)){ + err.x = c(mod.x,rev(mod.x)) + err.y = c(mod.ylow,rev(mod.yhigh)) + polygon(x=err.x,y=err.y,col=errcolours[1],angle=angle[1],density=dens[1] + ,lty="solid",lwd=shwd[1]) + }else{ + err.x = c(mod.x,rev(mod.x),NA,obs.x,rev(obs.x)) + err.y = c(mod.ylow,rev(mod.yhigh),NA,obs.ylow,rev(obs.yhigh)) + polygon(x=err.x,y=err.y,col=errcolours,angle=angle,density=dens + ,lty="solid",lwd=shwd) + }#end if }#end if points(x=montmont,y=thismean,col=lcolours[1],lwd=llwd[1],type=ltype ,pch=16,cex=1.0) @@ -2578,6 +2624,20 @@ for (place in myplaces){ #---------------------------------------------------------------------------------# + + #---------------------------------------------------------------------------------# + # Some observations do not have enough measurements to make a full year. If # + # this is the case, then we must split the observations into smaller intervals so # + # the polygon works. In case no observation is available, make the vectors NULL # + # so we will not plot observations at all. # + #---------------------------------------------------------------------------------# + obs.ylow = obsmean - obssdev + obs.yhigh = obsmean + obssdev + #---------------------------------------------------------------------------------# + + + + #---------------------------------------------------------------------------------# # Check whether the time series directory exists. If not, create it. # #---------------------------------------------------------------------------------# @@ -2590,33 +2650,28 @@ for (place in myplaces){ - #----- Define the number of layers. ----------------------------------------------# - thismean = dcyc12mn[[vname]] - thissdev = dcyc12sd[[vname]] - #---------------------------------------------------------------------------------# - - #---------------------------------------------------------------------------------# - # Some variables have no standard deviation in the model. Make them 0 if this # - # is the case. # + # Define the number of layers. Some variables have no standard deviation in # + # the model, so Make them 0 if this is the case. We also append the last hour # + # before the first one so 00 UTC appears in the left. # #---------------------------------------------------------------------------------# - if (length(thissdev) == 0){ + thismean = dcyc12mn[[vname]] + thismean = cbind(thismean[,ndcycle],thismean) + if (length(dcyc12sd[[vname]]) == 0){ thissdev = 0. * thismean + }else{ + thissdev = dcyc12sd[[vname]] + thissdev = cbind(thissdev[,ndcycle],thissdev) }#end if - #---------------------------------------------------------------------------------# - - - #----- Append the last hour before the first one. --------------------------------# - thismean = cbind(thismean[,ndcycle],thismean) - thissdev = cbind(thissdev[,ndcycle],thissdev) + mod.ylow = thismean - thissdev + mod.yhigh = thismean + thissdev #---------------------------------------------------------------------------------# #----- Find the plot range. ------------------------------------------------------# if (plotsd){ - ylimit = range(c(thismean + thissdev ,thismean - thissdev - ,obsmean + obssdev ,obsmean - obssdev ),na.rm=TRUE) + ylimit = range(c(mod.ylow,mod.yhigh,obs.ylow,obs.yhigh),na.rm=TRUE) }else{ ylimit = range(c(thismean,obsmean),na.rm=TRUE) }#end if @@ -2663,14 +2718,34 @@ for (place in myplaces){ abline(v=dcycplot$levels,h=axTicks(side=2),col="gray52",lty="solid") }#end if if (plotsd){ - err.x = c(thisday,rev(thisday),NA,thisday,rev(thisday)) - err.y = c(thismean[pmon,] + thissdev[pmon,] - ,rev(thismean[pmon,]) - rev(thissdev[pmon,]) - ,NA - ,obsmean[pmon,] + obssdev[pmon,] - ,rev(obsmean[pmon,]) - rev(obssdev[pmon,])) - polygon(x=err.x,y=err.y,col=errcolours,angle=angle,density=dens - ,lty="solid",lwd=shwd) + mod.x.now = thisday + mod.ylow.now = mod.ylow [pmon,] + mod.yhigh.now = mod.yhigh[pmon,] + #------ Find the periods with continous data. ---------------------------# + ok = is.finite(obs.ylow[pmon,]) & is.finite(obs.yhigh[pmon,]) + if (any(ok)){ + obs.x.now = thisday [ok] + obs.ylow.now = obs.ylow [pmon,ok] + obs.yhigh.now = obs.yhigh[pmon,ok] + }else{ + obs.x.now = NULL + obs.ylow.now = NULL + obs.yhigh.now = NULL + }#end if + #------------------------------------------------------------------------# + + if (is.null(obs.x.now)){ + err.x = c(mod.x.now,rev(mod.x.now)) + err.y = c(mod.ylow.now,rev(mod.yhigh.now)) + polygon(x=err.x,y=err.y,col=errcolours[1],angle=angle[1] + ,density=dens[1],lty="solid",lwd=shwd[1]) + }else{ + err.x = c(mod.x.now,rev(mod.x.now),NA,obs.x.now,rev(obs.x.now)) + err.y = c(mod.ylow.now,rev(mod.yhigh.now),NA + ,obs.ylow.now,rev(obs.yhigh.now)) + polygon(x=err.x,y=err.y,col=errcolours,angle=angle,density=dens + ,lty="solid",lwd=shwd) + }#end if }#end if points(x=thisday,y=thismean[pmon,],col=lcolours[1] ,lwd=llwd[1],type=ltype,pch=16,cex=1.0) diff --git a/ED/Template/epost.sh b/ED/Template/epost.sh index e71ca7a1a..374330498 100755 --- a/ED/Template/epost.sh +++ b/ED/Template/epost.sh @@ -504,7 +504,8 @@ do sed -i s@mymetcycz@${metcycz}@g ${here}/${polyname}/${script} #----- Run R to get the plots. ------------------------------------------------------# - comm="R CMD BATCH ${here}/${polyname}/${script} ${here}/${polyname}/${epostout}" + rbin="R CMD BATCH --no-save --no-restore" + comm="${rbin} ${here}/${polyname}/${script} ${here}/${polyname}/${epostout}" #------------------------------------------------------------------------------------# # plot_eval_ed won't run all at once due to the sheer number of HDF5 files. # diff --git a/ED/Template/spawn_poly.sh b/ED/Template/spawn_poly.sh index 85041e42a..e80d1697f 100755 --- a/ED/Template/spawn_poly.sh +++ b/ED/Template/spawn_poly.sh @@ -564,7 +564,7 @@ do #------------------------------------------------------------------------------------# case ${metdriver} in Bananal) - metdriverdb=${sitemet}'/Bananal_Island/Bananal_HEADER' + metdriverdb=${sitemet}'/Bananal/Bananal_HEADER' metcyc1=2004 metcycf=2006 imetavg=1 diff --git a/ED/run/ED2IN b/ED/run/ED2IN index 9a92c280a..cf1744928 100644 --- a/ED/run/ED2IN +++ b/ED/run/ED2IN @@ -1039,16 +1039,23 @@ $ED_NL ! often will never allow fires. ! ! 2. Fire will be triggered with enough biomass and the total soil ! ! water at the top 75 cm falls below a threshold. ! + ! 3. Fire will be triggered with enough biomass and accumulated ! + ! 30-day water deficit exceeds the threshold given by SM_FIRE. ! + ! This is soil independent. ! ! FIRE_PARAMETER -- If fire happens, this will control the intensity of the disturbance ! ! given the amount of fuel (currently the total above-ground ! ! biomass). ! - ! SM_FIRE -- This is used only when INCLUDE_FIRE = 2. The sign here matters. ! - ! >= 0. - Minimum relative soil moisture above dry air of the top 1m ! - ! that will prevent fires to happen. ! - ! < 0. - Minimum mean soil moisture potential in MPa of the top 1m ! - ! that will prevent fires to happen. The dry air soil ! - ! potential is defined as -3.1 MPa, so make sure SM_FIRE is ! - ! greater than this value. ! + ! SM_FIRE -- This is used only when INCLUDE_FIRE = 2 or 3, and it has different ! + ! meanings. The sign here matters. ! + ! When INCLUDE_FIRE = 2: ! + ! >= 0. - Minimum relative soil moisture above dry air of the top ! + ! 1m that will prevent fires to happen. ! + ! < 0. - Minimum mean soil moisture potential in MPa of the top ! + ! 1m that will prevent fires to happen. The dry air soil ! + ! potential is defined as -3.1 MPa, so make sure SM_FIRE ! + ! is greater than this value. ! + ! When INCLUDE_FIRE = 3, only positive values are allowed. This is ! + ! the minimum water deficit, in kg/m2/30 days, to trigger fires. ! !---------------------------------------------------------------------------------------! NL%INCLUDE_FIRE = 0 NL%FIRE_PARAMETER = 0.5 diff --git a/ED/src/dynamics/fire.f90 b/ED/src/dynamics/fire.f90 index f61b8db04..8506a01df 100644 --- a/ED/src/dynamics/fire.f90 +++ b/ED/src/dynamics/fire.f90 @@ -17,6 +17,7 @@ subroutine fire_frequency(cgrid) , dslz & ! intent(in) , dslzi ! ! intent(in) use disturb_coms , only : include_fire & ! intent(in) + , sm_fire & ! intent(in) , fire_dryness_threshold & ! intent(in) , fire_smoist_depth & ! intent(in) , k_fire_first & ! intent(in) @@ -79,6 +80,12 @@ subroutine fire_frequency(cgrid) csite%avg_monthly_gndwater(ipa) = csite%avg_monthly_gndwater(ipa) * normfac !------------------------------------------------------------------------------! + + !----- Normalise the monthly mean ground water. -------------------------------! + csite%avg_monthly_waterdef(ipa) = max(0.0,csite%avg_monthly_waterdef(ipa)) + !------------------------------------------------------------------------------! + + !----- Initialize patch fuel. -------------------------------------------------! fuel = 0.0 !------------------------------------------------------------------------------! @@ -145,51 +152,20 @@ subroutine fire_frequency(cgrid) fire_intensity = 0.0 end if !---------------------------------------------------------------------------! + case (3) !---------------------------------------------------------------------------! - ! The threshold not only determines whether fires will happen, it will ! - ! also control the fire intensity. ! + ! The threshold is independent on soil moisture. We use climatological ! + ! water deficit instead. ! !---------------------------------------------------------------------------! - fire_wmass_threshold = 0. - avg_slpot = 0. - do k = k_fire_first, nzg - nsoil = cpoly%ntext_soil(k,isi) - fire_wmass_threshold = fire_wmass_threshold & - + soil(nsoil)%soilfr * dslz(k) * wdns - end do - - if (csite%avg_monthly_gndwater(ipa) < fire_wmass_threshold) then - nsoil = cpoly%ntext_soil(nzg,isi) - !----- Find the equivalent soil moisture. -------------------------------! - avg_slmst = max( soil(nsoil)%soilcp & - , min( soil(nsoil)%slmsts & - , csite%avg_monthly_gndwater(ipa) & - / ( wdns * abs(slz(k_fire_first)) ) ) ) - !------------------------------------------------------------------------! - - - !----- Find the equivalent soil potential. ------------------------------! - avg_slpot = soil(nsoil)%slpots & - / ( avg_slmst / soil(nsoil)%slmsts ) ** soil(nsoil)%slbs - !------------------------------------------------------------------------! - - - !----- Find the scale to reduce or amplify fires. -----------------------! - fire_scale = log( avg_slpot / soil(nsoil)%slpotwp) & - / log(soil(nsoil)%slpotfr / soil(nsoil)%slpotwp) - fire_intensity = max(0.0, fire_parameter * (1.0 - fire_scale) ) - !------------------------------------------------------------------------! - + if (csite%avg_monthly_waterdef(ipa) > sm_fire) then + fire_intensity = fire_parameter + mean_fire_intensity = mean_fire_intensity & + + fire_intensity * csite%area(ipa) else - fire_intensity = 0.0 + fire_intensity = 0.0 end if !---------------------------------------------------------------------------! - - !---------------------------------------------------------------------------! - ! Find the contribution of this patch to fires. ! - !---------------------------------------------------------------------------! - mean_fire_intensity = mean_fire_intensity + fire_intensity * csite%area(ipa) - !---------------------------------------------------------------------------! end select !------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/hybrid_driver.f90 b/ED/src/dynamics/hybrid_driver.f90 index ee9445146..6d35094e4 100644 --- a/ED/src/dynamics/hybrid_driver.f90 +++ b/ED/src/dynamics/hybrid_driver.f90 @@ -700,7 +700,7 @@ subroutine copy_fb_patch(sourcep, targetp, cpatch) integer :: k !---------------------------------------------------------------------------------------! - targetp%can_enthalpy = sourcep%can_enthalpy + targetp%can_enthalpy = sourcep%can_enthalpy targetp%can_theta = sourcep%can_theta targetp%can_temp = sourcep%can_temp targetp%can_shv = sourcep%can_shv @@ -757,6 +757,9 @@ subroutine copy_fb_patch(sourcep, targetp, cpatch) targetp%cwd_rh = sourcep%cwd_rh targetp%rh = sourcep%rh + targetp%water_deficit = sourcep%water_deficit + + do k=rk4site%lsl,nzg targetp%soil_water (k) = sourcep%soil_water (k) targetp%soil_energy (k) = sourcep%soil_energy (k) @@ -878,6 +881,7 @@ subroutine copy_fb_patch(sourcep, targetp, cpatch) targetp%avg_drainage = sourcep%avg_drainage targetp%avg_drainage_heat = sourcep%avg_drainage_heat + do k=rk4site%lsl,nzg targetp%avg_sensible_gg(k) = sourcep%avg_sensible_gg(k) targetp%avg_smoist_gg(k) = sourcep%avg_smoist_gg(k) @@ -1059,6 +1063,8 @@ subroutine inc_fwd_patch(rkp, inc, fac, cpatch) rkp%qpwp = rkp%qpwp + fac * inc%qpwp rkp%cpwp = rkp%cpwp + fac * inc%cpwp + rkp%water_deficit = rkp%water_deficit + fac * inc%water_deficit + do ico = 1,cpatch%ncohorts rkp%leaf_water (ico) = rkp%leaf_water (ico) + fac * inc%leaf_water (ico) rkp%leaf_energy(ico) = rkp%leaf_energy(ico) + fac * inc%leaf_energy(ico) @@ -1118,6 +1124,7 @@ subroutine inc_fwd_patch(rkp, inc, fac, cpatch) rkp%avg_sensible_gc = rkp%avg_sensible_gc + fac * inc%avg_sensible_gc rkp%avg_sensible_ac = rkp%avg_sensible_ac + fac * inc%avg_sensible_ac + do k=rk4site%lsl,nzg rkp%avg_sensible_gg(k) = rkp%avg_sensible_gg(k) + fac * inc%avg_sensible_gg(k) rkp%avg_smoist_gg(k) = rkp%avg_smoist_gg(k) + fac * inc%avg_smoist_gg(k) diff --git a/ED/src/dynamics/rk4_derivs.F90 b/ED/src/dynamics/rk4_derivs.F90 index 7390caeed..ff6882bfa 100644 --- a/ED/src/dynamics/rk4_derivs.F90 +++ b/ED/src/dynamics/rk4_derivs.F90 @@ -1797,6 +1797,15 @@ subroutine canopy_derivs_two(mzg,initp,dinitp,csite,ipa,hflxgc,wflxgc,qwflxgc,de initp%wflxac = wflxac + + + !---------------------------------------------------------------------------------------! + ! Water deficit. ! + !---------------------------------------------------------------------------------------! + dinitp%water_deficit = - (wflxac + rk4site%pcpg) + !---------------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------------! ! Integrate diagnostic variables - These are not activated unless fast file-type ! ! outputs are selected. This will speed up the integrator. ! diff --git a/ED/src/dynamics/rk4_driver.F90 b/ED/src/dynamics/rk4_driver.F90 index 394681f51..770927e3a 100644 --- a/ED/src/dynamics/rk4_driver.F90 +++ b/ED/src/dynamics/rk4_driver.F90 @@ -399,9 +399,11 @@ subroutine initp2modelp(hdid,initp,csite,ipa,wbudget_loss2atm,ebudget_netrad real(kind=8) :: gnd_water real(kind=8) :: psiplusz real(kind=8) :: mcheight + real(kind=4) :: step_waterdef real(kind=4) :: can_rvap !----- Local contants ---------------------------------------------------------------! - real , parameter :: tendays_sec=10.*day_sec + real , parameter :: tendays_sec = 10. * day_sec + real , parameter :: thirtydays_sec = 30. * day_sec !----- External function ------------------------------------------------------------! real , external :: sngloff !------------------------------------------------------------------------------------! @@ -546,6 +548,15 @@ subroutine initp2modelp(hdid,initp,csite,ipa,wbudget_loss2atm,ebudget_netrad + !------------------------------------------------------------------------------------! + ! Update the water deficit. This is done as a 30-day running average. ! + !------------------------------------------------------------------------------------! + step_waterdef = sngloff(initp%water_deficit,tiny_offset) + csite%avg_monthly_waterdef(ipa) = csite%avg_monthly_waterdef(ipa) + step_waterdef + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! ! This variable is the monthly mean ground water that will be used to control ! ! fire disturbance. ! diff --git a/ED/src/dynamics/rk4_integ_utils.f90 b/ED/src/dynamics/rk4_integ_utils.f90 index 3f1193de7..7b70dd9bc 100644 --- a/ED/src/dynamics/rk4_integ_utils.f90 +++ b/ED/src/dynamics/rk4_integ_utils.f90 @@ -390,7 +390,8 @@ subroutine inc_rk4_patch(rkp, inc, fac, cpatch) rkp%virtual_water = rkp%virtual_water + fac * inc%virtual_water rkp%virtual_depth = rkp%virtual_depth + fac * inc%virtual_depth - + rkp%water_deficit = rkp%water_deficit + fac * inc%water_deficit + rkp%upwp = rkp%upwp + fac * inc%upwp rkp%wpwp = rkp%wpwp + fac * inc%wpwp rkp%tpwp = rkp%tpwp + fac * inc%tpwp @@ -1251,6 +1252,8 @@ subroutine copy_rk4_patch(sourcep, targetp, cpatch) targetp%cwd_rh = sourcep%cwd_rh targetp%rh = sourcep%rh + targetp%water_deficit = sourcep%water_deficit + do k=rk4site%lsl,nzg targetp%soil_water (k) = sourcep%soil_water (k) targetp%soil_energy (k) = sourcep%soil_energy (k) diff --git a/ED/src/dynamics/rk4_misc.f90 b/ED/src/dynamics/rk4_misc.f90 index acd4e242c..643899b8b 100644 --- a/ED/src/dynamics/rk4_misc.f90 +++ b/ED/src/dynamics/rk4_misc.f90 @@ -429,6 +429,12 @@ subroutine copy_patch_init(sourcesite,ipa,targetp) targetp%wbudget_loss2drainage = 0.d0 targetp%wbudget_loss2runoff = 0.d0 end if + !---------------------------------------------------------------------------------------! + + + !----- Water deficit, always start with zero. ------------------------------------------! + targetp%water_deficit = 0.d0 + !---------------------------------------------------------------------------------------! if (print_detailed) call reset_rk4_fluxes(targetp) diff --git a/ED/src/init/ed_params.f90 b/ED/src/init/ed_params.f90 index 8f9fb4ac7..1018b981f 100644 --- a/ED/src/init/ed_params.f90 +++ b/ED/src/init/ed_params.f90 @@ -1843,9 +1843,9 @@ subroutine init_decomp_params() rh_decay_high = 0.60 rh_low_temp = 18.0 + t00 rh_high_temp = 45.0 + t00 - rh_decay_dry = 12.0 - rh_decay_wet = 72.0 - rh_dry_smoist = 0.48 + rh_decay_dry = 10.0 + rh_decay_wet = 60.0 + rh_dry_smoist = 0.36 rh_wet_smoist = 0.96 !---------------------------------------------------------------------------------------! diff --git a/ED/src/init/ed_type_init.f90 b/ED/src/init/ed_type_init.f90 index 176761a27..8cc7860f5 100644 --- a/ED/src/init/ed_type_init.f90 +++ b/ED/src/init/ed_type_init.f90 @@ -396,6 +396,7 @@ subroutine init_ed_patch_vars(csite,ip1,ip2,lsl) csite%avg_daily_temp (ip1:ip2) = 0.0 csite%avg_monthly_gndwater(ip1:ip2) = 0.0 + csite%avg_monthly_waterdef(ip1:ip2) = 0.0 csite%mean_rh(ip1:ip2) = 0.0 csite%mean_cwd_rh(ip1:ip2) = 0.0 diff --git a/ED/src/io/ed_init_full_history.F90 b/ED/src/io/ed_init_full_history.F90 index 400fbd327..57c91a193 100644 --- a/ED/src/io/ed_init_full_history.F90 +++ b/ED/src/io/ed_init_full_history.F90 @@ -2151,6 +2151,7 @@ end subroutine hdf_getslab_i call hdf_getslab_r(csite%rough,'ROUGH ',dsetrank,iparallel,.true.,foundvar) call hdf_getslab_r(csite%avg_daily_temp,'AVG_DAILY_TEMP ',dsetrank,iparallel,.true.,foundvar) call hdf_getslab_r(csite%avg_monthly_gndwater,'AVG_MONTHLY_GNDWATER ',dsetrank,iparallel,.true.,foundvar) + call hdf_getslab_r(csite%avg_monthly_waterdef,'AVG_MONTHLY_WATERDEF ',dsetrank,iparallel,.true.,foundvar) call hdf_getslab_r(csite%mean_rh,'MEAN_RH ',dsetrank,iparallel,.true.,foundvar) call hdf_getslab_r(csite%mean_cwd_rh,'MEAN_CWD_RH ',dsetrank,iparallel,.true.,foundvar) diff --git a/ED/src/io/ed_opspec.F90 b/ED/src/io/ed_opspec.F90 index 7349c6fa6..5c2416e56 100644 --- a/ED/src/io/ed_opspec.F90 +++ b/ED/src/io/ed_opspec.F90 @@ -1895,14 +1895,26 @@ subroutine ed_opspec_misc end if end if - - if (sm_fire < -3.1 .or. sm_fire > 1.) then - write (reason,fmt='(a,1x,es12.5,a)') & - 'Invalid SM_FIRE, it must be between -3.1 and 1. Yours is set to' & - ,sm_fire,'...' - call opspec_fatal(reason,'opspec_misc') - ifaterr = ifaterr +1 - end if + select case (include_fire) + case(3) + if (sm_fire < 0.) then + write (reason,fmt='(a,1x,a,1x,i4,a,1x,es12.5,a)') & + 'Invalid SM_FIRE, it must be non-negative' & + ,'when INLCUDE_FIRE is ',include_fire,'. Your SM_FIRE is set to' & + ,sm_fire,'...' + call opspec_fatal(reason,'opspec_misc') + ifaterr = ifaterr +1 + end if + case(0:2) + if (sm_fire < -3.1 .or. sm_fire > 1.) then + write (reason,fmt='(a,1x,a,1x,i4,a,1x,es12.5,a)') & + 'Invalid SM_FIRE, it must be between -3.1 and 1.0' & + ,'when INLCUDE_FIRE is ',include_fire,'. Your SM_FIRE is set to' & + ,sm_fire,'...' + call opspec_fatal(reason,'opspec_misc') + ifaterr = ifaterr +1 + end if + end select if (ianth_disturb < 0 .or. ianth_disturb > 1) then write (reason,fmt='(a,1x,i4,a)') & diff --git a/ED/src/memory/disturb_coms.f90 b/ED/src/memory/disturb_coms.f90 index 754c208be..d06600879 100644 --- a/ED/src/memory/disturb_coms.f90 +++ b/ED/src/memory/disturb_coms.f90 @@ -160,9 +160,17 @@ module disturb_coms real :: fire_dryness_threshold !---------------------------------------------------------------------------------------! - ! Fire may occur when the total water (ground + underground) converted to ! - ! equivalent average soil moisture is below this threshold and include_fire is 2. ! - ! Units: relative fraction. ! + ! SM_FIRE -- This is used only when INCLUDE_FIRE = 2 or 3, and it has different ! + ! meanings. The sign here matters. ! + ! When INCLUDE_FIRE = 2: ! + ! >= 0. - Minimum relative soil moisture above dry air of the top ! + ! 1m that will prevent fires to happen. ! + ! < 0. - Minimum mean soil moisture potential in MPa of the top ! + ! 1m that will prevent fires to happen. The dry air soil ! + ! potential is defined as -3.1 MPa, so make sure SM_FIRE ! + ! is greater than this value. ! + ! When INCLUDE_FIRE = 3, only positive values are allowed. This is ! + ! the minimum water deficit, in kg/m2/30 days, to trigger fires. ! !---------------------------------------------------------------------------------------! real :: sm_fire diff --git a/ED/src/memory/ed_state_vars.f90 b/ED/src/memory/ed_state_vars.f90 index f47afbcd8..852335796 100644 --- a/ED/src/memory/ed_state_vars.f90 +++ b/ED/src/memory/ed_state_vars.f90 @@ -727,6 +727,9 @@ module ed_state_vars ! Average monthly ground water [kg/m2], used for fire ignition real , pointer,dimension(:) :: avg_monthly_gndwater + ! Average monthly water deficit [kg/m2], used for fire ignition + real , pointer,dimension(:) :: avg_monthly_waterdef + ! average of rh and cwd_rh [umol/m2/s] over FRQSTATE real , pointer,dimension(:) :: mean_rh real , pointer,dimension(:) :: mean_cwd_rh @@ -3134,6 +3137,7 @@ subroutine allocate_sitetype(csite,npatches) allocate(csite%avg_daily_temp(npatches)) allocate(csite%avg_monthly_gndwater(npatches)) + allocate(csite%avg_monthly_waterdef(npatches)) allocate(csite%mean_rh(npatches)) allocate(csite%mean_cwd_rh(npatches)) allocate(csite%mean_nep(npatches)) @@ -4356,6 +4360,7 @@ subroutine nullify_sitetype(csite) nullify(csite%A_c_max) nullify(csite%avg_daily_temp) nullify(csite%avg_monthly_gndwater) + nullify(csite%avg_monthly_waterdef) nullify(csite%mean_rh) nullify(csite%dmean_rh) nullify(csite%mmean_rh) @@ -5583,6 +5588,7 @@ subroutine deallocate_sitetype(csite) if(associated(csite%avg_daily_temp )) deallocate(csite%avg_daily_temp ) if(associated(csite%avg_monthly_gndwater )) deallocate(csite%avg_monthly_gndwater ) + if(associated(csite%avg_monthly_waterdef )) deallocate(csite%avg_monthly_waterdef ) if(associated(csite%mean_rh )) deallocate(csite%mean_rh ) if(associated(csite%dmean_rh )) deallocate(csite%dmean_rh ) if(associated(csite%qmean_rh )) deallocate(csite%qmean_rh ) @@ -6078,6 +6084,7 @@ subroutine copy_sitetype(isite,osite,ipaa,ipaz,opaa,opaz) osite%wai(opa) = isite%wai(ipa) osite%avg_daily_temp(opa) = isite%avg_daily_temp(ipa) osite%avg_monthly_gndwater(opa) = isite%avg_monthly_gndwater(ipa) + osite%avg_monthly_waterdef(opa) = isite%avg_monthly_waterdef(ipa) osite%mean_rh(opa) = isite%mean_rh(ipa) osite%mean_cwd_rh(opa) = isite%mean_cwd_rh(ipa) osite%mean_nep(opa) = isite%mean_nep(ipa) @@ -6386,6 +6393,7 @@ subroutine copy_sitetype_mask(sitein,siteout,logmask,masksz,newsz) siteout%wai(1:inc) = pack(sitein%wai,logmask) siteout%avg_daily_temp(1:inc) = pack(sitein%avg_daily_temp,logmask) siteout%avg_monthly_gndwater(1:inc) = pack(sitein%avg_monthly_gndwater,logmask) + siteout%avg_monthly_waterdef(1:inc) = pack(sitein%avg_monthly_waterdef,logmask) siteout%mean_rh(1:inc) = pack(sitein%mean_rh,logmask) siteout%mean_cwd_rh(1:inc) = pack(sitein%mean_cwd_rh,logmask) siteout%mean_nep(1:inc) = pack(sitein%mean_nep,logmask) @@ -12713,6 +12721,13 @@ subroutine filltab_sitetype(igr,ipy,isi,init) call metadata_edio(nvar,igr,'No metadata available','[NA]','NA') end if + if (associated(csite%avg_monthly_waterdef)) then + nvar=nvar+1 + call vtable_edio_r(npts,csite%avg_monthly_waterdef,nvar,igr,init,csite%paglob_id, & + var_len,var_len_global,max_ptrs,'AVG_MONTHLY_WATERDEF :31:hist:anal:dail:mont:dcyc') + call metadata_edio(nvar,igr,'Running average of water deficit','[kg/m2/30days]','NA') + end if + if (associated(csite%co2budget_plresp)) then nvar=nvar+1 call vtable_edio_r(npts,csite%co2budget_plresp,nvar,igr,init,csite%paglob_id, & diff --git a/ED/src/memory/rk4_coms.f90 b/ED/src/memory/rk4_coms.f90 index 93cc7df83..61dc12fcf 100644 --- a/ED/src/memory/rk4_coms.f90 +++ b/ED/src/memory/rk4_coms.f90 @@ -272,6 +272,8 @@ module rk4_coms real(kind=8),pointer,dimension(:) :: avg_sensible_gg ! Soil heat flux between layers real(kind=8) :: avg_drainage ! Drainage at the bottom. real(kind=8) :: avg_drainage_heat ! Drainage at the bottom. + !----- Water deficit. ---------------------------------------------------------------! + real(kind=8) :: water_deficit ! Step water deficit !------------------------------------------------------------------------------------! ! Fast time flux variables for each time step. These variables will be defined ! ! only when the user is debugging. ! @@ -1060,6 +1062,8 @@ subroutine zero_rk4_patch(y) y%avg_drainage = 0.d0 y%avg_drainage_heat = 0.d0 + y%water_deficit = 0.d0 + y%flx_carbon_ac = 0.d0 y%flx_carbon_st = 0.d0 y%flx_vapor_lc = 0.d0 @@ -1521,7 +1525,6 @@ subroutine reset_rk4_fluxes(y) type(rk4patchtype), target :: y !------------------------------------------------------------------------------------! - y%flx_carbon_ac = 0.d0 y%flx_vapor_lc = 0.d0 y%flx_vapor_wc = 0.d0 @@ -1587,7 +1590,6 @@ subroutine norm_rk4_fluxes(y,hdid) !----- Find the inverse of the time step. -------------------------------------------! hdidi = 1.d0 / hdid - y%flx_carbon_ac = y%flx_carbon_ac * hdidi y%flx_vapor_lc = y%flx_vapor_lc * hdidi y%flx_vapor_wc = y%flx_vapor_wc * hdidi diff --git a/ED/src/utils/fuse_fiss_utils.f90 b/ED/src/utils/fuse_fiss_utils.f90 index fbdfd26f0..61fb4d1c7 100644 --- a/ED/src/utils/fuse_fiss_utils.f90 +++ b/ED/src/utils/fuse_fiss_utils.f90 @@ -3484,6 +3484,12 @@ subroutine fuse_2_patches(csite,donp,recp,mzg,mzs,prss,lsl,ntext_soil,green_leaf * csite%area (donp) & + csite%avg_monthly_gndwater (recp) & * csite%area (recp) ) + + csite%avg_monthly_waterdef (recp) = newareai & + * ( csite%avg_monthly_waterdef (donp) & + * csite%area (donp) & + + csite%avg_monthly_waterdef (recp) & + * csite%area (recp) ) !------------------------------------------------------------------------------------!