From 014868e79635688300eb6351ea687ce7eb2fae83 Mon Sep 17 00:00:00 2001 From: Sean Swenson Date: Fri, 15 Jul 2022 13:18:46 -0600 Subject: [PATCH] minor cleanup --- src/biogeophys/HillslopeHydrologyMod.F90 | 51 +++++++++--------------- src/biogeophys/IrrigationMod.F90 | 9 ++++- src/biogeophys/SoilHydrologyMod.F90 | 10 ++--- src/biogeophys/SurfaceAlbedoType.F90 | 2 +- src/biogeophys/SurfaceRadiationMod.F90 | 12 +++--- src/biogeophys/WaterFluxType.F90 | 24 +++++------ src/main/TopoMod.F90 | 23 +++++++---- src/main/lnd2atmMod.F90 | 2 +- 8 files changed, 67 insertions(+), 66 deletions(-) diff --git a/src/biogeophys/HillslopeHydrologyMod.F90 b/src/biogeophys/HillslopeHydrologyMod.F90 index 270252adfd..0f0a3e70d3 100644 --- a/src/biogeophys/HillslopeHydrologyMod.F90 +++ b/src/biogeophys/HillslopeHydrologyMod.F90 @@ -558,27 +558,14 @@ subroutine InitHillslope(bounds,fsurdat,glc_behavior) lun%stream_channel_length(l) = 0.5_r8 * lun%stream_channel_length(l) endif - ! if missing hillslope information on surface dataset, fill data - ! and recalculate hillslope_area + ! if missing hillslope information on surface dataset, + ! call endrun if (sum(hillslope_area) == 0._r8) then - do c = lun%coli(l), lun%colf(l) - nh = col%hillslope_ndx(c) - col%hill_area(c) = (grc%area(g)/real(lun%ncolumns(l),r8))*1.e6_r8 ! km2 to m2 - col%hill_width(c) = sqrt(col%hill_area(c)) - col%hill_slope(c) = tan((rpi/180.)*col%topo_slope(c)) - col%hill_aspect(c) = (rpi/2.) ! east (arbitrarily chosen) - if (nh > 0) then - col%hill_elev(c) = col%topo_std(c) & - *((c-lun%coli(l))/ncol_per_hillslope(nh)) - col%hill_distance(c) = sqrt(col%hill_area(c)) & - *((c-lun%coli(l))/ncol_per_hillslope(nh)) - pct_hillslope(l,nh) = 100/nhillslope - else - col%hill_elev(c) = col%topo_std(c) - col%hill_distance(c) = sqrt(col%hill_area(c)) - endif - enddo - + if (masterproc) then + write(iulog,*) 'Problem with input data: nhillcolumns is non-zero, but hillslope area is zero' + write(iulog,*) 'Check surface data for gridcell at (lon/lat): ', grc%londeg(g),grc%latdeg(g) + call endrun( 'ERROR:: sum of hillslope areas is zero.'//errmsg(sourcefile, __LINE__) ) + end if endif ! Recalculate column weights using input areas @@ -1032,14 +1019,14 @@ subroutine HillslopeStreamOutflow(bounds, & !----------------------------------------------------------------------- associate( & stream_water_volume => waterstatebulk_inst%stream_water_volume_lun , & ! Input: [real(r8) (:) ] stream water volume (m3) - qstreamflow => waterfluxbulk_inst%qstreamflow_lun & ! Input: [real(r8) (:) ] stream water discharge (m3/s) + volumetric_streamflow => waterfluxbulk_inst%volumetric_streamflow_lun & ! Input: [real(r8) (:) ] stream water discharge (m3/s) ) ! Get time step dtime = get_step_size_real() do l = bounds%begl,bounds%endl - qstreamflow(l) = 0._r8 + volumetric_streamflow(l) = 0._r8 if(lun%itype(l) == istsoil .and. lun%active(l)) then ! Streamflow calculated from Manning equation if(streamflow_method == streamflow_manning) then @@ -1051,7 +1038,7 @@ subroutine HillslopeStreamOutflow(bounds, & /(lun%stream_channel_width(l) + 2*stream_depth) if(hydraulic_radius <= 0._r8) then - qstreamflow(l) = 0._r8 + volumetric_streamflow(l) = 0._r8 else flow_velocity = (hydraulic_radius)**manning_exponent & * sqrt(lun%stream_channel_slope(l)) & @@ -1060,15 +1047,15 @@ subroutine HillslopeStreamOutflow(bounds, & if (stream_depth > lun%stream_channel_depth(l)) then if (overbank_method == 1) then ! try increasing dynamic slope - qstreamflow(l) = cross_sectional_area * flow_velocity & + volumetric_streamflow(l) = cross_sectional_area * flow_velocity & *(stream_depth/lun%stream_channel_depth(l)) else if (overbank_method == 2) then ! try increasing flow area cross section overbank_area = (stream_depth -lun%stream_channel_depth(l)) * 30._r8 * lun%stream_channel_width(l) - qstreamflow(l) = (cross_sectional_area + overbank_area) * flow_velocity + volumetric_streamflow(l) = (cross_sectional_area + overbank_area) * flow_velocity else if (overbank_method == 3) then ! try removing all overbank flow instantly - qstreamflow(l) = cross_sectional_area * flow_velocity & + volumetric_streamflow(l) = cross_sectional_area * flow_velocity & + (stream_depth-lun%stream_channel_depth(l)) & *lun%stream_channel_width(l)*lun%stream_channel_length(l)/dtime else @@ -1078,13 +1065,13 @@ subroutine HillslopeStreamOutflow(bounds, & endif else - qstreamflow(l) = cross_sectional_area * flow_velocity + volumetric_streamflow(l) = cross_sectional_area * flow_velocity endif ! scale streamflow by number of channel reaches - qstreamflow(l) = qstreamflow(l) * lun%stream_channel_number(l) + volumetric_streamflow(l) = volumetric_streamflow(l) * lun%stream_channel_number(l) - qstreamflow(l) = max(0._r8,min(qstreamflow(l),stream_water_volume(l)/dtime)) + volumetric_streamflow(l) = max(0._r8,min(volumetric_streamflow(l),stream_water_volume(l)/dtime)) endif else if (masterproc) then @@ -1135,7 +1122,7 @@ subroutine HillslopeUpdateStreamWater(bounds, waterstatebulk_inst, & !----------------------------------------------------------------------- associate( & stream_water_volume => waterstatebulk_inst%stream_water_volume_lun, & ! Input/Output: [real(r8) (:) ] stream water volume (m3) - qstreamflow => waterfluxbulk_inst%qstreamflow_lun , & ! Input: [real(r8) (:) ] stream water discharge (m3/s) + volumetric_streamflow => waterfluxbulk_inst%volumetric_streamflow_lun , & ! Input: [real(r8) (:) ] stream water discharge (m3/s) qflx_drain => waterfluxbulk_inst%qflx_drain_col, & ! Input: [real(r8) (:) ] column level sub-surface runoff (mm H2O /s) qflx_drain_perched => waterfluxbulk_inst%qflx_drain_perched_col, & ! Input: [real(r8) (:) ] column level sub-surface runoff (mm H2O /s) qflx_surf => waterfluxbulk_inst%qflx_surf_col , & ! Input: [real(r8) (:) ] total surface runoff (mm H2O /s) @@ -1166,11 +1153,11 @@ subroutine HillslopeUpdateStreamWater(bounds, waterstatebulk_inst, & endif enddo stream_water_volume(l) = stream_water_volume(l) & - - qstreamflow(l) * dtime + - volumetric_streamflow(l) * dtime ! account for negative drainage (via searchforwater in soilhydrology) if(stream_water_volume(l) < 0._r8) then - qstreamflow(l) = qstreamflow(l) + stream_water_volume(l)/dtime + volumetric_streamflow(l) = volumetric_streamflow(l) + stream_water_volume(l)/dtime stream_water_volume(l) = 0._r8 endif diff --git a/src/biogeophys/IrrigationMod.F90 b/src/biogeophys/IrrigationMod.F90 index 27cf050dd3..c78572f0e8 100644 --- a/src/biogeophys/IrrigationMod.F90 +++ b/src/biogeophys/IrrigationMod.F90 @@ -78,7 +78,7 @@ module IrrigationMod use abortutils , only : endrun use clm_instur , only : irrig_method use pftconMod , only : pftcon - use clm_varctl , only : iulog + use clm_varctl , only : iulog, use_hillslope use clm_varcon , only : isecspday, denh2o, spval, ispval use clm_varpar , only : nlevsoi, nlevgrnd use clm_time_manager , only : get_step_size @@ -362,6 +362,7 @@ subroutine ReadNamelist(this, NLFilename, use_aquifer_layer) use spmdMod , only : masterproc, mpicom use shr_mpi_mod , only : shr_mpi_bcast use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) + ! ! !ARGUMENTS: class(irrigation_type) , intent(inout) :: this @@ -583,6 +584,12 @@ subroutine CheckNamelistValidity(this, use_aquifer_layer) errMsg(sourcefile, __LINE__)) end if + if (use_aquifer_layer .and. use_hillslope) then + write(iulog,*) ' ERROR: use_hillslope and use_aquifer_layer may not be used simultaneously' + call endrun(msg=' ERROR: use_hillslope and use_aquifer_layer cannot both be set to true' // & + errMsg(sourcefile, __LINE__)) + end if + end associate end subroutine CheckNamelistValidity diff --git a/src/biogeophys/SoilHydrologyMod.F90 b/src/biogeophys/SoilHydrologyMod.F90 index f52da8cd1a..0c3e748055 100644 --- a/src/biogeophys/SoilHydrologyMod.F90 +++ b/src/biogeophys/SoilHydrologyMod.F90 @@ -2178,7 +2178,7 @@ subroutine SubsurfaceLateralFlow(bounds, & hk_l => soilstate_inst%hk_l_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity (mm/s) qflx_latflow_out => waterfluxbulk_inst%qflx_latflow_out_col, & ! Output: [real(r8) (:) ] lateral saturated outflow (mm/s) qflx_latflow_in => waterfluxbulk_inst%qflx_latflow_in_col, & ! Output: [real(r8) (:) ] lateral saturated inflow (mm/s) - qdischarge => waterfluxbulk_inst%qdischarge_col , & ! Output: [real(r8) (:) ] discharge from column (m3/s) + volumetric_discharge => waterfluxbulk_inst%volumetric_discharge_col , & ! Output: [real(r8) (:) ] discharge from column (m3/s) tdepth => wateratm2lndbulk_inst%tdepth_grc , & ! Input: [real(r8) (:) ] depth of water in tributary channels (m) tdepth_bankfull => wateratm2lndbulk_inst%tdepthmax_grc , & ! Input: [real(r8) (:) ] bankfull depth of tributary channels (m) @@ -2224,7 +2224,7 @@ subroutine SubsurfaceLateralFlow(bounds, & qflx_latflow_in(c) = 0._r8 qflx_latflow_out(c) = 0._r8 qflx_net_latflow(c) = 0._r8 - qdischarge(c) = 0._r8 + volumetric_discharge(c) = 0._r8 qflx_latflow_out_vol(c) = 0._r8 end do @@ -2363,10 +2363,10 @@ subroutine SubsurfaceLateralFlow(bounds, & ! include ice impedance in transmissivity qflx_latflow_out_vol(c) = transmis*col%hill_width(c)*head_gradient - ! qdischarge from lowest column is qflx_latflow_out_vol + ! volumetric_discharge from lowest column is qflx_latflow_out_vol ! scaled by total area of column in gridcell divided by column area if (col%cold(c) == ispval) then - qdischarge(c) = qflx_latflow_out_vol(c) & + volumetric_discharge(c) = qflx_latflow_out_vol(c) & *(grc%area(g)*1.e6_r8*col%wtgcell(c)/col%hill_area(c)) endif @@ -2394,7 +2394,7 @@ subroutine SubsurfaceLateralFlow(bounds, & endif ! convert flux to volumetric flow qflx_latflow_out_vol(c) = 1.e-3_r8*qflx_latflow_out(c)*(grc%area(g)*1.e6_r8*col%wtgcell(c)) - qdischarge(c) = qflx_latflow_out_vol(c) + volumetric_discharge(c) = qflx_latflow_out_vol(c) endif enddo diff --git a/src/biogeophys/SurfaceAlbedoType.F90 b/src/biogeophys/SurfaceAlbedoType.F90 index 05cec3ae1a..110df6bfb5 100644 --- a/src/biogeophys/SurfaceAlbedoType.F90 +++ b/src/biogeophys/SurfaceAlbedoType.F90 @@ -203,7 +203,7 @@ subroutine InitHistory(this, bounds) this%coszen_col(begc:endc) = spval call hist_addfld1d (fname='COSZEN', units='none', & - avgflag='A', long_name='cosine of solar zenith angle', & + avgflag='A', long_name='cosine of solar zenith angle (downscaled if downscaling is activated)', & ptr_col=this%coszen_col, default='inactive') this%albgri_col(begc:endc,:) = spval diff --git a/src/biogeophys/SurfaceRadiationMod.F90 b/src/biogeophys/SurfaceRadiationMod.F90 index a34ac3c951..558b2bf2ee 100644 --- a/src/biogeophys/SurfaceRadiationMod.F90 +++ b/src/biogeophys/SurfaceRadiationMod.F90 @@ -899,8 +899,8 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & rnir = albd(p,2)*forc_solad_col(c,2) + albi(p,2)*forc_solai(g,2) fsr(p) = rvis + rnir if (use_SSRE) then - rvisSF = albdSF(p,1)*forc_solad(g,1) + albiSF(p,1)*forc_solai(g,1) - rnirSF = albdSF(p,2)*forc_solad(g,2) + albiSF(p,2)*forc_solai(g,2) + rvisSF = albdSF(p,1)*forc_solad_col(c,1) + albiSF(p,1)*forc_solai(g,1) + rnirSF = albdSF(p,2)*forc_solad_col(c,2) + albiSF(p,2)*forc_solai(g,2) fsrSF(p) = rvisSF + rnirSF ssre_fsr(p) = fsr(p)-fsrSF(p) end if @@ -913,8 +913,8 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsr_vis_i(p) = albi(p,1)*forc_solai(g,1) fsr_nir_i(p) = albi(p,2)*forc_solai(g,2) if (use_SSRE) then - fsrSF_vis_d(p) = albdSF(p,1)*forc_solad(g,1) - fsrSF_nir_d(p) = albdSF(p,2)*forc_solad(g,2) + fsrSF_vis_d(p) = albdSF(p,1)*forc_solad_col(c,1) + fsrSF_nir_d(p) = albdSF(p,2)*forc_solad_col(c,2) fsrSF_vis_i(p) = albiSF(p,1)*forc_solai(g,1) fsrSF_nir_i(p) = albiSF(p,2)*forc_solai(g,2) @@ -940,8 +940,8 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & end if if (use_SSRE) then if ( is_near_local_noon( grc%londeg(g), deltasec=nint(dtime)/2 ) )then - fsrSF_vis_d_ln(p) = albdSF(p,1)*forc_solad(g,1) - fsrSF_nir_d_ln(p) = albdSF(p,2)*forc_solad(g,2) + fsrSF_vis_d_ln(p) = albdSF(p,1)*forc_solad_col(c,1) + fsrSF_nir_d_ln(p) = albdSF(p,2)*forc_solad_col(c,2) else fsrSF_vis_d_ln(p) = spval fsrSF_nir_d_ln(p) = spval diff --git a/src/biogeophys/WaterFluxType.F90 b/src/biogeophys/WaterFluxType.F90 index 2ab92fead2..5ac717ca59 100644 --- a/src/biogeophys/WaterFluxType.F90 +++ b/src/biogeophys/WaterFluxType.F90 @@ -74,8 +74,8 @@ module WaterFluxType real(r8), pointer :: qflx_drain_col (:) ! col sub-surface runoff (mm H2O /s) real(r8), pointer :: qflx_latflow_in_col (:) ! col hillslope lateral flow input (mm/s) real(r8), pointer :: qflx_latflow_out_col (:) ! col hillslope lateral flow output (mm/s) - real(r8), pointer :: qdischarge_col (:) ! col hillslope discharge (m3/s) - real(r8), pointer :: qstreamflow_lun (:) ! lun stream discharge (m3/s) + real(r8), pointer :: volumetric_discharge_col (:) ! col hillslope discharge (m3/s) + real(r8), pointer :: volumetric_streamflow_lun(:) ! lun stream discharge (m3/s) real(r8), pointer :: qflx_drain_perched_col (:) ! col sub-surface runoff from perched wt (mm H2O /s) real(r8), pointer :: qflx_top_soil_col (:) ! col net water input into soil from top (mm/s) real(r8), pointer :: qflx_floodc_col (:) ! col flood water flux at column level @@ -288,10 +288,10 @@ subroutine InitAllocate(this, bounds, tracer_vars) call AllocateVar1d(var = this%qflx_latflow_out_col, name = 'qflx_latflow_out_col', & container = tracer_vars, & bounds = bounds, subgrid_level = subgrid_level_column) - call AllocateVar1d(var = this%qdischarge_col, name = 'qdischarge_col', & + call AllocateVar1d(var = this%volumetric_discharge_col, name = 'volumetric_discharge_col', & container = tracer_vars, & bounds = bounds, subgrid_level = subgrid_level_column) - call AllocateVar1d(var = this%qstreamflow_lun, name = 'qstreamflow_lun', & + call AllocateVar1d(var = this%volumetric_streamflow_lun, name = 'volumetric_streamflow_lun', & container = tracer_vars, & bounds = bounds, subgrid_level = subgrid_level_landunit) call AllocateVar1d(var = this%qflx_top_soil_col, name = 'qflx_top_soil_col', & @@ -512,24 +512,24 @@ subroutine InitHistory(this, bounds) l2g_scale_type='natveg', c2l_scale_type='urbanf', & ptr_col=this%qflx_latflow_out_col) - this%qdischarge_col(begc:endc) = spval + this%volumetric_discharge_col(begc:endc) = spval call hist_addfld1d ( & - fname=this%info%fname('QDISCHARGE'), & + fname=this%info%fname('VOLUMETRIC_DISCHARGE'), & units='m3/s', & avgflag='A', & long_name=this%info%lname('hillslope discharge from column'), & l2g_scale_type='natveg', c2l_scale_type='urbanf', & - ptr_col=this%qdischarge_col) + ptr_col=this%volumetric_discharge_col) if (use_hillslope_routing) then - this%qstreamflow_lun(begl:endl) = spval + this%volumetric_streamflow_lun(begl:endl) = spval call hist_addfld1d ( & - fname=this%info%fname('QSTREAMFLOW'), & + fname=this%info%fname('VOLUMETRIC_STREAMFLOW'), & units='m3/s', & avgflag='A', & long_name=this%info%lname('streamflow discharge'), & l2g_scale_type='natveg', & - ptr_lunit=this%qstreamflow_lun) + ptr_lunit=this%volumetric_streamflow_lun) endif endif @@ -915,13 +915,13 @@ subroutine InitCold(this, bounds) this%qflx_surf_col(c) = 0._r8 this%qflx_latflow_in_col(c) = 0._r8 this%qflx_latflow_out_col(c) = 0._r8 - this%qdischarge_col(c) = 0._r8 + this%volumetric_discharge_col(c) = 0._r8 end if end do if (use_hillslope_routing) then do l = bounds%begl, bounds%endl if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then - this%qstreamflow_lun(l) = 0._r8 + this%volumetric_streamflow_lun(l) = 0._r8 end if end do endif diff --git a/src/main/TopoMod.F90 b/src/main/TopoMod.F90 index 5e054ef00c..4e09d4886f 100644 --- a/src/main/TopoMod.F90 +++ b/src/main/TopoMod.F90 @@ -275,19 +275,26 @@ subroutine UpdateTopo(this, bounds, num_icec, filter_icec, & ! This could operate over a filter like 'allc' in order to just operate over active ! points, but I'm not sure that would speed things up much, and would require passing ! in this additional filter. + do c = bounds%begc, bounds%endc if (.not. this%needs_downscaling_col(c)) then + ! For any point that isn't already set to be downscaled, set its topo value to the + ! atmosphere's topographic height. This is important for the hillslope block + ! below. For non-hillslope columns, this shouldn't matter, but is useful if + ! topo_col is written to the history file. g = col%gridcell(c) - l = col%landunit(c) - this%topo_col(c) = atm_topo(g) - - if (col%is_hillslope_column(c) .and. downscale_hillslope_meteorology) then - this%topo_col(c) = this%topo_col(c) & - + (col%hill_elev(c) - mean_hillslope_elevation(l)) - this%needs_downscaling_col(c) = .true. - endif end if + ! If needs_downscaling_col was already set, then that implies + ! that topo_col was previously set by update_glc2lnd_topo. + ! In that case, topo_col should be used as a starting point, + ! rather than the atmosphere's topo value. + if (col%is_hillslope_column(c) .and. downscale_hillslope_meteorology) then + l = col%landunit(c) + this%topo_col(c) = this%topo_col(c) & + + (col%hill_elev(c) - mean_hillslope_elevation(l)) + this%needs_downscaling_col(c) = .true. + endif end do call glc_behavior%update_glc_classes(bounds, this%topo_col(begc:endc)) diff --git a/src/main/lnd2atmMod.F90 b/src/main/lnd2atmMod.F90 index 3ed2fa8c7e..8fb8a5fe2a 100644 --- a/src/main/lnd2atmMod.F90 +++ b/src/main/lnd2atmMod.F90 @@ -347,7 +347,7 @@ subroutine lnd2atm(bounds, & g = lun%gridcell(l) water_inst%waterlnd2atmbulk_inst%qflx_rofliq_stream_grc(g) = & water_inst%waterlnd2atmbulk_inst%qflx_rofliq_stream_grc(g) & - + water_inst%waterfluxbulk_inst%qstreamflow_lun(l) & + + water_inst%waterfluxbulk_inst%volumetric_streamflow_lun(l) & *1e3_r8/(grc%area(g)*1.e6_r8) enddo