Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Radiation changes #54

Open
wants to merge 26 commits into
base: master
Choose a base branch
from
Open

Radiation changes #54

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
26 commits
Select commit Hold shift + click to select a range
bb68d39
Various changes ref. land and update to Isca master at Execlim
Dec 21, 2017
5ede95e
Merge remote-tracking branch 'upstream/master'
Jun 16, 2018
55db070
Add T21 spectral resolution to experiment.py
Jun 17, 2018
4e7caed
fix T21 in experiment.py
Jun 17, 2018
b1ecfed
eccentricity working in rrtm
Jun 21, 2018
682ecb4
added interp_end for h20 interpolation
Jun 28, 2018
980c033
Add option to input prescribed specific humidity to radiation
Jun 29, 2018
1d741fd
Missed endif
Jun 29, 2018
4fcece5
Added Nakajima 1992 radiation scheme
ntlewis Jun 29, 2018
e4240cd
Included comment containing units for Nakajima radiation scheme const…
ntlewis Jun 29, 2018
60a7253
eccentricity now working in two_stream_gray_rad
ntlewis Jun 29, 2018
2f2e702
Removed #experiment.py#
ntlewis Jun 29, 2018
519db28
Fix additional T21 resolution option
ntlewis Jul 2, 2018
c520a8b
Merge branch 'master' into radiation_changes
jamesp Sep 4, 2018
6ec026d
Merge remote-tracking branch 'upstream/master'
ntlewis Jan 11, 2019
78c0e38
eccentricity working in rrtm
Jun 21, 2018
762bab2
added interp_end for h20 interpolation
Jun 28, 2018
5418487
Add option to input prescribed specific humidity to radiation
Jun 29, 2018
360db40
Missed endif
Jun 29, 2018
aad628a
Added Nakajima 1992 radiation scheme
ntlewis Jun 29, 2018
245b701
Included comment containing units for Nakajima radiation scheme const…
ntlewis Jun 29, 2018
18819e6
eccentricity now working in two_stream_gray_rad
ntlewis Jun 29, 2018
33df7ba
Removed #experiment.py#
ntlewis Jun 29, 2018
e4c680a
remove garbage files
ntlewis Feb 27, 2019
834b028
resolve conflicts in merge
ntlewis Feb 27, 2019
a1ae71e
remove artefact from merge commit
ntlewis Feb 28, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 7 additions & 4 deletions src/atmos_param/rrtm_radiation/rrtm_radiation.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -581,7 +581,7 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt,
real(kind=rb),dimension(size(q,1),size(q,2)) :: fracsun

integer :: year_in_s
real :: r_seconds, r_days, r_total_seconds, frac_of_day, frac_of_year, gmt, time_since_ae, rrsun, dt_rad_radians, day_in_s, r_solday, r_dt_rad_avg
real :: r_seconds, r_days, r_total_seconds, frac_of_day, frac_of_year, gmt, time_since_ae, rrsun, dt_rad_radians, day_in_s, r_solday, r_dt_rad_avg, insolation



Expand Down Expand Up @@ -656,6 +656,7 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt,
! Seasonal Cycle: Use astronomical parameters to calculate insolation
call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun)
end if
insolation = solr_cnst * rrsun
ntlewis marked this conversation as resolved.
Show resolved Hide resolved

! input files: only deal with case where we don't need to call radiation at all
if(do_read_radiation .and. do_read_sw_flux .and. do_read_lw_flux) then
Expand Down Expand Up @@ -787,7 +788,7 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt,
pfull , phalf , tfull , thalf , tsrf , &
h2o , o3 , co2 , ch4_val*ones , n2o_val*ones , o2_val*ones , &
albedo_rr , albedo_rr, albedo_rr, albedo_rr, &
cosz_rr , solrad , dyofyr , solr_cnst, &
cosz_rr , solrad , dyofyr , insolation, &
inflglw , iceflglw , liqflglw , &
! cloud parameters
zeros , taucld , sw_zro , sw_zro , sw_zro , &
Expand All @@ -801,7 +802,7 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt,
pfull , phalf , tfull , thalf , tsrf , &
h2o , o3 , co2 , zeros , zeros, zeros, &
albedo_rr , albedo_rr, albedo_rr, albedo_rr, &
cosz_rr , solrad , dyofyr , solr_cnst, &
cosz_rr , solrad , dyofyr , insolation, &
inflglw , iceflglw , liqflglw , &
! cloud parameters
zeros , taucld , sw_zro , sw_zro , sw_zro , &
Expand Down Expand Up @@ -1057,12 +1058,14 @@ end subroutine write_diag_rrtm
!*****************************************************************************************

