search

ESCOMP / CAM

Community Atmosphere Model
71 stars 133 forks source link

Turning on spectralflux namelist option doesn't return spectral fluxes #1032

Open GilbertCloud opened 1 month ago

GilbertCloud commented 1 month ago

What happened?

I was running SCAM on my personal computer (cesm2.1.3, cam_cesm2_1_rel41). I found the bug when I set the spectralflux namelist option to true. This should have returned the spectral longwave and shortwave fluxes in addition to broadband fluxes and it didn't. When I added the corresponding fields to the fincl1 namelist option (LD, LU, SD, SU), SCAM stopped with and error because it couldn't find those fields.

What are the steps to reproduce the bug?

To reproduce this bug, set 'spectralflux' to true and include the fields LD, LU, SD, SU in the fincl1 namelist option.

What CAM tag were you using?

cam_cesm2_1_rel41

What machine were you running CAM on?

Personal Computer

What compiler were you using?

Other (please specify below)

Path to a case directory, if applicable

No response

Will you be addressing this bug yourself?

No

Extra info

Notes

I am running CAM through SCAM on my personal computer and I am not sure what compiler I am using. I have checked the most recent tag of CAM (cam6_3_160) and this bug is still in the CAM code. The buggy files radconstants.F90 and radiation.F90 are located at /cam/src/physics/rrtmg.

Fixes

I have produced code that fixes the bug. The files are at /glade/u/home/glydia/mods/fixedspectral/.

In the radconstants file, I added the longwave and shortwave band midpoints as variables, as well as functions to retrieve them, similar to existing functions in the file to retrieve the band boundaries. I included the midpoints so that the history file has them and the boundaries as dimensions for the spectral fluxes.

The additions I made to the radconstants file are below:

real(r8),parameter :: wavenum_mid(nbndsw) = & ! in cm^-1
   (/2906.89_r8, 3605.55_r8, 4312.77_r8, 4893.62_r8, 5627.83_r8, 6881.50_r8, 7873.06_r8, &
   10170.67_r8,14338.76_r8,19036.81_r8,25629.08_r8,33196.39_r8,43588.99_r8, 1460.14_r8/)
real(r8), parameter :: wavenumber3_longwave(nlwbands) = &! Longwave spectral band mid-points (cm-1)
    (/  59.16_r8,  418.33_r8,  561.25_r8,  664.08_r8,  757.63_r8,  896.44_r8, 1028.79_r8, 1128.89_r8, &
       1280.70_r8, 1434.29_r8, 1632.18_r8, 1934.94_r8, 2163.33_r8, 2318.94_r8, 2492.79_r8, 2906.89_r8 /)
public :: get_number_sw_bands, &
          get_sw_spectral_boundaries, &
          get_sw_spectral_midpts, &
          get_lw_spectral_boundaries, &
          get_lw_spectral_midpts, &
          get_ref_solar_band_irrad, &
          get_ref_total_solar_irrad, &
          get_solar_band_fraction_irrad
