module radiation 10,11
!---------------------------------------------------------------------------------
! Purpose:
!
! Provides the CAM interface to the radiation code
!
! Revision history:
! May 2004, D. B. Coleman, Initial version of interface module.
! July 2004, B. Eaton, Use interfaces from new shortwave, longwave, and ozone modules.
! Feb 2005, B. Eaton, Add namelist variables and control of when calcs are done.
! Mar 2008, J. Kay, Add FLDS (downwelling LW rad at surface) as an outfld.
! May 2008, J. Kay, ADD FSUTOA (upwelling SW radiation at TOA) as an outfld.
! Sep 2009, R. Neale, Add FLDSC (clear sky downwelling LW rad at surface) as an outfld.
!---------------------------------------------------------------------------------
use shr_kind_mod
, only: r8=>shr_kind_r8
use spmd_utils, only: masterproc
use ppgrid, only: pcols, pver, pverp
use physics_types, only: physics_state, physics_ptend
use physconst, only: cpair, cappa
use time_manager, only: get_nstep
use abortutils, only: endrun
use error_messages, only: handle_err
use cam_control_mod, only: lambm0, obliqr, mvelpp, eccen
use scamMod, only: scm_crm_mode, single_column,have_cld,cldobs,&
have_clwp,clwpobs,have_tg,tground
use perf_mod, only: t_startf, t_stopf
use cam_logfile, only: iulog
implicit none
private
save
public :: &
radiation_register, &! registers radiation physics buffer fields
radiation_defaultopts, &! set default values of namelist variables in runtime_opts
radiation_setopts, &! set namelist values from runtime_opts
radiation_printopts, &! print namelist values to log
radiation_get, &! provide read access to private module data
radiation_nextsw_cday, &! calendar day of next radiation calculation
radiation_do, &! query which radiation calcs are done this timestep
radiation_init, &! calls radini
radiation_tend ! moved from radctl.F90
! Private module data
! Default values for namelist variables
integer :: iradsw = -1 ! freq. of shortwave radiation calc in time steps (positive)
! or hours (negative).
integer :: iradlw = -1 ! frequency of longwave rad. calc. in time steps (positive)
! or hours (negative).
integer :: iradae = -12 ! frequency of absorp/emis calc in time steps (positive)
! or hours (negative).
integer :: irad_always = 0 ! Specifies length of time in timesteps (positive)
! or hours (negative) SW/LW radiation will be
! run continuously from the start of an
! initial or restart run
!===============================================================================
contains
!===============================================================================
subroutine radiation_register 1,3
!-----------------------------------------------------------------------
!
! Register radiation fields in the physics buffer
!
!-----------------------------------------------------------------------
use phys_buffer, only: pbuf_times, pbuf_add
integer idx
call pbuf_add('QRS' , 'global', 1,pver,1, idx) ! shortwave radiative heating rate
call pbuf_add('QRL' , 'global', 1,pver,1, idx) ! longwave radiative heating rate
end subroutine radiation_register
!================================================================================================
subroutine radiation_defaultopts(iradsw_out, iradlw_out, iradae_out, irad_always_out) 1
!-----------------------------------------------------------------------
! Purpose: Return default runtime options
!-----------------------------------------------------------------------
integer, intent(out), optional :: iradsw_out
integer, intent(out), optional :: iradlw_out
integer, intent(out), optional :: iradae_out
integer, intent(out), optional :: irad_always_out
!-----------------------------------------------------------------------
if ( present(iradsw_out) ) iradsw_out = iradsw
if ( present(iradlw_out) ) iradlw_out = iradlw
if ( present(iradae_out) ) iradae_out = iradae
if ( present(irad_always_out) ) irad_always_out = irad_always
end subroutine radiation_defaultopts
!================================================================================================
subroutine radiation_setopts(dtime, nhtfrq, iradsw_in, iradlw_in, iradae_in, & 1,1
irad_always_in)
!-----------------------------------------------------------------------
! Purpose: Set runtime options
! *** NOTE *** This routine needs information about dtime (init by dycore)
! and nhtfrq (init by history) to do its checking. Being called
! from runtime_opts provides these values possibly before they
! have been set in the modules responsible for them.
!-----------------------------------------------------------------------
integer, intent(in) :: dtime ! timestep size (s)
integer, intent(in) :: nhtfrq ! output frequency of primary history file
integer, intent(in), optional :: iradsw_in
integer, intent(in), optional :: iradlw_in
integer, intent(in), optional :: iradae_in
integer, intent(in), optional :: irad_always_in
! Local
integer :: ntspdy ! no. timesteps per day
integer :: nhtfrq1 ! local copy of input arg nhtfrq
!-----------------------------------------------------------------------
if ( present(iradsw_in) ) iradsw = iradsw_in
if ( present(iradlw_in) ) iradlw = iradlw_in
if ( present(iradae_in) ) iradae = iradae_in
if ( present(irad_always_in) ) irad_always = irad_always_in
! Convert iradsw, iradlw and irad_always from hours to timesteps if necessary
if (iradsw < 0) iradsw = nint((-iradsw *3600._r8)/dtime)
if (iradlw < 0) iradlw = nint((-iradlw *3600._r8)/dtime)
if (irad_always < 0) irad_always = nint((-irad_always*3600._r8)/dtime)
! Convert iradae from hours to timesteps if necessary and check that
! iradae must be an even multiple of iradlw
if (iradae < 0) iradae = nint((-iradae*3600._r8)/dtime)
if (mod(iradae,iradlw)/=0) then
write(iulog,*)'radiation_setopts: iradae must be an even multiple of iradlw.'
write(iulog,*)' iradae = ',iradae,', iradlw = ',iradlw
call endrun('radiation_setopts: iradae must be an even multiple of iradlw.')
end if
! Do absorptivities/emissivities have to go on a restart dataset?
! The value of nhtfrq from the namelist may need conversion.
nhtfrq1 = nhtfrq
if (nhtfrq1 < 0) nhtfrq1 = nint((-nhtfrq1*3600._r8)/dtime)
ntspdy = nint(86400._r8/dtime)
if (nhtfrq1 /= 0) then
if (masterproc .and. mod(nhtfrq1,iradae)/=0) then
write(iulog,*)'radiation_setopts: *** NOTE: Extra overhead invoked putting', &
' a/e numbers on restart dataset. *** ', &
' To avoid, make mod(nhtfrq,iradae) = 0'
end if
else
if (masterproc) then
if (mod(ntspdy,iradae) /= 0 .or. iradae > ntspdy) then
write(iulog,*)'radiation_setopts: *** NOTE: Extra overhead invoked', &
' putting a/e numbers on restart dataset. ***'
write(iulog,*)' To avoid, make mod(timesteps per day,iradae)= 0'
end if
end if
end if
end subroutine radiation_setopts
!===============================================================================
subroutine radiation_get(iradsw_out, iradlw_out, iradae_out, irad_always_out)
!-----------------------------------------------------------------------
! Purpose: Provide access to private module data. (This should be eliminated.)
!-----------------------------------------------------------------------
integer, intent(out), optional :: iradsw_out
integer, intent(out), optional :: iradlw_out
integer, intent(out), optional :: iradae_out
integer, intent(out), optional :: irad_always_out
!-----------------------------------------------------------------------
if ( present(iradsw_out) ) iradsw_out = iradsw
if ( present(iradlw_out) ) iradlw_out = iradlw
if ( present(iradae_out) ) iradae_out = iradae
if ( present(irad_always_out) ) irad_always_out = irad_always
end subroutine radiation_get
!================================================================================================
subroutine radiation_printopts 1
!-----------------------------------------------------------------------
! Purpose: Print runtime options to log.
!-----------------------------------------------------------------------
if(irad_always /= 0) write(iulog,10) irad_always
write(iulog,20) iradsw,iradlw,iradae
10 format(' Execute SW/LW radiation continuously for the first ',i5, ' timestep(s) of this run')
20 format(' Frequency of Shortwave Radiation calc. (IRADSW) ',i5/, &
' Frequency of Longwave Radiation calc. (IRADLW) ',i5/, &
' Frequency of Absorptivity/Emissivity calc. (IRADAE) ',i5)
end subroutine radiation_printopts
!================================================================================================
function radiation_do(op, timestep) 8,2
!-----------------------------------------------------------------------
! Purpose: Returns true if the specified operation is done this timestep.
!-----------------------------------------------------------------------
character(len=*), intent(in) :: op ! name of operation
integer, intent(in), optional:: timestep
logical :: radiation_do ! return value
! Local variables
integer :: nstep ! current timestep number
!-----------------------------------------------------------------------
if (present(timestep)) then
nstep = timestep
else
nstep = get_nstep()
end if
select case (op)
case ('sw') ! do a shortwave heating calc this timestep?
radiation_do = nstep == 0 .or. iradsw == 1 &
.or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) &
.or. nstep <= irad_always
case ('lw') ! do a longwave heating calc this timestep?
radiation_do = nstep == 0 .or. iradlw == 1 &
.or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) &
.or. nstep <= irad_always
case ('absems') ! do an absorptivity/emissivity calculation this timestep?
radiation_do = nstep == 0 .or. iradae == 1 &
.or. (mod(nstep-1,iradae) == 0 .and. nstep /= 1)
case ('aeres') ! write absorptivity/emissivity to restart file this timestep?
radiation_do = mod(nstep,iradae) /= 0
case default
call endrun('radiation_do: unknown operation:'//op)
end select
end function radiation_do
!================================================================================================
real(r8) function radiation_nextsw_cday() 4,6
!-----------------------------------------------------------------------
! Purpose: Returns calendar day of next sw radiation calculation
!-----------------------------------------------------------------------
use time_manager, only: get_curr_calday, get_nstep, get_step_size
! Local variables
integer :: nstep ! timestep counter
logical :: dosw ! true => do shosrtwave calc
integer :: offset ! offset for calendar day calculation
integer :: dTime ! integer timestep size
real(r8):: calday ! calendar day of
!-----------------------------------------------------------------------
radiation_nextsw_cday = -1._r8
dosw = .false.
nstep = get_nstep()
dtime = get_step_size()
offset = 0
do while (.not. dosw)
nstep = nstep + 1
offset = offset + dtime
if (radiation_do('sw', nstep)) then
radiation_nextsw_cday = get_curr_calday(offset=offset)
dosw = .true.
end if
end do
if(radiation_nextsw_cday == -1._r8) then
call endrun('error in radiation_nextsw_cday')
end if
end function radiation_nextsw_cday
!================================================================================================
subroutine radiation_init() 1,106
!-----------------------------------------------------------------------
!
! Initialize the radiation parameterization, add fields to the history buffer
!
!-----------------------------------------------------------------------
use cam_history, only: addfld, add_default, phys_decomp
use physconst, only: gravit, cpair, epsilo, stebol, &
pstd, mwdry, mwco2, mwo3
use param_cldoptics, only: param_cldoptics_init
use radsw, only: radsw_init
use radlw, only: radlw_init
use radae, only: radae_init
use cloudsimulator, only: doisccp, cloudsimulator_init
use radconstants, only: radconstants_init
use radiation_data, only: init_rad_data
integer :: nstep ! current timestep number
!-----------------------------------------------------------------------
call radconstants_init()
call init_rad_data()
call radsw_init(gravit)
call radlw_init(gravit, stebol)
call radae_init(gravit, epsilo, stebol, pstd, mwdry, mwco2, mwo3)
call param_cldoptics_init
! "irad_always" is number of time steps to execute radiation continuously from start of
! initial OR restart run
nstep = get_nstep()
if ( irad_always > 0) then
nstep = get_nstep()
irad_always = irad_always + nstep
end if
if (doisccp) call cloudsimulator_init
! Shortwave radiation
call addfld ('SOLIN ','W/m2 ',1, 'A','Solar insolation',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('SOLL ','W/m2 ',1, 'A','Solar downward near infrared direct to surface',phys_decomp, &
sampling_seq='rad_lwsw')
call addfld ('SOLS ','W/m2 ',1, 'A','Solar downward visible direct to surface',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('SOLLD ','W/m2 ',1, 'A','Solar downward near infrared diffuse to surface',phys_decomp, &
sampling_seq='rad_lwsw')
call addfld ('SOLSD ','W/m2 ',1, 'A','Solar downward visible diffuse to surface',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('QRS ','K/s ',pver, 'A','Solar heating rate',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('QRSC ','K/s ',pver, 'A','Clearsky solar heating rate',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSNS ','W/m2 ',1, 'A','Net solar flux at surface',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSNT ','W/m2 ',1, 'A','Net solar flux at top of model',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSNTOA ','W/m2 ',1, 'A','Net solar flux at top of atmosphere',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSUTOA ','W/m2 ',1, 'A','Upwelling solar flux at top of atmosphere',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSNTOAC ','W/m2 ',1, 'A','Clearsky net solar flux at top of atmosphere',phys_decomp, &
sampling_seq='rad_lwsw')
call addfld ('FSDTOA ','W/m2 ',1, 'A','Downwelling solar flux at top of atmosphere',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSN200 ','W/m2 ',1, 'A','Net shortwave flux at 200 mb',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSN200C ','W/m2 ',1, 'A','Clearsky net shortwave flux at 200 mb',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSNTC ','W/m2 ',1, 'A','Clearsky net solar flux at top of model',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSNSC ','W/m2 ',1, 'A','Clearsky net solar flux at surface',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSDSC ','W/m2 ',1, 'A','Clearsky downwelling solar flux at surface',phys_decomp, &
sampling_seq='rad_lwsw')
call addfld ('FSDS ','W/m2 ',1, 'A','Downwelling solar flux at surface',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FUS ','W/m2 ',pverp,'I','Shortwave upward flux',phys_decomp)
call addfld ('FDS ','W/m2 ',pverp,'I','Shortwave downward flux',phys_decomp)
call addfld ('FUSC ','W/m2 ',pverp,'I','Shortwave clear-sky upward flux',phys_decomp)
call addfld ('FDSC ','W/m2 ',pverp,'I','Shortwave clear-sky downward flux',phys_decomp)
call addfld ('FSNIRTOA','W/m2 ',1, 'A','Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', &
phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSNRTOAC','W/m2 ',1, 'A','Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', &
phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FSNRTOAS','W/m2 ',1, 'A','Net near-infrared flux (>= 0.7 microns) at top of atmosphere',phys_decomp, &
sampling_seq='rad_lwsw')
call addfld ('SWCF ','W/m2 ',1, 'A','Shortwave cloud forcing',phys_decomp, sampling_seq='rad_lwsw')
call add_default ('SOLIN ', 1, ' ')
call add_default ('QRS ', 1, ' ')
call add_default ('FSNS ', 1, ' ')
call add_default ('FSNT ', 1, ' ')
call add_default ('FSDTOA ', 1, ' ')
call add_default ('FSNTOA ', 1, ' ')
call add_default ('FSUTOA ', 1, ' ')
call add_default ('FSNTOAC ', 1, ' ')
call add_default ('FSNTC ', 1, ' ')
call add_default ('FSNSC ', 1, ' ')
call add_default ('FSDSC ', 1, ' ')
call add_default ('FSDS ', 1, ' ')
call add_default ('SWCF ', 1, ' ')
if (single_column.and.scm_crm_mode) then
call add_default ('FUS ', 1, ' ')
call add_default ('FUSC ', 1, ' ')
call add_default ('FDS ', 1, ' ')
call add_default ('FDSC ', 1, ' ')
endif
! Longwave radiation
call addfld ('QRL ','K/s ',pver, 'A','Longwave heating rate',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('QRLC ','K/s ',pver, 'A','Clearsky longwave heating rate',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FLNS ','W/m2 ',1, 'A','Net longwave flux at surface',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FLDS ','W/m2 ',1, 'A','Downwelling longwave flux at surface',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FLNT ','W/m2 ',1, 'A','Net longwave flux at top of model',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FLUT ','W/m2 ',1, 'A','Upwelling longwave flux at top of model',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FLUTC ','W/m2 ',1, 'A','Clearsky upwelling longwave flux at top of model',phys_decomp, &
sampling_seq='rad_lwsw')
call addfld ('FLNTC ','W/m2 ',1, 'A','Clearsky net longwave flux at top of model',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FLN200 ','W/m2 ',1, 'A','Net longwave flux at 200 mb',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FLN200C ','W/m2 ',1, 'A','Clearsky net longwave flux at 200 mb',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FLNSC ','W/m2 ',1, 'A','Clearsky net longwave flux at surface',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FLDSC ','W/m2 ',1, 'A','Clearsky downwelling longwave flux at surface',phys_decomp, &
sampling_seq='rad_lwsw')
call addfld ('LWCF ','W/m2 ',1, 'A','Longwave cloud forcing',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('FUL ','W/m2 ',pverp,'I','Longwave upward flux',phys_decomp)
call addfld ('FDL ','W/m2 ',pverp,'I','Longwave downward flux',phys_decomp)
call addfld ('FULC ','W/m2 ',pverp,'I','Longwave clear-sky upward flux',phys_decomp)
call addfld ('FDLC ','W/m2 ',pverp,'I','Longwave clear-sky downward flux',phys_decomp)
call add_default ('QRL ', 1, ' ')
call add_default ('FLNS ', 1, ' ')
call add_default ('FLDS ', 1, ' ')
call add_default ('FLNT ', 1, ' ')
call add_default ('FLUT ', 1, ' ')
call add_default ('FLUTC ', 1, ' ')
call add_default ('FLNTC ', 1, ' ')
call add_default ('FLNSC ', 1, ' ')
call add_default ('FLDSC ', 1, ' ')
call add_default ('LWCF ', 1, ' ')
if (single_column.and.scm_crm_mode) then
call add_default ('FUL ', 1, ' ')
call add_default ('FULC ', 1, ' ')
call add_default ('FDL ', 1, ' ')
call add_default ('FDLC ', 1, ' ')
endif
! Heating rate needed for d(theta)/dt computation
call addfld ('HR ','K/s ',pver, 'A','Heating rate needed for d(theta)/dt computation',phys_decomp)
! (Almost) net radiative flux at surface, does not have lwup.
call addfld ('SRFRAD ','W/m2 ',1, 'A','Net radiative flux at surface',phys_decomp)
call add_default ('SRFRAD ', 1, ' ')
! Cloud variables
call addfld ('CLOUD ','fraction',pver, 'A','Cloud fraction',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('CLDTOT ','fraction',1, 'A','Vertically-integrated total cloud',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('CLDLOW ','fraction',1, 'A','Vertically-integrated low cloud',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('CLDMED ','fraction',1, 'A','Vertically-integrated mid-level cloud',phys_decomp, sampling_seq='rad_lwsw')
call addfld ('CLDHGH ','fraction',1, 'A','Vertically-integrated high cloud',phys_decomp, sampling_seq='rad_lwsw')
call add_default ('CLOUD ', 1, ' ')
call add_default ('CLDTOT ', 1, ' ')
call add_default ('CLDLOW ', 1, ' ')
call add_default ('CLDMED ', 1, ' ')
call add_default ('CLDHGH ', 1, ' ')
end subroutine radiation_init
!===============================================================================
subroutine radiation_tend(state,ptend,pbuf, & 1,101
cam_out, cam_in, &
landfrac,landm,icefrac,snowh, &
fsns, fsnt, flns, flnt, &
fsds, net_flx)
!-----------------------------------------------------------------------
!
! Purpose:
! Driver for radiation computation.
!
! Method:
! Radiation uses cgs units, so conversions must be done from
! model fields to radiation fields.
!
! Revision history:
! May 2004 D.B. Coleman Merge of code from radctl.F90 and parts of tphysbc.F90.
! 2004-08-09 B. Eaton Add pointer variables for constituents.
! 2004-08-24 B. Eaton Access O3 and GHG constituents from chem_get_cnst.
! 2004-08-30 B. Eaton Replace chem_get_cnst by rad_constituent_get.
!-----------------------------------------------------------------------
use phys_buffer, only: pbuf_size_max, pbuf_fld, pbuf_old_tim_idx, pbuf_get_fld_idx
use phys_grid, only: get_rlat_all_p, get_rlon_all_p
use param_cldoptics, only: param_cldoptics_calc
use physics_types, only: physics_state, physics_ptend
use time_manager, only: get_curr_calday
use camsrfexch_types,only: surface_state, srfflx_state
use cam_history, only: outfld
use cloudsimulator, only: doisccp, cloudsimulator_run
use radheat, only: radheat_tend
use ppgrid
use pspect
use physconst, only: cpair, stebol
use radconstants, only: nlwbands, nswbands
use radsw, only: radcswmx
use radlw, only: radclwmx
use rad_constituents, only: rad_cnst_get_clim_gas, rad_cnst_out
use aer_rad_props, only: aer_rad_props_get_clim_sw, aer_rad_props_get_clim_lw
use interpolate_data, only: vertinterp
use radiation_data, only: output_rad_data
! Arguments
real(r8), intent(in) :: landfrac(pcols) ! land fraction
real(r8), intent(in) :: landm(pcols) ! land fraction ramp
real(r8), intent(in) :: icefrac(pcols) ! land fraction
real(r8), intent(in) :: snowh(pcols) ! Snow depth (liquid water equivalent)
real(r8), intent(inout) :: fsns(pcols) ! Surface solar absorbed flux
real(r8), intent(inout) :: fsnt(pcols) ! Net column abs solar flux at model top
real(r8), intent(inout) :: flns(pcols) ! Srf longwave cooling (up-down) flux
real(r8), intent(inout) :: flnt(pcols) ! Net outgoing lw flux at model top
real(r8), intent(inout) :: fsds(pcols) ! Surface solar down flux
real(r8), intent(inout) :: net_flx(pcols)
type(physics_state), intent(in), target :: state
type(physics_ptend), intent(out) :: ptend
type(pbuf_fld), intent(inout) :: pbuf(pbuf_size_max)
type(surface_state), intent(inout) :: cam_out
type(srfflx_state), intent(in) :: cam_in
! Local variables
logical :: dosw, dolw, doabsems
integer nmxrgn(pcols) ! Number of maximally overlapped regions
real(r8) pmxrgn(pcols,pverp) ! Maximum values of pressure for each
! maximally overlapped region.
! 0->pmxrgn(i,1) is range of pressure for
! 1st region,pmxrgn(i,1)->pmxrgn(i,2) for
! 2nd region, etc
real(r8) emis(pcols,pver) ! Cloud longwave emissivity
real(r8) :: cldtau(pcols,pver) ! Cloud longwave optical depth
real(r8) :: cicewp(pcols,pver) ! in-cloud cloud ice water path
real(r8) :: cliqwp(pcols,pver) ! in-cloud cloud liquid water path
real(r8) cltot(pcols) ! Diagnostic total cloud cover
real(r8) cllow(pcols) ! " low cloud cover
real(r8) clmed(pcols) ! " mid cloud cover
real(r8) clhgh(pcols) ! " hgh cloud cover
real(r8) :: ftem(pcols,pver) ! Temporary workspace for outfld variables
integer itim, ifld
real(r8), pointer, dimension(:,:) :: rel ! liquid effective drop radius (microns)
real(r8), pointer, dimension(:,:) :: rei ! ice effective drop size (microns)
real(r8), pointer, dimension(:,:) :: cld ! cloud fraction
real(r8), pointer, dimension(:,:) :: concld ! convective cloud fraction
real(r8), pointer, dimension(:,:) :: qrs ! shortwave radiative heating rate
real(r8), pointer, dimension(:,:) :: qrl ! longwave radiative heating rate
real(r8) :: qrsc(pcols,pver) ! clearsky shortwave radiative heating rate
real(r8) :: qrlc(pcols,pver) ! clearsky longwave radiative heating rate
integer lchnk, ncol
real(r8) :: calday ! current calendar day
real(r8) :: clat(pcols) ! current latitudes(radians)
real(r8) :: clon(pcols) ! current longitudes(radians)
real(r8) coszrs(pcols) ! Cosine solar zenith angle
logical :: conserve_energy = .true. ! flag to carry (QRS,QRL)*dp across time steps
integer i, k ! index
integer :: istat
real(r8) solin(pcols) ! Solar incident flux
real(r8) fsntoa(pcols) ! Net solar flux at TOA
real(r8) fsutoa(pcols) ! Upwelling solar flux at TOA
real(r8) fsntoac(pcols) ! Clear sky net solar flux at TOA
real(r8) fsnirt(pcols) ! Near-IR flux absorbed at toa
real(r8) fsnrtc(pcols) ! Clear sky near-IR flux absorbed at toa
real(r8) fsnirtsq(pcols) ! Near-IR flux absorbed at toa >= 0.7 microns
real(r8) fsntc(pcols) ! Clear sky total column abs solar flux
real(r8) fsdtoa(pcols) ! Solar input = Flux Solar Downward Top of Atmosphere
real(r8) fsnsc(pcols) ! Clear sky surface abs solar flux
real(r8) fsdsc(pcols) ! Clear sky surface downwelling solar flux
real(r8) flut(pcols) ! Upward flux at top of model
real(r8) lwcf(pcols) ! longwave cloud forcing
real(r8) swcf(pcols) ! shortwave cloud forcing
real(r8) flutc(pcols) ! Upward Clear Sky flux at top of model
real(r8) flntc(pcols) ! Clear sky lw flux at model top
real(r8) fldsc(pcols) ! Clear sky down lw flux at srf
real(r8) flnsc(pcols) ! Clear sky lw flux at srf (up-down)
real(r8) fln200(pcols) ! net longwave flux interpolated to 200 mb
real(r8) fln200c(pcols) ! net clearsky longwave flux interpolated to 200 mb
real(r8) flds(pcols) ! Srf downwelling longwave flux
real(r8) fns(pcols,pverp) ! net shortwave flux
real(r8) fcns(pcols,pverp) ! net clear-sky shortwave flux
real(r8) fsn200(pcols) ! fns interpolated to 200 mb
real(r8) fsn200c(pcols) ! fcns interpolated to 200 mb
real(r8) fnl(pcols,pverp) ! net longwave flux
real(r8) fcnl(pcols,pverp) ! net clear-sky longwave flux
real(r8) pbr(pcols,pver) ! Model mid-level pressures (dynes/cm2)
real(r8) pnm(pcols,pverp) ! Model interface pressures (dynes/cm2)
real(r8) eccf ! Earth/sun distance factor
real(r8) lwupcgs(pcols) ! Upward longwave flux in cgs units
real(r8), parameter :: cgs2mks = 1.e-3_r8
real(r8), pointer, dimension(:,:) :: n2o ! nitrous oxide mass mixing ratio
real(r8), pointer, dimension(:,:) :: ch4 ! methane mass mixing ratio
real(r8), pointer, dimension(:,:) :: cfc11 ! cfc11 mass mixing ratio
real(r8), pointer, dimension(:,:) :: cfc12 ! cfc12 mass mixing ratio
real(r8), pointer, dimension(:,:) :: o3 ! Ozone mass mixing ratio
real(r8), pointer, dimension(:,:) :: o2 ! Oxygen mass mixing ratio
real(r8), dimension(pcols) :: o2_col ! column oxygen mmr
real(r8), pointer, dimension(:,:) :: co2 ! co2 mass mixing ratio
real(r8), dimension(pcols) :: co2_col_mean ! co2 column mean mmr
real(r8), pointer, dimension(:,:) :: sp_hum ! specific humidity
! Aerosol shortwave radiative properties
real(r8) :: aer_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth
real(r8) :: aer_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau
real(r8) :: aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol assymetry parameter * w * tau
real(r8) :: aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau
! Aerosol longwave absorption optical depth
real(r8) :: odap_aer(pcols,pver,nlwbands)
! Gathered indicies of day and night columns
! chunk_column_index = IdxDay(daylight_column_index)
integer :: Nday ! Number of daylight columns
integer :: Nnite ! Number of night columns
integer, dimension(pcols) :: IdxDay ! Indicies of daylight coumns
integer, dimension(pcols) :: IdxNite ! Indicies of night coumns
character(*), parameter :: name = 'radiation_tend'
!----------------------------------------------------------------------
lchnk = state%lchnk
ncol = state%ncol
calday = get_curr_calday()
itim = pbuf_old_tim_idx()
ifld = pbuf_get_fld_idx('CLD')
cld => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,itim)
ifld = pbuf_get_fld_idx('CONCLD')
concld => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,itim)
ifld = pbuf_get_fld_idx('QRS')
qrs => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
ifld = pbuf_get_fld_idx('QRL')
qrl => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
ifld = pbuf_get_fld_idx('REL')
rel => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk, 1)
ifld = pbuf_get_fld_idx('REI')
rei => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,lchnk, 1)
! For CRM, make cloud equal to input observations:
if (single_column.and.scm_crm_mode.and.have_cld) then
do k = 1,pver
cld(:ncol,k)= cldobs(k)
enddo
endif
! Cosine solar zenith angle for current time step
call get_rlat_all_p(lchnk, ncol, clat)
call get_rlon_all_p(lchnk, ncol, clon)
call zenith (calday, clat, clon, coszrs, ncol)
call output_rad_data( pbuf, state, cam_in, landm, coszrs )
! Gather night/day column indices.
Nday = 0
Nnite = 0
do i = 1, ncol
if ( coszrs(i) > 0.0_r8 ) then
Nday = Nday + 1
IdxDay(Nday) = i
else
Nnite = Nnite + 1
IdxNite(Nnite) = i
end if
end do
dosw = radiation_do('sw') ! do shortwave heating calc this timestep?
dolw = radiation_do('lw') ! do longwave heating calc this timestep?
doabsems = radiation_do('absems') ! do absorptivity/emissivity calc this timestep?
if (dosw .or. dolw) then
! Compute cloud water/ice paths and optical properties for input to radiation
call t_startf('cldoptics')
call param_cldoptics_calc(state, cld, landfrac, landm,icefrac, &
cicewp, cliqwp, emis, cldtau, rel, rei, pmxrgn, nmxrgn, snowh, pbuf)
call t_stopf('cldoptics')
! For CRM, make cloud liquid water path equal to input observations
if(single_column.and.scm_crm_mode.and.have_clwp)then
do k=1,pver
cliqwp(:ncol,k) = clwpobs(k)
end do
endif
! Get specific humidity
call rad_cnst_get_clim_gas('H2O', state, pbuf, sp_hum)
! Get ozone mass mixing ratio.
call rad_cnst_get_clim_gas('O3', state, pbuf, o3)
! Get CO2 mass mixing ratio and compute column mean values
call rad_cnst_get_clim_gas('CO2', state, pbuf, co2)
call calc_col_mean(state, co2, co2_col_mean)
! construct cgs unit reps of pmid and pint and get "eccf" - earthsundistancefactor
call radinp(ncol, state%pmid, state%pint, pbr, pnm, eccf)
! Solar radiation computation
if (dosw) then
call t_startf('rad_sw')
! Get Oxygen mass mixing ratio.
call rad_cnst_get_clim_gas('O2', state, pbuf, o2)
call calc_col_mean(state, o2, o2_col)
! Get aerosol radiative properties.
call t_startf('aero_optics_sw')
call aer_rad_props_get_clim_sw(state, pbuf, nnite, idxnite, &
aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f)
call t_stopf('aero_optics_sw')
call radcswmx(lchnk, &
ncol, pnm, pbr, sp_hum, o3, &
o2_col, cld, cicewp, cliqwp, rel, &
rei, eccf, coszrs, solin, &
cam_in%asdir, cam_in%asdif, cam_in%aldir, cam_in%aldif, nmxrgn, &
pmxrgn, qrs, qrsc, fsnt, fsntc, fsdtoa, &
fsntoa, fsutoa, fsntoac, fsnirt, fsnrtc, fsnirtsq, &
fsns, fsnsc, fsdsc, fsds, cam_out%sols, &
cam_out%soll, cam_out%solsd, cam_out%solld, fns, fcns, &
Nday, Nnite, IdxDay, IdxNite, co2_col_mean, &
aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f )
call t_stopf('rad_sw')
! Output net fluxes at 200 mb
call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, fsn200c)
call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns, fsn200)
!
! Convert units of shortwave fields needed by rest of model from CGS to MKS
!
do i=1,ncol
solin(i) = solin(i)*cgs2mks
fsds(i) = fsds(i)*cgs2mks
fsnirt(i)= fsnirt(i)*cgs2mks
fsnrtc(i)= fsnrtc(i)*cgs2mks
fsnirtsq(i)= fsnirtsq(i)*cgs2mks
fsnt(i) = fsnt(i) *cgs2mks
fsdtoa(i)= fsdtoa(i) *cgs2mks
fsns(i) = fsns(i) *cgs2mks
fsntc(i) = fsntc(i)*cgs2mks
fsnsc(i) = fsnsc(i)*cgs2mks
fsdsc(i) = fsdsc(i)*cgs2mks
fsntoa(i)=fsntoa(i)*cgs2mks
fsutoa(i)=fsutoa(i)*cgs2mks
fsntoac(i)=fsntoac(i)*cgs2mks
fsn200(i) = fsn200(i)*cgs2mks
fsn200c(i) = fsn200c(i)*cgs2mks
swcf(i)=fsntoa(i) - fsntoac(i)
end do
ftem(:ncol,:pver) = qrs(:ncol,:pver)/cpair
! Dump shortwave radiation information to history buffer (diagnostics)
call outfld('QRS ',ftem ,pcols,lchnk)
ftem(:ncol,:pver) = qrsc(:ncol,:pver)/cpair
call outfld('QRSC ',ftem ,pcols,lchnk)
call outfld('SOLIN ',solin ,pcols,lchnk)
call outfld('FSDS ',fsds ,pcols,lchnk)
call outfld('FSNIRTOA',fsnirt,pcols,lchnk)
call outfld('FSNRTOAC',fsnrtc,pcols,lchnk)
call outfld('FSNRTOAS',fsnirtsq,pcols,lchnk)
call outfld('FSNT ',fsnt ,pcols,lchnk)
call outfld('FSDTOA ',fsdtoa,pcols,lchnk)
call outfld('FSNS ',fsns ,pcols,lchnk)
call outfld('FSNTC ',fsntc ,pcols,lchnk)
call outfld('FSNSC ',fsnsc ,pcols,lchnk)
call outfld('FSDSC ',fsdsc ,pcols,lchnk)
call outfld('FSNTOA ',fsntoa,pcols,lchnk)
call outfld('FSUTOA ',fsutoa,pcols,lchnk)
call outfld('FSNTOAC ',fsntoac,pcols,lchnk)
call outfld('SOLS ',cam_out%sols ,pcols,lchnk)
call outfld('SOLL ',cam_out%soll ,pcols,lchnk)
call outfld('SOLSD ',cam_out%solsd ,pcols,lchnk)
call outfld('SOLLD ',cam_out%solld ,pcols,lchnk)
call outfld('FSN200 ',fsn200,pcols,lchnk)
call outfld('FSN200C ',fsn200c,pcols,lchnk)
call outfld('SWCF ',swcf ,pcols,lchnk)
end if ! dosw
! Longwave radiation computation
if (dolw) then
call t_startf("rad_lw")
! Convert upward longwave flux units to CGS
do i=1,ncol
lwupcgs(i) = cam_in%lwup(i)*1000._r8
if(single_column.and.scm_crm_mode.and.have_tg) &
lwupcgs(i) = 1000*stebol*tground(1)**4
end do
! Get gas phase constituents.
call rad_cnst_get_clim_gas('N2O', state, pbuf, n2o)
call rad_cnst_get_clim_gas('CH4', state, pbuf, ch4)
call rad_cnst_get_clim_gas('CFC11', state, pbuf, cfc11)
call rad_cnst_get_clim_gas('CFC12', state, pbuf, cfc12)
! absems requires lw absorption optical depth and transmission through aerosols
call t_startf('aero_optics_lw')
if (doabsems) call aer_rad_props_get_clim_lw(state, pbuf, odap_aer)
call t_stopf('aero_optics_lw')
call radclwmx(lchnk, ncol, doabsems, &
lwupcgs, state%t, sp_hum, o3, pbr, &
pnm, state%lnpmid, state%lnpint, n2o, ch4, &
cfc11, cfc12, cld, emis, pmxrgn, &
nmxrgn, qrl, qrlc, flns, flnt, flnsc, &
flntc, cam_out%flwds, fldsc, flut, flutc, &
fnl, fcnl, co2_col_mean, odap_aer)
call t_stopf("rad_lw")
! Output fluxes at 200 mb
call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, fln200)
call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, fln200c)
!
! Convert units of longwave fields needed by rest of model from CGS to MKS
!
do i=1,ncol
flnt(i) = flnt(i)*cgs2mks
flut(i) = flut(i)*cgs2mks
flutc(i) = flutc(i)*cgs2mks
flns(i) = flns(i)*cgs2mks
flds(i) = cam_out%flwds(i)*cgs2mks
fldsc(i) = fldsc(i)*cgs2mks
flntc(i) = flntc(i)*cgs2mks
fln200(i) = fln200(i)*cgs2mks
fln200c(i) = fln200c(i)*cgs2mks
flnsc(i) = flnsc(i)*cgs2mks
cam_out%flwds(i) = cam_out%flwds(i)*cgs2mks
lwcf(i)=flutc(i) - flut(i)
end do
! Dump longwave radiation information to history tape buffer (diagnostics)
call outfld('QRL ',qrl (:ncol,:)/cpair,ncol,lchnk)
call outfld('QRLC ',qrlc(:ncol,:)/cpair,ncol,lchnk)
call outfld('FLNT ',flnt ,pcols,lchnk)
call outfld('FLUT ',flut ,pcols,lchnk)
call outfld('FLUTC ',flutc ,pcols,lchnk)
call outfld('FLNTC ',flntc ,pcols,lchnk)
call outfld('FLNS ',flns ,pcols,lchnk)
call outfld('FLDS ',flds ,pcols,lchnk)
call outfld('FLNSC ',flnsc ,pcols,lchnk)
call outfld('FLDSC ',fldsc ,pcols,lchnk)
call outfld('LWCF ',lwcf ,pcols,lchnk)
call outfld('FLN200 ',fln200,pcols,lchnk)
call outfld('FLN200C ',fln200c,pcols,lchnk)
end if !dolw
! Output aerosol mmr
call rad_cnst_out(state, pbuf)
! Cloud cover diagnostics
! radctl can change pmxrgn and nmxrgn so cldsav needs to follow
! radctl.
!
call cldsav (lchnk, ncol, cld, state%pmid, cltot, &
cllow, clmed, clhgh, nmxrgn, pmxrgn)
!
! Dump cloud field information to history tape buffer (diagnostics)
!
call outfld('CLDTOT ',cltot ,pcols,lchnk)
call outfld('CLDLOW ',cllow ,pcols,lchnk)
call outfld('CLDMED ',clmed ,pcols,lchnk)
call outfld('CLDHGH ',clhgh ,pcols,lchnk)
call outfld('CLOUD ',cld ,pcols,lchnk) !This should be written out by stratiform package
if (doisccp) then
call cloudsimulator_run(state, cam_in%ts, concld, cld, cliqwp, &
cicewp, rel, rei, emis, coszrs )
end if
else ! if (dosw .or. dolw) then
! convert radiative heating rates from Q*dp to Q for energy conservation
if (conserve_energy) then
!DIR$ CONCURRENT
do k =1 , pver
!DIR$ CONCURRENT
do i = 1, ncol
qrs(i,k) = qrs(i,k)/state%pdel(i,k)
qrl(i,k) = qrl(i,k)/state%pdel(i,k)
end do
end do
end if
end if ! if (dosw .or. dolw) then
! Compute net radiative heating tendency
call radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, &
fsnt, flns, flnt, cam_in%asdir, net_flx)
! Compute heating rate for dtheta/dt
do k=1,pver
do i=1,ncol
ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5_r8/state%pmid(i,k))**cappa
end do
end do
call outfld('HR ',ftem ,pcols ,lchnk )
! convert radiative heating rates to Q*dp for energy conservation
if (conserve_energy) then
!DIR$ CONCURRENT
do k =1 , pver
!DIR$ CONCURRENT
do i = 1, ncol
qrs(i,k) = qrs(i,k)*state%pdel(i,k)
qrl(i,k) = qrl(i,k)*state%pdel(i,k)
end do
end do
end if
! Compute net surface radiative flux for use by surface temperature code.
! Note that units have already been converted to mks in RADCTL. Since
! fsns and flwds are in the buffer, array values will be carried across
! timesteps when the radiation code is not invoked.
cam_out%srfrad(:ncol) = fsns(:ncol) + cam_out%flwds(:ncol)
call outfld('SRFRAD ',cam_out%srfrad,pcols,lchnk)
end subroutine radiation_tend
!===============================================================================
subroutine radinp(ncol, pmid, pint, pmidrd, pintrd, eccf) 1,4
!-----------------------------------------------------------------------
!
! Purpose:
! Set latitude and time dependent arrays for input to solar
! and longwave radiation.
! Convert model pressures to cgs.
!
! Author: CCM1, CMS Contact J. Kiehl
!-----------------------------------------------------------------------
use shr_orb_mod
use time_manager, only: get_curr_calday
!------------------------------Arguments--------------------------------
!
! Input arguments
!
integer, intent(in) :: ncol ! number of atmospheric columns
real(r8), intent(in) :: pmid(pcols,pver) ! Pressure at model mid-levels (pascals)
real(r8), intent(in) :: pint(pcols,pverp) ! Pressure at model interfaces (pascals)
!
! Output arguments
!
real(r8), intent(out) :: pmidrd(pcols,pver) ! Pressure at mid-levels (dynes/cm*2)
real(r8), intent(out) :: pintrd(pcols,pverp) ! Pressure at interfaces (dynes/cm*2)
real(r8), intent(out) :: eccf ! Earth-sun distance factor
!
!---------------------------Local variables-----------------------------
!
integer i ! Longitude loop index
integer k ! Vertical loop index
real(r8) :: calday ! current calendar day
real(r8) :: delta ! Solar declination angle
!-----------------------------------------------------------------------
!
calday = get_curr_calday()
call shr_orb_decl (calday ,eccen ,mvelpp ,lambm0 ,obliqr , &
delta ,eccf)
!
! Convert pressure from pascals to dynes/cm2
!
do k=1,pver
do i=1,ncol
pmidrd(i,k) = pmid(i,k)*10.0_r8
pintrd(i,k) = pint(i,k)*10.0_r8
end do
end do
do i=1,ncol
pintrd(i,pverp) = pint(i,pverp)*10.0_r8
end do
end subroutine radinp
!===============================================================================
subroutine calc_col_mean(state, mmr_pointer, mean_value) 2,1
!-----------------------------------------------------------------------
!
! Radiation only knows how to work with the column-mean co2 value.
! Compute the column mean.
!
!-----------------------------------------------------------------------
use cam_logfile, only: iulog
type(physics_state), intent(in) :: state
real(r8), dimension(:,:), pointer :: mmr_pointer ! mass mixing ratio (lev)
real(r8), dimension(pcols), intent(out) :: mean_value ! column mean mmr
integer :: i, k, ncol
real(r8) :: ptot(pcols)
!-----------------------------------------------------------------------
ncol = state%ncol
mean_value = 0.0_r8
ptot = 0.0_r8
do k=1,pver
do i=1,ncol
mean_value(i) = mean_value(i) + mmr_pointer(i,k)*state%pdeldry(i,k)
ptot(i) = ptot(i) + state%pdeldry(i,k)
end do
end do
do i=1,ncol
mean_value(i) = mean_value(i) / ptot(i)
end do
end subroutine calc_col_mean
!===============================================================================
end module radiation