subroutine rrtm_radiation_end
use rrtm_vars, only: do_read_ozone,o3_interp, do_read_co2, co2_interp
use rrtm_vars, only: do_read_ozone,o3_interp, do_read_co2, co2_interp, &
do_read_h2o, h2o_interp
use interpolator_mod, only: interpolator_end
implicit none

if(do_read_ozone)call interpolator_end(o3_interp)
if(do_read_co2)call interpolator_end(co2_interp)
if(do_read_h2o)call interpolator_end(h2o_interp) ! NTL change here, interp_end for h2o

end subroutine rrtm_radiation_end

Expand Down
86 changes: 71 additions & 15 deletions src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@ module two_stream_gray_rad_mod
mpp_pe, close_file, error_mesg, &
NOTE, FATAL, uppercase

use constants_mod, only: stefan, cp_air, grav, pstd_mks, pstd_mks_earth, seconds_per_sol, orbital_period
use constants_mod, only: stefan, cp_air, grav, pstd_mks, pstd_mks_earth, seconds_per_sol, orbital_period, &
rdgas, rvgas

use diag_manager_mod, only: register_diag_field, send_data

Expand Down Expand Up @@ -88,7 +89,7 @@ module two_stream_gray_rad_mod
character(len=32) :: rad_scheme = 'frierson'

integer, parameter :: B_GEEN = 1, B_FRIERSON = 2, &
B_BYRNE = 3, B_SCHNEIDER_LIU=4
B_BYRNE = 3, B_SCHNEIDER_LIU=4, B_NAKAJIMA=5
integer, private :: sw_scheme = B_FRIERSON
integer, private :: lw_scheme = B_FRIERSON

Expand Down Expand Up @@ -117,6 +118,14 @@ module two_stream_gray_rad_mod
real :: bog_b = 1997.9
real :: bog_mu = 1.0

! parameters for Nakajima 1992 radiation scheme
real :: naka_kv = 0.01 !vapour absorption coefficient, units: m2 kg-1
real :: naka_kn = 0.0001 !non-condensable absorption coefficient, units: m2 kg-1

! constants for Nakajima 1992 radiation scheme
real :: naka_mv = 18e-3 ! vapour molecular weight, units: kg mol-1
real :: naka_mn = 18e-3 ! non-condensable molecular weight, units: kg mol-1

real, allocatable, dimension(:,:) :: insolation, p2, lw_tau_0, sw_tau_0 !s albedo now defined in mixed_layer_init
real, allocatable, dimension(:,:) :: b_surf, b_surf_gp
real, allocatable, dimension(:,:,:) :: b, tdt_rad, tdt_solar
Expand All @@ -136,6 +145,11 @@ module two_stream_gray_rad_mod
type(interpolate_type),save :: co2_interp ! use external file for co2
character(len=256) :: co2_file='co2' ! file name of co2 file to read
character(len=256) :: co2_variable_name='co2' ! file name of co2 file to read
!extras for reading in h2o concentration
logical :: do_read_h2o=.false.
type(interpolate_type),save :: h2o_interp ! same as for co2
character(len=256) :: h2o_file='h2o'



