Skip to content

Commit

Permalink
Syntax cleaning, mostly related to using defined constants, and defin…
Browse files Browse the repository at this point in the history 
…ed indexes related to radiation scattering and tree_lai.
  • Loading branch information
@rgknox
rgknox committed Apr 17, 2018
1 parent edbc92c commit eb2d66e
Show file tree
Hide file tree
Showing 3 changed files with 71 additions and 44 deletions.
16 changes: 6 additions & 10 deletions biogeochem/FatesAllometryMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,9 @@ module FatesAllometryMod
use EDPFTvarcon , only : EDPftvarcon_inst
use FatesConstantsMod, only : r8 => fates_r8
use FatesConstantsMod, only : i4 => fates_int
use FatesConstantsMod, only : g_per_kg
use FatesConstantsMod, only : cm2_per_m2
use FatesConstantsMod, only : kg_per_Megag
use shr_log_mod , only : errMsg => shr_log_errMsg
use FatesGlobals , only : fates_log
use FatesGlobals , only : endrun => fates_endrun
Expand Down Expand Up @@ -520,7 +523,7 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n )
write(fates_log(),*) 'problem in treelai',bl,pft
endif

slat = 1000.0_r8 * EDPftvarcon_inst%slatop(pft) ! m2/g to m2/kg
slat = g_per_kg * EDPftvarcon_inst%slatop(pft) ! m2/g to m2/kg
leafc_per_unitarea = bl/(c_area/n) !KgC/m2
if(leafc_per_unitarea > 0.0_r8)then
tree_lai = leafc_per_unitarea * slat !kg/m2 * m2/kg = unitless LAI
Expand Down Expand Up @@ -560,7 +563,7 @@ real(r8) function tree_sai( dbh, pft, canopy_trim, c_area, n )
real(r8) :: sai_scaler
real(r8) :: b_leaf

sai_scaler = 1000. * EDPftvarcon_inst%allom_sai_scaler(pft) ! m2/g to m2/kg
sai_scaler = g_per_kg * EDPftvarcon_inst%allom_sai_scaler(pft) ! m2/g to m2/kg

call bleaf(dbh,pft,canopy_trim,b_leaf)

Expand Down Expand Up @@ -881,9 +884,7 @@ end subroutine bbgw_const

subroutine bsap_deprecated(d,h,dhdd,bleaf,dbleafdd,ipft,bsap,dbsapdd)

use FatesConstantsMod, only : g_per_kg
use FatesConstantsMod, only : cm2_per_m2
use FatesConstantsMod, only : kg_per_Megag


! -------------------------------------------------------------------------
! -------------------------------------------------------------------------
Expand Down Expand Up @@ -930,10 +931,6 @@ end subroutine bsap_deprecated

subroutine bsap_dlinear(d,h,dhdd,bleaf,dbleafdd,ipft,bsap,dbsapdd)

use FatesConstantsMod, only : g_per_kg
use FatesConstantsMod, only : cm2_per_m2
use FatesConstantsMod, only : kg_per_Megag

! -------------------------------------------------------------------------
! Calculate sapwood biomass based on leaf area to sapwood area
! proportionality. In this function, the leaftosapwood area is a function
Expand Down Expand Up @@ -1762,7 +1759,6 @@ subroutine StructureResetOfDH( bdead, ipft, canopy_trim, d, h )
! T
! ============================================================================


use FatesConstantsMod , only : calloc_abs_error
! Arguments

Expand Down
56 changes: 30 additions & 26 deletions biogeophys/EDSurfaceAlbedoMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,12 @@ module EDSurfaceRadiationMod
use EDTypesMod , only : maxSWb
use EDTypesMod , only : nclmax
use EDTypesMod , only : nlevleaf
use EDTypesMod , only : n_rad_stream_types
use EDTypesMod , only : idiffuse
use EDTypesMod , only : idirect
use EDTypesMod , only : ivis
use EDTypesMod , only : inir
use EDTypesMod , only : ipar
use EDCanopyStructureMod, only: calc_areaindex
use FatesGlobals , only : fates_log

Expand All @@ -39,10 +45,6 @@ module EDSurfaceRadiationMod
real(r8), public :: albice(maxSWb) = & ! albedo land ice by waveband (1=vis, 2=nir)
(/ 0.80_r8, 0.55_r8 /)

