Skip to content

Commit

Permalink
add diagnostic output for fluxes on the RRTMGP grid
Browse files Browse the repository at this point in the history
  • Loading branch information
@brian-eaton
brian-eaton committed Aug 14, 2023
1 parent 8d31535 commit d46aa4b
Showing 1 changed file with 85 additions and 26 deletions.
111 changes: 85 additions & 26 deletions src/physics/rrtmgp/radiation.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ module radiation
use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp

use cam_history, only: addfld, add_default, horiz_only, outfld, hist_fld_active
use cam_history_support, only: fillvalue
use cam_history_support, only: fillvalue, add_vert_coord

use ioFileMod, only: getfil
use cam_pio_utils, only: cam_pio_openfile
Expand Down Expand Up @@ -126,11 +126,21 @@ module radiation
real(r8) :: flux_sw_dn(pcols,pverp) ! downward flux
real(r8) :: flux_sw_clr_dn(pcols,pverp) ! downward clearsky flux

real(r8), allocatable :: fsdn(:,:) ! Downward SW flux on rrtmgp grid
real(r8), allocatable :: fsdnc(:,:) ! Downward SW clear sky flux on rrtmgp grid
real(r8), allocatable :: fsup(:,:) ! Upward SW flux on rrtmgp grid
real(r8), allocatable :: fsupc(:,:) ! Upward SW clear sky flux on rrtmgp grid

real(r8) :: flux_lw_up(pcols,pverp) ! upward shortwave flux on interfaces
real(r8) :: flux_lw_clr_up(pcols,pverp) ! upward shortwave clearsky flux
real(r8) :: flux_lw_dn(pcols,pverp) ! downward flux
real(r8) :: flux_lw_clr_dn(pcols,pverp) ! downward clearsky flux

real(r8), allocatable :: fldn(:,:) ! Downward LW flux on rrtmgp grid
real(r8), allocatable :: fldnc(:,:) ! Downward LW clear sky flux on rrtmgp grid
real(r8), allocatable :: flup(:,:) ! Upward LW flux on rrtmgp grid
real(r8), allocatable :: flupc(:,:) ! Upward LW clear sky flux on rrtmgp grid

real(r8) :: qrlc(pcols,pver)

real(r8) :: flntc(pcols) ! Clear sky lw flux at model top
Expand Down Expand Up @@ -209,6 +219,9 @@ module radiation
integer :: ktopradm ! index in RRTMGP arrays of layer corresponding to CAM top layer
integer :: ktopradi ! index in RRTMGP arrays of interface corresponding to CAM top interface

! vertical coordinate for output of fluxes on radiation grid
real(r8), allocatable, target :: plev_rad(:)

! LW coefficients
type(ty_gas_optics_rrtmgp) :: kdist_lw ! bpm changed here
integer :: ngpt_lw
Expand Down Expand Up @@ -349,11 +362,11 @@ end subroutine radiation_readnl

subroutine radiation_register

! Register radiation fields in the physics buffer

use physics_buffer, only: pbuf_add_field, dtype_r8
use radiation_data, only: rad_data_register

! Register radiation fields in the physics buffer

call pbuf_add_field('QRS' , 'global',dtype_r8,(/pcols,pver/), qrs_idx) ! shortwave radiative heating rate
call pbuf_add_field('QRL' , 'global',dtype_r8,(/pcols,pver/), qrl_idx) ! longwave radiative heating rate

Expand Down Expand Up @@ -502,6 +515,7 @@ subroutine radiation_init(pbuf2d)
! below 1 Pa then an extra layer is added to the top of the model for
! the purpose of the radiation calculation.
nlay = count( pref_edge(:) > 1._r8 ) ! pascals (0.01 mbar)
allocate(plev_rad(nlay+1))

if (nlay == pverp) then
! Top model interface is below 1 Pa. RRTMGP is active in all model layers plus
Expand All @@ -510,14 +524,21 @@ subroutine radiation_init(pbuf2d)
ktopcami = 1
ktopradm = 2
ktopradi = 2
plev_rad(1) = 1.01_r8 ! Top of extra layer, Pa.
plev_rad(2:) = pref_edge
else
! nlay < pverp. nlay layers are set by radiation
ktopcamm = pverp - nlay + 1
ktopcami = pverp - nlay + 1
ktopradm = 1
ktopradi = 1
plev_rad = pref_edge(ktopcami:)
end if

! Define a pressure coordinate to allow output of data on the radiation grid.
call add_vert_coord('plev_rad', nlay+1, 'Pressures at radiation flux calculations', &
'Pa', plev_rad)

call set_available_gases(active_gases, available_gases) ! gases needed to initialize spectral info

