From 26cc6b1340dd7bec3bf85f79c30ded2faf6a024a Mon Sep 17 00:00:00 2001 From: Dustin Swales Date: Fri, 6 Dec 2019 11:00:33 -0700 Subject: [PATCH] Cleaned up daytime masking in SW calculation --- physics/GFS_rrtmgp_pre.F90 | 7 +- physics/GFS_rrtmgp_sw_post.F90 | 337 ++++++++++++++--------------- physics/GFS_rrtmgp_sw_post.meta | 2 +- physics/GFS_rrtmgp_sw_pre.F90 | 7 +- physics/rrtmgp_sw_cloud_optics.F90 | 208 +++++++++--------- physics/rrtmgp_sw_gas_optics.F90 | 60 +++-- physics/rrtmgp_sw_gas_optics.meta | 16 ++ physics/rrtmgp_sw_rte.F90 | 129 +++++------ physics/rrtmgp_sw_rte.meta | 2 +- 9 files changed, 388 insertions(+), 380 deletions(-) diff --git a/physics/GFS_rrtmgp_pre.F90 b/physics/GFS_rrtmgp_pre.F90 index 4584fe2ac..a80c4b5a6 100644 --- a/physics/GFS_rrtmgp_pre.F90 +++ b/physics/GFS_rrtmgp_pre.F90 @@ -216,7 +216,8 @@ subroutine GFS_rrtmgp_pre_run (Model, Grid, Statein, Coupling, Radtend, Sfcprop, ! ####################################################################################### ! Water-vapor mixing-ratio - q_lay(1:ncol,:) = max( 1.e-6, Statein%qgrs(1:NCOL,:,1)) + q_lay(1:ncol,:) = Statein%qgrs(1:NCOL,:,1) + where(q_lay .lt. 1.e-6) q_lay = 1.e-6 ! Pressure at layer-interface p_lev(1:NCOL,:) = Statein%prsi(1:NCOL,:) @@ -227,7 +228,6 @@ subroutine GFS_rrtmgp_pre_run (Model, Grid, Statein, Coupling, Radtend, Sfcprop, ! Temperature at layer-center t_lay(1:NCOL,:) = Statein%tgrs(1:NCOL,:) - ! ! Temperature at layer-interfaces if (top_at_1) then t_lev(1:NCOL,1) = t_lay(1:NCOL,iTOA) @@ -260,7 +260,8 @@ subroutine GFS_rrtmgp_pre_run (Model, Grid, Statein, Coupling, Radtend, Sfcprop, ! ####################################################################################### ! First recast remaining all tracers (except sphum) forcing them all to be positive do j = 2, model%NTRAC - tracer(1:NCOL,:,j) = max(0.0, Statein%qgrs(1:NCOL,:,j)) + tracer(1:NCOL,:,j) = Statein%qgrs(1:NCOL,:,j) + where(tracer(:,:,j) .lt. 0.0) tracer(:,:,j) = 0._kind_phys enddo if (Model%ntoz > 0) then diff --git a/physics/GFS_rrtmgp_sw_post.F90 b/physics/GFS_rrtmgp_sw_post.F90 index 27edd06b7..72e841b07 100644 --- a/physics/GFS_rrtmgp_sw_post.F90 +++ b/physics/GFS_rrtmgp_sw_post.F90 @@ -1,17 +1,10 @@ -!>\file GFS_rrtmgp_sw_post -!!This file contains module GFS_rrtmgp_sw_post - use machine, only: kind_phys - use GFS_typedefs, only: GFS_coupling_type, & - GFS_control_type, & - GFS_grid_type, & - GFS_radtend_type, & - GFS_diag_type, & - GFS_statein_type, & - GFS_interstitial_type + use machine, only: kind_phys + use GFS_typedefs, only: GFS_coupling_type, GFS_control_type, GFS_grid_type, & + GFS_radtend_type, GFS_diag_type, GFS_statein_type, & + GFS_interstitial_type use module_radiation_aerosols, only: NSPC1 use module_radsw_parameters, only: topfsw_type, sfcfsw_type, profsw_type, cmpfsw_type - ! RRTMGP DDT's use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp use mo_fluxes_byband, only: ty_fluxes_byband use mo_heating_rates, only: compute_heating_rate @@ -35,7 +28,7 @@ end subroutine GFS_rrtmgp_sw_post_init !! \htmlinclude GFS_rrtmgp_sw_post.html !! subroutine GFS_rrtmgp_sw_post_run (Model, Interstitial, Grid, Diag, Radtend, Coupling, & - Statein, scmpsw, im, p_lev, sw_gas_props, nday, idxday, fluxswUP_allsky, & + Statein, scmpsw, nCol, p_lev, sw_gas_props, nday, idxday, fluxswUP_allsky, & fluxswDOWN_allsky, fluxswUP_clrsky, fluxswDOWN_clrsky, raddt, aerodp, cldsa, mbota, & mtopa, cld_frac, cldtausw, flxprf_sw, hsw0, errmsg, errflg) @@ -55,32 +48,32 @@ subroutine GFS_rrtmgp_sw_post_run (Model, Interstitial, Grid, Diag, Radtend, Cou type(GFS_statein_type), intent(in) :: & Statein ! Fortran DDT: FV3-GFS prognostic state data in from dycore integer, intent(in) :: & - im, & ! Horizontal loop extent + nCol, & ! Horizontal loop extent nDay ! Number of daylit columns integer, intent(in), dimension(nday) :: & idxday ! Index array for daytime points type(ty_gas_optics_rrtmgp),intent(in) :: & sw_gas_props ! DDT containing SW spectral information - real(kind_phys), dimension(size(Grid%xlon,1), Model%levs+1), intent(in) :: & + real(kind_phys), dimension(nCol, Model%levs+1), intent(in) :: & p_lev ! Pressure @ model layer-interfaces (hPa) - real(kind_phys), dimension(size(Grid%xlon,1), Model%levs+1), intent(in) :: & + real(kind_phys), dimension(nCol, Model%levs+1), intent(in) :: & fluxswUP_allsky, & ! SW All-sky flux (W/m2) fluxswDOWN_allsky, & ! SW All-sky flux (W/m2) fluxswUP_clrsky, & ! SW Clear-sky flux (W/m2) fluxswDOWN_clrsky ! SW All-sky flux (W/m2) real(kind_phys), intent(in) :: & raddt ! Radiation time step - real(kind_phys), dimension(im,NSPC1), intent(in) :: & + real(kind_phys), dimension(nCol,NSPC1), intent(in) :: & aerodp ! Vertical integrated optical depth for various aerosol species - real(kind_phys), dimension(im,5), intent(in) :: & + real(kind_phys), dimension(nCol,5), intent(in) :: & cldsa ! Fraction of clouds for low, middle, high, total and BL - integer, dimension(im,3), intent(in) ::& + integer, dimension(nCol,3), intent(in) ::& mbota, & ! vertical indices for low, middle and high cloud tops mtopa ! vertical indices for low, middle and high cloud bases - real(kind_phys), dimension(im,Model%levs), intent(in) :: & + real(kind_phys), dimension(nCol,Model%levs), intent(in) :: & cld_frac, & ! Total cloud fraction in each layer cldtausw ! approx .55mu band layer cloud optical depth - real(kind_phys),dimension(size(Grid%xlon,1), Model%levs) :: & + real(kind_phys),dimension(nCol, Model%levs) :: & hswc ! All-sky heating rates (K/s) ! Outputs (mandatory) @@ -90,15 +83,15 @@ subroutine GFS_rrtmgp_sw_post_run (Model, Interstitial, Grid, Diag, Radtend, Cou errflg ! Outputs (optional) - real(kind_phys), dimension(size(Grid%xlon,1), Model%levs), optional, intent(inout) :: & + real(kind_phys), dimension(nCol, Model%levs), optional, intent(inout) :: & hsw0 ! Shortwave clear-sky heating-rate (K/sec) - type(profsw_type), dimension(size(Grid%xlon,1), Model%levs+1), intent(inout), optional :: & + type(profsw_type), dimension(nCol, Model%levs+1), intent(inout), optional :: & flxprf_sw ! 2D radiative fluxes, components: ! upfxc - total sky upward flux (W/m2) ! dnfxc - total sky dnward flux (W/m2) ! upfx0 - clear sky upward flux (W/m2) ! dnfx0 - clear sky dnward flux (W/m2) - type(cmpfsw_type), dimension(size(Grid%xlon,1)), intent(inout), optional :: & + type(cmpfsw_type), dimension(nCol), intent(inout), optional :: & scmpsw ! 2D surface fluxes, components: ! uvbfc - total sky downward uv-b flux at (W/m2) ! uvbf0 - clear sky downward uv-b flux at (W/m2) @@ -112,61 +105,50 @@ subroutine GFS_rrtmgp_sw_post_run (Model, Interstitial, Grid, Diag, Radtend, Cou real(kind_phys), dimension(nDay, Model%levs) :: thetaTendClrSky, thetaTendAllSky logical :: l_clrskysw_hr, l_fluxessw2d, top_at_1, l_sfcFluxessw1D - ! Initialize CCPP error handling variables + ! Initialize CCPP error handling variables errmsg = '' errflg = 0 - ! Are any optional outputs requested? - l_clrskysw_hr = present(hsw0) - l_fluxessw2d = present(flxprf_sw) - l_sfcfluxessw1D = present(scmpsw) + if (.not. Model%lsswr) return + if (nDay .gt. 0) then - ! ####################################################################################### - ! What is vertical ordering? - ! ####################################################################################### - top_at_1 = (p_lev(1,1) .lt. p_lev(1, Model%levs)) - if (top_at_1) then - iSFC = Model%levs+1 - iTOA = 1 - else - iSFC = 1 - iTOA = Model%levs+1 - endif - - ! ####################################################################################### - ! Compute SW heating-rates - ! ####################################################################################### - ! Initialize - hswc = 0 - Diag%topfsw = topfsw_type ( 0., 0., 0. ) - Radtend%sfcfsw = sfcfsw_type ( 0., 0., 0., 0. ) - if (l_clrskysw_hr) then - hsw0(:,:) = 0. - endif - if (l_fluxessw2D) then - flxprf_sw = profsw_type ( 0., 0., 0., 0. ) - endif - if (l_sfcfluxessw1D) then - scmpsw = cmpfsw_type (0.,0.,0.,0.,0.,0.) - endif + ! Are any optional outputs requested? + l_clrskysw_hr = present(hsw0) + l_fluxessw2d = present(flxprf_sw) + l_sfcfluxessw1D = present(scmpsw) - if (Model%lsswr .and. nDay .gt. 0) then + ! ####################################################################################### + ! What is vertical ordering? + ! ####################################################################################### + top_at_1 = (p_lev(1,1) .lt. p_lev(1, Model%levs)) + if (top_at_1) then + iSFC = Model%levs+1 + iTOA = 1 + else + iSFC = 1 + iTOA = Model%levs+1 + endif + + ! ####################################################################################### + ! Compute SW heating-rates + ! ####################################################################################### ! Clear-sky heating-rate (optional) if (l_clrskysw_HR) then call check_error_msg('GFS_rrtmgp_post',compute_heating_rate( & - fluxswUP_clrsky(idxday,:), & - fluxswDOWN_clrsky(idxday,:), & - p_lev(idxday,:), & - thetaTendClrSky)) - hsw0(idxday,:)=thetaTendClrSky + fluxswUP_clrsky(idxday(1:nDay),:), & ! IN - Shortwave upward clear-sky flux profiles (W/m2) + fluxswDOWN_clrsky(idxday(1:nDay),:), & ! IN - Shortwave downward clear-sky flux profiles (W/m2) + p_lev(idxday(1:nDay),:), & ! IN - Pressure at model-interface (Pa) + thetaTendClrSky)) ! OUT - Clear-sky heating-rate (K/sec) + hsw0(idxday(1:nDay),:)=thetaTendClrSky endif + ! All-sky heating-rate (mandatory) call check_error_msg('GFS_rrtmgp_post',compute_heating_rate( & - fluxswUP_allsky(idxday,:), & - fluxswDOWN_allsky(idxday,:), & - p_lev(idxday,:), & - thetaTendAllSky)) - hswc(idxday,:) = thetaTendAllSky + fluxswUP_allsky(idxday(1:nDay),:), & ! IN - Shortwave upward all-sky flux profiles (W/m2) + fluxswDOWN_allsky(idxday(1:nDay),:), & ! IN - Shortwave downward all-sky flux profiles (W/m2) + p_lev(idxday(1:nDay),:), & ! IN - Pressure at model-interface (Pa) + thetaTendAllSky)) ! OUT - All-sky heating-rate (K/sec) + hswc(idxday(1:nDay),:) = thetaTendAllSky ! Copy fluxes from RRTGMP types into model radiation types. ! Mandatory outputs @@ -185,65 +167,65 @@ subroutine GFS_rrtmgp_sw_post_run (Model, Interstitial, Grid, Diag, Radtend, Cou flxprf_sw(:,:)%upfx0 = fluxswUP_clrsky(:,:) flxprf_sw(:,:)%dnfx0 = fluxswDOWN_clrsky(:,:) endif - endif - ! ####################################################################################### - ! Save SW outputs - ! ####################################################################################### - if (Model%lsswr) then - if (nday > 0) then - ! All-sky heating rate + + ! ####################################################################################### + ! Save SW outputs + ! ####################################################################################### + ! All-sky heating rate + do k = 1, Model%levs + Radtend%htrsw(1:nCol,k) = hswc(1:nCol,k) + enddo + ! Clear-sky heating rate + if (Model%swhtr) then do k = 1, Model%levs - Radtend%htrsw(1:im,k) = hswc(1:im,k) - enddo - ! Clear-sky heating rate - if (Model%swhtr) then - do k = 1, Model%levs - Radtend%swhc(1:im,k) = hsw0(1:im,k) - enddo - endif - - ! Surface down and up spectral component fluxes - ! - Save two spectral bands' surface downward and upward fluxes for output. - do i=1,im - Coupling%nirbmdi(i) = scmpsw(i)%nirbm - Coupling%nirdfdi(i) = scmpsw(i)%nirdf - Coupling%visbmdi(i) = scmpsw(i)%visbm - Coupling%visdfdi(i) = scmpsw(i)%visdf - - Coupling%nirbmui(i) = scmpsw(i)%nirbm * Interstitial%sfc_alb_nir_dir(1,i) - Coupling%nirdfui(i) = scmpsw(i)%nirdf * Interstitial%sfc_alb_nir_dif(1,i) - Coupling%visbmui(i) = scmpsw(i)%visbm * Interstitial%sfc_alb_uvvis_dir(1,i) - Coupling%visdfui(i) = scmpsw(i)%visdf * Interstitial%sfc_alb_uvvis_dif(1,i) + Radtend%swhc(1:nCol,k) = hsw0(1:nCol,k) enddo - else ! if_nday_block - Radtend%htrsw(:,:) = 0.0 - Radtend%sfcfsw = sfcfsw_type( 0.0, 0.0, 0.0, 0.0 ) - Diag%topfsw = topfsw_type( 0.0, 0.0, 0.0 ) - scmpsw = cmpfsw_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ) + endif + + ! Surface down and up spectral component fluxes + ! - Save two spectral bands' surface downward and upward fluxes for output. + do i=1,nCol + Coupling%nirbmdi(i) = scmpsw(i)%nirbm + Coupling%nirdfdi(i) = scmpsw(i)%nirdf + Coupling%visbmdi(i) = scmpsw(i)%visbm + Coupling%visdfdi(i) = scmpsw(i)%visdf - do i=1,im - Coupling%nirbmdi(i) = 0.0 - Coupling%nirdfdi(i) = 0.0 - Coupling%visbmdi(i) = 0.0 - Coupling%visdfdi(i) = 0.0 - - Coupling%nirbmui(i) = 0.0 - Coupling%nirdfui(i) = 0.0 - Coupling%visbmui(i) = 0.0 - Coupling%visdfui(i) = 0.0 - enddo + Coupling%nirbmui(i) = scmpsw(i)%nirbm * Interstitial%sfc_alb_nir_dir(1,i) + Coupling%nirdfui(i) = scmpsw(i)%nirdf * Interstitial%sfc_alb_nir_dif(1,i) + Coupling%visbmui(i) = scmpsw(i)%visbm * Interstitial%sfc_alb_uvvis_dir(1,i) + Coupling%visdfui(i) = scmpsw(i)%visdf * Interstitial%sfc_alb_uvvis_dif(1,i) + enddo + else ! if_nday_block + ! ####################################################################################### + ! Save SW outputs + ! ####################################################################################### + Radtend%htrsw(:,:) = 0.0 + Radtend%sfcfsw = sfcfsw_type( 0.0, 0.0, 0.0, 0.0 ) + Diag%topfsw = topfsw_type( 0.0, 0.0, 0.0 ) + scmpsw = cmpfsw_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ) + + do i=1,nCol + Coupling%nirbmdi(i) = 0.0 + Coupling%nirdfdi(i) = 0.0 + Coupling%visbmdi(i) = 0.0 + Coupling%visdfdi(i) = 0.0 - if (Model%swhtr) then - Radtend%swhc(:,:) = 0 - endif - endif ! end_if_nday + Coupling%nirbmui(i) = 0.0 + Coupling%nirdfui(i) = 0.0 + Coupling%visbmui(i) = 0.0 + Coupling%visdfui(i) = 0.0 + enddo - ! Radiation fluxes for other physics processes - do i=1,im - Coupling%sfcnsw(i) = Radtend%sfcfsw(i)%dnfxc - Radtend%sfcfsw(i)%upfxc - Coupling%sfcdsw(i) = Radtend%sfcfsw(i)%dnfxc - enddo - endif ! end_if_lsswr + if (Model%swhtr) then + Radtend%swhc(:,:) = 0 + endif + endif ! end_if_nday + + ! Radiation fluxes for other physics processes + do i=1,nCol + Coupling%sfcnsw(i) = Radtend%sfcfsw(i)%dnfxc - Radtend%sfcfsw(i)%upfxc + Coupling%sfcdsw(i) = Radtend%sfcfsw(i)%dnfxc + enddo ! ####################################################################################### ! Save SW diagnostics @@ -254,68 +236,63 @@ subroutine GFS_rrtmgp_sw_post_run (Model, Interstitial, Grid, Diag, Radtend, Cou ! - Collect the fluxr data for wrtsfc ! ####################################################################################### if (Model%lssav) then - if (Model%lsswr) then - do i=1,im - Diag%fluxr(i,34) = Diag%fluxr(i,34) + Model%fhswr*aerodp(i,1) ! total aod at 550nm - Diag%fluxr(i,35) = Diag%fluxr(i,35) + Model%fhswr*aerodp(i,2) ! DU aod at 550nm - Diag%fluxr(i,36) = Diag%fluxr(i,36) + Model%fhswr*aerodp(i,3) ! BC aod at 550nm - Diag%fluxr(i,37) = Diag%fluxr(i,37) + Model%fhswr*aerodp(i,4) ! OC aod at 550nm - Diag%fluxr(i,38) = Diag%fluxr(i,38) + Model%fhswr*aerodp(i,5) ! SU aod at 550nm - Diag%fluxr(i,39) = Diag%fluxr(i,39) + Model%fhswr*aerodp(i,6) ! SS aod at 550nm - if (Radtend%coszen(i) > 0.) then - ! SW all-sky fluxes - tem0d = Model%fhswr * Radtend%coszdg(i) / Radtend%coszen(i) - Diag%fluxr(i,2 ) = Diag%fluxr(i,2) + Diag%topfsw(i)%upfxc * tem0d ! total sky top sw up - Diag%fluxr(i,3 ) = Diag%fluxr(i,3) + Radtend%sfcfsw(i)%upfxc * tem0d - Diag%fluxr(i,4 ) = Diag%fluxr(i,4) + Radtend%sfcfsw(i)%dnfxc * tem0d ! total sky sfc sw dn - ! SW uv-b fluxes - Diag%fluxr(i,21) = Diag%fluxr(i,21) + scmpsw(i)%uvbfc * tem0d ! total sky uv-b sw dn - Diag%fluxr(i,22) = Diag%fluxr(i,22) + scmpsw(i)%uvbf0 * tem0d ! clear sky uv-b sw dn - ! SW TOA incoming fluxes - Diag%fluxr(i,23) = Diag%fluxr(i,23) + Diag%topfsw(i)%dnfxc * tem0d ! top sw dn - ! SW SFC flux components - Diag%fluxr(i,24) = Diag%fluxr(i,24) + scmpsw(i)%visbm * tem0d ! uv/vis beam sw dn - Diag%fluxr(i,25) = Diag%fluxr(i,25) + scmpsw(i)%visdf * tem0d ! uv/vis diff sw dn - Diag%fluxr(i,26) = Diag%fluxr(i,26) + scmpsw(i)%nirbm * tem0d ! nir beam sw dn - Diag%fluxr(i,27) = Diag%fluxr(i,27) + scmpsw(i)%nirdf * tem0d ! nir diff sw dn - ! SW clear-sky fluxes - Diag%fluxr(i,29) = Diag%fluxr(i,29) + Diag%topfsw(i)%upfx0 * tem0d - Diag%fluxr(i,31) = Diag%fluxr(i,31) + Radtend%sfcfsw(i)%upfx0 * tem0d - Diag%fluxr(i,32) = Diag%fluxr(i,32) + Radtend%sfcfsw(i)%dnfx0 * tem0d - - endif - enddo - - ! Save total and boundary-layer clouds - do i=1,im - Diag%fluxr(i,17) = Diag%fluxr(i,17) + raddt * cldsa(i,4) - Diag%fluxr(i,18) = Diag%fluxr(i,18) + raddt * cldsa(i,5) - enddo - - ! Save cld frac,toplyr,botlyr and top temp, note that the order of h,m,l cloud - ! is reversed for the fluxr output. save interface pressure (pa) of top/bot - do j = 1, 3 - do i = 1, im - tem0d = raddt * cldsa(i,j) - itop = mtopa(i,j) - ibtc = mbota(i,j) - Diag%fluxr(i, 8-j) = Diag%fluxr(i, 8-j) + tem0d - Diag%fluxr(i,11-j) = Diag%fluxr(i,11-j) + tem0d * Statein%prsi(i,itop) - Diag%fluxr(i,14-j) = Diag%fluxr(i,14-j) + tem0d * Statein%prsi(i,ibtc) - Diag%fluxr(i,17-j) = Diag%fluxr(i,17-j) + tem0d * Statein%tgrs(i,itop) - - ! Add optical depth and emissivity output - tem1 = 0. - do k=ibtc,itop - tem1 = tem1 + cldtausw(i,k) ! approx .55 mu channel - enddo - Diag%fluxr(i,43-j) = Diag%fluxr(i,43-j) + tem0d * tem1 + do i=1,nCol + Diag%fluxr(i,34) = Diag%fluxr(i,34) + Model%fhswr*aerodp(i,1) ! total aod at 550nm + Diag%fluxr(i,35) = Diag%fluxr(i,35) + Model%fhswr*aerodp(i,2) ! DU aod at 550nm + Diag%fluxr(i,36) = Diag%fluxr(i,36) + Model%fhswr*aerodp(i,3) ! BC aod at 550nm + Diag%fluxr(i,37) = Diag%fluxr(i,37) + Model%fhswr*aerodp(i,4) ! OC aod at 550nm + Diag%fluxr(i,38) = Diag%fluxr(i,38) + Model%fhswr*aerodp(i,5) ! SU aod at 550nm + Diag%fluxr(i,39) = Diag%fluxr(i,39) + Model%fhswr*aerodp(i,6) ! SS aod at 550nm + if (Radtend%coszen(i) > 0.) then + ! SW all-sky fluxes + tem0d = Model%fhswr * Radtend%coszdg(i) / Radtend%coszen(i) + Diag%fluxr(i,2 ) = Diag%fluxr(i,2) + Diag%topfsw(i)%upfxc * tem0d ! total sky top sw up + Diag%fluxr(i,3 ) = Diag%fluxr(i,3) + Radtend%sfcfsw(i)%upfxc * tem0d + Diag%fluxr(i,4 ) = Diag%fluxr(i,4) + Radtend%sfcfsw(i)%dnfxc * tem0d ! total sky sfc sw dn + ! SW uv-b fluxes + Diag%fluxr(i,21) = Diag%fluxr(i,21) + scmpsw(i)%uvbfc * tem0d ! total sky uv-b sw dn + Diag%fluxr(i,22) = Diag%fluxr(i,22) + scmpsw(i)%uvbf0 * tem0d ! clear sky uv-b sw dn + ! SW TOA incoming fluxes + Diag%fluxr(i,23) = Diag%fluxr(i,23) + Diag%topfsw(i)%dnfxc * tem0d ! top sw dn + ! SW SFC flux components + Diag%fluxr(i,24) = Diag%fluxr(i,24) + scmpsw(i)%visbm * tem0d ! uv/vis beam sw dn + Diag%fluxr(i,25) = Diag%fluxr(i,25) + scmpsw(i)%visdf * tem0d ! uv/vis diff sw dn + Diag%fluxr(i,26) = Diag%fluxr(i,26) + scmpsw(i)%nirbm * tem0d ! nir beam sw dn + Diag%fluxr(i,27) = Diag%fluxr(i,27) + scmpsw(i)%nirdf * tem0d ! nir diff sw dn + ! SW clear-sky fluxes + Diag%fluxr(i,29) = Diag%fluxr(i,29) + Diag%topfsw(i)%upfx0 * tem0d + Diag%fluxr(i,31) = Diag%fluxr(i,31) + Radtend%sfcfsw(i)%upfx0 * tem0d + Diag%fluxr(i,32) = Diag%fluxr(i,32) + Radtend%sfcfsw(i)%dnfx0 * tem0d + endif + enddo + + ! Save total and boundary-layer clouds + do i=1,nCol + Diag%fluxr(i,17) = Diag%fluxr(i,17) + raddt * cldsa(i,4) + Diag%fluxr(i,18) = Diag%fluxr(i,18) + raddt * cldsa(i,5) + enddo + + ! Save cld frac,toplyr,botlyr and top temp, note that the order of h,m,l cloud + ! is reversed for the fluxr output. save interface pressure (pa) of top/bot + do j = 1, 3 + do i = 1, nCol + tem0d = raddt * cldsa(i,j) + itop = mtopa(i,j) + ibtc = mbota(i,j) + Diag%fluxr(i, 8-j) = Diag%fluxr(i, 8-j) + tem0d + Diag%fluxr(i,11-j) = Diag%fluxr(i,11-j) + tem0d * Statein%prsi(i,itop) + Diag%fluxr(i,14-j) = Diag%fluxr(i,14-j) + tem0d * Statein%prsi(i,ibtc) + Diag%fluxr(i,17-j) = Diag%fluxr(i,17-j) + tem0d * Statein%tgrs(i,itop) + + ! Add optical depth and emissivity output + tem1 = 0. + do k=ibtc,itop + tem1 = tem1 + cldtausw(i,k) ! approx .55 mu channel enddo + Diag%fluxr(i,43-j) = Diag%fluxr(i,43-j) + tem0d * tem1 enddo - endif + enddo endif - - end subroutine GFS_rrtmgp_sw_post_run ! ######################################################################################### diff --git a/physics/GFS_rrtmgp_sw_post.meta b/physics/GFS_rrtmgp_sw_post.meta index bff311cb7..bf71798bd 100644 --- a/physics/GFS_rrtmgp_sw_post.meta +++ b/physics/GFS_rrtmgp_sw_post.meta @@ -65,7 +65,7 @@ type = cmpfsw_type intent = inout optional = T -[im] +[ncol] standard_name = horizontal_loop_extent long_name = horizontal loop extent units = count diff --git a/physics/GFS_rrtmgp_sw_pre.F90 b/physics/GFS_rrtmgp_sw_pre.F90 index 8c4b5b7fd..e20f8224f 100644 --- a/physics/GFS_rrtmgp_sw_pre.F90 +++ b/physics/GFS_rrtmgp_sw_pre.F90 @@ -111,10 +111,9 @@ subroutine GFS_rrtmgp_sw_pre_run(Model, Interstitial, Grid, Sfcprop, Statein, nc Radtend%coszen, Radtend%coszdg) ! ####################################################################################### - ! For SW, gather daylit points, compute surface albedo in each band, + ! For SW gather daylit points ! ####################################################################################### - ! Check for daytime points for SW radiation. - nday = 0 + nday = 0 idxday = 0 do i = 1, NCOL if (Radtend%coszen(i) >= 0.0001) then @@ -137,7 +136,9 @@ subroutine GFS_rrtmgp_sw_pre_run(Model, Interstitial, Grid, Sfcprop, Statein, nc endif endif + ! ####################################################################################### ! Call module_radiation_surface::setalb() to setup surface albedo. + ! ####################################################################################### call setalb (Sfcprop%slmsk, Sfcprop%snowd, Sfcprop%sncovr, Sfcprop%snoalb, Sfcprop%zorl, & Radtend%coszen, Sfcprop%tsfc, Sfcprop%tsfc, Sfcprop%hprime(:,1), Sfcprop%alvsf, & Sfcprop%alnsf, Sfcprop%alvwf, Sfcprop%alnwf, Sfcprop%facsf, Sfcprop%facwf, & diff --git a/physics/rrtmgp_sw_cloud_optics.F90 b/physics/rrtmgp_sw_cloud_optics.F90 index 48dde613d..c29bb2142 100644 --- a/physics/rrtmgp_sw_cloud_optics.F90 +++ b/physics/rrtmgp_sw_cloud_optics.F90 @@ -14,12 +14,12 @@ module rrtmgp_sw_cloud_optics contains -!! \section arg_table_rrtmgp_sw_cloud_optics_init -!! \htmlinclude rrtmgp_lw_cloud_optics.html -!! ! ######################################################################################### ! SUBROUTINE sw_cloud_optics_init ! ######################################################################################### +!! \section arg_table_rrtmgp_sw_cloud_optics_init +!! \htmlinclude rrtmgp_lw_cloud_optics.html +!! subroutine rrtmgp_sw_cloud_optics_init(Model,mpicomm, mpirank, mpiroot, sw_cloud_props, & errmsg, errflg) use netcdf @@ -337,12 +337,12 @@ subroutine rrtmgp_sw_cloud_optics_init(Model,mpicomm, mpirank, mpiroot, sw_cloud endif end subroutine rrtmgp_sw_cloud_optics_init -!! \section arg_table_rrtmgp_sw_cloud_optics_run -!! \htmlinclude rrtmgp_sw_cloud_optics.html -!! ! ######################################################################################### ! SUBROTUINE rrtmgp_sw_cloud_optics_run() ! ######################################################################################### +!! \section arg_table_rrtmgp_sw_cloud_optics_run +!! \htmlinclude rrtmgp_sw_cloud_optics.html +!! subroutine rrtmgp_sw_cloud_optics_run(Model, ncol, icseed_sw, cld_frac, cld_lwp, & cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, cld_rerain, aerosolssw, & sw_cloud_props, sw_gas_props, ipsdsw0, nday, idxday, sw_optical_props_clouds, & @@ -355,7 +355,7 @@ subroutine rrtmgp_sw_cloud_optics_run(Model, ncol, icseed_sw, cld_frac, cld_lwp, ncol, & ! Number of horizontal gridpoints nday, & ! Number of daylit points. ipsdsw0 ! Initial permutation seed for McICA - integer,intent(in),dimension(nday) :: & + integer,intent(in),dimension(ncol) :: & idxday ! Indices for daylit points. integer,intent(in),dimension(ncol) :: & icseed_sw ! auxiliary special cloud related array when module @@ -393,7 +393,7 @@ subroutine rrtmgp_sw_cloud_optics_run(Model, ncol, icseed_sw, cld_frac, cld_lwp, ! Local variables integer :: iCol integer,dimension(ncol) :: ipseed_sw - logical,dimension(ncol,model%levs) :: liqmask, icemask + logical,dimension(nday,model%levs) :: liqmask, icemask type(ty_optical_props_2str) :: sw_optical_props_cloudsByBand type(random_stat) :: rng_stat real(kind_phys), dimension(sw_gas_props%get_ngpt(),model%levs,ncol) :: rng3D @@ -407,102 +407,108 @@ subroutine rrtmgp_sw_cloud_optics_run(Model, ncol, icseed_sw, cld_frac, cld_lwp, errflg = 0 if (.not. Model%lsswr) return + if (nDay .gt. 0) then + + ! ####################################################################################### + ! Compute ice/liquid cloud masks, needed by rrtmgp_cloud_optics + ! ####################################################################################### + liqmask = (cld_frac(idxday(1:nday),:) .gt. 0 .and. cld_lwp(idxday(1:nday),:) .gt. 0) + icemask = (cld_frac(idxday(1:nday),:) .gt. 0 .and. cld_iwp(idxday(1:nday),:) .gt. 0) + + ! ####################################################################################### + ! Allocate space for RRTMGP DDTs containing cloud and aerosol radiative properties + ! ####################################################################################### + ! Cloud optics [nday,model%levs,nBands] + call check_error_msg('rrtmgp_sw_cloud_optics_run',sw_optical_props_cloudsByBand%alloc_2str(& + nday, model%levs, sw_gas_props%get_band_lims_wavenumber())) + ! Aerosol optics [nday,model%levs,nBands] + call check_error_msg('rrtmgp_sw_cloud_optics_run',sw_optical_props_aerosol%alloc_2str( & + nday, model%levs, sw_gas_props%get_band_lims_wavenumber())) + ! Cloud optics [nday,model%levs,nGpt] + call check_error_msg('rrtmgp_sw_cloud_optics_run',sw_optical_props_clouds%alloc_2str( & + nday, model%levs, sw_gas_props)) + + ! ####################################################################################### + ! Copy aerosol optical information to RRTMGP DDT + ! ####################################################################################### + sw_optical_props_aerosol%tau = aerosolssw(idxday(1:nday),:,:,1) + sw_optical_props_aerosol%ssa = aerosolssw(idxday(1:nday),:,:,2) + sw_optical_props_aerosol%g = aerosolssw(idxday(1:nday),:,:,3) + + ! ####################################################################################### + ! Compute cloud-optics for RTE. + ! ####################################################################################### + if (Model%rrtmgp_cld_optics .gt. 0) then + ! RRTMGP cloud-optics. + call check_error_msg('rrtmgp_sw_cloud_optics_run',sw_cloud_props%cloud_optics(& + nday, & ! IN - Number of daylit gridpoints + model%levs, & ! IN - Number of vertical layers + sw_cloud_props%get_nband(), & ! IN - Number of SW bands + Model%rrtmgp_nrghice, & ! IN - Number of ice-roughness categories + liqmask, & ! IN - Liquid-cloud mask + icemask, & ! IN - Ice-cloud mask + cld_lwp(idxday(1:nday),:), & ! IN - Cloud liquid water path + cld_iwp(idxday(1:nday),:), & ! IN - Cloud ice water path + cld_reliq(idxday(1:nday),:), & ! IN - Cloud liquid effective radius + cld_reice(idxday(1:nday),:), & ! IN - Cloud ice effective radius + sw_optical_props_cloudsByBand)) ! OUT - RRTMGP DDT: Shortwave optical properties, + ! in each band (tau,ssa,g) + else + ! RRTMG cloud-optics + if (any(cld_frac .gt. 0)) then + sw_optical_props_cloudsByBand%tau(:,:,:) = 0._kind_phys + sw_optical_props_cloudsByBand%ssa(:,:,:) = 0._kind_phys + sw_optical_props_cloudsByBand%g(:,:,:) = 0._kind_phys + call rrtmg_sw_cloud_optics(nday, model%levs, sw_gas_props%get_nband(), & + cld_lwp(idxday(1:nday),:), cld_reliq(idxday(1:nday),:), & + cld_iwp(idxday(1:nday),:), cld_reice(idxday(1:nday),:), & + cld_rwp(idxday(1:nday),:), cld_rerain(idxday(1:nday),:), & + cld_swp(idxday(1:nday),:), cld_resnow(idxday(1:nday),:), & + cld_frac(idxday(1:nday),:), tau_cld, ssa_cld, asy_cld) + sw_optical_props_cloudsByBand%tau(:,:,:) = tau_cld + sw_optical_props_cloudsByBand%ssa(:,:,:) = ssa_cld + sw_optical_props_cloudsByBand%g(:,:,:) = asy_cld + endif + endif - ! ####################################################################################### - ! Change random number seed value for each radiation invocation (isubcsw =1 or 2). - ! ####################################################################################### - if(isubcsw == 1) then ! advance prescribed permutation seed - do iCol = 1, ncol - ipseed_sw(iCol) = ipsdsw0 + iCol - enddo - elseif (isubcsw == 2) then ! use input array of permutaion seeds - do iCol = 1, ncol - ipseed_sw(iCol) = icseed_sw(iCol) - enddo - endif - - ! ####################################################################################### - ! Compute ice/liquid cloud masks, needed by rrtmgp_cloud_optics - ! ####################################################################################### - liqmask = (cld_frac .gt. 0 .and. cld_lwp .gt. 0) - icemask = (cld_frac .gt. 0 .and. cld_iwp .gt. 0) - - ! ####################################################################################### - ! Allocate space for RRTMGP DDTs containing cloud and aerosol radiative properties - ! ####################################################################################### - ! Cloud optics [ncol,model%levs,nBands] - call check_error_msg('rrtmgp_sw_cloud_optics_run',sw_optical_props_cloudsByBand%alloc_2str(& - ncol, model%levs, sw_gas_props%get_band_lims_wavenumber())) - ! Aerosol optics [ncol,model%levs,nBands] - call check_error_msg('rrtmgp_sw_cloud_optics_run',sw_optical_props_aerosol%alloc_2str( & - ncol, model%levs, sw_gas_props%get_band_lims_wavenumber())) - ! Cloud optics [ncol,model%levs,nGpt] - call check_error_msg('rrtmgp_sw_cloud_optics_run',sw_optical_props_clouds%alloc_2str( & - ncol, model%levs, sw_gas_props)) - - ! ####################################################################################### - ! Copy aerosol optical information to RRTMGP DDT - ! ####################################################################################### - sw_optical_props_aerosol%tau = aerosolssw(:,:,:,1) - sw_optical_props_aerosol%ssa = aerosolssw(:,:,:,2) - sw_optical_props_aerosol%g = aerosolssw(:,:,:,3) - - ! ####################################################################################### - ! Compute cloud-optics for RTE. - ! ####################################################################################### - if (Model%rrtmgp_cld_optics .gt. 0) then - ! RRTMGP cloud-optics. - call check_error_msg('rrtmgp_sw_cloud_optics_run',sw_cloud_props%cloud_optics(& - ncol, & ! IN - Number of daylit gridpoints - model%levs, & ! IN - Number of vertical layers - sw_cloud_props%get_nband(), & ! IN - Number of SW bands - Model%rrtmgp_nrghice, & ! IN - Number of ice-roughness categories - liqmask, & ! IN - Liquid-cloud mask - icemask, & ! IN - Ice-cloud mask - cld_lwp, & ! IN - Cloud liquid water path - cld_iwp, & ! IN - Cloud ice water path - cld_reliq, & ! IN - Cloud liquid effective radius - cld_reice, & ! IN - Cloud ice effective radius - sw_optical_props_cloudsByBand)) ! OUT - RRTMGP DDT: Shortwave optical properties, - ! in each band (tau,ssa,g) - else - ! RRTMG cloud-optics - if (any(cld_frac .gt. 0)) then - sw_optical_props_cloudsByBand%tau(:,:,:) = 0._kind_phys - sw_optical_props_cloudsByBand%ssa(:,:,:) = 0._kind_phys - sw_optical_props_cloudsByBand%g(:,:,:) = 0._kind_phys - call rrtmg_sw_cloud_optics(nday, model%levs, sw_gas_props%get_nband(), cld_lwp(idxday,:), & - cld_reliq(idxday,:), cld_iwp(idxday,:), cld_reice(idxday,:), cld_rwp(idxday,:), & - cld_rerain(idxday,:), cld_swp(idxday,:), cld_resnow(idxday,:), cld_frac(idxday,:),& - tau_cld, ssa_cld, asy_cld) - sw_optical_props_cloudsByBand%tau(idxday,:,:) = tau_cld - sw_optical_props_cloudsByBand%ssa(idxday,:,:) = ssa_cld - sw_optical_props_cloudsByBand%g(idxday,:,:) = asy_cld + ! ####################################################################################### + ! Change random number seed value for each radiation invocation (isubcsw =1 or 2). + ! ####################################################################################### + if(isubcsw == 1) then ! advance prescribed permutation seed + do iCol = 1, ncol + ipseed_sw(iCol) = ipsdsw0 + iCol + enddo + elseif (isubcsw == 2) then ! use input array of permutaion seeds + do iCol = 1, ncol + ipseed_sw(iCol) = icseed_sw(iCol) + enddo endif - endif - ! ####################################################################################### - ! Call McICA to generate subcolumns. - ! ####################################################################################### - ! Call RNG. Mersennse Twister accepts 1D array, so loop over columns and collapse along G-points - ! and layers. ([nGpts,model%levs,nColumn]-> [nGpts*model%levs]*nColumn) - do iCol=1,ncol - call random_setseed(ipseed_sw(icol),rng_stat) - call random_number(rng1D,rng_stat) - rng3D(:,:,iCol) = reshape(source = rng1D,shape=[sw_gas_props%get_ngpt(),model%levs]) - enddo - - ! Call McICA - select case ( iovrsw ) - ! Maximumn-random - case(1) - call check_error_msg('rrtmgp_sw_cloud_optics_run',sampled_mask_max_ran(rng3D,cld_frac,cldfracMCICA)) - end select - - ! Map band optical depth to each g-point using McICA - call check_error_msg('rrtmgp_sw_cloud_optics_run',draw_samples(cldfracMCICA,sw_optical_props_cloudsByBand,sw_optical_props_clouds)) - ! GFS_RRTMGP_POST_RUN() requires the SW optical depth ~0.55microns - cldtausw = sw_optical_props_cloudsByBand%tau(:,:,11) + ! ####################################################################################### + ! Call McICA to generate subcolumns. + ! ####################################################################################### + ! Call RNG. Mersennse Twister accepts 1D array, so loop over columns and collapse along G-points + ! and layers. ([nGpts,model%levs,nColumn]-> [nGpts*model%levs]*nColumn) + do iCol=1,ncol + call random_setseed(ipseed_sw(icol),rng_stat) + call random_number(rng1D,rng_stat) + rng3D(:,:,iCol) = reshape(source = rng1D,shape=[sw_gas_props%get_ngpt(),model%levs]) + enddo + + ! Call McICA + select case ( iovrsw ) + ! Maximumn-random + case(1) + call check_error_msg('rrtmgp_sw_cloud_optics_run',sampled_mask_max_ran(rng3D,cld_frac,cldfracMCICA)) + end select + + ! Map band optical depth to each g-point using McICA + call check_error_msg('rrtmgp_sw_cloud_optics_run',draw_samples(& + cldfracMCICA(idxday(1:nDay),:,:),sw_optical_props_cloudsByBand,sw_optical_props_clouds)) + + ! GFS_RRTMGP_POST_RUN() requires the SW optical depth ~0.55microns + cldtausw(idxday(1:nDay),:) = sw_optical_props_cloudsByBand%tau(:,:,11) + endif end subroutine rrtmgp_sw_cloud_optics_run diff --git a/physics/rrtmgp_sw_gas_optics.F90 b/physics/rrtmgp_sw_gas_optics.F90 index f235b8d02..de1bfecb1 100644 --- a/physics/rrtmgp_sw_gas_optics.F90 +++ b/physics/rrtmgp_sw_gas_optics.F90 @@ -398,8 +398,9 @@ end subroutine rrtmgp_sw_gas_optics_init !! \section arg_table_rrtmgp_sw_gas_optics_run !! \htmlinclude rrtmgp_sw_gas_optics.html !! - subroutine rrtmgp_sw_gas_optics_run(Model, Interstitial, sw_gas_props, ncol, p_lay, p_lev, & - t_lay, t_lev, gas_concentrations, lsswr, solcon, sw_optical_props_clrsky, errmsg, errflg) + subroutine rrtmgp_sw_gas_optics_run(Model, Interstitial, sw_gas_props, ncol, nday, idxday,& + p_lay, p_lev, t_lay, t_lev, gas_concentrations, lsswr, solcon, & + sw_optical_props_clrsky, errmsg, errflg) ! Inputs type(GFS_control_type), intent(in) :: & @@ -409,7 +410,10 @@ subroutine rrtmgp_sw_gas_optics_run(Model, Interstitial, sw_gas_props, ncol, p_l type(ty_gas_optics_rrtmgp),intent(in) :: & sw_gas_props ! RRTMGP DDT: spectral information for RRTMGP SW radiation scheme integer,intent(in) :: & + nday, & ! Number of daylit points. ncol ! Number of horizontal points + integer,intent(in),dimension(ncol) :: & + idxday ! Indices for daylit points. real(kind_phys), dimension(ncol,model%levs), intent(in) :: & p_lay, & ! Pressure @ model layer-centers (hPa) t_lay ! Temperature (K) @@ -432,7 +436,11 @@ subroutine rrtmgp_sw_gas_optics_run(Model, Interstitial, sw_gas_props, ncol, p_l sw_optical_props_clrsky ! RRTMGP DDT: clear-sky shortwave optical properties, spectral (tau,ssa,g) ! Local variables - integer :: ij + integer :: ij,iGas + real(kind_phys), dimension(ncol,Model%levs) :: vmrTemp + real(kind_phys), dimension(nday,sw_gas_props%get_ngpt()) :: toa_src_sw_temp + type(ty_gas_concs) :: & + gas_concentrations_daylit ! RRTMGP DDT: trace gas concentrations (vmr) ! Initialize CCPP error handling variables errmsg = '' @@ -440,23 +448,35 @@ subroutine rrtmgp_sw_gas_optics_run(Model, Interstitial, sw_gas_props, ncol, p_l if (.not. Model%lsswr) return - ! Allocate space - call check_error_msg('rrtmgp_sw_gas_optics_run',sw_optical_props_clrsky%alloc_2str(ncol, model%levs, sw_gas_props)) - - ! Gas-optics - call check_error_msg('rrtmgp_sw_gas_optics_run',sw_gas_props%gas_optics(& - p_lay, & ! IN - Pressure @ layer-centers (Pa) - p_lev, & ! IN - Pressure @ layer-interfaces (Pa) - t_lay, & ! IN - Temperature @ layer-centers (K) - gas_concentrations, & ! IN - RRTMGP DDT: trace gas volumne mixing-ratios - sw_optical_props_clrsky, & ! OUT - RRTMGP DDT: Shortwave optical properties, by - ! spectral point (tau,ssa,g) - Interstitial%toa_src_sw)) ! OUT - TOA incident shortwave radiation (spectral) - - ! Scale incident flux - do ij=1,ncol - Interstitial%toa_src_sw(ij,:) = Interstitial%toa_src_sw(ij,:)*solcon/sum(Interstitial%toa_src_sw(ij,:)) - enddo + if (nDay .gt. 0) then + ! Allocate space + call check_error_msg('rrtmgp_sw_gas_optics_run',sw_optical_props_clrsky%alloc_2str(nday, model%levs, sw_gas_props)) + + ! Subset the gas concentrations, only need daylit points. + do iGas=1,Model%nGases + call check_error_msg('rrtmgp_sw_rte_run',& + gas_concentrations%get_vmr(trim(Model%active_gases_array(iGas)),vmrTemp)) + call check_error_msg('rrtmgp_sw_rte_run',& + gas_concentrations_daylit%set_vmr(trim(Model%active_gases_array(iGas)),vmrTemp(idxday(1:nday),:))) + enddo + + ! Gas-optics + call check_error_msg('rrtmgp_sw_gas_optics_run',sw_gas_props%gas_optics(& + p_lay(idxday(1:nday),:), & ! IN - Pressure @ layer-centers (Pa) + p_lev(idxday(1:nday),:), & ! IN - Pressure @ layer-interfaces (Pa) + t_lay(idxday(1:nday),:), & ! IN - Temperature @ layer-centers (K) + gas_concentrations_daylit, & ! IN - RRTMGP DDT: trace gas volumne mixing-ratios + sw_optical_props_clrsky, & ! OUT - RRTMGP DDT: Shortwave optical properties, by + ! spectral point (tau,ssa,g) + toa_src_sw_temp)) ! OUT - TOA incident shortwave radiation (spectral) + Interstitial%toa_src_sw(idxday(1:nday),:) = toa_src_sw_temp + ! Scale incident flux + do ij=1,nday + Interstitial%toa_src_sw(idxday(ij),:) = Interstitial%toa_src_sw(idxday(ij),:)*solcon/ & + sum(Interstitial%toa_src_sw(idxday(ij),:)) + enddo + endif + end subroutine rrtmgp_sw_gas_optics_run diff --git a/physics/rrtmgp_sw_gas_optics.meta b/physics/rrtmgp_sw_gas_optics.meta index 129f278e6..1ec914d63 100644 --- a/physics/rrtmgp_sw_gas_optics.meta +++ b/physics/rrtmgp_sw_gas_optics.meta @@ -103,6 +103,22 @@ type = integer intent = in optional = F +[nday] + standard_name = daytime_points_dimension + long_name = daytime points dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[idxday] + standard_name = daytime_points + long_name = daytime points + units = index + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F [p_lay] standard_name = air_pressure_at_layer_for_RRTMGP_in_hPa long_name = air pressure layer diff --git a/physics/rrtmgp_sw_rte.F90 b/physics/rrtmgp_sw_rte.F90 index ca2c0248a..426826fb7 100644 --- a/physics/rrtmgp_sw_rte.F90 +++ b/physics/rrtmgp_sw_rte.F90 @@ -45,7 +45,7 @@ subroutine rrtmgp_sw_rte_run(Model, Interstitial, Radtend, Statein, ncol, sw_gas integer, intent(in) :: & ncol, & ! Number of horizontal gridpoints nday ! Number of daytime points - integer, intent(in), dimension(nday) :: & + integer, intent(in), dimension(ncol) :: & idxday ! Index array for daytime points real(kind_phys), dimension(ncol,Model%levs), intent(in) :: & p_lay, & ! Pressure @ model layer-centers (Pa) @@ -54,7 +54,7 @@ subroutine rrtmgp_sw_rte_run(Model, Interstitial, Radtend, Statein, ncol, sw_gas p_lev ! Pressure @ model layer-interfaces (Pa) type(ty_gas_optics_rrtmgp),intent(in) :: & sw_gas_props ! RRTMGP DDT: SW spectral information - type(ty_optical_props_2str),intent(in) :: & + type(ty_optical_props_2str),intent(inout) :: & sw_optical_props_clrsky, & ! RRTMGP DDT: longwave clear-sky radiative properties sw_optical_props_clouds, & ! RRTMGP DDT: longwave cloud radiative properties sw_optical_props_aerosol ! RRTMGP DDT: longwave aerosol radiative properties @@ -108,52 +108,38 @@ subroutine rrtmgp_sw_rte_run(Model, Interstitial, Radtend, Statein, ncol, sw_gas errflg = 0 if (.not. lsswr) return - - ! Vertical ordering? - top_at_1 = (p_lev(1,1) .lt. p_lev(1, Model%levs)) - if (top_at_1) then - iSFC = Model%levs+1 - iTOA = 1 - else - iSFC = 1 - iTOA = Model%levs+1 - endif - - ! Are any optional outputs requested? Need to know now to compute correct fluxes. - l_ClrSky_HR = present(hsw0) - l_AllSky_HR_byband = present(hswb) - l_scmpsw = present(scmpsw) - if ( l_scmpsw ) then - scmpsw = cmpfsw_type (0., 0., 0., 0., 0., 0.) - endif - fluxswUP_allsky(:,:) = 0._kind_phys - fluxswDOWN_allsky(:,:) = 0._kind_phys - fluxswUP_clrsky(:,:) = 0._kind_phys - fluxswDOWN_clrsky(:,:) = 0._kind_phys - if (nDay .gt. 0) then - ! Subset the cloud and aerosol radiative properties over daylit points. - ! Cloud optics [nDay,Model%levs,nGpts] - call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_clouds_daylit%alloc_2str(nday, Model%levs, sw_gas_props)) - sw_optical_props_clouds_daylit%tau = sw_optical_props_clouds%tau(idxday,:,:) - sw_optical_props_clouds_daylit%ssa = sw_optical_props_clouds%ssa(idxday,:,:) - sw_optical_props_clouds_daylit%g = sw_optical_props_clouds%g(idxday,:,:) - ! Aerosol optics [nDay,Model%levs,nBands] - call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_aerosol_daylit%alloc_2str(nday, Model%levs, sw_gas_props%get_band_lims_wavenumber())) - sw_optical_props_aerosol_daylit%tau = sw_optical_props_aerosol%tau(idxday,:,:) - sw_optical_props_aerosol_daylit%ssa = sw_optical_props_aerosol%ssa(idxday,:,:) - sw_optical_props_aerosol_daylit%g = sw_optical_props_aerosol%g(idxday,:,:) - ! Clear-sky optics [nDay,Model%levs,nGpts] - call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_clrsky_daylit%alloc_2str(nday, Model%levs, sw_gas_props)) - sw_optical_props_clrsky_daylit%tau = sw_optical_props_clrsky%tau(idxday,:,:) - sw_optical_props_clrsky_daylit%ssa = sw_optical_props_clrsky%ssa(idxday,:,:) - sw_optical_props_clrsky_daylit%g = sw_optical_props_clrsky%g(idxday,:,:) - - ! Similarly, subset the gas concentrations. + ! Vertical ordering? + top_at_1 = (p_lev(1,1) .lt. p_lev(1, Model%levs)) + if (top_at_1) then + iSFC = Model%levs+1 + iTOA = 1 + else + iSFC = 1 + iTOA = Model%levs+1 + endif + + ! Are any optional outputs requested? Need to know now to compute correct fluxes. + l_ClrSky_HR = present(hsw0) + l_AllSky_HR_byband = present(hswb) + l_scmpsw = present(scmpsw) + if ( l_scmpsw ) then + scmpsw = cmpfsw_type (0., 0., 0., 0., 0., 0.) + endif + + ! Initialize fluxes + fluxswUP_allsky(:,:) = 0._kind_phys + fluxswDOWN_allsky(:,:) = 0._kind_phys + fluxswUP_clrsky(:,:) = 0._kind_phys + fluxswDOWN_clrsky(:,:) = 0._kind_phys + + ! Subset the gas concentrations, only need daylit points. do iGas=1,Model%nGases - call check_error_msg('rrtmgp_sw_rte_run',gas_concentrations%get_vmr(trim(Model%active_gases_array(iGas)),vmrTemp)) - call check_error_msg('rrtmgp_sw_rte_run',gas_concentrations_daylit%set_vmr(trim(Model%active_gases_array(iGas)),vmrTemp(idxday,:))) + call check_error_msg('rrtmgp_sw_rte_run',& + gas_concentrations%get_vmr(trim(Model%active_gases_array(iGas)),vmrTemp)) + call check_error_msg('rrtmgp_sw_rte_run',& + gas_concentrations_daylit%set_vmr(trim(Model%active_gases_array(iGas)),vmrTemp(idxday(1:nday),:))) enddo ! Initialize RRTMGP DDT containing 2D(3D) fluxes @@ -165,42 +151,43 @@ subroutine rrtmgp_sw_rte_run(Model, Interstitial, Radtend, Statein, ncol, sw_gas ! Compute clear-sky fluxes (if requested) ! Clear-sky fluxes (gas+aerosol) - call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_aerosol_daylit%increment(sw_optical_props_clrsky_daylit)) + call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_aerosol%increment(sw_optical_props_clrsky)) ! Delta-scale optical properties - call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_clrsky_daylit%delta_scale()) + call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_clrsky%delta_scale()) if (l_ClrSky_HR) then - call check_error_msg('rrtmgp_sw_rte_run',rte_sw( & - sw_optical_props_clrsky_daylit, & ! IN - optical-properties - top_at_1, & ! IN - veritcal ordering flag - Radtend%coszen(idxday), & ! IN - Cosine of solar zenith angle - Interstitial%toa_src_sw(idxday,:), & ! IN - incident solar flux at TOA - Interstitial%sfc_alb_nir_dir(:,idxday), & ! IN - Shortwave surface albedo (direct) - Interstitial%sfc_alb_nir_dif(:,idxday), & ! IN - Shortwave surface albedo (diffuse) - flux_clrsky)) ! OUT - Fluxes, clear-sky, 3D (nCol,Model%levs,nBand) + call check_error_msg('rrtmgp_sw_rte_run',rte_sw( & + sw_optical_props_clrsky, & ! IN - optical-properties + top_at_1, & ! IN - veritcal ordering flag + Radtend%coszen(idxday(1:nday)), & ! IN - Cosine of solar zenith angle + Interstitial%toa_src_sw(idxday(1:nday),:), & ! IN - incident solar flux at TOA + Interstitial%sfc_alb_nir_dir(:,idxday(1:nday)), & ! IN - Shortwave surface albedo (direct) + Interstitial%sfc_alb_nir_dif(:,idxday(1:nday)), & ! IN - Shortwave surface albedo (diffuse) + flux_clrsky)) ! OUT - Fluxes, clear-sky, 3D (nCol,Model%levs,nBand) ! Store fluxes - fluxswUP_clrsky(idxday,:) = sum(flux_clrsky%bnd_flux_up,dim=3) - fluxswDOWN_clrsky(idxday,:) = sum(flux_clrsky%bnd_flux_dn,dim=3) + fluxswUP_clrsky(idxday(1:nday),:) = sum(flux_clrsky%bnd_flux_up,dim=3) + fluxswDOWN_clrsky(idxday(1:nday),:) = sum(flux_clrsky%bnd_flux_dn,dim=3) endif ! Compute all-sky fluxes ! All-sky fluxes (clear-sky + clouds) - call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_clouds_daylit%increment(sw_optical_props_clrsky_daylit)) + call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_clouds%increment(sw_optical_props_clrsky)) ! Delta-scale optical properties - call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_clouds_daylit%delta_scale()) - call check_error_msg('rrtmgp_sw_rte_run',rte_sw( & - sw_optical_props_clrsky_daylit, & ! IN - optical-properties - top_at_1, & ! IN - veritcal ordering flag - Radtend%coszen(idxday), & ! IN - Cosine of solar zenith angle - Interstitial%toa_src_sw(idxday,:), & ! IN - incident solar flux at TOA - Interstitial%sfc_alb_nir_dir(:,idxday), & ! IN - Shortwave surface albedo (direct) - Interstitial%sfc_alb_nir_dif(:,idxday), & ! IN - Shortwave surface albedo (diffuse) - flux_allsky)) ! OUT - Fluxes, clear-sky, 3D (nCol,Model%levs,nBand) + call check_error_msg('rrtmgp_sw_rte_run',sw_optical_props_clouds%delta_scale()) + call check_error_msg('rrtmgp_sw_rte_run',rte_sw( & + sw_optical_props_clrsky, & ! IN - optical-properties + top_at_1, & ! IN - veritcal ordering flag + Radtend%coszen(idxday(1:nday)), & ! IN - Cosine of solar zenith angle + Interstitial%toa_src_sw(idxday(1:nday),:), & ! IN - incident solar flux at TOA + Interstitial%sfc_alb_nir_dir(:,idxday(1:nday)), & ! IN - Shortwave surface albedo (direct) + Interstitial%sfc_alb_nir_dif(:,idxday(1:nday)), & ! IN - Shortwave surface albedo (diffuse) + flux_allsky)) ! OUT - Fluxes, clear-sky, 3D (nCol,Model%levs,nBand) ! Store fluxes - fluxswUP_allsky(idxday,:) = sum(flux_allsky%bnd_flux_up,dim=3) - fluxswDOWN_allsky(idxday,:) = sum(flux_allsky%bnd_flux_dn,dim=3) + fluxswUP_allsky(idxday(1:nday),:) = sum(flux_allsky%bnd_flux_up,dim=3) + fluxswDOWN_allsky(idxday(1:nday),:) = sum(flux_allsky%bnd_flux_dn,dim=3) if ( l_scmpsw ) then - scmpsw(idxday)%nirbm = sum(flux_allsky%bnd_flux_dn_dir(idxday,iSFC,:),dim=2) - scmpsw(idxday)%nirdf = sum(flux_allsky%bnd_flux_dn(idxday,iSFC,:),dim=2) - sum(flux_allsky%bnd_flux_dn_dir(idxday,iSFC,:),dim=2) + scmpsw(idxday(1:nday))%nirbm = sum(flux_allsky%bnd_flux_dn_dir(idxday(1:nday),iSFC,:),dim=2) + scmpsw(idxday(1:nday))%nirdf = sum(flux_allsky%bnd_flux_dn(idxday(1:nday),iSFC,:),dim=2) - & + sum(flux_allsky%bnd_flux_dn_dir(idxday(1:nday),iSFC,:),dim=2) endif endif end subroutine rrtmgp_sw_rte_run diff --git a/physics/rrtmgp_sw_rte.meta b/physics/rrtmgp_sw_rte.meta index 197824495..b677f15a6 100644 --- a/physics/rrtmgp_sw_rte.meta +++ b/physics/rrtmgp_sw_rte.meta @@ -82,7 +82,7 @@ units = DDT dimensions = () type = ty_optical_props_2str - intent = in + intent = inout optional = F [sw_optical_props_clouds] standard_name = shortwave_optical_properties_for_cloudy_atmosphere