! INTERF-TODO: THIS NEEDS SOME CONSISTENCY AND SHOULD BE SET IN THE INTERFACE
! WITH OTHER DIMENSIONS
integer, parameter :: ipar = 1 ! The band index for PAR

contains

subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
Expand Down Expand Up @@ -214,9 +216,10 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
end do !ft
end do !L

do radtype = 1,2 !do this once for one unit of diffuse, and once for one unit of direct radiation
!do this once for one unit of diffuse, and once for one unit of direct radiation
do radtype = 1, n_rad_stream_types
do ib = 1,hlm_numSWb
if (radtype == 1) then
if (radtype == idirect) then
! Set the hypothetical driving radiation. We do this once for a single unit of direct and
! once for a single unit of diffuse radiation.
forc_dir(ifp,ib) = 1.00_r8
Expand Down Expand Up @@ -706,9 +709,9 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
! Absorbed radiation, shaded and sunlit portions of leaf layers
!here we get one unit of diffuse radiation... how much of
!it is absorbed?
if (ib == 1) then ! only set the absorbed PAR for the visible light band.
if (ib == ivis) then ! only set the absorbed PAR for the visible light band.
do iv = 1, currentPatch%nrad(L,ft)
if (radtype==1) then
if (radtype==idirect) then
if ( DEBUG ) then
write(fates_log(),*) 'EDsurfAlb 730 ',Abs_dif_z(ft,iv),currentPatch%f_sun(L,ft,iv)
write(fates_log(),*) 'EDsurfAlb 731 ', currentPatch%fabd_sha_z(L,ft,iv), &
Expand Down Expand Up @@ -742,8 +745,9 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
end if
! Solar radiation absorbed by vegetation and sunlit/shaded leaves
do iv = 1,currentPatch%nrad(L,ft)
if (radtype == 1)then
currentPatch%fabd(ib) = currentPatch%fabd(ib) + Abs_dir_z(ft,iv)+Abs_dif_z(ft,iv)
if (radtype == idirect)then
currentPatch%fabd(ib) = currentPatch%fabd(ib) + &
Abs_dir_z(ft,iv)+Abs_dif_z(ft,iv)
! bc_out(s)%fabd_parb(ifp,ib) = currentPatch%fabd(ib)
else
currentPatch%fabi(ib) = currentPatch%fabi(ib) + Abs_dif_z(ft,iv)
Expand All @@ -752,7 +756,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
end do
! Albefor
if (L==1)then !top canopy layer.
if (radtype == 1)then
if (radtype == idirect)then
bc_out(s)%albd_parb(ifp,ib) = bc_out(s)%albd_parb(ifp,ib) + &
Dif_up(L,ft,1) * ftweight(L,ft,1)
else
Expand All @@ -762,7 +766,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
end if

! pass normalized PAR profiles for use in diagnostic averaging for history fields
if (ib == 1) then ! only diagnose PAR profiles for the visible band
if (ib == ivis) then ! only diagnose PAR profiles for the visible band
do iv = 1, currentPatch%nrad(L,ft)
currentPatch%nrmlzd_parprof_pft_dir_z(radtype,L,ft,iv) = &
forc_dir(ifp,ib) * tr_dir_z(L,ft,iv)
Expand All @@ -781,7 +785,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
end if ! ib = visible
end if ! present
end do !ft
if (radtype == 1)then
if (radtype == idirect)then
bc_out(s)%fabd_parb(ifp,ib) = currentPatch%fabd(ib)
else
bc_out(s)%fabi_parb(ifp,ib) = currentPatch%fabi(ib)
Expand All @@ -798,7 +802,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
tr_soild = tr_soild + forc_dir(ifp,ib) * weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1)))
endif

if (radtype == 1)then
if (radtype == idirect)then
currentPatch%tr_soil_dir(ib) = tr_soild
currentPatch%tr_soil_dir_dif(ib) = tr_soili
currentPatch%sabs_dir(ib) = abs_rad(ib)
Expand All @@ -818,7 +822,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
!==============================================================================!
! Total radiation balance: absorbed = incoming - outgoing