call coefs_init(coefs_lw_file, kdist_lw, available_gases, band2gpt_lw)
Expand Down Expand Up @@ -690,6 +711,12 @@ subroutine radiation_init(pbuf2d)
call addfld('FUSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave clear-sky upward flux')
call addfld('FDSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave clear-sky downward flux')

! Fluxes on rrtmgp grid
call addfld('FSDN'//diag(icall), (/ 'plev_rad' /), 'I', 'W/m2', 'SW downward flux on rrtmgp grid')
call addfld('FSDNC'//diag(icall), (/ 'plev_rad' /), 'I', 'W/m2', 'SW downward clear sky flux on rrtmgp grid')
call addfld('FSUP'//diag(icall), (/ 'plev_rad' /), 'I', 'W/m2', 'SW upward flux on rrtmgp grid')
call addfld('FSUPC'//diag(icall), (/ 'plev_rad' /), 'I', 'W/m2', 'SW upward clear sky flux on rrtmgp grid')

if (history_amwg) then
call add_default('SOLIN'//diag(icall), 1, ' ')
call add_default('QRS'//diag(icall), 1, ' ')
Expand Down Expand Up @@ -746,6 +773,12 @@ subroutine radiation_init(pbuf2d)
call addfld('FULC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Longwave clear-sky upward flux')
call addfld('FDLC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Longwave clear-sky downward flux')

! Fluxes on rrtmgp grid
call addfld('FLDN'//diag(icall), (/ 'plev_rad' /), 'I', 'W/m2', 'LW downward flux on rrtmgp grid')
call addfld('FLDNC'//diag(icall), (/ 'plev_rad' /), 'I', 'W/m2', 'LW downward clear sky flux on rrtmgp grid')
call addfld('FLUP'//diag(icall), (/ 'plev_rad' /), 'I', 'W/m2', 'LW upward flux on rrtmgp grid')
call addfld('FLUPC'//diag(icall), (/ 'plev_rad' /), 'I', 'W/m2', 'LW upward clear sky flux on rrtmgp grid')

if (history_amwg) then
call add_default('QRL'//diag(icall), 1, ' ')
call add_default('FLNT'//diag(icall), 1, ' ')
Expand Down Expand Up @@ -1095,6 +1128,11 @@ subroutine radiation_tend( &
write_output = .false.
else
allocate(rd)
! allocate some elements of rd
if (.not. allocated(rd%fsdn)) then
allocate(rd%fsdn(pcols,nlay+1), rd%fsdnc(pcols,nlay+1), rd%fsup(pcols,nlay+1), rd%fsupc(pcols,nlay+1), &
rd%fldn(pcols,nlay+1), rd%fldnc(pcols,nlay+1), rd%flup(pcols,nlay+1), rd%flupc(pcols,nlay+1) )
end if
write_output = .true.
end if

Expand Down Expand Up @@ -1831,34 +1869,39 @@ subroutine set_sw_diags()
size(fsw%bnd_flux_dn,2), &
size(fsw%bnd_flux_dn,3)) :: flux_dn_diffuse
!-------------------------------------------------------------------------
fns = 0._r8 ! net sw flux
fcns = 0._r8 ! net sw clearsky flux
fsds = 0._r8 ! downward sw flux at surface
rd%fsdsc = 0._r8 ! downward sw clearsky flux at surface
rd%fsutoa = 0._r8 ! upward sw flux at TOA
rd%fsntoa = 0._r8 ! net sw at TOA
rd%fsntoac = 0._r8 ! net sw clearsky flux at TOA
rd%solin = 0._r8 ! solar irradiance at TOA

! Initializing these arrays to 0.0 provides fill in the night columns:
fns = 0._r8 ! net sw flux
fcns = 0._r8 ! net sw clearsky flux
fsds = 0._r8 ! downward sw flux at surface
rd%fsdsc = 0._r8 ! downward sw clearsky flux at surface
rd%fsutoa = 0._r8 ! upward sw flux at TOA
rd%fsntoa = 0._r8 ! net sw at TOA
rd%fsntoac = 0._r8 ! net sw clearsky flux at TOA
rd%solin = 0._r8 ! solar irradiance at TOA
rd%fsdn = 0._r8
rd%fsdnc = 0._r8
rd%fsup = 0._r8
rd%fsupc = 0._r8

! fns, fcns, rd are on CAM grid (do not have "extra layer" when it is present.)
! fill in the daylit columns:
do i = 1, nday
fns(idxday(i),ktopcami:) = fsw%flux_net(i, ktopradi:)
fcns(idxday(i),ktopcami:) = fswc%flux_net(i,ktopradi:)
rd%flux_sw_up(idxday(i),ktopcami:) = &
fsw%flux_up(i,ktopradi:)
rd%flux_sw_dn(idxday(i),ktopcami:) = &
fsw%flux_dn(i,ktopradi:)
rd%flux_sw_clr_up(idxday(i),ktopcami:) = &
fswc%flux_up(i,ktopradi:)
rd%flux_sw_clr_dn(idxday(i),ktopcami:) = &
fswc%flux_dn(i,ktopradi:)
fsds(idxday(i)) = fsw%flux_dn(i, nlay+1)
rd%fsdsc(idxday(i)) = fswc%flux_dn(i, nlay+1)
rd%fsutoa(idxday(i)) = fsw%flux_up(i, 1)
rd%fsntoa(idxday(i)) = fsw%flux_net(i, 1) ! net sw flux at TOA (*NOT* the same as fsnt)
rd%fsntoac(idxday(i)) = fswc%flux_net(i, 1) ! net sw clearsky flux at TOA (*NOT* the same as fsntc)
rd%solin(idxday(i)) = fswc%flux_dn(i, 1)
fsds(idxday(i)) = fsw%flux_dn(i, nlay+1)
rd%fsdsc(idxday(i)) = fswc%flux_dn(i, nlay+1)
rd%fsutoa(idxday(i)) = fsw%flux_up(i, 1)
rd%fsntoa(idxday(i)) = fsw%flux_net(i, 1) ! net sw flux at TOA (*NOT* the same as fsnt)
rd%fsntoac(idxday(i)) = fswc%flux_net(i, 1) ! net sw clearsky flux at TOA (*NOT* the same as fsntc)
rd%solin(idxday(i)) = fswc%flux_dn(i, 1)
rd%flux_sw_up(idxday(i),ktopcami:) = fsw%flux_up(i,ktopradi:)
rd%flux_sw_dn(idxday(i),ktopcami:) = fsw%flux_dn(i,ktopradi:)
rd%flux_sw_clr_up(idxday(i),ktopcami:) = fswc%flux_up(i,ktopradi:)
rd%flux_sw_clr_dn(idxday(i),ktopcami:) = fswc%flux_dn(i,ktopradi:)
rd%fsdn(idxday(i),:) = fsw%flux_dn(i,:)
rd%fsdnc(idxday(i),:) = fswc%flux_dn(i,:)
rd%fsup(idxday(i),:) = fsw%flux_up(i,:)
rd%fsupc(idxday(i),:) = fswc%flux_up(i,:)
end do

call heating_rate('SW', ncol, fns, qrs)
Expand Down Expand Up @@ -1962,6 +2005,11 @@ subroutine set_lw_diags()
rd%flut(:ncol) = flw%flux_up(:, ktopradi)
rd%flutc(:ncol) = flwc%flux_up(:, ktopradi)

rd%fldn(:ncol,:) = flw%flux_dn
rd%fldnc(:ncol,:) = flwc%flux_dn
rd%flup(:ncol,:) = flw%flux_up
rd%flupc(:ncol,:) = flwc%flux_up

! Output fluxes at 200 mb
call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, rd%fln200)
call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c)
Expand Down Expand Up @@ -2090,6 +2138,11 @@ subroutine radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
call outfld('FDS'//diag(icall), rd%flux_sw_dn, pcols, lchnk)
call outfld('FDSC'//diag(icall), rd%flux_sw_clr_dn, pcols, lchnk)

call outfld('FSDN'//diag(icall), rd%fsdn, pcols, lchnk)
call outfld('FSDNC'//diag(icall), rd%fsdnc, pcols, lchnk)
call outfld('FSUP'//diag(icall), rd%fsup, pcols, lchnk)
call outfld('FSUPC'//diag(icall), rd%fsupc, pcols, lchnk)

end subroutine radiation_output_sw


Expand Down Expand Up @@ -2169,6 +2222,12 @@ subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out)
call outfld('FDLC'//diag(icall), rd%flux_lw_clr_dn, pcols, lchnk)
call outfld('FUL'//diag(icall), rd%flux_lw_up, pcols, lchnk)
call outfld('FULC'//diag(icall), rd%flux_lw_clr_up, pcols, lchnk)

call outfld('FLDN'//diag(icall), rd%fldn, pcols, lchnk)
call outfld('FLDNC'//diag(icall), rd%fldnc, pcols, lchnk)
call outfld('FLUP'//diag(icall), rd%flup, pcols, lchnk)
call outfld('FLUPC'//diag(icall), rd%flupc, pcols, lchnk)

end subroutine radiation_output_lw

!===============================================================================
Expand Down

0 comments on commit d46aa4b

Please sign in to comment.