subroutine get_lw_spectral_midpts(midpts, units)
   ! provide spectral midpts of each longwave band

   real(r8), intent(out) :: midpts(nlwbands)
   character(*), intent(in) :: units ! requested units

   select case (units)
   case ('inv_cm','cm^-1','cm-1')
      midpts  = wavenumber3_longwave
   case('m','meter','meters')
      midpts  = 1.e-2_r8/wavenumber3_longwave
   case('nm','nanometer','nanometers')
      midpts  = 1.e7_r8/wavenumber3_longwave
   case('um','micrometer','micrometers','micron','microns')
      midpts  = 1.e4_r8/wavenumber3_longwave
   case('cm','centimeter','centimeters')
      midpts  = 1._r8/wavenumber3_longwave
   case default
      call endrun('get_lw_spectral_midpts: spectral units not acceptable'//units)
   end select

end subroutine get_lw_spectral_midpts
subroutine get_sw_spectral_midpts(midpts, units)
   ! provide spectral midpts of each shortwave band

   real(r8), intent(out) :: midpts(nswbands)
   character(*), intent(in) :: units ! requested units

   select case (units)
   case ('inv_cm','cm^-1','cm-1')
      midpts = wavenum_mid
   case('m','meter','meters')
      midpts = 1.e-2_r8/wavenum_mid
   case('nm','nanometer','nanometers')
      midpts = 1.e7_r8/wavenum_mid
   case('um','micrometer','micrometers','micron','microns')
      midpts = 1.e4_r8/wavenum_mid
   case('cm','centimeter','centimeters')
      midpts  = 1._r8/wavenum_mid
   case default
      call endrun('rad_constants.F90: spectral units not acceptable'//units)
   end select

end subroutine get_sw_spectral_midpts

In the radiation file, I added the longwave and shortwave boundaries and midpoints as variables and added the LD, LU, SD, SU fields to the history master field list. I added the longwave and shortwave bands to the history file coordinates and added calls to the pbuf in the longwave and shortwave output functions to output the spectral fluxes to the history file.

The additions/changes I made to the radiation file are below:

use radconstants,        only: nswbands, nlwbands, rrtmg_sw_cloudsim_band, rrtmg_lw_cloudsim_band, &
                               idx_sw_diag, get_lw_spectral_boundaries, get_sw_spectral_boundaries, &
                               get_lw_spectral_midpts, get_sw_spectral_midpts
use cam_history_support, only: fillvalue, add_hist_coord
! LW/SW band variables
real(r8) :: lw_bounds(2,nlwbands)
real(r8) :: sw_bounds(2,nswbands)
real(r8) :: lw_vals(nlwbands)
real(r8) :: sw_vals(nswbands)
subroutine radiation_register
   use cam_history_support, only: add_hist_coord
   ...
   if (spectralflux) then
      ! Add history coordinates for lw/sw bands
      call get_sw_spectral_midpts(sw_vals,'cm^-1')
      call get_lw_spectral_midpts(lw_vals,'cm^-1')
      call get_sw_spectral_boundaries(sw_bounds(1,:),sw_bounds(2,:),'cm^-1')
      call get_lw_spectral_boundaries(lw_bounds(1,:),lw_bounds(2,:),'cm^-1')

      call add_hist_coord('sw_band', nswbands, 'Shortwave bands',  &
               'cm^-1', sw_vals, bounds_name='sw_band_bds', bounds=sw_bounds)
      call add_hist_coord('lw_band', nlwbands, 'Longwave bands',  &
               'cm^-1', lw_vals, bounds_name='lw_band_bds', bounds=lw_bounds)
subroutine radiation_init(pbuf2d)
  ...
   ! Add shortwave radiation fields to history master field list.
   do icall = 0, N_DIAG

      if (active_calls(icall)) then
         ...
         call addfld('SU'//diag(icall),      (/ character(len=8) :: 'ilev','sw_band' /),   'A', 'W/m2', 'Shortwave spectral flux up',        &
                                                                                 sampling_seq='rad_lwsw')
         call addfld('SD'//diag(icall),      (/ character(len=8) :: 'ilev','sw_band' /),   'A', 'W/m2', 'Shortwave spectral flux down',   &
                                                                                 sampling_seq='rad_lwsw')
   ...
   ! Add longwave radiation fields to history master field list.

   do icall = 0, N_DIAG

      if (active_calls(icall)) then
         ...
         call addfld('LU'//diag(icall),      (/  character(len=8) :: 'ilev','lw_band' /),   'A', 'W/m2', 'Longwave spectral flux up',         &
                                                                           sampling_seq='rad_lwsw')      
         call addfld('LD'//diag(icall),      (/  character(len=8) :: 'ilev','lw_band' /),   'A', 'W/m2', 'Longwave spectral flux down',    &
                                                                           sampling_seq='rad_lwsw') 
subroutine radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out)
   real(r8), pointer, dimension(:,:,:) :: su
   real(r8), pointer, dimension(:,:,:) :: sd

   call pbuf_get_field(pbuf, su_idx, su)
   call pbuf_get_field(pbuf, sd_idx, sd)

   call outfld('SU'//diag(icall),       su,         ncol, lchnk)
   call outfld('SD'//diag(icall),       sd,         ncol, lchnk)
subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr)
   real(r8), pointer, dimension(:,:,:) :: lu
   real(r8), pointer, dimension(:,:,:) :: ld

   call pbuf_get_field(pbuf, lu_idx, lu)
   call pbuf_get_field(pbuf, ld_idx, ld)

   call outfld('LU'//diag(icall),       lu,         ncol, lchnk)
   call outfld('LD'//diag(icall),       ld,         ncol, lchnk)
brian-eaton commented 1 month ago

Hi @GilbertCloud. The spectralflux option was originally added to CAM in order to make the spectrally resolved fluxes available to the CARMA aerosol code. That is why those fluxes are added to the physics buffer. But the modifications for CARMA never included the addfld/outfld calls necessary to output these fields to the history files. I'm not sure why. It seems reasonable to do so. @brianpm, what do you think?

I have looked through your code mods to output the spectrally resolved fluxes and they look well done. Thanks.

brianpm commented 4 weeks ago

I think this is a good idea. There are probably lots of uses for spectrally-resolved output.