if (radtype == 1)then
if (radtype == idirect)then
error = abs(currentPatch%sabs_dir(ib) - (currentPatch%tr_soil_dir(ib) * &
(1.0_r8-bc_in(s)%albgr_dir_rb(ib)) + &
currentPatch%tr_soil_dir_dif(ib) * (1.0_r8-bc_in(s)%albgr_dif_rb(ib))))
Expand All @@ -844,7 +848,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
endif
endif

if (radtype == 1)then
if (radtype == idirect)then
error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - &
(bc_out(s)%fabd_parb(ifp,ib) + bc_out(s)%albd_parb(ifp,ib) + currentPatch%sabs_dir(ib))
else
Expand All @@ -864,7 +868,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
enddo
enddo

if (radtype == 1)then
if (radtype == idirect)then
!here we are adding a within-ED radiation scheme tolerance, and then adding the diffrence onto the albedo
!it is important that the lower boundary for this is ~1000 times smaller than the tolerance in surface albedo.
if (abs(error) > 1.e-9_r8 .and. abs(error) < 0.15_r8)then
Expand Down Expand Up @@ -913,7 +917,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + error
end if

if (radtype == 1)then
if (radtype == idirect)then
error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - &
(bc_out(s)%fabd_parb(ifp,ib) + bc_out(s)%albd_parb(ifp,ib) + currentPatch%sabs_dir(ib))
else
Expand Down Expand Up @@ -1087,27 +1091,27 @@ subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out)
do FT = 1,numpft
do iv = 1, cpatch%nrad(CL,ft)
cpatch%parprof_pft_dir_z(CL,FT,iv) = (bc_in(s)%solad_parb(ifp,ipar) * &
cpatch%nrmlzd_parprof_pft_dir_z(1,CL,FT,iv)) + &
cpatch%nrmlzd_parprof_pft_dir_z(idirect,CL,FT,iv)) + &
(bc_in(s)%solai_parb(ifp,ipar) * &
cpatch%nrmlzd_parprof_pft_dir_z(2,CL,FT,iv))
cpatch%nrmlzd_parprof_pft_dir_z(idiffuse,CL,FT,iv))
cpatch%parprof_pft_dif_z(CL,FT,iv) = (bc_in(s)%solad_parb(ifp,ipar) * &
cpatch%nrmlzd_parprof_pft_dif_z(1,CL,FT,iv)) + &
cpatch%nrmlzd_parprof_pft_dif_z(idirect,CL,FT,iv)) + &
(bc_in(s)%solai_parb(ifp,ipar) * &
cpatch%nrmlzd_parprof_pft_dif_z(2,CL,FT,iv))
cpatch%nrmlzd_parprof_pft_dif_z(idiffuse,CL,FT,iv))
end do ! iv
end do ! FT
end do ! CL

do CL = 1, cpatch%NCL_p
do iv = 1, maxval(cpatch%nrad(CL,:))
cpatch%parprof_dir_z(CL,iv) = (bc_in(s)%solad_parb(ifp,ipar) * &
cpatch%nrmlzd_parprof_dir_z(1,CL,iv)) + &
cpatch%nrmlzd_parprof_dir_z(idirect,CL,iv)) + &
(bc_in(s)%solai_parb(ifp,ipar) * &
cpatch%nrmlzd_parprof_dir_z(2,CL,iv))
cpatch%nrmlzd_parprof_dir_z(idiffuse,CL,iv))
cpatch%parprof_dif_z(CL,iv) = (bc_in(s)%solad_parb(ifp,ipar) * &
cpatch%nrmlzd_parprof_dif_z(1,CL,iv)) + &
cpatch%nrmlzd_parprof_dif_z(idirect,CL,iv)) + &
(bc_in(s)%solai_parb(ifp,ipar) * &
cpatch%nrmlzd_parprof_dif_z(2,CL,iv))
cpatch%nrmlzd_parprof_dif_z(idiffuse,CL,iv))
end do ! iv
end do ! CL

Expand Down
43 changes: 35 additions & 8 deletions main/EDTypesMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,11 +22,25 @@ module EDTypesMod
! the parameter file may determine that fewer
! are used, but this helps allocate scratch
! space and output arrays.