namelist/two_stream_gray_rad_nml/ solar_constant, del_sol, &
Expand All @@ -147,7 +161,9 @@ module two_stream_gray_rad_mod
window, carbon_conc, rad_scheme, &
do_read_co2, co2_file, co2_variable_name, solday, equinox_day, bog_a, bog_b, bog_mu, &
use_time_average_coszen, dt_rad_avg,&
diabatic_acce !Schneider Liu values
diabatic_acce, & !Schneider Liu values
do_read_h2o, h2o_file, & ! for reading h2o
naka_kn, naka_kv ! nakajima rad absorption coefficients

!==================================================================================
!-------------------- diagnostics fields -------------------------------
Expand Down Expand Up @@ -209,6 +225,10 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb
call interpolator_init (co2_interp, trim(co2_file)//'.nc', lonb, latb, data_out_of_bounds=(/ZERO/))
endif

if(do_read_h2o)then
call interpolator_init(h2o_interp, trim(h2o_file)//'.nc', lonb, latb, data_out_of_bounds=(/ZERO/))
endif


if(uppercase(trim(rad_scheme)) == 'GEEN') then
lw_scheme = B_GEEN
Expand All @@ -223,6 +243,9 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb
else if(uppercase(trim(rad_scheme)) == 'SCHNEIDER') then
lw_scheme = B_SCHNEIDER_LIU
call error_mesg('two_stream_gray_rad','Using Schneider & Liu (2009) radiation scheme for GIANT PLANETS.', NOTE)
else if (uppercase(trim(rad_scheme)) == 'NAKAJIMA') then
lw_scheme = B_NAKAJIMA
call error_mesg('two_stream_gray_rad','Using Nakajima (1992) radiation scheme.', NOTE)
else
call error_mesg('two_stream_gray_rad','"'//trim(rad_scheme)//'"'//' is not a valid radiation scheme.', FATAL)
endif
Expand Down Expand Up @@ -286,7 +309,7 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb
allocate (del_sol_tau (ie-is+1, je-js+1))
allocate (sw_tau_k (ie-is+1, je-js+1))

case(B_BYRNE)
case(B_BYRNE, B_NAKAJIMA)
allocate (lw_del_tau (ie-is+1, je-js+1))

case(B_SCHNEIDER_LIU)
Expand Down Expand Up @@ -383,7 +406,7 @@ end subroutine two_stream_gray_rad_init
! ==================================================================================

subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, &
net_surf_sw_down, surf_lw_down, albedo, q)
net_surf_sw_down, surf_lw_down, albedo, q_in)

! Begin the radiation calculation by computing downward fluxes.
! This part of the calculation does not depend on the surface temperature.
Expand All @@ -393,22 +416,36 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
real, intent(in), dimension(:,:) :: lat, lon, albedo
real, intent(out), dimension(:,:) :: net_surf_sw_down
real, intent(out), dimension(:,:) :: surf_lw_down
real, intent(in), dimension(:,:,:) :: t, q, p_half
real, intent(in), dimension(:,:,:) :: t, q_in, p_half
integer :: i, j, k, n, dyofyr

integer :: seconds, year_in_s, days
real :: r_seconds, frac_of_day, frac_of_year, gmt, time_since_ae, rrsun, day_in_s, r_solday, r_total_seconds, r_days, r_dt_rad_avg, dt_rad_radians
logical :: used


real ,dimension(size(q,1),size(q,2),size(q,3)) :: co2f
real ,dimension(size(q_in,1),size(q_in,2),size(q_in,3)) :: co2f, q ! q to be passed to radiation

n = size(t,3)



! albedo(:,:) = albedo_value !s albedo now set in mixed_layer_init.


! If prescribing specific humidity from input file, read in here
if(do_read_h2o)then
call interpolator( h2o_interp, Time_diag, p_half, q, trim(h2o_file))
else
q = q_in ! otherwise set to input model specific humidity
endif

! If prescribing co2 from input file, read in here
if(do_read_co2)then
call interpolator( co2_interp, Time_diag, p_half, co2f, trim(co2_variable_name))
carbon_conc = maxval(co2f) !Needs maxval just because co2f is a 3d array of a constant, so maxval is just a way to pick out one number
endif

! =================================================================================
! SHORTWAVE RADIATION

Expand Down Expand Up @@ -444,7 +481,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun)
end if

