From eb2d66ec8cbbecd62dc7e34e60fd2854d06afb29 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 17 Apr 2018 11:53:35 -0700 Subject: [PATCH] Syntax cleaning, mostly related to using defined constants, and defined indexes related to radiation scattering and tree_lai. --- biogeochem/FatesAllometryMod.F90 | 16 ++++----- biogeophys/EDSurfaceAlbedoMod.F90 | 56 +++++++++++++++++-------------- main/EDTypesMod.F90 | 43 +++++++++++++++++++----- 3 files changed, 71 insertions(+), 44 deletions(-) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 702764bfae..97b6cd8db6 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -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 @@ -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 @@ -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) @@ -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 + ! ------------------------------------------------------------------------- ! ------------------------------------------------------------------------- @@ -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 @@ -1762,7 +1759,6 @@ subroutine StructureResetOfDH( bdead, ipft, canopy_trim, d, h ) ! T ! ============================================================================ - use FatesConstantsMod , only : calloc_abs_error ! Arguments diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index d71b9715cd..95c37de148 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -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 @@ -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 ) @@ -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 @@ -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), & @@ -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) @@ -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 @@ -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) @@ -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) @@ -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) @@ -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)))) @@ -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 @@ -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 @@ -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 @@ -1087,13 +1091,13 @@ 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 @@ -1101,13 +1105,13 @@ subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out) 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 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 8ecea6054a..89dd9794a8 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -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 @@ -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 @@ -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)