! -------------------------------------------------------------------------------------
! Radiation parameters
! These should be part of the radiation module, but since we only have one option
! this is ok for now. (RGK 04-2018)
! -------------------------------------------------------------------------------------


integer, parameter :: n_rad_stream_types = 2 ! The number of radiation streams used (direct/diffuse)

integer, parameter :: idirect = 1 ! This is the array index for direct radiation
integer, parameter :: idiffuse = 2 ! This is the array index for diffuse radiation


! TODO: we use this cp_maxSWb only because we have a static array q(size=2) of
! land-ice abledo for vis and nir. This should be a parameter, which would
! get us on track to start using multi-spectral or hyper-spectral (RGK 02-2017)

integer, parameter :: maxSWb = 2 ! maximum number of broad-bands in the
! shortwave spectrum cp_numSWb <= cp_maxSWb
! this is just for scratch-array purposes
Expand All @@ -44,6 +58,10 @@ module EDTypesMod
! files. This will be compared with
! the HLM's expectation in FatesInterfaceMod

integer, parameter :: ipar = ivis ! The photosynthetically active band
! can be approximated to be equal to the visible band


! Switches that turn on/off ED dynamics process (names are self explanatory)
! IMPORTANT NOTE!!! THESE SWITCHES ARE EXPERIMENTAL.
! THEY SHOULD CORRECTLY TURN OFF OR ON THE PROCESS, BUT.. THERE ARE VARIOUS
Expand Down Expand Up @@ -342,14 +360,23 @@ module EDTypesMod
real(r8) :: f_sun(nclmax,maxpft,nlevleaf) ! fraction of leaves in the sun in each canopy layer, pft,

! radiation profiles for comparison against observations
real(r8) :: nrmlzd_parprof_pft_dir_z(2,nclmax,maxpft,nlevleaf) ! normalized direct photosynthetically active radiation profiles by
! incident type (direct/diffuse at top of canopy),leaf,pft,leaf (unitless)
real(r8) :: nrmlzd_parprof_pft_dif_z(2,nclmax,maxpft,nlevleaf) ! normalized diffuse photosynthetically active radiation profiles by
! incident type (direct/diffuse at top of canopy),leaf,pft,leaf (unitless)
real(r8) :: nrmlzd_parprof_dir_z(2,nclmax,nlevleaf) ! normalized direct photosynthetically active radiation profiles by
! incident type (direct/diffuse at top of canopy),leaf,leaf (unitless)
real(r8) :: nrmlzd_parprof_dif_z(2,nclmax,nlevleaf) ! normalized diffuse photosynthetically active radiation profiles by
! incident type (direct/diffuse at top of canopy),leaf,leaf (unitless)

! normalized direct photosynthetically active radiation profiles by
! incident type (direct/diffuse at top of canopy),leaf,pft,leaf (unitless)
real(r8) :: nrmlzd_parprof_pft_dir_z(n_rad_stream_types,nclmax,maxpft,nlevleaf)

! normalized diffuse photosynthetically active radiation profiles by
! incident type (direct/diffuse at top of canopy),leaf,pft,leaf (unitless)
real(r8) :: nrmlzd_parprof_pft_dif_z(n_rad_stream_types,nclmax,maxpft,nlevleaf)

! normalized direct photosynthetically active radiation profiles by
! incident type (direct/diffuse at top of canopy),leaf,leaf (unitless)
real(r8) :: nrmlzd_parprof_dir_z(n_rad_stream_types,nclmax,nlevleaf)

! normalized diffuse photosynthetically active radiation profiles by
! incident type (direct/diffuse at top of canopy),leaf,leaf (unitless)
real(r8) :: nrmlzd_parprof_dif_z(n_rad_stream_types,nclmax,nlevleaf)

real(r8) :: parprof_pft_dir_z(nclmax,maxpft,nlevleaf) ! direct-beam PAR profile through canopy, by canopy,PFT,leaf level (w/m2)
real(r8) :: parprof_pft_dif_z(nclmax,maxpft,nlevleaf) ! diffuse PAR profile through canopy, by canopy,PFT,leaf level (w/m2)
real(r8) :: parprof_dir_z(nclmax,nlevleaf) ! direct-beam PAR profile through canopy, by canopy,leaf level (w/m2)
Expand Down

0 comments on commit eb2d66e

Please sign in to comment.