insolation = solar_constant * coszen
insolation = solar_constant * coszen * rrsun

else if (sw_scheme==B_SCHNEIDER_LIU) then
insolation = (solar_constant/pi)*cos(lat)
Expand Down Expand Up @@ -475,7 +512,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
sw_down(:,:,k+1) = sw_down(:,:,k) * sw_dtrans(:,:,k)
end do

case(B_FRIERSON, B_BYRNE)
case(B_FRIERSON, B_BYRNE, B_NAKAJIMA)
! Default: Frierson handling of SW radiation
! SW optical thickness
sw_tau_0 = (1.0 - sw_diff*sin(lat)**2)*atm_abs
Expand Down Expand Up @@ -515,10 +552,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
! longwave source function
b = stefan*t**4

if(do_read_co2)then
call interpolator( co2_interp, Time_diag, p_half, co2f, trim(co2_variable_name))
carbon_conc = maxval(co2f) !Needs maxval just because co2f is a 3d array of a constant, so maxval is just a way to pick out one number
endif





select case(lw_scheme)
Expand Down Expand Up @@ -569,7 +605,26 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
do k = 1, n
lw_down(:,:,k+1) = lw_down(:,:,k)*lw_dtrans(:,:,k) + b(:,:,k)*(1. - lw_dtrans(:,:,k))
end do

case(B_NAKAJIMA)
! ref: Nakajima, S., Hayashi, Y.-Y., Abe, Y.
! A Study on the "Runaway Greenhouse Effect" with a One-Dimensional
! Radiative-Convective Equilibrium Model
! J. Atmos. Sci. 49, 2256-2266 (1992).

do k = 1, n
lw_del_tau = 3.0/2.0 * (naka_kv*q(:,:,k)*naka_mv + naka_kn*(rdgas/rvgas - q(:,:,k))*naka_mn) &
/ (q(:,:,k)*naka_mv + (rdgas/rvgas - q(:,:,k))*naka_mn) * (p_half(:,:,k+1)-p_half(:,:,k)) / grav
lw_dtrans(:,:,k) = exp( - lw_del_tau )

end do

! compute downward longwave flux by integrating downward
lw_down(:,:,1) = 0.
do k = 1, n
lw_down(:,:,k+1) = lw_down(:,:,k)*lw_dtrans(:,:,k) + b(:,:,k)*(1. - lw_dtrans(:,:,k))
end do

case(B_FRIERSON)
! longwave optical thickness function of latitude and pressure
lw_tau_0 = ir_tau_eq + (ir_tau_pole - ir_tau_eq)*sin(lat)**2
Expand Down Expand Up @@ -685,7 +740,7 @@ subroutine two_stream_gray_rad_up (is, js, Time_diag, lat, p_half, t_surf, t, td
end do
lw_up = lw_up + lw_up_win

case(B_FRIERSON, B_BYRNE)
case(B_FRIERSON, B_BYRNE, B_NAKAJIMA)
! compute upward longwave flux by integrating upward
lw_up(:,:,n+1) = b_surf
do k = n, 1, -1
Expand Down Expand Up @@ -791,14 +846,15 @@ subroutine two_stream_gray_rad_end
if (sw_scheme.eq.B_GEEN) then
deallocate (sw_dtrans, sw_wv, del_sol_tau, sw_tau_k)
endif
if (lw_scheme.eq.B_BYRNE) then
if ((lw_scheme.eq.B_BYRNE).or.(lw_scheme.eq.B_NAKAJIMA)) then
deallocate (lw_del_tau)
endif
if(lw_scheme.eq.B_SCHNEIDER_LIU) then
deallocate (b_surf_gp)
endif

if(do_read_co2)call interpolator_end(co2_interp)
if(do_read_h2o)call interpolator_end(h2o_interp)

end subroutine two_stream_gray_rad_end

Expand Down