From f57e02726a0736902d2de7aeb6489f656a29effe Mon Sep 17 00:00:00 2001 From: Brian Eaton Date: Fri, 1 Sep 2023 19:58:06 -0400 Subject: [PATCH] cleanup solar forcing and gas_concs init; bugfix in rrtmgp_inputs --- src/physics/cam/aer_rad_props.F90 | 8 +- src/physics/cam/phys_prop.F90 | 11 +- src/physics/cam/rad_constituents.F90 | 4 +- src/physics/camrt/radconstants.F90 | 17 -- src/physics/rrtmg/cloud_rad_props.F90 | 2 +- src/physics/rrtmg/ebert_curry.F90 | 2 +- src/physics/rrtmg/oldcloud.F90 | 2 +- src/physics/rrtmg/radconstants.F90 | 16 -- src/physics/rrtmg/slingo.F90 | 2 +- src/physics/rrtmgp/cloud_rad_props.F90 | 2 +- src/physics/rrtmgp/ebert_curry.F90 | 15 +- src/physics/rrtmgp/oldcloud.F90 | 8 +- src/physics/rrtmgp/rad_solar_var.F90 | 145 ---------- src/physics/rrtmgp/radconstants.F90 | 350 ++++++++----------------- src/physics/rrtmgp/radiation.F90 | 249 +++++++----------- src/physics/rrtmgp/rrtmgp_inputs.F90 | 150 ++++------- src/physics/rrtmgp/slingo.F90 | 22 +- src/physics/simple/radconstants.F90 | 2 - 18 files changed, 262 insertions(+), 745 deletions(-) delete mode 100644 src/physics/rrtmgp/rad_solar_var.F90 diff --git a/src/physics/cam/aer_rad_props.F90 b/src/physics/cam/aer_rad_props.F90 index 058f53f784..be8f0708a6 100644 --- a/src/physics/cam/aer_rad_props.F90 +++ b/src/physics/cam/aer_rad_props.F90 @@ -11,7 +11,8 @@ module aer_rad_props use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc -use radconstants, only: nrh, nswbands, nlwbands, idx_sw_diag, ot_length +use radconstants, only: nswbands, nlwbands, idx_sw_diag +use phys_prop, only: nrh, ot_length use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, & rad_cnst_get_aer_props use wv_saturation, only: qsat @@ -304,9 +305,6 @@ end subroutine aer_rad_props_sw subroutine aer_rad_props_lw(list_idx, state, pbuf, odap_aer) - use radconstants, only: ot_length - - use physics_buffer, only : pbuf_get_field, pbuf_get_index, physics_buffer_desc ! Purpose: Compute aerosol transmissions needed in absorptivity/ ! emissivity calculations @@ -314,6 +312,8 @@ subroutine aer_rad_props_lw(list_idx, state, pbuf, odap_aer) ! species. If this changes, this routine will need to do something ! similar to the sw with routines like get_hygro_lw_abs + use physics_buffer, only : pbuf_get_field, pbuf_get_index, physics_buffer_desc + ! Arguments integer, intent(in) :: list_idx ! index of the climate or a diagnostic list type(physics_state), intent(in), target :: state diff --git a/src/physics/cam/phys_prop.F90 b/src/physics/cam/phys_prop.F90 index 568427e44e..ecbf6f85e0 100644 --- a/src/physics/cam/phys_prop.F90 +++ b/src/physics/cam/phys_prop.F90 @@ -11,7 +11,7 @@ module phys_prop use shr_kind_mod, only: r8 => shr_kind_r8 use spmd_utils, only: masterproc -use radconstants, only: nrh, nlwbands, nswbands, idx_sw_diag +use radconstants, only: nlwbands, nswbands, idx_sw_diag use ioFileMod, only: getfil use cam_pio_utils, only: cam_pio_openfile use pio, only: file_desc_t, var_desc_t, pio_get_var, pio_inq_varid, & @@ -26,6 +26,7 @@ module phys_prop save integer, parameter, public :: ot_length = 32 + public :: & physprop_accum_unique_files, &! Make a list of the unique set of files that contain properties ! This is an initialization step that must be done before calling physprop_init @@ -105,6 +106,10 @@ module phys_prop ! array. character(len=256), allocatable :: uniquefilenames(:) +! Number of evenly spaced intervals in rh used in this module and in the aer_rad_props module +! for calculations of aerosol hygroscopic growth. +integer, parameter, public :: nrh = 1000 + !================================================================================================ contains !================================================================================================ @@ -1106,6 +1111,8 @@ subroutine bulk_props_init(physprop, nc_id) type(var_desc_T) :: vid + ! ***N.B.*** RRTMGP hasn't set the value of idx_sw_diag when this routine is + ! called. The debug option will need to be modified for RRTMGP. logical :: debug = .true. character(len=*), parameter :: subname = 'bulk_props_init' @@ -1134,7 +1141,7 @@ subroutine bulk_props_init(physprop, nc_id) ierr = pio_get_var(nc_id, vid, physprop%num_to_mass_aer) ! Output select data to log file - if (debug .and. masterproc) then + if (debug .and. masterproc .and. idx_sw_diag > 0) then if (trim(physprop%aername) == 'SULFATE') then write(iulog, '(2x, a)') '_______ hygroscopic growth in visible band _______' call aer_optics_log_rh('SO4', physprop%sw_hygro_ext(:,idx_sw_diag), & diff --git a/src/physics/cam/rad_constituents.F90 b/src/physics/cam/rad_constituents.F90 index ced2c35cfa..42c978cc72 100644 --- a/src/physics/cam/rad_constituents.F90 +++ b/src/physics/cam/rad_constituents.F90 @@ -17,9 +17,9 @@ module rad_constituents use physics_types, only: physics_state use phys_control, only: use_simple_phys use constituents, only: cnst_get_ind -use radconstants, only: nradgas, rad_gas_index, ot_length +use radconstants, only: nradgas, rad_gas_index use phys_prop, only: physprop_accum_unique_files, physprop_init, & - physprop_get_id + physprop_get_id, ot_length use cam_history, only: addfld, fieldname_len, outfld, horiz_only use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_get_index diff --git a/src/physics/camrt/radconstants.F90 b/src/physics/camrt/radconstants.F90 index 89503fd0f5..c95c8d2154 100644 --- a/src/physics/camrt/radconstants.F90 +++ b/src/physics/camrt/radconstants.F90 @@ -21,9 +21,6 @@ module radconstants public :: radconstants_init public :: rad_gas_index -! optics files specify a type. What length is it? -integer, parameter, public :: ot_length = 32 - ! SHORTWAVE DATA ! number of shorwave spectral intervals @@ -40,20 +37,6 @@ module radconstants integer, parameter, public :: idx_lw_diag = 2 ! index to (H20 window) LW band - -! Number of evenly spaced intervals in rh -! The globality of this mesh may not be necessary -! Perhaps it could be specific to the aerosol -! But it is difficult to see how refined it must be -! for lookup. This value was found to be sufficient -! for Sulfate and probably necessary to resolve the -! high variation near rh = 1. Alternative methods -! were found to be too slow. -! Optimal approach would be for cam to specify size of aerosol -! based on each aerosol's characteristics. Radiation -! should know nothing about hygroscopic growth! -integer, parameter, public :: nrh = 1000 - ! LONGWAVE DATA ! number of lw bands diff --git a/src/physics/rrtmg/cloud_rad_props.F90 b/src/physics/rrtmg/cloud_rad_props.F90 index 2911e0ac21..c629c38e4b 100644 --- a/src/physics/rrtmg/cloud_rad_props.F90 +++ b/src/physics/rrtmg/cloud_rad_props.F90 @@ -7,7 +7,7 @@ module cloud_rad_props use ppgrid, only: pcols, pver, pverp use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx -use radconstants, only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag +use radconstants, only: nswbands, nlwbands, idx_sw_diag, idx_lw_diag use cam_abortutils, only: endrun use rad_constituents, only: iceopticsfile, liqopticsfile use oldcloud, only: oldcloud_lw, old_liq_get_rad_props_lw, old_ice_get_rad_props_lw, oldcloud_init diff --git a/src/physics/rrtmg/ebert_curry.F90 b/src/physics/rrtmg/ebert_curry.F90 index a1e1c031b1..7bca4ce257 100644 --- a/src/physics/rrtmg/ebert_curry.F90 +++ b/src/physics/rrtmg/ebert_curry.F90 @@ -7,7 +7,7 @@ module ebert_curry use ppgrid, only: pcols, pver, pverp use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx -use radconstants, only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag, get_sw_spectral_boundaries +use radconstants, only: nswbands, nlwbands, idx_sw_diag, idx_lw_diag, get_sw_spectral_boundaries use cam_abortutils, only: endrun use cam_history, only: outfld diff --git a/src/physics/rrtmg/oldcloud.F90 b/src/physics/rrtmg/oldcloud.F90 index 609c6b4668..fb0ae4d80e 100644 --- a/src/physics/rrtmg/oldcloud.F90 +++ b/src/physics/rrtmg/oldcloud.F90 @@ -7,7 +7,7 @@ module oldcloud use ppgrid, only: pcols, pver, pverp use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field -use radconstants, only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag, get_sw_spectral_boundaries +use radconstants, only: nswbands, nlwbands, idx_sw_diag, idx_lw_diag, get_sw_spectral_boundaries use cam_abortutils, only: endrun use cam_history, only: outfld use rad_constituents, only: iceopticsfile, liqopticsfile diff --git a/src/physics/rrtmg/radconstants.F90 b/src/physics/rrtmg/radconstants.F90 index f4f8c76b9c..601bcd3cf6 100644 --- a/src/physics/rrtmg/radconstants.F90 +++ b/src/physics/rrtmg/radconstants.F90 @@ -63,19 +63,6 @@ module radconstants integer, parameter, public :: rrtmg_sw_cloudsim_band = 9 ! rrtmg band for .67 micron -! Number of evenly spaced intervals in rh -! The globality of this mesh may not be necessary -! Perhaps it could be specific to the aerosol -! But it is difficult to see how refined it must be -! for lookup. This value was found to be sufficient -! for Sulfate and probably necessary to resolve the -! high variation near rh = 1. Alternative methods -! were found to be too slow. -! Optimal approach would be for cam to specify size of aerosol -! based on each aerosol's characteristics. Radiation -! should know nothing about hygroscopic growth! -integer, parameter, public :: nrh = 1000 - ! LONGWAVE DATA ! These are indices to the band for diagnostic output @@ -123,9 +110,6 @@ module radconstants real(r8), public, parameter :: minmmr(nradgas) & = epsilon(1._r8) -! Length of "optics type" string specified in optics files. -integer, parameter, public :: ot_length = 32 - public :: rad_gas_index public :: get_number_sw_bands, & diff --git a/src/physics/rrtmg/slingo.F90 b/src/physics/rrtmg/slingo.F90 index aedb44bcee..b9d68565ec 100644 --- a/src/physics/rrtmg/slingo.F90 +++ b/src/physics/rrtmg/slingo.F90 @@ -9,7 +9,7 @@ module slingo use ppgrid, only: pcols, pver, pverp use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx -use radconstants, only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag, get_sw_spectral_boundaries +use radconstants, only: nswbands, nlwbands, idx_sw_diag, idx_lw_diag, get_sw_spectral_boundaries use cam_abortutils, only: endrun use cam_history, only: outfld diff --git a/src/physics/rrtmgp/cloud_rad_props.F90 b/src/physics/rrtmgp/cloud_rad_props.F90 index 1099fb714a..1581e04d9a 100644 --- a/src/physics/rrtmgp/cloud_rad_props.F90 +++ b/src/physics/rrtmgp/cloud_rad_props.F90 @@ -7,7 +7,7 @@ module cloud_rad_props use ppgrid, only: pcols, pver, pverp use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx -use radconstants, only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag +use radconstants, only: nswbands, nlwbands, idx_sw_diag use cam_abortutils, only: endrun use rad_constituents, only: iceopticsfile, liqopticsfile use oldcloud, only: oldcloud_lw, old_liq_get_rad_props_lw, old_ice_get_rad_props_lw, oldcloud_init diff --git a/src/physics/rrtmgp/ebert_curry.F90 b/src/physics/rrtmgp/ebert_curry.F90 index a1e1c031b1..c04a864ef0 100644 --- a/src/physics/rrtmgp/ebert_curry.F90 +++ b/src/physics/rrtmgp/ebert_curry.F90 @@ -7,7 +7,7 @@ module ebert_curry use ppgrid, only: pcols, pver, pverp use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx -use radconstants, only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag, get_sw_spectral_boundaries +use radconstants, only: nswbands, nlwbands, get_sw_spectral_boundaries use cam_abortutils, only: endrun use cam_history, only: outfld @@ -143,10 +143,7 @@ subroutine cloud_rad_props_get_sw(state, pbuf, & tau_w_g(:,1:ncol,:) = 0._r8 tau_w_f(:,1:ncol,:) = 0._r8 - call ec_ice_optics_sw (state, pbuf, tau, tau_w, tau_w_g, tau_w_f, oldicewp=.true.) -! call outfld ('CI_OD_SW_OLD', ice_tau(idx_sw_diag,:,:), pcols, lchnk) - end subroutine cloud_rad_props_get_sw !============================================================================== @@ -182,7 +179,6 @@ subroutine cloud_rad_props_get_lw(state, pbuf, cld_abs_od, diagnosticindex, oldl cld_abs_od = 0._r8 call ec_ice_get_rad_props_lw(state, pbuf, cld_abs_od, oldicewp=.true.) - !call outfld('CI_OD_LW_OLD', ice_tau_abs_od(idx_lw_diag ,:,:), pcols, lchnk) end subroutine cloud_rad_props_get_lw @@ -390,18 +386,11 @@ subroutine ec_ice_get_rad_props_lw(state, pbuf, abs_od, oldicewp) cldtau(i,k) = kabs*cwp(i,k) end do end do -! + do lwband = 1,nlwbands abs_od(lwband,1:ncol,1:pver)=cldtau(1:ncol,1:pver) enddo - !if(oldicewp) then - ! call outfld('CIWPTH_OLD',cicewp(:,:)/1000,pcols,lchnk) - !else - ! call outfld('CIWPTH_OLD',iciwpth(:,:),pcols,lchnk) - !endif - !call outfld('CI_OD_LW_OLD',cldtau(:,:),pcols,lchnk) - end subroutine ec_ice_get_rad_props_lw !============================================================================== diff --git a/src/physics/rrtmgp/oldcloud.F90 b/src/physics/rrtmgp/oldcloud.F90 index 609c6b4668..06a91b232e 100644 --- a/src/physics/rrtmgp/oldcloud.F90 +++ b/src/physics/rrtmgp/oldcloud.F90 @@ -7,7 +7,7 @@ module oldcloud use ppgrid, only: pcols, pver, pverp use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field -use radconstants, only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag, get_sw_spectral_boundaries +use radconstants, only: nswbands, nlwbands, get_sw_spectral_boundaries use cam_abortutils, only: endrun use cam_history, only: outfld use rad_constituents, only: iceopticsfile, liqopticsfile @@ -226,12 +226,6 @@ subroutine old_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, li end do ! End do k=1,pver end do ! nswbands - !call outfld('CL_OD_SW_OLD',liq_tau(idx_sw_diag,:,:), pcols, lchnk) - !call outfld('REL_OLD',rel(:,:), pcols, lchnk) - !call outfld('CLWPTH_OLD',cliqwp(:,:), pcols, lchnk) - !call outfld('KEXT_OLD',kext(:,:), pcols, lchnk) - - end subroutine old_liquid_optics_sw !============================================================================== diff --git a/src/physics/rrtmgp/rad_solar_var.F90 b/src/physics/rrtmgp/rad_solar_var.F90 deleted file mode 100644 index 0cf996e901..0000000000 --- a/src/physics/rrtmgp/rad_solar_var.F90 +++ /dev/null @@ -1,145 +0,0 @@ -!------------------------------------------------------------------------------- -! This module uses the Lean solar irradiance data to provide a solar cycle -! scaling factor used in heating rate calculations -!------------------------------------------------------------------------------- -module rad_solar_var - - use radconstants, only : nswbands - use shr_kind_mod , only : r8 => shr_kind_r8 - use solar_irrad_data, only : sol_irrad, we, nbins, has_spectrum, sol_tsi - use solar_irrad_data, only : do_spctrl_scaling - use cam_abortutils, only : endrun - - implicit none - save - - private - public :: rad_solar_var_init - public :: get_variability - - real(r8), allocatable :: ref_band_irrad(:) ! scaling will be relative to ref_band_irrad in each band - real(r8), allocatable :: irrad(:) ! solar irradiance at model timestep in each band - real(r8) :: tsi_ref ! total solar irradiance assumed by RRTMGP - - real(r8), allocatable :: radbinmax(:) - real(r8), allocatable :: radbinmin(:) - -!------------------------------------------------------------------------------- -contains -!------------------------------------------------------------------------------- - - subroutine rad_solar_var_init( ) - use radconstants, only : get_sw_spectral_boundaries - use radconstants, only : get_ref_solar_band_irrad - use radconstants, only : get_ref_total_solar_irrad - - integer :: i - integer :: ierr - integer :: yr, mon, tod - integer :: radmax_loc - - - if ( do_spctrl_scaling ) then - - if ( .not.has_spectrum ) then - call endrun('rad_solar_var_init: solar input file must have irradiance spectrum') - endif - - allocate (radbinmax(nswbands),stat=ierr) - if (ierr /= 0) then - call endrun('rad_solar_var_init: Error allocating space for radbinmax') - end if - - allocate (radbinmin(nswbands),stat=ierr) - if (ierr /= 0) then - call endrun('rad_solar_var_init: Error allocating space for radbinmin') - end if - - allocate (ref_band_irrad(nswbands), stat=ierr) - if (ierr /= 0) then - call endrun('rad_solar_var_init: Error allocating space for ref_band_irrad') - end if - - allocate (irrad(nswbands), stat=ierr) - if (ierr /= 0) then - call endrun('rad_solar_var_init: Error allocating space for irrad') - end if - - call get_sw_spectral_boundaries(radbinmin, radbinmax, 'nm') - - ! Make sure that the far-IR is included, even if RRTMG does not - ! extend that far down. 10^5 nm corresponds to a wavenumber of - ! 100 cm^-1. - radmax_loc = maxloc(radbinmax,1) - radbinmax(radmax_loc) = max(100000._r8,radbinmax(radmax_loc)) - - ! for rrtmg, reference spectrum from rrtmg - call get_ref_solar_band_irrad( ref_band_irrad ) - - else - - call get_ref_total_solar_irrad(tsi_ref) - - endif - - end subroutine rad_solar_var_init - -!------------------------------------------------------------------------------- -!------------------------------------------------------------------------------- - subroutine get_variability( sfac ) - - real(r8), intent(out) :: sfac(nswbands) ! scaling factors for CAM heating - - integer :: yr, mon, day, tod - - if ( do_spctrl_scaling ) then - call integrate_spectrum( nbins, nswbands, we, radbinmin, radbinmax, sol_irrad, irrad) - sfac(:nswbands) = irrad(:nswbands)/ref_band_irrad(:nswbands) - else - sfac(:nswbands) = sol_tsi/tsi_ref - endif - - end subroutine get_variability - -!------------------------------------------------------------------------------- -! private method......... -!------------------------------------------------------------------------------- - - subroutine integrate_spectrum( nsrc, ntrg, src_x, min_trg, max_trg, src, trg ) - - use mo_util, only : rebin - - implicit none - - !--------------------------------------------------------------- - ! ... dummy arguments - !--------------------------------------------------------------- - integer, intent(in) :: nsrc ! dimension source array - integer, intent(in) :: ntrg ! dimension target array - real(r8), intent(in) :: src_x(nsrc+1) ! source coordinates - real(r8), intent(in) :: max_trg(ntrg) ! target coordinates - real(r8), intent(in) :: min_trg(ntrg) ! target coordinates - real(r8), intent(in) :: src(nsrc) ! source array - real(r8), intent(out) :: trg(ntrg) ! target array - - !--------------------------------------------------------------- - ! ... local variables - !--------------------------------------------------------------- - real(r8) :: trg_x(2), targ(1) ! target coordinates - integer :: i - - do i = 1, ntrg - - trg_x(1) = min_trg(i) - trg_x(2) = max_trg(i) - - call rebin( nsrc, 1, src_x, trg_x, src(1:nsrc), targ(:) ) - ! W/m2/nm --> W/m2 - trg( i ) = targ(1)*(trg_x(2)-trg_x(1)) - - enddo - - - end subroutine integrate_spectrum - -end module rad_solar_var diff --git a/src/physics/rrtmgp/radconstants.F90 b/src/physics/rrtmgp/radconstants.F90 index e573bfb792..d086d1ce16 100644 --- a/src/physics/rrtmgp/radconstants.F90 +++ b/src/physics/rrtmgp/radconstants.F90 @@ -3,16 +3,6 @@ module radconstants ! This module contains constants that are specific to the radiative transfer ! code used in the RRTMGP model. -! This comment from E3SM implementation, and is entirely relevant here: -! TODO: Should this data be handled in a more robust way? Much of this contains -! explicit mappings to indices, which would probably be better handled with get_ -! functions. I.e., get_nswbands() could query the kdist objects in case of -! RRTMGP, and the diag indices could look up the actual bands used in the kdist -! objects as well. On that note, this module should probably go away if -! possible in the future, and we should provide more robust access to the -! radiation interface. - - use shr_kind_mod, only: r8 => shr_kind_r8 use cam_abortutils, only: endrun @@ -20,228 +10,93 @@ module radconstants private save -! Number of bands in SW and LW (these will be checked when RRTMGP initializes) +! Number of bands in SW and LW. These values must match data in the RRTMGP coefficients datasets. +! But they are needed to allocate space in the physics buffer and need to be available before the +! RRTMGP datasets are read. So they are set as parameters here and checked in radiation_init after +! the datasets are read. integer, parameter, public :: nswbands = 14 integer, parameter, public :: nlwbands = 16 -! Band limits (these get also get set at initialization) -real(r8), public, allocatable :: wavenumber_low_shortwave(:) -real(r8), public, allocatable :: wavenumber_high_shortwave(:) -real(r8), public, allocatable :: wavenumber_low_longwave(:) -real(r8), public, allocatable :: wavenumber_high_longwave(:) -! Reference irradiance per band -real(r8), public, allocatable :: solar_ref_band_irradiance(:) -real(r8), public, protected :: ref_tsi - -! SHORTWAVE DATA - - -! Wavenumbers of band boundaries -! -! Note: Currently rad_solar_var extends the lowest band down to -! 100 cm^-1 if it is too high to cover the far-IR. Any changes meant -! to affect IR solar variability should take note of this. - -! NOTE: these follow the non-monotonic ordering used for RRTMG -! - This is necessary because the optical properties files made for RRTMG use this order too. - -! NOTE: aside from order, as noted, these values match the ones in -! RRTMGP coefficients files. But I think we should be *setting* these -! values based on what is in that file, rather than hard-coding it here. - -! BPM: comment this data structure --> set it from radiation_init -! real(r8),parameter :: wavenumber_low_shortwave(nswbands) = & ! in cm^-1 -! (/2600._r8, 3250._r8, 4000._r8, 4650._r8, 5150._r8, 6150._r8, 7700._r8, & -! 8050._r8,12850._r8,16000._r8,22650._r8,29000._r8,38000._r8, 820._r8/) -! real(r8),parameter :: wavenumber_high_shortwave(nswbands) = & ! in cm^-1 -! (/3250._r8, 4000._r8, 4650._r8, 5150._r8, 6150._r8, 7700._r8, 8050._r8, & -! 12850._r8,16000._r8,22650._r8,29000._r8,38000._r8,50000._r8, 2600._r8/) - -! Mapping from RRTMG shortwave bands to RRTMGP -integer, parameter, dimension(14), public :: rrtmg_to_rrtmgp_swbands = & - (/ & - 14, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 & - /) - -! BPM <-- commented this block. Replaced by allocatable, get values by calling set_irrad_by_band --> -! Solar irradiance at 1 A.U. in W/m^2 assumed by radiation code -! Rescaled so that sum is precisely 1368.22 and fractional amounts sum to 1.0 -! real(r8), parameter :: solar_ref_band_irradiance(nswbands) = & -! (/ & -! 12.11_r8, 20.3600000000001_r8, 23.73_r8, & -! 22.43_r8, 55.63_r8, 102.93_r8, 24.29_r8, & -! 345.74_r8, 218.19_r8, 347.20_r8, & -! 129.49_r8, 50.15_r8, 3.08_r8, 12.89_r8 & -! /) - -! These are indices to the band for diagnostic output -! CHANGE: rather than make these parameters, provide subroutines that set them -! using the function get_band_index_by_value (which should be called on initializing radiation) -! integer, parameter, public :: idx_sw_diag = 10 ! index to sw visible band (441 - 625 nm) -! integer, parameter, public :: idx_nir_diag = 8 ! index to sw near infrared (778-1240 nm) band -! integer, parameter, public :: idx_uv_diag = 11 ! index to sw uv (345-441 nm) band - -! integer, parameter, public :: rrtmg_sw_cloudsim_band = 9 ! rrtmgp band for .67 micron -! integer, parameter, public :: rrtmgp_sw_cloudsim_band = 10 ! b/c one band moves to beginning - -integer, public :: idx_sw_diag ! index to sw visible band (441 - 625 nm) -integer, public :: idx_nir_diag! index to sw near infrared (778-1240 nm) band -integer, public :: idx_uv_diag ! index to sw uv (345-441 nm) band - -! CHANGE: instead of setting rrtmg[p]_sw_cloudsim_band in radconstants, just make it in radiation -! rrtmgp_sw_cloudsim_band = get_band_index_by_value('sw', 0.67_r8, 'micron') ! rrtmgp band for .67 micron -! same for lw: -! rrtmgp_lw_cloudsim_band = get_band_index_by_value('lw', 10.5_r8, 'micron') - -! Number of evenly spaced intervals in rh -! The globality of this mesh may not be necessary -! Perhaps it could be specific to the aerosol -! But it is difficult to see how refined it must be -! for lookup. This value was found to be sufficient -! for Sulfate and probably necessary to resolve the -! high variation near rh = 1. Alternative methods -! were found to be too slow. -! Optimal approach would be for cam to specify size of aerosol -! based on each aerosol's characteristics. Radiation -! should know nothing about hygroscopic growth! -integer, parameter, public :: nrh = 1000 - -! LONGWAVE DATA - -! These are indices to the band for diagnostic output (see comment above about change) -! integer, parameter, public :: idx_lw_diag = 7 ! index to (H20 window) LW band -integer, public :: idx_lw_diag - - -! These are commented, and intended to be replaced by reading the RRTMGP optics object -! real(r8), parameter :: wavenumber_low_longwave(nlwbands) = &! Longwave spectral band limits (cm-1) -! (/ 10._r8, 350._r8, 500._r8, 630._r8, 700._r8, 820._r8, 980._r8, 1080._r8, & -! 1180._r8, 1390._r8, 1480._r8, 1800._r8, 2080._r8, 2250._r8, 2380._r8, 2600._r8 /) - -! real(r8), parameter :: wavenumber_high_longwave(nlwbands) = &! Longwave spectral band limits (cm-1) -! (/ 350._r8, 500._r8, 630._r8, 700._r8, 820._r8, 980._r8, 1080._r8, 1180._r8, & -! 1390._r8, 1480._r8, 1800._r8, 2080._r8, 2250._r8, 2380._r8, 2600._r8, 3250._r8 /) +! Band limits (set from data in RRTMGP coefficient datasets) +real(r8), allocatable, target :: wavenumber_low_shortwave(:) +real(r8), allocatable, target :: wavenumber_high_shortwave(:) +real(r8), allocatable, target :: wavenumber_low_longwave(:) +real(r8), allocatable, target :: wavenumber_high_longwave(:) + +! These are indices to specific bands for diagnostic output and COSP input. +integer, public, protected :: idx_sw_diag = -1 ! band contains 500-nm wave +integer, public, protected :: idx_nir_diag = -1 ! band contains 1000-nm wave +integer, public, protected :: idx_uv_diag = -1 ! band contains 400-nm wave +integer, public, protected :: idx_lw_diag = -1 ! band contains 1000 cm-1 wave (H20 window) +integer, public, protected :: idx_sw_cloudsim = -1 ! band contains 670-nm wave (for COSP) +integer, public, protected :: idx_lw_cloudsim = -1 ! band contains 10.5 micron wave (for COSP) ! GASES TREATED BY RADIATION (line spectrae) +! These names are recognized by RRTMGP. They are in the coefficients files as +! lower case strings. These upper case names are used by CAM's namelist and can +! be used to initialize the ty_gas_conc object because the name matching is case +! insensitive. integer, public, parameter :: gasnamelength = 5 integer, public, parameter :: nradgas = 8 character(len=gasnamelength), public, parameter :: gaslist(nradgas) & = (/'H2O ','O3 ', 'O2 ', 'CO2 ', 'N2O ', 'CH4 ', 'CFC11', 'CFC12'/) ! what is the minimum mass mixing ratio that can be supported by radiation implementation? -real(r8), public, parameter :: minmmr(nradgas) & - = epsilon(1._r8) - -! Length of "optics type" string specified in optics files. -integer, parameter, public :: ot_length = 32 - -public :: rad_gas_index - -public :: get_sw_spectral_boundaries, & - get_lw_spectral_boundaries, & - get_ref_solar_band_irrad, & - get_ref_total_solar_irrad, & - get_idx_sw_diag, & - get_idx_nir_diag, & - get_idx_uv_diag, & - get_idx_lw_diag, & - get_band_index_by_value, & - set_wavenumber_bands,& - set_irrad_by_band, & - set_reference_tsi +real(r8), public, parameter :: minmmr(nradgas) = epsilon(1._r8) + +public :: & + set_wavenumber_bands, & + get_sw_spectral_boundaries, & + get_lw_spectral_boundaries, & + get_band_index_by_value, & + rad_gas_index !=============================================================================== contains !=============================================================================== -subroutine get_ref_total_solar_irrad(tsi) - ! provide Total Solar Irradiance assumed by RRTMGP - - real(r8), intent(out) :: tsi - - ! tsi = sum(solar_ref_band_irradiance) - tsi = ref_tsi - -end subroutine get_ref_total_solar_irrad -!------------------------------------------------------------------------------ -subroutine set_reference_tsi(tsi) - ! set ref_tsi to provide total solar irradiance - ! this usually comes from reading a file - ! provided by the radiation scheme developers - real(r8), intent(in) :: tsi - ref_tsi = tsi -end subroutine set_reference_tsi -!------------------------------------------------------------------------------ -subroutine get_ref_solar_band_irrad( band_irrad ) - ! note: this shouldn't be used. - ! Instead, just use radconstants, only: solar_ref_band_irradiance - ! to access the data directly - ! solar irradiance in each band (W/m^2) - real(r8), intent(out) :: band_irrad(nswbands) - - if (allocated(solar_ref_band_irradiance)) then - band_irrad = solar_ref_band_irradiance - else - ! what to do - end if +subroutine set_wavenumber_bands(swlw, nbands, values) -end subroutine get_ref_solar_band_irrad + ! Set the low and high limits of the wavenumber grid for sw or lw. + ! Values comes from RRTMGP coefficients datasets. + ! Also set band indices for bands containing specific wavelengths. -!------------------------------------------------------------------------------ + character(*), intent(in) :: swlw ! which bands to set ['sw', 'lw'] + integer, intent(in) :: nbands + real(r8), intent(in) :: values(2,nbands) ! cm-1 -subroutine set_wavenumber_bands(swlw, nbands, values) - ! set the low and high limits of the wavenumber grid for sw or lw - ! expect that values comes from RRTMGP method get_band_lims_wavenumber - character(*), intent(in) :: swlw ! which set of bands to set ['sw', 'lw'] - integer, intent(in) :: nbands - real(r8), intent(in) :: values(2,nbands) select case(swlw) case ('sw') allocate(wavenumber_low_shortwave(nbands)) allocate(wavenumber_high_shortwave(nbands)) wavenumber_low_shortwave = values(1,:) wavenumber_high_shortwave = values(2,:) + + idx_sw_diag = get_band_index_by_value('sw', 500.0_r8, 'nm') + idx_nir_diag = get_band_index_by_value('sw', 1000.0_r8, 'nm') + idx_uv_diag = get_band_index_by_value('sw', 400._r8, 'nm') + idx_sw_cloudsim = get_band_index_by_value('sw', 0.67_r8, 'micron') + case ('lw') allocate(wavenumber_low_longwave(nbands)) allocate(wavenumber_high_longwave(nbands)) wavenumber_low_longwave = values(1,:) wavenumber_high_longwave = values(2,:) - end select -end subroutine set_wavenumber_bands -!------------------------------------------------------------------------------ -subroutine get_lw_spectral_boundaries(low_boundaries, high_boundaries, units) - ! provide spectral boundaries of each longwave band - real(r8), intent(out) :: low_boundaries(nlwbands), high_boundaries(nlwbands) - character(*), intent(in) :: units ! requested units + idx_lw_diag = get_band_index_by_value('lw', 1000.0_r8, 'cm^-1') + idx_lw_cloudsim = get_band_index_by_value('lw', 10.5_r8, 'micron') - select case (units) - case ('inv_cm','cm^-1','cm-1') - low_boundaries = wavenumber_low_longwave - high_boundaries = wavenumber_high_longwave - case('m','meter','meters') - low_boundaries = 1.e-2_r8/wavenumber_high_longwave - high_boundaries = 1.e-2_r8/wavenumber_low_longwave - case('nm','nanometer','nanometers') - low_boundaries = 1.e7_r8/wavenumber_high_longwave - high_boundaries = 1.e7_r8/wavenumber_low_longwave - case('um','micrometer','micrometers','micron','microns') - low_boundaries = 1.e4_r8/wavenumber_high_longwave - high_boundaries = 1.e4_r8/wavenumber_low_longwave - case('cm','centimeter','centimeters') - low_boundaries = 1._r8/wavenumber_high_longwave - high_boundaries = 1._r8/wavenumber_low_longwave - case default - call endrun('get_lw_spectral_boundaries: spectral units not acceptable'//units) end select -end subroutine get_lw_spectral_boundaries +end subroutine set_wavenumber_bands !------------------------------------------------------------------------------ + subroutine get_sw_spectral_boundaries(low_boundaries, high_boundaries, units) + ! provide spectral boundaries of each shortwave band - real(r8), intent(out) :: low_boundaries(nswbands), high_boundaries(nswbands) + real(r8), intent(out) :: low_boundaries(nswbands), high_boundaries(nswbands) character(*), intent(in) :: units ! requested units select case (units) @@ -261,12 +116,44 @@ subroutine get_sw_spectral_boundaries(low_boundaries, high_boundaries, units) low_boundaries = 1._r8/wavenumber_high_shortwave high_boundaries = 1._r8/wavenumber_low_shortwave case default - call endrun('rad_constants.F90: requested spectral units not acceptable: '//units) + call endrun('rad_constants.F90: requested spectral units not recognized: '//units) end select end subroutine get_sw_spectral_boundaries !------------------------------------------------------------------------------ + +subroutine get_lw_spectral_boundaries(low_boundaries, high_boundaries, units) + + ! provide spectral boundaries of each longwave band + + real(r8), intent(out) :: low_boundaries(nlwbands), high_boundaries(nlwbands) + character(*), intent(in) :: units ! requested units + + select case (units) + case ('inv_cm','cm^-1','cm-1') + low_boundaries = wavenumber_low_longwave + high_boundaries = wavenumber_high_longwave + case('m','meter','meters') + low_boundaries = 1.e-2_r8/wavenumber_high_longwave + high_boundaries = 1.e-2_r8/wavenumber_low_longwave + case('nm','nanometer','nanometers') + low_boundaries = 1.e7_r8/wavenumber_high_longwave + high_boundaries = 1.e7_r8/wavenumber_low_longwave + case('um','micrometer','micrometers','micron','microns') + low_boundaries = 1.e4_r8/wavenumber_high_longwave + high_boundaries = 1.e4_r8/wavenumber_low_longwave + case('cm','centimeter','centimeters') + low_boundaries = 1._r8/wavenumber_high_longwave + high_boundaries = 1._r8/wavenumber_low_longwave + case default + call endrun('get_lw_spectral_boundaries: spectral units not recognized: '//units) + end select + +end subroutine get_lw_spectral_boundaries + +!------------------------------------------------------------------------------ + integer function rad_gas_index(gasname) ! return the index in the gaslist array of the specified gasname @@ -283,48 +170,36 @@ integer function rad_gas_index(gasname) enddo call endrun ("rad_gas_index: can not find gas with name "//gasname) end function rad_gas_index -!------------------------------------------------------------------------------ -subroutine get_idx_sw_diag() - idx_sw_diag = get_band_index_by_value('sw', 500.0_r8, 'nm') -end subroutine -subroutine get_idx_nir_diag() - idx_nir_diag = get_band_index_by_value('sw', 1000.0_r8, 'nm') -end subroutine +!------------------------------------------------------------------------------ -subroutine get_idx_uv_diag() - idx_uv_diag = get_band_index_by_value('sw', 400._r8, 'nm') -end subroutine +function get_band_index_by_value(swlw, targetvalue, units) result(ans) -subroutine get_idx_lw_diag() - idx_lw_diag = get_band_index_by_value('lw', 1000.0_r8, 'cm^-1') - ! value chosen to match the band used in CESM1/CESM2 -end subroutine + ! Find band index for requested wavelength/wavenumber. -function get_band_index_by_value(swlw, targetvalue, units) result(ans) - character(len=*),intent(in) :: swlw ! sw or lw bands - real(r8),intent(in) :: targetvalue - character(len=*),intent(in) :: units ! units of targetvalue + character(len=*), intent(in) :: swlw ! sw or lw bands + real(r8), intent(in) :: targetvalue + character(len=*), intent(in) :: units ! units of targetvalue integer :: ans + ! local - real(r8), allocatable, dimension(:) :: lowboundaries, highboundaries + real(r8), pointer, dimension(:) :: lowboundaries, highboundaries real(r8) :: tgt integer :: nbnds, i select case (swlw) case ('sw','SW','shortwave') nbnds = nswbands - allocate(lowboundaries(nbnds), highboundaries(nbnds)) - lowboundaries = wavenumber_low_shortwave - highboundaries = wavenumber_high_shortwave + lowboundaries => wavenumber_low_shortwave + highboundaries => wavenumber_high_shortwave case ('lw', 'LW', 'longwave') nbnds = nlwbands - allocate(lowboundaries(nbnds), highboundaries(nbnds)) - lowboundaries = wavenumber_low_longwave - highboundaries = wavenumber_high_longwave + lowboundaries => wavenumber_low_longwave + highboundaries => wavenumber_high_longwave case default - call endrun('rad_constants.F90: get_band_index_by_value: type of bands not accepted '//swlw) + call endrun('radconstants.F90: get_band_index_by_value: type of bands not recognized: '//swlw) end select + ! band info is in cm^-1 but target value may be other units, ! so convert targetvalue to cm^-1 select case (units) @@ -339,43 +214,24 @@ function get_band_index_by_value(swlw, targetvalue, units) result(ans) case('cm','centimeter','centimeters') tgt = 1._r8/targetvalue case default - call endrun('rad_constants.F90: get_band_index_by_value: units not acceptable'//units) + call endrun('radconstants.F90: get_band_index_by_value: units not recognized: '//units) end select + ! now just loop through the array + ans = 0 do i = 1,nbnds if ((tgt > lowboundaries(i)) .and. (tgt <= highboundaries(i))) then ans = i exit end if end do - ! Do something if the answer is not found? -end function get_band_index_by_value - -subroutine set_irrad_by_band(solar_source, g2b) - ! Sets the solar irradiance in each shortwave band by summing the irradiance from gpoints. - ! solar_source = kdist_sw%solar_source <-- private TRY solar_source = kdist_sw%solar_source_quiet - ! g2b = kdist_sw%get_gpoint_bands() - real(r8), intent(in) :: solar_source(:) ! size ngpoints: irradiance per gpoint - integer, intent(in) :: g2b(:) ! size ngpoints: mapping from gpoint to band - integer :: i - allocate(solar_ref_band_irradiance(nswbands)) - solar_ref_band_irradiance(:) = 0.0_r8 - do i = 1,size(g2b) - solar_ref_band_irradiance(g2b(i)) = solar_ref_band_irradiance(g2b(i)) + solar_source(i) - end do -end subroutine set_irrad_by_band - -function get_irrad_by_band(solar_source, g2b) result(ans) - real(r8) :: solar_source(:) - integer :: g2b(:) - real(r8), allocatable :: ans(:) - if (.not. allocated(solar_ref_band_irradiance)) then - call set_irrad_by_band(solar_source, g2b) + if (ans == 0) then + call endrun('radconstants.F90: get_band_index_by_value: band not found: ') end if - allocate(ans(size(solar_ref_band_irradiance))) - ans = solar_ref_band_irradiance -end function get_irrad_by_band + +end function get_band_index_by_value +!------------------------------------------------------------------------------ end module radconstants diff --git a/src/physics/rrtmgp/radiation.F90 b/src/physics/rrtmgp/radiation.F90 index baf9620389..12955ae4ed 100644 --- a/src/physics/rrtmgp/radiation.F90 +++ b/src/physics/rrtmgp/radiation.F90 @@ -15,6 +15,7 @@ module radiation use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx use camsrfexch, only: cam_out_t, cam_in_t use physconst, only: cappa, cpair, gravit +use solar_irrad_data, only: sol_tsi use time_manager, only: get_nstep, is_first_restart_step, & get_curr_calday, get_step_size @@ -27,20 +28,9 @@ module radiation liqcldoptics, & icecldoptics -use radconstants, only: nswbands, nlwbands, & ! number of bands - idx_sw_diag, & ! indices for diagnostics - idx_nir_diag, & - idx_uv_diag, & - idx_lw_diag, & - get_idx_sw_diag, & ! sets the idx_*_diag in radconstants module - get_idx_nir_diag, & - get_idx_uv_diag, & - get_idx_lw_diag, & - rrtmg_to_rrtmgp_swbands, & ! maps bands between rrtmg and rrtmgp - get_band_index_by_value, & ! function that figures out band for a wavelength - gasnamelength, & - nradgas, & - gaslist +use radconstants, only: nswbands, nlwbands, idx_sw_diag, idx_nir_diag, idx_uv_diag, & + idx_lw_diag, idx_sw_cloudsim, idx_lw_cloudsim, & + nradgas, gasnamelength, gaslist use mo_gas_concentrations, only: ty_gas_concs use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp @@ -67,14 +57,16 @@ module radiation PIO_NOWRITE, & pio_closefile -use cam_abortutils, only: endrun -use error_messages, only: handle_err -use cam_logfile, only: iulog use scamMod, only: scm_crm_mode, single_column, have_cld, cldobs use cospsimulator_intr, only: docosp, cospsimulator_intr_init, & cospsimulator_intr_run, cosp_nradsteps +use string_utils, only: to_lower +use cam_abortutils, only: endrun +use error_messages, only: handle_err +use cam_logfile, only: iulog + implicit none private @@ -94,7 +86,6 @@ module radiation integer,public, allocatable :: cosp_cnt(:) ! counter for cosp integer,public :: cosp_cnt_init = 0 !initial value for cosp counter -integer, public :: sw_cloudsim_band, lw_cloudsim_band ! radiation bands that COSP uses real(r8), public, protected :: nextsw_cday ! future radiation calday for surface models @@ -181,6 +172,10 @@ module radiation logical :: graupel_in_rad = .false. ! graupel in radiation code logical :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the coszrs calculation +! active_calls is set by a rad_constituents method after parsing namelist input +! for the rad_climate and rad_diag_N entries. +logical :: active_calls(0:N_DIAG) + ! Physics buffer indices integer :: qrs_idx = 0 integer :: qrl_idx = 0 @@ -222,29 +217,22 @@ module radiation ! 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 - -! SW coefficients -type(ty_gas_optics_rrtmgp) :: kdist_sw ! bpm changed here -integer :: ngpt_sw +! Gas optics objects contain the data read from the coefficients files. +type(ty_gas_optics_rrtmgp) :: kdist_lw +type(ty_gas_optics_rrtmgp) :: kdist_sw -! data to go from bands to gpoints (bpm) -integer, allocatable :: band2gpt_sw(:,:) ! n[s,l]wbands come from radconstants for now +! data to go from bands to gpoints +integer, allocatable :: band2gpt_sw(:,:) integer, allocatable :: band2gpt_lw(:,:) +! Mapping from RRTMG shortwave bands to RRTMGP. Currently needed to continue using +! the SW optics datasets from RRTMG (even thought there is a slight mismatch in the +! band boundaries of the 2 bands that overlap with the LW bands). +integer, parameter, dimension(14) :: rrtmg_to_rrtmgp_swbands = & + [ 14, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 ] -! Gases to use in the radiative calculations. -! RRTMGP kdist initialization needs to know the names of the -! gases before these are available via the rad_cnst interface. -! TODO: Move this to namelist or somewhere appropriate. -! NOTE: This list is not the same as `gaslist` in radconstants; is that a problem? Implication for diagnostic calls? -! character(len=5), dimension(10) :: active_gases = (/ & -! 'H2O ', 'CO2 ', 'O3 ', 'N2O ', & -! 'CO ', 'CH4 ', 'O2 ', 'N2 ', & -! 'CFC11', 'CFC12' /) -! BPM: use radconstants to define the active gases: -character(len=gasnamelength), dimension(nradgas) :: active_gases = gaslist +! lower case version of gaslist for RRTMGP +character(len=gasnamelength) :: gaslist_lc(nradgas) type(var_desc_t) :: cospcnt_desc ! cosp type(var_desc_t) :: nextsw_cday_desc @@ -474,13 +462,12 @@ subroutine radiation_init(pbuf2d) use physics_buffer, only: pbuf_get_index, pbuf_set_field use phys_control, only: phys_getopts - use rad_solar_var, only: rad_solar_var_init ! This initializes total solar irradiance use radiation_data, only: rad_data_init use cloud_rad_props, only: cloud_rad_props_init use modal_aer_opt, only: modal_aer_opt_init use rrtmgp_inputs, only: rrtmgp_inputs_init use time_manager, only: is_first_step - use radconstants, only: set_wavenumber_bands, set_irrad_by_band, set_reference_tsi + use radconstants, only: set_wavenumber_bands ! arguments type(physics_buffer_desc), pointer :: pbuf2d(:,:) @@ -492,8 +479,7 @@ subroutine radiation_init(pbuf2d) ! -- needed for the kdist initialization routines type(ty_gas_concs) :: available_gases - integer :: icall, nmodes - logical :: active_calls(0:N_DIAG) + integer :: i, icall, nmodes integer :: nstep ! current timestep number logical :: history_amwg ! output the variables used by the AMWG diag package logical :: history_vdiag ! output the variables used by the AMWG variability diag package @@ -504,7 +490,6 @@ subroutine radiation_init(pbuf2d) integer :: ierr integer :: dtime - real(r8) :: ref_tsi character(len=*), parameter :: sub = 'radiation_init' !----------------------------------------------------------------------- @@ -538,11 +523,20 @@ subroutine radiation_init(pbuf2d) 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 + ! Create lowercase version of the gaslist for RRTMGP. The ty_gas_concs objects + ! work with CAM's uppercase names, but other objects that get input from the gas + ! concs objects don't work. + do i = 1, nradgas + gaslist_lc(i) = to_lower(gaslist(i)) + end do + + errmsg = available_gases%init(gaslist_lc) + if (len_trim(errmsg) > 0) then + call endrun(sub//': ERROR: available_gases%init: '//trim(errmsg)) + end if call coefs_init(coefs_lw_file, kdist_lw, available_gases, band2gpt_lw) - call coefs_init(coefs_sw_file, kdist_sw, available_gases, band2gpt_sw, ref_tsi) ! bpm : these now provide band2gpt which should be global - call set_reference_tsi(ref_tsi) + call coefs_init(coefs_sw_file, kdist_sw, available_gases, band2gpt_sw) ! check number of sw/lw bands in gas optics files if (kdist_sw%get_nband() /= nswbands) then @@ -563,25 +557,11 @@ subroutine radiation_init(pbuf2d) call set_wavenumber_bands('sw', kdist_sw%get_nband(), kdist_sw%get_band_lims_wavenumber()) call set_wavenumber_bands('lw', kdist_lw%get_nband(), kdist_lw%get_band_lims_wavenumber()) - call rad_solar_var_init() ! sets the total solar irradiance (I wonder whether this should use kdist information instead of radconstants; alternative use kdist%set_tsi to ensure consistency?) call rrtmgp_inputs_init(ktopcamm, ktopradm, ktopcami, ktopradi) ! this sets these values as module data in rrtmgp_inputs call rad_data_init(pbuf2d) ! initialize output fields for offline driver call cloud_rad_props_init() - ngpt_sw = kdist_sw%get_ngpt() - - ! bpm: set the indices used for diagnostics using specific band: - call get_idx_sw_diag() ! index to sw visible band (441 - 625 nm) - call get_idx_nir_diag() ! index to sw near infrared (778-1240 nm) band - call get_idx_uv_diag() ! index to sw uv (345-441 nm) band - if (docosp) then - sw_cloudsim_band = get_band_index_by_value('sw', 0.67_r8, 'micron') ! rrtmgp band for .67 micron - lw_cloudsim_band = get_band_index_by_value('lw', 10.5_r8, 'micron') - end if - call get_idx_lw_diag() - - if (is_first_step()) then call pbuf_set_field(pbuf2d, qrl_idx, 0._r8) end if @@ -906,7 +886,7 @@ subroutine radiation_tend( & !----------------------------------------------------------------------- ! - ! Driver for radiation computation. + ! CAM driver for radiation computation. ! !----------------------------------------------------------------------- @@ -915,7 +895,6 @@ subroutine radiation_tend( & use cam_control_mod, only: eccen, mvelpp, lambm0, obliqr use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz - use mo_gas_concentrations, only: ty_gas_concs use rrtmgp_inputs, only: rrtmgp_set_state, rrtmgp_set_gases_lw, rrtmgp_set_cloud_lw, & rrtmgp_set_aer_lw, rrtmgp_set_gases_sw, rrtmgp_set_cloud_sw, & rrtmgp_set_aer_sw @@ -936,8 +915,8 @@ subroutine radiation_tend( & use mo_fluxes_byband, only: ty_fluxes_byband - ! use mo_rrtmgp_clr_all_sky, only: rte_lw, rte_sw - use rrtmgp_driver, only: rte_lw, rte_sw + ! RRTMGP drivers for flux calculations. + use rrtmgp_driver, only: rte_lw, rte_sw use radheat, only: radheat_tend @@ -979,8 +958,8 @@ subroutine radiation_tend( & ! chunk_column_index = IdxDay(daylight_column_index) integer :: Nday ! Number of daylight columns integer :: Nnite ! Number of night columns - integer :: IdxDay(pcols) ! Indices of daylight columns -- Dimension is pcols, and is filled from beginning, so idxday(1:nday) are the indices of daylit columns. - integer :: IdxNite(pcols) ! Indices of night columns + integer :: IdxDay(pcols) ! chunk indices of daylight columns + integer :: IdxNite(pcols) ! chunk indices of night columns integer :: itim_old @@ -1016,7 +995,6 @@ subroutine radiation_tend( & real(r8), allocatable :: coszrs_day(:) real(r8), allocatable :: alb_dir(:,:) real(r8), allocatable :: alb_dif(:,:) - real(r8) :: tsi ! cloud radiative parameters are "in cloud" not "in cell" @@ -1074,13 +1052,12 @@ subroutine radiation_tend( & type(ty_optical_props_1scl) :: cloud_lw type(ty_optical_props_2str) :: cloud_sw - ! Irradiance integer :: icall ! index through climate/diagnostic radiation calls - logical :: active_calls(0:N_DIAG) - ! gas vmr + ! gas vmr. Separate objects because SW only does calculations for daylight columns. type(ty_gas_concs) :: gas_concs_lw type(ty_gas_concs) :: gas_concs_sw + ! RRTMGP aerosol objects type(ty_optical_props_1scl) :: aer_lw type(ty_optical_props_2str) :: aer_sw @@ -1264,13 +1241,10 @@ subroutine radiation_tend( & cam_in, & ! input (%lwup, %aldir, %asdir, %aldif, %asdif) ncol, & ! input nlay, & ! input - nlwbands, & ! input - nswbands, & ! input - ngpt_sw, & ! input nday, & ! input idxday, & ! input, [would prefer to truncate as 1:ncol] coszrs, & ! input - kdist_sw, & ! input (from init) ! removed: eccf, & ! input + kdist_sw, & ! input (from init) band2gpt_sw, & ! input (from init), gpoints by band t_sfc, & ! output emis_sfc, & ! output @@ -1282,14 +1256,10 @@ subroutine radiation_tend( & pint_day, & ! output coszrs_day, & ! output alb_dir, & ! output - alb_dif, & ! output - tsi & ! output, total solar irradiance (not scaled) - ) + alb_dif) ! output - !!--> Set TSI used in radiation to the value in the solar forcing file. - !!--> This replaces get_variability() and does same thing. - !!--> The Earth-Sun distance (eccf) provides another scaling, applied later. - errmsg = kdist_sw%set_tsi(tsi) ! scales the TSI but does not change spectral distribution + ! Set TSI used in rrtmgp to the value from CAM's solar forcing file. + errmsg = kdist_sw%set_tsi(sol_tsi) if (len_trim(errmsg) > 0) then call endrun(sub//': ERROR: kdist_sw%set_tsi: '//trim(errmsg)) end if @@ -1468,16 +1438,20 @@ subroutine radiation_tend( & if (write_output) then call radiation_output_cld(lchnk, ncol, rd) end if - ! - ! SHORTWAVE CALCULATION(S) - ! - ! Get the active climate/diagnostic shortwave calculations - call rad_cnst_get_call_list(active_calls) + + !=============================! + ! SHORTWAVE flux calculations ! + !=============================! + + ! initialize object for gas concentrations + errmsg = gas_concs_sw%init(gaslist_lc) + if (len_trim(errmsg) > 0) then + call endrun(sub//': ERROR: gas_concs_sw%init: '//trim(errmsg)) + end if ! The climate (icall==0) calculation must occur last. do icall = N_DIAG, 0, -1 if (active_calls(icall)) then - call set_available_gases(active_gases, gas_concs_sw) ! set gas concentrations call rrtmgp_set_gases_sw( & ! Put gas volume mixing ratio into gas_concs_sw icall, & ! input @@ -1576,9 +1550,10 @@ subroutine radiation_tend( & ! This happens between SW and LW (Why?) call rad_cnst_out(0, state, pbuf) - ! - ! -- LONGWAVE -- - ! + !============================! + ! LONGWAVE flux calculations ! + !============================! + if (dolw) then if (oldcldoptics) then call cloud_rad_props_get_lw(state, pbuf, cld_lw_abs, oldcloud=.true.) @@ -1653,23 +1628,21 @@ subroutine radiation_tend( & nlay, & kdist_lw%get_band_lims_wavenumber(), & name='longwave aerosol optics') - if (len_trim(errmsg) > 0) then - call endrun(sub//': ERROR: aer_lw%init_1scalar: '//trim(errmsg)) + call endrun(sub//': ERROR: aer_lw%alloc_1scalar: '//trim(errmsg)) end if - call rad_cnst_get_call_list(active_calls) ! get list of diagnostic calls + ! initialize object for gas concentrations + errmsg = gas_concs_lw%init(gaslist_lc) + if (len_trim(errmsg) > 0) then + call endrun(sub//': ERROR, gas_concs_lw%init: '//trim(errmsg)) + end if ! The climate (icall==0) calculation must occur last. do icall = N_DIAG, 0, -1 if (active_calls(icall)) then - ! initialize the gas concentrations - call set_available_gases(active_gases, gas_concs_lw) -! errmsg = gas_concs_lw%init(active_gases) -! if (len_trim(errmsg) > 0) then -! call endrun(sub//': ERROR code returned by gas_concs_lw%init: '//trim(errmsg)) -! end if + call rrtmgp_set_gases_lw(icall, state, pbuf, nlay, gas_concs_lw) call aer_rad_props_lw( & ! get absorption optical depth @@ -1703,14 +1676,7 @@ subroutine radiation_tend( & aer_props=aer_lw & ! optional input, (rrtmgp_set_aer_lw) ) ! note inc_flux is an optional input, but as defined in set_rrtmgp_state, it is only for shortwave if (len_trim(errmsg) > 0) then - ! - ! DEBUG -- if we die here, find out why - ! - write(iulog,*) '** [radiation_tend] DIAGNOSE LW CRASH **' - do i = 1,ncol - write(iulog,*) 'ncol = ',ncol,' t_sfc = ',t_sfc(i),' AT LOCATION lat = ', clat(i), ' lon = ', clon(i) - end do - call endrun(sub//': ERROR code returned by rte_lw: '//trim(errmsg)) + call endrun(sub//': ERROR: rte_lw: '//trim(errmsg)) end if ! ! -- longwave output -- @@ -1743,7 +1709,7 @@ subroutine radiation_tend( & if (docosp) then ! initialize and calculate emis emis(:,:) = 0._r8 - emis(:ncol,:) = 1._r8 - exp(-cld_lw_abs(lw_cloudsim_band,:ncol,:)) + emis(:ncol,:) = 1._r8 - exp(-cld_lw_abs(idx_lw_cloudsim,:ncol,:)) call outfld('EMIS', emis, pcols, lchnk) ! compute grid-box mean SW and LW snow optical depth for use by COSP @@ -1756,13 +1722,13 @@ subroutine radiation_tend( & ! Add graupel to snow tau for cosp if (cldfgrau_idx > 0 .and. graupel_in_rad) then - gb_snow_tau(i,k) = snow_tau(sw_cloudsim_band,i,k)*cldfsnow(i,k) + & - grau_tau(sw_cloudsim_band,i,k)*cldfgrau(i,k) - gb_snow_lw(i,k) = snow_lw_abs(lw_cloudsim_band,i,k)*cldfsnow(i,k) + & - grau_lw_abs(lw_cloudsim_band,i,k)*cldfgrau(i,k) + gb_snow_tau(i,k) = snow_tau(idx_sw_cloudsim,i,k)*cldfsnow(i,k) + & + grau_tau(idx_sw_cloudsim,i,k)*cldfgrau(i,k) + gb_snow_lw(i,k) = snow_lw_abs(idx_lw_cloudsim,i,k)*cldfsnow(i,k) + & + grau_lw_abs(idx_lw_cloudsim,i,k)*cldfgrau(i,k) else - gb_snow_tau(i,k) = snow_tau(sw_cloudsim_band,i,k)*cldfsnow(i,k) - gb_snow_lw(i,k) = snow_lw_abs(lw_cloudsim_band,i,k)*cldfsnow(i,k) + gb_snow_tau(i,k) = snow_tau(idx_sw_cloudsim,i,k)*cldfsnow(i,k) + gb_snow_lw(i,k) = snow_lw_abs(idx_lw_cloudsim,i,k)*cldfsnow(i,k) end if end if end do @@ -1778,14 +1744,14 @@ subroutine radiation_tend( & ! N.B.: For snow optical properties, the GRID-BOX MEAN shortwave and longwave ! optical depths are passed. call cospsimulator_intr_run(state, pbuf, cam_in, emis, coszrs, & - cld_swtau_in=cld_tau(sw_cloudsim_band,:,:),& + cld_swtau_in=cld_tau(idx_sw_cloudsim,:,:),& snow_tau_in=gb_snow_tau, snow_emis_in=gb_snow_lw) cosp_cnt(lchnk) = 0 end if end if !!! *** END COSP *** - else ! if (dosw .or. dolw) --> no radiation being done. + else ! --> radiative flux calculations not updated ! convert radiative heating rates from Q*dp to Q for energy conservation ! qrs and qrl are whatever are in pbuf ! since those might have been multiplied by pdel, we actually need to divide by pdel @@ -1848,9 +1814,9 @@ subroutine radiation_tend( & call free_fluxes(flw) call free_fluxes(flwc) -!------------------------------------------------------------------------------- -contains -!------------------------------------------------------------------------------- + !------------------------------------------------------------------------------- + contains + !------------------------------------------------------------------------------- subroutine set_sw_diags() @@ -2255,18 +2221,17 @@ subroutine calc_col_mean(state, mmr_pointer, mean_value) end subroutine calc_col_mean -!=============================================================================== +!========================================================================================= -subroutine coefs_init(coefs_file, kdist, available_gases, band2gpt, tsi_default) +subroutine coefs_init(coefs_file, kdist, available_gases, band2gpt) ! Read data from coefficients file. Initialize the kdist object. + ! available_gases object provides the gas names that CAM provides. ! arguments character(len=*), intent(in) :: coefs_file class(ty_gas_optics_rrtmgp), intent(out) :: kdist - class(ty_gas_concs), intent(in) :: available_gases ! Which gases does the host model have available? - - real(r8), intent(out), optional :: tsi_default ! RRTMGP reference TSI + class(ty_gas_concs), intent(in) :: available_gases ! local variables type(file_desc_t) :: fh ! pio file handle @@ -2302,6 +2267,7 @@ subroutine coefs_init(coefs_file, kdist, available_gases, band2gpt, tsi_default) real(r8), dimension(:,:), allocatable :: totplnk real(r8), dimension(:,:,:,:), allocatable :: planck_frac real(r8), dimension(:), allocatable :: solar_src_quiet, solar_src_facular, solar_src_sunspot ! updated from solar_src + real(r8) :: tsi_default real(r8), dimension(:,:,:), allocatable :: rayl_lower, rayl_upper character(len=32), dimension(:), allocatable :: gas_minor, & identifier_minor, & @@ -2540,15 +2506,6 @@ subroutine coefs_init(coefs_file, kdist, available_gases, band2gpt, tsi_default) if (ierr /= PIO_NOERR) call endrun(sub//': error reading optimal_angle_fit') end if - ! solar_src - ! !bpm -- solar_source is not in file, there are solar_source_[facular, sunspot, quiet] - ! There's a method that adds them together to get solar_source. - ! ierr = pio_inq_varid(fh, 'solar_source', vid) - ! if (ierr == PIO_NOERR) then - ! allocate(solar_src(gpt)) - ! ierr = pio_get_var(fh, vid, solar_src) - ! if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source') - ! end if ierr = pio_inq_varid(fh, 'solar_source_quiet', vid) if (ierr == PIO_NOERR) then allocate(solar_src_quiet(gpt)) @@ -2568,7 +2525,6 @@ subroutine coefs_init(coefs_file, kdist, available_gases, band2gpt, tsi_default) if (ierr /= PIO_NOERR) call endrun(sub//': error reading solar_source_sunspot') end if - ! +bpm also need to have tsi_default, mg_default, and sb_default ierr = pio_inq_varid(fh, 'tsi_default', vid) if (ierr == PIO_NOERR) then ierr = pio_get_var(fh, vid, tsi_default) @@ -2836,32 +2792,7 @@ subroutine coefs_init(coefs_file, kdist, available_gases, band2gpt, tsi_default) if (allocated(rayl_upper)) deallocate(rayl_upper) end subroutine coefs_init - - -subroutine set_available_gases(gases, gas_concentrations) - ! This subroutine is based on the E3SM implementation. -bpm - ! For each gas name in gases, initialize that gas in gas_concentrations. - use mo_gas_concentrations, only: ty_gas_concs - use mo_rrtmgp_util_string, only: lower_case - ! Arguments - type(ty_gas_concs), intent(inout) :: gas_concentrations - character(len=*), intent(in) :: gases(:) - ! Local - character(len=32), dimension(size(gases)) :: gases_lowercase - integer :: igas - character(len=128) :: error_msg - ! Initialize with lowercase gas names; we should work in lowercase - ! whenever possible because we cannot trust string comparisons in RRTMGP - ! to be case insensitive ... it *should* work regardless of case. - do igas = 1,size(gases) - gases_lowercase(igas) = trim(lower_case(gases(igas))) - end do - error_msg = gas_concentrations%init(gases_lowercase) - if (len_trim(error_msg) > 0) then - call endrun('Setting available gases. ERROR: '//trim(error_msg)) - end if -end subroutine set_available_gases - +!========================================================================================= subroutine reset_fluxes(fluxes) diff --git a/src/physics/rrtmgp/rrtmgp_inputs.F90 b/src/physics/rrtmgp/rrtmgp_inputs.F90 index 116093add4..6823d5aaa0 100644 --- a/src/physics/rrtmgp/rrtmgp_inputs.F90 +++ b/src/physics/rrtmgp/rrtmgp_inputs.F90 @@ -5,11 +5,6 @@ module rrtmgp_inputs ! RRTMGP. Subset the number of model levels if CAM's top exceeds RRTMGP's ! valid domain. ! -! This code is currently set up to send RRTMGP vertical layers ordered bottom -! to top of model. Although the RRTMGP is supposed to be agnostic about the -! vertical ordering problems have arisen trying to use the top to bottom order -! as used by CAM's infrastructure. -! !-------------------------------------------------------------------------------- use shr_kind_mod, only: r8=>shr_kind_r8 @@ -21,24 +16,20 @@ module rrtmgp_inputs use physics_buffer, only: physics_buffer_desc use camsrfexch, only: cam_in_t -use radconstants, only: get_ref_solar_band_irrad, rad_gas_index -use radconstants, only: nradgas, gaslist, rrtmg_to_rrtmgp_swbands -use rad_solar_var, only: get_variability -use solar_irrad_data, only : do_spctrl_scaling, sol_tsi +use radconstants, only: nswbands, nlwbands, get_sw_spectral_boundaries +use radconstants, only: nradgas, gaslist + use rad_constituents, only: rad_cnst_get_gas use mcica_subcol_gen, only: mcica_subcol_sw, mcica_subcol_lw use mo_gas_concentrations, only: ty_gas_concs -use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp -use mo_optical_props, only: ty_optical_props, ty_optical_props_2str, ty_optical_props_1scl +use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp +use mo_optical_props, only: ty_optical_props_2str, ty_optical_props_1scl -! unneeded use mo_rrtmgp_util_string, only: lower_case -use cam_logfile, only: iulog +use cam_logfile, only: iulog use cam_abortutils, only: endrun -use cam_history, only: outfld ! just for getting ozone VMR above model top. - implicit none private save @@ -71,12 +62,18 @@ module rrtmgp_inputs integer :: ktopcami ! cam index of top interface integer :: ktopradi ! rrtmgp index of interface corresponding to ktopcami +! wavenumber (cm^-1) boundaries of shortwave bands +real(r8) :: sw_low_bounds(nswbands), sw_high_bounds(nswbands) + !================================================================================================== contains !================================================================================================== subroutine rrtmgp_inputs_init(ktcamm, ktradm, ktcami, ktradi) + ! Note that this routine must be called after the calls to set_wavenumber_bands which set + ! the sw/lw band boundaries in the radconstants module. + integer, intent(in) :: ktcamm integer, intent(in) :: ktradm integer, intent(in) :: ktcami @@ -87,27 +84,26 @@ subroutine rrtmgp_inputs_init(ktcamm, ktradm, ktcami, ktradi) ktopcami = ktcami ktopradi = ktradi + call get_sw_spectral_boundaries(sw_low_bounds, sw_high_bounds, 'cm^-1') + end subroutine rrtmgp_inputs_init !================================================================================================== subroutine rrtmgp_set_state( & - pstate, cam_in, ncol, nlay, nlwbands, & - nswbands, ngpt_sw, nday, idxday, coszrs, & - kdist_sw, & ! eccf, & !!! Removing eccf from arguments, as it is not needed here + pstate, cam_in, ncol, nlay, & + nday, idxday, coszrs, & + kdist_sw, & band2gpt_sw, & t_sfc, emis_sfc, t_rad, & pmid_rad, pint_rad, t_day, pmid_day, pint_day, & - coszrs_day, alb_dir, alb_dif, tsi) + coszrs_day, alb_dir, alb_dif) ! arguments type(physics_state), target, intent(in) :: pstate type(cam_in_t), intent(in) :: cam_in integer, intent(in) :: ncol integer, intent(in) :: nlay - integer, intent(in) :: nlwbands - integer, intent(in) :: nswbands - integer, intent(in) :: ngpt_sw integer, intent(in) :: nday integer, intent(in) :: idxday(:) real(r8), intent(in) :: coszrs(:) @@ -127,10 +123,6 @@ subroutine rrtmgp_set_state( & real(r8), intent(out) :: coszrs_day(nday) ! cosine of solar zenith angle real(r8), intent(out) :: alb_dir(nswbands,nday) ! surface albedo, direct radiation real(r8), intent(out) :: alb_dif(nswbands,nday) ! surface albedo, diffuse radiation - ! real(r8), intent(out) :: solin(ncol) ! incident flux at domain top [W/m2] - ! real(r8), intent(out) :: solar_irrad_gpt(nday,ngpt_sw) ! incident flux at domain top per gpoint [W/m2] AT DAYLIT POINTS - ! real(r8), intent(out) :: tsi_scaling_gpt(ngpt_sw) ! scale factor for irradiance by gpoint [fraction] - real(r8), intent(out) :: tsi ! total irradiance W/m2 ! local variables integer :: k, kk, i, iband @@ -139,10 +131,6 @@ subroutine rrtmgp_set_state( & real(r8) :: sfac(nswbands) ! time varying scaling factors due to Solar Spectral ! Irrad at 1 A.U. per band - real(r8) :: wavenumber_limits(2,nswbands) - - ! real(r8) :: toa_flx_by_band(nswbands) ! temporary array of incoming flux by band - ! real(r8) :: toa_flx_by_gpt(ngpt_sw) ! temporary array of incoming flux by gpt character(len=*), parameter :: sub='rrtmgp_set_state' character(len=512) :: errmsg @@ -192,65 +180,16 @@ subroutine rrtmgp_set_state( & coszrs_day(i) = coszrs(idxday(i)) end do - - ! total solar incident radiation - tsi = sol_tsi ! when using sol_tsi from solar_irrad_data, this is read from a file. - - ! TO BE REMOVED - ! We can get TSI from the solar forcing file (above). - ! We can't get the scaling here because we might not have access - ! to RRTMGP's reference irradiance on bands yet (without running kdist%gas_optics). - ! The scaling can be derived in rrtmgp_driver / rte_sw (after %gas_optics provides the toa_flux). - ! call get_ref_solar_band_irrad(solar_band_irrad) - ! call get_variability(sfac) - ! solar_band_irrad = solar_band_irrad(rrtmg_to_rrtmgp_swbands) - ! tsi = sum(solar_band_irrad(:)) ! total TSI integrated across bands, BUT NOT scaled for variability - ! ! convert from irradiance scale factor per band (sfac) to per gpoint - ! ! --> this can then be used in rrtmgp_driver module, rte_sw to scale TOA flux - ! tsi_scaling_gpt = 0.0 - - ! do iband = 1,nswbands - ! tsi_scaling_gpt(band2gpt_sw(1,iband):band2gpt_sw(2,iband)) = sfac(iband) - ! end do - - ! if we had a method to produce toa flux by gpoint, we could make that an output here. - - ! <-- begin: old way of setting albedo hard-wired to 14 SW bands --> - ! ! Surface albedo (band mapping is hardcoded for RRTMG(P) code) - ! ! This mapping assumes nswbands=14. - ! if (nswbands /= 14) & - ! call endrun(sub//': ERROR: albedo band mapping assumes nswbands=14') - - ! do i = 1, nday - ! ! Near-IR bands (1-9 and 14), 820-16000 cm-1, 0.625-12.195 microns - ! alb_dir(1:8,i) = cam_in%aldir(idxday(i)) - ! alb_dif(1:8,i) = cam_in%aldif(idxday(i)) - ! alb_dir(14,i) = cam_in%aldir(idxday(i)) - ! alb_dif(14,i) = cam_in%aldif(idxday(i)) - - ! ! Set band 24 (or, band 9 counting from 1) to use linear average of UV/visible - ! ! and near-IR values, since this band straddles 0.7 microns: - ! alb_dir(9,i) = 0.5_r8*(cam_in%aldir(idxday(i)) + cam_in%asdir(idxday(i))) - ! alb_dif(9,i) = 0.5_r8*(cam_in%aldif(idxday(i)) + cam_in%asdif(idxday(i))) - - ! ! UV/visible bands 25-28 (10-13), 16000-50000 cm-1, 0.200-0.625 micron - ! alb_dir(10:13,i) = cam_in%asdir(idxday(i)) - ! alb_dif(10:13,i) = cam_in%asdif(idxday(i)) - ! enddo - ! <-- end: old way of setting albedo hard-wired to 14 SW bands --> - - ! More flexible way to assign albedo (from E3SM implementation) - ! adapted here to loop over bands and cols b/c cam_in has all cols but albedos are daylit cols - ! We could remove cols loop if we just set albedos for all columns separate from rrtmgp_set_state. - ! Albedos are input as broadband (visible, and near-IR), and we need to map - ! these to appropriate bands. Bands are categorized broadly as "visible" or - ! "infrared" based on wavenumber, so we get the wavenumber limits here - wavenumber_limits = kdist_sw%get_band_lims_wavenumber() + ! Assign albedos to the daylight columns (from E3SM implementation) + ! Albedos are imported from the surface models as broadband (visible, and near-IR), + ! and we need to map these to appropriate narrower bands used in RRTMGP. Bands + ! are categorized broadly as "visible/UV" or "infrared" based on wavenumber. ! Loop over bands, and determine for each band whether it is broadly in the - ! visible or infrared part of the spectrum (visible or "not visible") + ! visible or infrared part of the spectrum based on a dividing line of + ! 0.7 micron, or 14286 cm^-1 do iband = 1,nswbands - if (is_visible(wavenumber_limits(1,iband)) .and. & - is_visible(wavenumber_limits(2,iband))) then + if (is_visible(sw_low_bounds(iband)) .and. & + is_visible(sw_high_bounds(iband))) then ! Entire band is in the visible do i = 1, nday @@ -258,8 +197,8 @@ subroutine rrtmgp_set_state( & alb_dif(iband,i) = cam_in%asdif(idxday(i)) end do - else if (.not.is_visible(wavenumber_limits(1,iband)) .and. & - .not.is_visible(wavenumber_limits(2,iband))) then + else if (.not.is_visible(sw_low_bounds(iband)) .and. & + .not.is_visible(sw_high_bounds(iband))) then ! Entire band is in the longwave (near-infrared) do i = 1, nday alb_dir(iband,i) = cam_in%aldir(idxday(i)) @@ -276,7 +215,6 @@ subroutine rrtmgp_set_state( & end if end do - ! Strictly enforce albedo bounds where (alb_dir < 0) alb_dir = 0.0_r8 @@ -292,19 +230,21 @@ subroutine rrtmgp_set_state( & end where end subroutine rrtmgp_set_state -! -! Function to check if a wavenumber is in the visible or IR +!================================================================================================== + logical function is_visible(wavenumber) + ! Wavenumber is in the visible if it is above the visible threshold + ! wavenumber, and in the infrared if it is below the threshold + ! This function doesn't distinquish between visible and UV. + ! wavenumber in inverse cm (cm^-1) real(r8), intent(in) :: wavenumber ! Threshold between visible and infrared is 0.7 micron, or 14286 cm^-1 real(r8), parameter :: visible_wavenumber_threshold = 14286._r8 ! cm^-1 - ! Wavenumber is in the visible if it is above the visible threshold - ! wavenumber, and in the infrared if it is below the threshold if (wavenumber > visible_wavenumber_threshold) then is_visible = .true. else @@ -313,29 +253,29 @@ logical function is_visible(wavenumber) end function is_visible - !================================================================================================== + function get_molar_mass_ratio(gas_name) result(massratio) ! return the molar mass ratio of dry air to gas based on gas_name character(len=*),intent(in) :: gas_name real(r8) :: massratio select case (trim(gas_name)) - case ('h2o', 'H2O') + case ('H2O') massratio = 1.607793_r8 - case ('co2', 'CO2') + case ('CO2') massratio = 0.658114_r8 - case ('o3', 'O3') + case ('O3') massratio = 0.603428_r8 - case ('ch4', 'CH4') + case ('CH4') massratio = 1.805423_r8 - case ('n2o', 'N2O') + case ('N2O') massratio = 0.658090_r8 - case ('o2', 'O2') + case ('O2') massratio = 0.905140_r8 - case ('cfc11', 'CFC11') + case ('CFC11') massratio = 0.210852_r8 - case ('cfc12', 'CFC12') + case ('CFC12') massratio = 0.239546_r8 case default call endrun("Invalid gas: "//trim(gas_name)) @@ -379,7 +319,7 @@ subroutine rad_gas_get_vmr(icall, gas_name, pstate, pbuf, nlay, numactivecols, g mmr = gas_mmr ! special case: H2O is specific humidity, not mixing ratio. Use r = q/(1-q): - if (gas_name == 'h2o') then + if (gas_name == 'H2O') then mmr = mmr / (1._r8 - mmr) end if @@ -457,7 +397,7 @@ subroutine rad_gas_get_vmr(icall, gas_name, pstate, pbuf, nlay, numactivecols, g errmsg = gas_concs%set_vmr(gas_name, gas_vmr) if (len_trim(errmsg) > 0) then - call endrun(sub//': error setting CO2: '//trim(errmsg)) + call endrun(sub//': ERROR, gas_concs%set_vmr: '//trim(errmsg)) end if deallocate(gas_vmr) diff --git a/src/physics/rrtmgp/slingo.F90 b/src/physics/rrtmgp/slingo.F90 index aedb44bcee..64d614365e 100644 --- a/src/physics/rrtmgp/slingo.F90 +++ b/src/physics/rrtmgp/slingo.F90 @@ -9,7 +9,7 @@ module slingo use ppgrid, only: pcols, pver, pverp use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx -use radconstants, only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag, get_sw_spectral_boundaries +use radconstants, only: nswbands, nlwbands, get_sw_spectral_boundaries use cam_abortutils, only: endrun use cam_history, only: outfld @@ -80,20 +80,6 @@ subroutine slingo_rad_props_init() call cnst_get_ind('CLDLIQ', ixcldliq) call cnst_get_ind('CLDICE', ixcldice) - !call addfld ('CLWPTH_OLD',(/ 'lev' /), 'I','Kg/m2','old In Cloud Liquid Water Path', sampling_seq='rad_lwsw') - !call addfld ('KEXT_OLD',(/ 'lev' /),'I','m^2/kg','old extinction') - !call addfld ('CLDOD_OLD',(/ 'lev' /),'I','1','old liquid OD') - !call addfld ('REL_OLD',(/ 'lev' /),'I','1','old liquid effective radius (liquid)') - - !call addfld ('CLWPTH_NEW',(/ 'lev' /), 'I','Kg/m2','In Cloud Liquid Water Path', sampling_seq='rad_lwsw') - !call addfld ('KEXT_NEW',(/ 'lev' /),'I','m^2/kg','extinction') - !call addfld ('CLDOD_NEW',(/ 'lev' /),'I','1','liquid OD') - - !call addfld('CIWPTH_NEW',(/ 'lev' /), 'I','Kg/m2','In Cloud Ice Water Path', sampling_seq='rad_lwsw') - !call addfld('CIWPTH_OLD',(/ 'lev' /), 'I','Kg/m2','In Cloud Ice Water Path (old)', sampling_seq='rad_lwsw') - - return - end subroutine slingo_rad_props_init !============================================================================== @@ -318,12 +304,6 @@ subroutine slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, li end do ! End do k=1,pver end do ! nswbands - !call outfld('CL_OD_SW_OLD',liq_tau(idx_sw_diag,:,:), pcols, lchnk) - !call outfld('REL_OLD',rel(:,:), pcols, lchnk) - !call outfld('CLWPTH_OLD',cliqwp(:,:), pcols, lchnk) - !call outfld('KEXT_OLD',kext(:,:), pcols, lchnk) - - end subroutine slingo_liq_optics_sw subroutine slingo_liq_get_rad_props_lw(state, pbuf, abs_od, oldliqwp) diff --git a/src/physics/simple/radconstants.F90 b/src/physics/simple/radconstants.F90 index b69fac1552..60585713d6 100644 --- a/src/physics/simple/radconstants.F90 +++ b/src/physics/simple/radconstants.F90 @@ -15,8 +15,6 @@ module radconstants integer, parameter, public :: idx_lw_diag = 1 integer, parameter, public :: idx_nir_diag = 1 integer, parameter, public :: idx_uv_diag = 1 -integer, parameter, public :: nrh = 1 -integer, parameter, public :: ot_length = 32 public :: rad_gas_index public :: get_lw_spectral_boundaries